# Pseudospectrum of de Sitter black holes

Kyriakos Destounis1,2, Valentin Boyanov3, Rodrigo Panosso Macedo4
1Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 5,
00185, Roma, Italy  2INFN, Sezione di Roma, Piazzale Aldo Moro 2, 00185, Roma,
Italy  3CENTRA, Departamento de Física, Instituto Superior Técnico – IST,
Universidade de Lisboa, Avenida Rovisco Pais 1, 1049 Lisboa, Portugal  4 Niels
Bohr International Academy, Niels Bohr Institute, Blegdamsvej 17, 2100
Copenhagen, Denmark

###### Abstract

Pseudospectral analyses have broadened our understanding of ringdown waveforms
from binary remnants, by providing insight into both the stability of their
characteristic frequencies under environmental perturbations, as well as the
underlying transient and non-modal phenomenology that a mode analysis may
miss. In this work we present the pseudospectrum of scalar perturbations on
spherically-symmetric black holes in de Sitter spacetimes. We expand upon
previous analyses in this setting by calculating the pseudospectrum of
Reissner-Nordström-de Sitter black holes, and revisit results regarding the
stability of quasinormal modes under perturbations in several cases. Of
particular note is the case of scalar quasinormal modes with angular parameter
$\ell=0$, which possess a zero mode related to the presence of a cosmological
horizon. We show that the non-trivial eigenfunction associated to this mode
has a vanishing energy norm which poses a challenge in quantifying the
magnitude of external perturbations to the wave equation’s potential, as well
as in calculating the pseudospectrum. Nonetheless, we present results which
suggest that the spectral instability manifestation of $\ell=0$ scalar
quasinormal modes is qualitatively the same as in other cases, in contrast to
recent claims. We also analyze the stability of the fundamental mode for
$\ell\geq 1$, finding it to be spectrally stable, except for certain
configurations in which a perturbation leads to a discontinuous overtaking of
the fundamental unperturbed purely-imaginary mode by a perturbed complex
quasinormal mode.

##  I Introduction

One of the most important tests of general relativity (GR) in the strong-field
regime is the detection of gravitational waves (GWs) from coalescing binary
compact objects [1]. In particular, the ringdown phase of the remnant has a
short burst during which it relaxes to stationarity through the emission of an
exponentially damped oscillatory signal. This signal is characteristic to the
remnant and corresponds to its resonant frequencies, i.e. its quasinormal
modes (QNMs) [2, 3, 4, 5].

In practice, a GW measurement is analyzed by comparing it to templates
generated from theoretical calculations in GR, or in other, modified theories
of gravity. Such calculations are typically performed with certain simplifying
assumptions, the most common one being that the coalescing objects form an
isolated gravitational system. In other words, the astrophysical environment
both at small scales (other gravitating objects in the vicinity of the binary)
and at large scales (the cosmological background) are neglected. At first
glance this assumption appears to be justified, as going beyond it does not
seem to affect the dominant part of the resulting signals, i.e. the part
resolved by current GW detectors. However, an analysis in the frequency domain
reveals that the characteristics of the underlying system can change
completely, altering the sub-dominant and late time aspects of the signals.

At large scales, it is known that changing the asymptotic structure of
spacetime, where the system resides, changes some key features of the QNM
spectrum. Particularly, in asymptotically flat spacetimes the discrete
spectrum is accompanied by a continuous part corresponding to a branch cut in
the Green’s function of the evolution operator, which translates into a late-
time inverse-polynomial decay of GWs [6, 7, 8, 9, 10, 11]. On the other hand,
if one includes a cosmological constant [12], even the small one which is
seemingly present in our universe [13, 14, 15], the presence of a cosmological
horizon changes the late-time decay to one which is exponentially damped [16,
17, 18, 19, 20], replacing the branch cut with a discrete set of non-
oscillatory, purely damped modes.

At small scales, changes to the background of the system seem to affect the
QNM spectrum even more significantly. It was first observed in Refs. [21, 22,
23] that slight small-scale changes to the effective potential in the wave
equation which governs the evolution of perturbations leads to a
disproportionately large displacement of the QNM frequencies in the complex
plane. This then led to the seminal work on the spectral instability of QNMs
in Schwarzschild black holes (BHs) [24]. A systematic analysis of the effect
that perturbations to the effective potential have on the QNM spectrum was
performed, accompanied by a calculation of the corresponding pseudospectrum
[25]. The instability was thus quantified by relating the migration of the QNM
frequencies to the energy scale of the background perturbation as measured by
the energy norm [26]. It was discovered that, among perturbations which
preserve the asymptotic structure, those concentrated on small scales affect
the spectrum the most, thus all overtones can be destabilized by them but the
fundamental mode that seems to remain stable. This picture seems to be in tune
with the BH spectroscopy interpretation of the dominant part of the ringdown
signals observed corresponding to the fundamental mode [27, 28], but it also
serves a warning for results which may rely on stable overtones [29, 30] (see
Ref. [31] for a review on BH spectroscopy and spectral instabilities). We also
note that it has been found that even the fundamental QNM can be
disproportionately destabilized [32, 33] when seemingly “small” envinronmental
bumps are added to the effective potential of the system at large distances
[34, 35, 36, 37], though the “smallness” of these perturbations was not
quantified in the same manner.

The overtone instability, as measured by the pseudospectrum in the energy
norm, was subsequently analyzed in a variety of spacetimes with different
asymptotic structures and boundaries, such as Reissner-Nordström (RN) [38],
Schwarzschild-de Sitter (SdS) [39], Schwarzschild-anti de Sitter black branes
[40] and BHs [41, 42], as well as horizonless exotic compact objects with a
reflective surface [43]. Each of these spacetimes, even before being
perturbed, have an overall QNM spectrum which is quite different from the
Schwarzschild case alone. In all these cases, the qualitative picture of
spectral instability is similar, though there are certain differences. For
instance, in the exotic compact object case, the stability extends past the
fundamental mode and includes the first overtones, when the object is
sufficiently compact. Another series of differences were observed in the SdS
case, which we aim to revisit in the present work, also extending part of the
analysis to the Reissner-Nordström-de Sitter (RNdS) family of BH solutions.

The RNdS case is of particular importance, given the richness of its spectral
content. Particularly, its QNMs are typically categorized into three distinct
families. Firstly, there is the family of _photon sphere_ (PS) QNMs,
associated to the angular frequency and instability timescale of null
geodesics at the PS [44, 45]. Second, there is a set of purely imaginary modes
(non-oscillatory damped modes) associated to the presence of the cosmological
horizon, referred to as _de Sitter_ (dS) modes [46, 47, 48, 49, 50]. Third,
there is another family of purely damped modes, for which the decay timescales
are typically very short (their imaginary parts are very large), except when
the BH is near extremality, i.e. when its charge is large enough for its inner
and outer horizons to be close to each other; this family is referred to as
the _near-extremal_ (NE) family of modes [51, 52, 48, 53, 54, 49, 55]. The
different behaviors of each of these families for different parameters of the
RNdS BH present an ideal testing ground for spectral and pseudospectral
analyses, since they represent different “large scale” settings and thus allow
for a more thorough analysis of environmental effects on the spectrum and its
(in)stability in different environments. Another important characteristic is
the lack of the aforementioned branch cut in the Green’s function, which, when
present, leads to numerically non-convergent eigenvalues on the imaginary
axis.

The goal of the present work is to revisit and expand upon previous work done
on the spectral instability of SdS and RNdS spacetimes [39]. We begin by
analyzing the case of $\ell=0$ scalar QNMs, where a surprisingly strong
spectral stability was reported in [39] for both SdS and RNdS geometries. We
show that this apparent stability is likely just an artefact of an issue with
the norm used to measure the size of perturbations. Particularly, the energy
norm of the zero-mode contained in the $\ell=0$ spectrum vanishes, implying
that the operator associated with the norm has a zero eigenvalue and thus its
inverse is ill-defined. As this inverse is used in calculations of the norms
of perturbations, the issue can lead to numerically non-convergent measures of
their size, and attempts to normalize them would in turn lead to practically
vanishing values. We show that under appropriately scaled perturbations,
overtones do indeed seem to be unstable. Additionally, we find that despite
the strict non-convergence of the norm, at a given finite numerical resolution
there are no divergences. The choice of a finite number of gridpoints
therefore serves as cutoff to the divergence, and can provide a qualitative
(but not quantitative) example of a well behaved norm. We find that in such a
norm, the pseudospectrum also suggests that the standard spectral instability
is present for $\ell=0$, though it does not serve to quantify it.111It is
worth noting that the pseudospectrum is in fact not used to directly quantify
the BH QNM instability, neither in this nor in previous works, due to
convergence issues which are more general in nature, as discussed in [42].

We then calculate the pseudospectra of scalar QNMs for various SdS and RNdS
BHs, for choices in the parameter space where different families of modes are
dominant, and for angular numbers $\ell\geq 1$, where the energy norm is well
behaved. Combined with calculations of the spectrum with particular
perturbations to the effective potential, we arrive at a picture of spectral
instability that is qualitatively similar in almost all cases. We find that
overtones are almost always unstable under relatively small perturbations, and
we quantify the stability of the fundamental mode. We find that the mode which
is initially fundamental, whether it be a PS mode or a purely imaginary one,
is by itself spectrally stable under perturbations of energy norm up to order
$10^{-1}$ (that is, its relative displacement is less than the size of the
perturbation), as hinted at in [39]. However, when the first overtone has an
imaginary part which is close enough to that of the fundamental mode that the
difference lies within the range of displacement by a small perturbation,
which can occur for particular values in the parameter space of dS BHs, an
_overtaking instability_ occurs, whereby the displacement of the overtone can
bring it below the the (previously) fundamental mode, making it the new
longest-lived one. Despite being dubbed an instability, we emphasize that it
can and does occur when both modes involved are individually stable, i.e. when
after the perturbation they each have a relative displacement less than the
magnitude of the perturbation. This brings us to the additional observation
that the only overtones which appear to be spectrally stable are the ones
which have an imaginary part very close to that of the (initially) fundamental
mode. We also note that the phenomenon of overtaking was also observed in Ref.
[32], where the effective potential was perturbed in a different manner.

##  II Asymptotically de-Sitter black holes

###  II.1 Schwarzschild coordinates

In the standard Schwarzschild coordinates $(t,r,\theta,\varphi)$, the
spherically-symmetric electrovacuum solution of the Einstein-Maxwell equations
with a positive cosmological constant, i.e. the RNdS geometry [56], has a line
element of the form

$ds^{2}=-f(r)dt^{2}+f(r)^{-1}dr^{2}+r^{2}d\Omega^{2},$ (1)

where $d\Omega^{2}=d\theta^{2}+\text{sin}^{2}\theta\,d\varphi^{2}$ is the line
element of the unit 2-sphere and the redshift function $f(r)$ reads

$f(r)=1-\frac{2M}{r}+\dfrac{Q^{2}}{r^{2}}-\frac{\Lambda r^{2}}{3}.\\\ $ (2)

Here, $M$ and $Q$ are the BH mass and electric charge, respectively, and
$\Lambda>0$ is the cosmological constant. By setting $Q=0$ we obtain the SdS
spacetime [56]. With the redshift function we define the tortoise coordinate
$r_{*}$ such that $f(r)=dr/dr_{*}$. Sub-extremal RNdS geometries possess three
horizons, which are located at the positive roots of $f(r)$. We will use
$r=r_{\Lambda}$ to designate the cosmological horizon, $r=r_{h}$ the event
horizon and $r=r_{c}$ the Cauchy horizon radius. For a given value of the
parameters $M$, $Q$ and $\Lambda$, these horizons are uniquely determined by
the relations

$\displaystyle M$ $\displaystyle=$
$\displaystyle\dfrac{(r_{h}+r_{\Lambda})(r_{h}+r_{c})(r_{c}+r_{\Lambda})}{2(r_{h}^{2}+r_{h}r_{\Lambda}+r_{h}r_{c}+\
r_{\Lambda}^{2}+r_{c}^{2}+r_{c}r_{\Lambda})},$ (3) $\displaystyle Q^{2}$
$\displaystyle=$
$\displaystyle\dfrac{r_{h}r_{\Lambda}r_{c}(r_{h}+r_{\Lambda}+r_{c})}{r_{h}^{2}+r_{h}r_{\Lambda}+r_{h}r_{c}+\
r_{\Lambda}^{2}+r_{c}^{2}+r_{c}r_{\Lambda}},$ (4) $\displaystyle\Lambda$
$\displaystyle=$
$\displaystyle\dfrac{3}{r_{h}^{2}+r_{h}r_{\Lambda}+r_{h}r_{c}+\
r_{\Lambda}^{2}+r_{c}^{2}+r_{c}r_{\Lambda}}.$ (5)

The redshift function $f(r)$ has one additional, negative root
$r_{o}=-(r_{h}+r_{c}+r_{\Lambda})$, which lies outside the range of $r$
corresponding to the spacetime description. Even though the most standard
representation of RNdS is given by Eq. (2), in what follows it is useful to
express it in its factorized form, with respect to the horizons and the
negative root, as

$f(r)=-r^{2}\dfrac{\Lambda}{3}\left(1-\dfrac{r_{h}}{r}\right)\left(1-\dfrac{r_{c}}{r}\right)\left(1-\dfrac{r_{\Lambda}}{r}\right)\left(1-\dfrac{r_{o}}{r}\right).$
(6)

This form is convenient for constructing the hyperboloidal slices [57], which
we will do in Sec. II.2. Furthermore, it also makes it easier to identify the
extremal regimes: extremal charge corresponds to $r_{h}=r_{c}$, and extremal
cosmological constant (i.e. the Nariai limit) to $r_{h}=r_{\Lambda}$.

Parametric studies on RNdS spacetimes may, therefore, follow two strategies:
(i) one fixes the parameters $M$, $Q$ and $\Lambda$ and them determines the
horizons $r_{h}$, $r_{c}$ and $r_{\Lambda}$ by solving Eqs. (3)-(5), or (ii)
one fixes the horizons $r_{h}$, $r_{c}$ and $r_{\Lambda}$ and then calculates
$M$, $Q$ and $\Lambda$ directly from Eqs. (3)-(5). The former is more
practical when expressing dimensional observables with respect to a typical
length scale such as the BH mass $M$. This is the more common approach in the
literature (see e.g. [48, 53, 54, 58, 49, 55, 59, 33]), which allows one to
prescribe dimensionless values such as $Q/M$ and $\Lambda M^{2}$. The
expressions for the hyperboloidal framework from Sec. II.2, on the other hand,
assume simpler forms if one employs the horizon $r_{h}$ as the typical length
scale. For that purpose, strategy (ii) is more convenient. Thus, while we
follow the typical convention of setting $Q/M$ and $\Lambda M^{2}$ when
choosing the parameters for the examples we use in Sec. III, in the numerical
results the QNM frequency space is represented in units of $r_{h}$ as a direct
consequence of the hyperboloidal framework which we will introduce.

###  II.2 Hyperboloidal framework

Following Ref. [24], together with the techniques reviewed in Ref. [57], we
introduce a hyperboloidal coordinate system that imposes geometrically the
boundary conditions associated to waves propagating into the BH and
cosmological horizons [60]. Specifically, the scri-fixing approach [61], with
definitions from Ref. [57] introduces hyperboloidal coordinates
$(\tau,\chi,\theta,\varphi)$ via

$t=\lambda\bigg{(}\tau-H(\chi)\bigg{)},\,\quad r=\dfrac{r_{h}}{\chi}.$ (7)

Here, $\lambda$ is some characteristic length scale of spacetime, and the
height function $H(\chi)$ is constructed by requiring that $\tau=$ constant
slices are spacelike hypersurfaces penetrating the future event and
cosmological horizons.

Following previous works [62, 63, 64, 57], we will choose a height function
corresponding to the so-called _minimal gauge_. The particular form of the
radial compactification already reflects a first step associated to this
choice. In terms of the compact radial coordinate $\chi$, the horizons are
located at $\chi_{h}=1$ (event horizon), $\chi_{c}=r_{h}/r_{c}$ (Cauchy
horizon) and $\chi_{\Lambda}=r_{h}/r_{\Lambda}$ (cosmological horizon). We
also identify the negative root of $f(r)$ as $\chi_{o}=r_{h}/r_{o}$. Due to
the inversion in the radial transformation, the horizons now satisfy the
relation $\chi_{\Lambda}\leq\chi_{h}\leq\chi_{c}$. We are interested in the
behavior of scalar fields propagating in the external region of a sub-extremal
BH, delimited by the domain of the observable universe, i.e.
$\chi\in[\chi_{\Lambda},\chi_{h}]$.

To fix $H(\chi)$ in the minimal gauge, it is convenient to first define the
dimensionless tortoise coordinate [57]

$\displaystyle x(\chi)$ $\displaystyle=$
$\displaystyle\dfrac{r_{*}(r(\chi))}{\lambda}$ (8) $\displaystyle=$
$\displaystyle\underbrace{x_{h}(\chi)+x_{c}(\chi)+x_{\Lambda}(\chi)}_{x_{\rm
sing}(\chi)}+x_{\rm reg}(\chi),$ (9)

with

$\displaystyle x_{\rm h}(\chi)$ $\displaystyle=$
$\displaystyle-\dfrac{3\,\chi_{h}}{\lambda
r_{h}\Lambda}\dfrac{\ln(\chi_{h}-\chi)}{(1-\chi_{h}/\chi_{c})(1-\chi_{h}/\chi_{\Lambda})(1-\chi_{h}/\chi_{o})},$
$\displaystyle x_{\rm c}(\chi)$ $\displaystyle=$
$\displaystyle-\dfrac{3\,\chi_{c}}{\lambda
r_{h}\Lambda}\dfrac{\ln(\chi_{c}-\chi)}{(1-\chi_{c}/\chi_{h})(1-\chi_{c}/\chi_{\Lambda})(1-\chi_{c}/\chi_{o})},$
$\displaystyle x_{\Lambda}(\chi)$ $\displaystyle=$
$\displaystyle-\dfrac{3\,\chi_{\Lambda}}{\lambda
r_{h}\Lambda}\dfrac{\ln(\chi-\chi_{\Lambda})}{(1-\chi_{\Lambda}/\chi_{h})(1-\chi_{\Lambda}/\chi_{\rm
C})(1-\chi_{\Lambda}/\chi_{o})},$ $\displaystyle x_{\rm reg}(\chi)$
$\displaystyle=$ $\displaystyle-\dfrac{3\chi_{o}}{\lambda
r_{h}\Lambda}\dfrac{\ln(\chi-\chi_{o})}{(1-\chi_{o}/\chi_{h})(1-\chi_{o}/\chi_{c})(1-\chi_{o}/\chi_{\Lambda})}.$

As pointed out in Eq. (9), the terms $x_{\rm h}(\chi)$, $x_{\rm c}(\chi)$ and
$x_{\Lambda}(\chi)$ contribute to the singular piece of $x(\chi)$ because they
diverge at $\chi_{h}$, $\chi_{c}$ and $\chi_{\Lambda}$, respectively. The
regular piece $x_{\rm reg}(\chi)$, resulting from the negative root
$\chi_{o}$, assumes a finite value in the entire domain, including all
horizons.

The height function is constructed straightforwardly [57] by simply changing
the signs of the $x_{\Lambda}(\chi)$ and $x_{\rm reg}(\chi)$, i.e.

$H(\chi)=x_{h}(\chi)+x_{c}(\chi)-x_{\Lambda}(\chi)-x_{\rm reg}(\chi).$ (14)

This practical approach is a direct consequence of the out-in strategy for
constructing hyperboloidal slices in the minimal gauge, and it ensures that
the surfaces $\tau=$ constant remain spacelike in the domain
$\chi\in[\chi_{\Lambda},\chi_{h}]$, for all parameters $Q/M$ and $\Lambda
M^{2}$ [57].

With the coordinate transformation (7), the line element (1) conformally
rescales as $ds^{2}=\Xi^{2}d\bar{s}^{2}$, with $\Xi=\sigma/\lambda$ and [57]

$d\bar{s}^{2}=-p(\chi)d\tau^{2}+2\gamma(\chi)d\tau
d\chi+w(\chi)d\chi^{2}+\dfrac{r_{h}^{2}}{\lambda^{2}}d\Omega^{2}.$ (15)

The metric functions relate to $x(\chi)$ and $H(\chi)$ via

$\displaystyle p=-\dfrac{1}{x^{\prime}},\quad\gamma=pH^{\prime},\quad
w=\dfrac{1-\gamma^{2}}{p}.$ (16)

In particular, $p(\chi)$ assumes the form

$p(\chi)=\dfrac{r_{h}^{2}\Lambda}{3}\left(1-\chi\right)\left(1-\dfrac{\chi}{\chi_{c}}\right)\left(\dfrac{\chi}{\chi_{\Lambda}}-1\right)\left(1-\dfrac{\chi}{\chi_{o}}\right),$
(17)

which will be explicitly used in Sec. III.1.1. The expressions for
$\gamma(\chi)$ and $w(\chi)$ are quite lengthy and thus we will omit them
here, though their explicit representations can be easily obtained.

###  II.3 Scalar perturbations and quasinormal modes

At linear order, the propagation of scalar perturbations $\Psi$ is governed by
the free Klein-Gordon equation

$\Box\Psi=\frac{1}{\sqrt{-g}}\left(g^{\mu\nu}\sqrt{-g}\Psi_{,\nu}\right)_{,\mu}=0,$
(18)

where $g$ is the determinant of the metric tensor of the background geometry
$g_{\mu\nu}$, and a comma in the subscript denotes partial differentiation. A
standard separation of variables of the form

$\Psi(t,r,\theta,\varphi)=\sum_{\ell,m}\frac{\psi_{\ell}(t,r)}{r}Y_{\ell,m}(\theta,\varphi),$
(19)

where $Y_{\ell,m}$ are the spherical harmonics, leads to the wave equation

$-\psi_{,tt}+\psi_{,r_{*}r_{*}}-V\psi=0.$ (20)

The effective potential $V$ for scalar perturbations has the form [4]

$\displaystyle
V=f(r)\left(\frac{\ell(\ell+1)}{r^{2}}+\frac{f^{\prime}(r)}{r}\right),$ (21)

where prime denotes differentiation with respect to the function’s variable.
While gravitational and electromagnetic perturbations can also be reduced to
solving an equation of the type (20), we find no qualitative difference in the
results presented below between these different types of perturbations, and we
thus limit the discussion to the scalar case, for which the spectrum is
already quite intricate [48] and leads to interesting phenomenology [39].

The region of interest for our analysis, as discussed, lies in the static
patch between the event horizon $r_{\rm h}$ and the cosmological horizon
$r_{\Lambda}$. This translates to solutions in the static region that satisfy
a purely ingoing boundary condition at the event horizon,

$\left(\psi_{,t}-\psi_{,r_{*}}\right)|_{r=r_{\rm h}}=0,$ (22)

and a purely outgoing boundary condition at the cosmological horizon,

$\left(\psi_{,t}+\psi_{,r_{*}}\right)|_{r=r_{\rm c}}=0.$ (23)

The vanishing of the effective potential $V$ at the horizons makes these
conditions compatible with the wave equation. The QNMs are then obtained by
considering a Fourier decomposition of the form
$\psi(t,r_{*})\sim\psi(r_{*})e^{i\omega t}$ and solving (20), along with the
aforementioned boundary conditions, as an eigenvalue problem. However, the
boundary value problem, as cast in the present coordinate system, becomes
singular at the horizons. As in previous works, we resolve this issue
geometrically by switching to the hyperboloidal slicing described in the
previous subsection.

###  II.4 Non-modal tools: pseudospectrum, energy norm and the eigenvalue
condition number

We follow Ref. [24] and cast the QNM eigenvalue problem in the hyperboloidal
frame in order to calculate the spectrum and pseudospectrum. The first step
consists of mapping the wave equation in $(t,r_{*})$ coordinates to its
equivalent form in the hyperboloidal $(\tau,\sigma)$ coordinates. For the sake
of completeness we rewrite the relevant expressions here (see Ref. [57] for a
systematic review). With a reduction of order in time
$\phi\equiv\psi_{,\tau}$, Eq. (20) can be rewritten as a system of two
equations,

$u_{,\tau}=iLu,\quad L=\frac{1}{i}\\!\left(\begin{array}[]{c c}0&1\\\
L_{1}&L_{2}\end{array}\right),\quad u=\left(\begin{array}[]{c}\psi\\\
\phi\end{array}\right),$ (24)

where we have defined the differential operators

$\begin{split}L_{1}&=\frac{1}{w(\chi)}\bigg{[}\partial_{\chi}\left(p(\chi)\partial_{\chi}\right)-q(\chi)\bigg{]},\\\
L_{2}&=\frac{1}{w(\chi)}\bigg{[}2\gamma(\chi)\partial_{\chi}+\partial_{\chi}\gamma(\chi)\bigg{]},\end{split}$
(25)

with the functions $p(\chi)$, $\gamma(\chi)$ and $w(\chi)$ given in Eq. (16).
In the hyperboloidal framework, the potential rescales to

$q=\lambda^{2}\dfrac{V}{p}.$ (26)

The discrete set of QNMs, $\omega_{n}$ with $n$ being the overtone number,
comprise the spectrum $\sigma(L)$ of the differential operator $L$ in (34).
Expanding this notion, the $\epsilon$-_pseudospectrum_
$\sigma^{(\epsilon)}(L)$ is defined as

$\sigma^{(\epsilon)}(L)=\\{\omega\in\mathbb{C}:\|R_{L}(\omega)\|>1/\epsilon\\},$
(27)

where $R_{L}(\omega)=(L-\omega\mathbb{I})^{-1}$ is the _resolvent operator_.
In the $\epsilon\to 0$ limit, the set of points $\omega$ reduces to the
spectrum $\omega_{n}$, while for larger values of $\epsilon$ the boundary of
the set corresponds to points in the complex plane which are “further away”
from being eigenvalues. The quantity $\epsilon$ provides a measure of this
“closeness” of points to the spectrum, which has a clear interpretation in
terms of perturbations to the underlying operator. If one assumes a
perturbation

$L\to L+\epsilon\delta L,$ (28)

where $\|\delta L\|=1$ and $\epsilon>0$, then the spectrum can migrate
anywhere inside the $\epsilon$-pseudospectral region of the same $\epsilon$.
Considering all possible perturbations in fact provides an alternative
definition of the pseudospectrum [25]. The shape and size of the
$\epsilon$-pseudospectral regions therefore quantifies the spectral stability
of $L$. Particularly, if their boundaries form concentric circles of radius
$\epsilon$, then the operator is spectrally stable, while other behaviors are
indicative of instability, as seen, e.g. in Refs. [24, 38].

We note that the pseudospectrum depends on the choice of norm with which one
computes $\|R_{L}(\omega)\|$, and $\|\delta L\|$. As in previous works [24,
38, 43], we will use the _energy norm_ [26], which provides a physical measure
of the magnitude of solutions of the wave equation, and by extension of
operators which can act on them. The definition of this norm for a
spherically-symmetric BH in dS spacetime, for a solution vector $u$ is given
by

$\|u\|^{2}_{{}_{E}}=\frac{1}{2}\int_{\chi_{{\rm\Lambda}}}^{\chi_{{\rm
h}}}\left(w(\chi)|\phi|^{2}+p(\chi)|\partial_{\chi}\psi|^{2}+q(\chi)|\psi|^{2}\right)d\chi,$
(29)

which allows one to introduce the scalar product between two vectors $u$ and
$v$ via

$\displaystyle\displaystyle|\langle
u,v\rangle|_{{}_{E}}=\frac{1}{2}\int_{\chi_{{\rm\Lambda}}}^{\chi_{{\rm
h}}}\Big{(}w(\chi)\phi_{v}^{*}\phi_{u}$
$\displaystyle+p(\chi)\partial_{\chi}\psi_{u}^{*}\,\partial_{\chi}\psi_{v}+q(\chi)\psi_{u}^{*}\,\psi_{v}\Big{)}d\chi.$
(30)

The pseudospectra of the $L$ operator which we calculate below are non-
trivial, that is, their contours do not form concentric circles around the
spectrum. This is a consequence of the non-selfadjoint (more generally, non-
normal) nature of $L$ in this norm. Another measure of this feature is the the
_eigenvalue condition number_ of the spectrum [25, 24]. For the QNM operator
$L$ and its eigenvalues $\omega_{n}$, the corresponding left and right
eigenvectors $u_{n}$ and $v_{n}$, respectively, are characterized by the
following eigenvalue expressions

$L^{\dagger}\,u_{n}=\tilde{\omega}_{n}\,u_{n},\qquad
L\,v_{n}=\omega_{n}\,u_{n},$ (31)

where $\tilde{\omega}_{n}$ is the complex conjugate of $\omega_{n}$. If we
perturb the operator, as in Eq. (28), then the perturbed eigenvalues
$\omega_{n}(\epsilon)$ satisfy [24]

$|\omega_{n}(\epsilon)-\omega_{n}|\leq\epsilon\,\kappa_{n},$ (32)

with the eigenvalue condition number being defined as

$\kappa(\omega_{n})=\kappa_{n}\equiv\frac{\|u_{n}\|_{{}_{E}}\,\|v_{n}\|_{{}_{E}}}{|\langle
u_{n},v_{n}\rangle|_{{}_{E}}}.$ (33)

When the operator $L$ is normal, i.e. when $[L^{\dagger},L]=0$, the left and
right eigenvectors can be diagonalized in the same basis. Therefore, from Eq.
(33), $\kappa_{n}=1$ for normal operators, which is indicative of spectral
stability. In these cases, generic perturbations of order $\epsilon$ to the
operator can only lead to migration of $\omega_{n}(\epsilon)$ up to a distance
$\epsilon$ in the complex plane. In contrast, when $L$ is non-selfadjoint,
then the left and right eigenvectors are not necessarily collinear and can
even become almost orthogonal, which leads to a very large condition number,
as can be seen from Eq. (33); a minuscule perturbation $\epsilon$ to $L$ can
then lead to a disproportionately large migration of $\omega_{n}(\epsilon)$ in
the complex plane. Therefore, $\kappa_{n}$ provides another measure of the
sensitivity of eigenvalues to perturbations, once again associated with our
choice of norm.

###  II.5 Numerical discretization

As in Ref. [24], the calculations presented below have been performed with a
Chebyshev spectral discretization method with $N$ grid points along the
compactified radial coordinate $\chi$. In particular, Appendix C from Ref.
[24] shows that the energy scalar product has a Chebyshev representation in
terms of the discrete vector $\ket{u}$ and the energy operator $\hat{H}$ as

$|\langle u,v\rangle|_{{}_{E}}=\bra{u}\hat{H}\ket{u}.$ (34)

Then, the discrete representation of the adjoint operator $\hat{L}^{\dagger}$
in the the energy norm follows from the discrete energy operator $\hat{H}$ via

$\hat{L}^{\dagger}=\hat{H}^{-1}\cdot\left(\hat{L}^{t}\right)^{*}\cdot\hat{H}.$
(35)

From the energy scalar product, one can also derive [24] an induced norm
$\|\hat{A}\|_{{}_{E}}$ for a given matrix $\hat{A}$ describing a discrete
representation of some physical quantity via

$\|\hat{A}\|_{{}_{E}}=\sqrt{\rho\left(A^{\dagger}A\right)},\quad\rho\left(M\right)=\max_{m\in\sigma(M)}\left\\{m\right\\},$
(36)

where the adjoint $A^{\dagger}$ is constructed through Eq. (35).

##  III Pseudospectrum of de Sitter black holes

In a recent work [39], the spectral stability of SdS and RNdS BHs was studied.
For SdS, the QNM spectrum of perturbed versions of the evolution operator $L$
was calculated, along with a topographical map of the pseudospectrum of some
particular cases, using the numerical implementation of the definition (27)
discussed above. For RNdS, only the perturbed QNM spectra were investigated,
while no pseudospectra calculations were performed. In the present work, we
aim to corroborate, reinterpret, and expand upon some of these results. We
focus on the case of a scalar field, as a representative example which covers
all behaviors of interest, as we will see in the following. We split our
presentation into two subsections covering the SdS and RNdS cases, though some
results will be common to both.

###  III.1 Schwarzschild-de Sitter black holes

####  III.1.1 The special case of $\ell=0$ scalar perturbations

We begin by reinterpreting one particular result of Ref. [39]: the apparent
spectral stability of the $\ell=0$ scalar QNMs. No migration of scalar QNMs
with $\ell=0$ was observed in that work, even when the effective potential was
perturbed with (seemingly) high-amplitude fluctuations. In what follows we
identify the particular features for the $\ell=0$ case which led to this
conclusion. We will focus on SdS geometries, though we note that the results
are independent of the BH charge and thus generalize to the broader family of
RNdS.

Linear-order scalar perturbations in dS spacetimes with angular number
$\ell=0$ have one particular feature which distinguished them from other types
of perturbations and other angular numbers. The eigenvalue problem related to
Eq. (34) in the frequency domain has a non-trivial zero mode, i.e. the problem
has a non-vanishing eigenvector $\phi_{0}(\chi)$ associated with the zero
eigenvalue $\omega=0$, that is regular at the event and cosmological horizon.
For the zero-mode function $\phi_{0}(\chi)$, Eq. (34) becomes

$\dfrac{d}{d\chi}\bigg{(}p(\chi)\,\phi_{0}^{\prime}(\chi)\bigg{)}-q_{0}(\chi)=0,$
(37)

with $p(\chi)$ given by Eq. (17) and $q_{0}(\chi)$ given by setting $\ell=0$
in Eqs. (21) and (26), i.e.

$q_{0}(\chi)=\dfrac{r_{h}^{2}\Lambda}{3\chi^{2}}\left[-2+\chi^{3}\left(1+\dfrac{r_{\Lambda}}{r_{h}}\right)\dfrac{r_{\Lambda}}{r_{h}}\right].$
(38)

A solution to Eq. (37), regular at $r=r_{h}$ and $r=r_{\Lambda}$, is

$\phi_{0}(\chi)=\dfrac{1}{\chi}.$ (39)

The zero mode, however, does not contribute to the energy of the system.
Indeed, Eq. (29) (with $\phi\sim\omega\psi=0$) gives a vanishing result,

$\|\phi_{0}\|^{2}_{{}_{E}}=\dfrac{r_{h}^{2}\Lambda}{3\chi^{3}}\left[1-\chi^{2}\left(1+\dfrac{r_{\Lambda}}{r_{h}}+\dfrac{r_{\Lambda}^{2}}{r_{h}^{2}}\right)\right]\Bigg{|}_{\chi_{h}}^{\chi_{\Lambda}}=0.$
(40)

Figure 1: Eigenvalues of the discretized representation of the energy operator
$\hat{H}$ (see Eqs. (II.4) and (34)) for scalar perturbations on a SdS
spacetime with $\Lambda M^{2}=0.01$ and $N=100$ grid points. Scalar
perturbations with $\ell=0$ (red circles) have a vanishing fundamental
eigenvalue $\uplambda_{0}\sim 10^{-65}$ (up to round-off error) associated
with the non-trivial zero mode. For $\ell>0$, i.e. $\ell=1$ (black squares),
$\ell=2$ (blue ‘x’ crosses) and $\ell=3$ (purple ‘+’ crosses), $\hat{H}$ is
positive-definite and $\uplambda_{n}>0$. Due to the scaling of the plot, all
these modes are bunched together at the top.  Figure 2: Numerical convergence
of the normalized perturbation amplitude $\epsilon$ from Eq. (42) for scalar
fields on SdS spacetime with $\Lambda M^{2}=0.01$ and frequency $k=10$ for a
deterministic perturbation added to the BH potential. For $\ell>0$, i.e.
$\ell=1$ (black squares), $\ell=2$ (blue ‘x’ crosses) and $\ell=3$ (purple ‘+’
crosses), $\|\delta L_{1}\|_{{}_{E}}\sim 1$ and the normalized amplitudes
converge to values with the same order of magnitude of the target scale (here
$10^{-1}$). For $\ell=0$ (red circles), $\epsilon$ tends to zero as $N$
increases, due to the norm $\|\delta L_{1}\|_{{}_{E}}$ being ill-defined.
These extremely small values of the normalized perturbation amplitude explain
why Ref. [39] did not observe the same spectral instabilities for $\ell=0$ as
those for $\ell>0$.

Consequently, the zero mode is also an eigenvector of the energy operator,
with vanishing eigenvalue.

We confirm this by calculating the eigenvalues of the discrete representation
for the energy operator $\hat{H}$ (see Sec. II.5). Fig. 1 displays the
eigenvalues of $\hat{H}$ for $\ell=0$ (red circles), $\ell=1$ (black squares),
$\ell=2$ (blue ‘x’ crosses) and $\ell=3$ (purple ‘+’ crosses). The cases with
$\ell>0$ have real, positive eigenvalues $\uplambda_{n}$, as expected from a
positive-definite operator. For $\ell=0$, we encounter $\uplambda_{0}=0$ up to
numerical round-off error, which in our case is $\uplambda_{0}\sim 10^{-65}$
for $N=100$ grid points.

Formally, the existence of the zero eigenvalue prevents the energy operator
$\hat{H}$ from being invertible, i.e. $\hat{H}$ is singular. In principle,
this property prevents the calculation of the adjoint operator $L^{\dagger}$
via Eq. (35). In practice, numerical techniques compute finite large numbers
for $\hat{H}^{-1}$ and $L^{\dagger}$ due to grid discretization with a finite
$N$, in contrast to the exact analytical results. However, conclusions
extracted from manipulating these discrete operators must be dealt with
extreme care.

Let us then consider the results from Ref. [39] on the spectral stability of
scalar $\ell=0$ QNMs under perturbations of the SdS potentials. A direct
assessment of the stability follows from calculating the QNMs of the operator
with a small perturbation, as indicated in Eq. (28). The particular type of
perturbation $\delta L$, in accordance with the setup of Ref. [24], is one
which only affect the effective potential $q(\chi)$ in Eq. (25),

$\delta L=\dfrac{\delta L_{1}}{\|\delta L_{1}\|_{{}_{E}}},\quad\delta
L_{1}=\left(\begin{array}[]{cc}0&0\\\ -\delta
q(\chi)/w(\chi)&0\end{array}\right).$ (41)

Note that the factor $\|\delta L_{1}\|_{{}_{E}}$ in the above expression
ensures that $\|\delta L\|_{{}_{E}}=1$ as required by Eq. (28). Effectively,
one can redefine Eq. (28) in the form

$L\rightarrow L+\epsilon\,\delta
L_{1},\quad\epsilon\rightarrow\dfrac{\epsilon}{\|\delta L_{1}\|_{{}_{E}}},$
(42)

where $\epsilon$ is redefined such that it represent an “energy-normalized”
perturbation scale222From this point on, $\epsilon$ will refer to the
normalized perturbation scale if not stated otherwise.. For $\ell=0$, the
normalized perturbation amplitude $\epsilon$ in Eq. (42) is ill-defined, since
the quantity $\|\delta L_{1}\|_{{}_{E}}$ relies on the existence of the
inverse operator $\hat{H}^{-1}$.

As a representative example, Fig. 2 shows the normalized perturbation scale
$\epsilon$ from Eq. (42) for the cases $\ell=0$ (red circles), $\ell=1$ (black
squares), $\ell=2$ (blue ‘x’ crosses) and $\ell=3$ (purple ‘+’ crosses) as a
function of the numerical resolution $N$. Following Ref. [39], we consider a
deterministic perturbation

$\delta q(\chi)=\cos(2\pi k\chi),$ (43)

with targeted amplitude $10^{-1}$. For the wavenumber, we set $k=10$, as
opposed to $k=60$ employed in Ref. [39]. The reason is that an accurate
discrete representation of $\delta q(\chi)$ for $k=60$ requires a rather high
grid resolution ($N\sim 10^{3}$), and the relevant conclusions can be drawn
with considerably less numerical expenditure. For $\ell>0$, one finds
approximately $\|\delta L_{1}\|_{{}_{E}}\sim 1$, and the normalized amplitude
$\epsilon$ already has the same order of magnitude as the targeted value.
Moreover, its numerical value converges to a fixed value as the numerical
resolution $N$ increases. The $\ell=0$ case, on the other hand, tends to zero
when we increase the numerical resolution, since $\|\delta L_{1}\|_{{}_{E}}$
is ill-defined. i.e. $\|\delta L_{1}\|_{{}_{E}}\rightarrow\infty$ as
$N\rightarrow\infty$.

Hence, one may be wrongly lead to the conclusions that the $\ell=0$ QNMs
possess stronger spectral stability properties than the cases where $\ell>0$,
if one uses this nearly-zero normalized amplitude $\epsilon$ and sees no
significant displacement in the lower part of the QNM spectrum (the
fundamental mode and first overtones). As we will show in Sec. III.2, the
$\ell=0$ spectra do suffer from the same instabilities as $\ell>0$ if one
promotes the normalized amplitude $\epsilon$ as a free parameter (not
normalized by the energy norm) to quantify the perturbation’s scale. This
approach is appropriate due to the fact that the perturbation’s targeted scale
and its normalized scale $\epsilon$ have typically the same order of magnitude
whenever $\|\delta L_{1}\|_{{}_{E}}$ is well-defined, as we discussed above.
We assume the validity of this property in a limiting approach
$\ell\rightarrow 0$, and work directly with $\epsilon$ as a free (non-
normalized) perturbation parameter when $\ell=0$.

Figure 3: Top panel: Scalar $\ell=0$ normal modes (red dots) and pseudospectra
(white contour curves) of a SdS BH with $\Lambda M^{2}=0.01$. The operator $L$
in this test case is self-adjoint, since $L_{2}=0$. The
$\epsilon$-pseudospectra contours range from $-1.8$ (innermost contour) to
$-0.4$ (outermost contour) with steps of $0.2$. On top of figure, the
condition number is shown for various modes. The question mark for the non-
trivial zero mode corresponds to a divergent condition number as the grid
points are increased. Bottom panel: Same as top panel for scalar $\ell=1$
normal modes. For both figures, $N=150$ grid points were used.

Figure 4: Left: Scalar $\ell=1$ QNMs (red dots) and pseudospectra (white
contour curves) of a SdS BH with $\Lambda M^{2}=0.01$. The
$\epsilon$-pseudospectra contours range from $-27.5$ (uppermost contour) to
$-2.5$ (bottom contour) with steps of $2.5$ and $N=150$. Right: Same as left
with $\Lambda M^{2}=0.001$. The $\epsilon$-pseudospectra contours range from
$-30$ (uppermost contour) to $-2.5$ (bottom contour) with steps of $2.5$ and
$N=250$.

####  III.1.2 $L_{2}=0$: The self-adjoint case

We have analytically and numerically established that calculations involving
the energy norm associated with the scalar $\ell=0$ QNM problem of BHs in dS
geometries are not to be taken at face value. However, the numerical
discretization of the operators we use leads to an effective cutoff to the
divergences involved. And while some quantities are manifestly unphysical even
with this cutoff, such as the normalized perturbation amplitudes discussed
above, others are seemingly better behaved. Particularly, a calculation of the
pseudospectrum at a finite numerical resolution seems to yield a result which
is similar to those for other values of $\ell$, both qualitatively and
quantitatively.

Motivated by Ref. [24], we perform a test to verify whether the numerically
evaluated pseudospectrum provides a qualitatively correct picture by first
performing a calculation for which we set $L_{2}=0$. This choice cancels the
dissipative boundary conditions at the event and cosmological horizon, thus
rendering the operator $L$ selfadjoint. The system becomes conservative and
its eigenvalues become normal modes. Figure 3 shows two pseudospectra for a
conservative SdS system for which we have set $L_{2}=0$. In the bottom panel,
we observe a case of a scalar field with $\ell=1$. The pseudospectrum is well-
behaved, as expected given the discussion in the previous sections. The
boundaries of the $\epsilon$-pseudospectral regions form concentric circles of
radius $\epsilon$ around the modes, indicating spectral stability. The
condition numbers of the modes tell the same story; since $L$ is selfadjoint,
the left and right eigenvectors of each eigenvalue coincide, leading to
$\kappa_{n}=1$.

In the top panel of Fig. 3, we show the pseudospectrum of the $L_{2}=0$
construction corresponding to an $\ell=0$ scalar field. The picture is similar
to the $\ell=1$ case, though with one particular oddity. Although the
concentric circles around modes $\omega\neq 0$ have radii of approximately
$\epsilon$, the circles around the zero mode are about twice as large,
modulating the function as they spread out. This indicates that, while the
operator is still manifestly selfadjoint, the divergence of the resolvent
around the zero mode has a greater “width”, suggesting that it captures not
only the standard singularity present at an eigenvalue, but also a divergence
of another nature, namely, the one due to the norm. What is perhaps surprising
is that at small enough values of $\epsilon$ the effect seems to be contained
to a vicinity of the zero mode itself, suggesting that the pseudospectrum
still provides a qualitatively valid description elsewhere. We will make use
of this fact in interpreting the $\ell=0$ pseudospectra we compute for the
RNdS case in the next subsection. We leave this analysis for the RNdS case for
the sake of brevity, and generality.

We also note that the special status of the $\omega=0$ mode is also observed
when we calculate the condition number via Eq. (33). Our results give
$\kappa_{n}=1$ for all modes with $\omega\neq 0$, as expected (see inset in
the top panel of Fig. 3). However, Eq. (33) is not directly applicable for
$\omega=0$, because the zero mode has a vanishing energy norm, i.e.,
$\|u_{n}\|_{{}_{E}}=\|v_{n}\|_{{}_{E}}=|\langle
u_{n},v_{n}\rangle|_{{}_{E}}=0$. Numerically, the ill-defined condition number
manifests itself with values of $\kappa_{0}$ that tend to a divergence as the
resolution $N$ increases (designated with a red question mark on the top panel
of Fig. 3). While formally it may be possible that $\kappa_{0}\rightarrow 1$
if one considers the limit $\ell\rightarrow 0$ in a suitable manner,
developing this hypotheses goes beyond the scope of the present work.

####  III.1.3 Schwarzschild-de Sitter pseudospectra

For completeness, in Fig. 4 we present two plots of SdS scalar-field
pseudospectra, corresponding to two representative examples in the parameter
space. Particularly, we set the angular number $\ell=1$, and we use the values
of the cosmological constant $\Lambda M^{2}=10^{-2}$ and $\Lambda
M^{2}=10^{-3}$. We observe the expected picture of the spectral instability of
QNMs [24, 39], with contour branches that open up in the complex plane.

We note that besides the PS QNMs, purely-imaginary modes appear as well. These
are not the numerically non-convergent “branch-cut” modes, such as those found
in asymptotically-flat compact objects [24, 38, 43], but rather convergent
eigenvalues of the QNM spectrum that belong to the dS family. In both cases
depicted in Fig. 4, the dS family is the dominant one, i.e. the least damped
mode is purely imaginary. Particularly, for the $\Lambda M^{2}=10^{-2}$ case,
the dominant dS mode is followed by a PS oscillatory mode, followed by two
more dS modes, and so on and so forth. For the case $\Lambda M^{2}=10^{-3}$,
the fundamental mode and first two overtones belong to the dS family, with the
next overtone belonging to the PS family, after which nine dS overtones take
over before a PS mode reappears (that lies outside the range of the complex
plane that we have plotted).

Figure 5: Left: Scalar $\ell=2$ QNMs (red dots) and pseudospectra (white
contour curves) of a RNdS BH with $\Lambda M^{2}=0.06$ and $Q/M=0.5$. The
$\epsilon$-pseudospectra contours range from $-28$ (uppermost contour) to $0$
(bottom contour) with steps of $2$ and $N=150$. Right: Same as left for scalar
$\ell=1$ QNMs of a RNdS BH with $\Lambda M^{2}=0.001$ and $Q/M=0.3$. The
$\epsilon$-pseudospectra contours range from $-34$ (uppermost contour) to $0$
(bottom contour) with steps of $2$ and $N=250$.

Figure 6: Left: Scalar $\ell=0$ QNMs (red dots) and pseudospectra (white
contour curves) of a RNdS BH with $\Lambda M^{2}=0.01$ and $Q/M=0.9$. The
$\epsilon$-pseudospectra contours range from $-34$ (uppermost contour) to $0$
(bottom contour) with steps of $2$ and $N=150$. Right: Same as left panel with
$\Lambda M^{2}=0.03$ and $Q/Q_{\rm{max}}=0.999$. The $\epsilon$-pseudospectra
contours range from $-35$ (uppermost contour) to $-2.5$ (bottom contour) with
steps of $2.5$ and $N=200$.

An important aspect of the pseudospectrum, both here and in the RNdS
calculations presented below, is that the contour lines generally have their
smallest imaginary part on the imaginary axis, from where they go up and
presumably approach the logarithmic behavior described in [38], but they are
closer together in that region. This implies that the dS modes (as well as the
NE ones in RNdS) are affected by lower-energy perturbations, but that their
relative displacement may be less than that of PS modes when the latter are
also destabilized. This may contribute to the likelihood of an “overtaking”
instability of the fundamental mode after a perturbation, as we will see in
the next section for particular examples in the parameter space.

###  III.2 Reissner-Nordström-de Sitter black holes

####  III.2.1 Reissner-Nordström-de Sitter pseudospectra

We now move on to the RNdS case, presenting results involving the
pseudospectrum and spectrum of a perturbed evolution operator, for a scalar
field. As discussed above, as well as in [48], RNdS geometries possess three
families of QNMs: the complex PS modes, the purely imaginary dS, and purely
imaginary NE modes. Since the NE modes are in the vicinity of dominant
(longest-lived) modes only for RNdS geometries with charge which is quite
close to the extremal one $Q_{\text{max}}$ (for which the inner and outer BH
horizons almost overlap), and are particularly long-lived for the complicated
case of $\ell=0$, it is interesting to test how the pseudospectra behave in
this previously uncharted case.

Figure 5 shows the pseudospectrum for two representative cases with different
angular numbers $\ell\geq 1$, cosmological constants, and BH charges. In both
cases the charge $Q$ is less than $50\%$ of the extremal, and $\ell\neq 0$,
making the NE family of modes subdominant (in fact they lie outside the
plotted region). On the one hand, we observe the effect of the cosmological
constant and charge on the spectrum, which has been thoroughly analyzed in
[48, 53, 49]. Note that in the first case (on the left), the dominant mode
belongs to the PS family, while the three purely imaginary modes that appear
belong to the dS family, though on the right, the picture is similar to that
of SdS BHs with $\Lambda M^{2}=10^{-3}$, where a plethora of dS modes dominate
the spectrum in the part of the complex plane shown. On the other hand, we
observe the standard QNM spectral instability, and we note that the minima of
the contour lines once again generally lie on the imaginary axis, except for
the lines in the vicinity of the fundamental mode and first overtone in the
left panel, suggesting that this is a general feature for PS modes which are
below any purely imaginary ones.

Figure 7: Scalar unperturbed $\ell=0$ QNMs (red dots) and perturbed $\ell=0$
QNMs (blue boxes) of a RNdS BH with $\Lambda M^{2}=0.01$ and
$Q/Q_{\rm{max}}=0.999$. The perturbed QNMs result from a random-vector
perturbation added to the effective potential with non-normalized perturbation
scale $\epsilon=10^{-3}$. The number of grid points used for both spectra is
$N=200$.

In Fig. 6 we show the pseudospectra of two RNdS BHs for $\ell=0$, for which
the spectrum contains the non-trivial zero mode. As discussed in previous
sections, this zero mode has an associated eigenfunction with a vanishing
energy norm, which makes the calculation of operator norms generally ill-
defined. However, as discussed above, qualitatively reasonable results can be
obtained for the pseudospectrum for a given numerical resolution $N$, which
acts as a cutoff to the divergences involved. Particularly, the results from
Sec. III.1.2 provide strong indication that the numerical cutoff leads to
meaningful results for the pseudospectrum in the vicinity of the overtones.
Since we only use the pseudospectrum as a qualitative measure of the QNM
instability, as discussed in [42], this result is sufficient for the present
analysis.

Particularly, the left panel of Fig. 6 depicts the pseudospectrum associated
with the results shown in Fig. 13 of Ref. [39] ($Q/M=0.9$ and $\Lambda
M^{2}=10^{-2}$), for which a stronger spectral stability was suspected.
However, the qualitative aspects of the pseudospectrum suggest the opposite,
which is in line with our discussion in Sec. III.1.1. The opening of branches
indicate a spectral instability of both complex PS and purely imaginary dS
modes. As a side note, it is worth observing that in this left panel, the
purely imaginary modes (which here belong to the dS family) seem almost
equidistant as the overtone number increases, if one excludes the zero mode,
in contrast with the figure in the right panel.

The right panel of Fig. 6 depicts a case of a near-extremally-charged RNdS BH
($Q/Q_{\text{max}}=0.999$ and $\Lambda\,M^{2}=0.03$). In this case the part of
the spectrum depicted contains all three families of QNMs: oscillatory PS
modes, dS modes, and NE modes. Notice the interplay occurring in the imaginary
axis, where now the equidistant dS modes are scrambled with NE modes, though
both families can be clearly tracked with the use of the approximations
provided in [48]. For the specified BH parameters, the NE family possesses the
least-damped fundamental mode that dominates the late-time tail until the zero
mode takes over to relax the ringdown signal to a non-zero constant, as
discussed in [65]. Again, the pseudosepctra contour lines indicate a spectral
instability of all families of QNMs, as in all other cases studied here.

The presence of spectral instability in all cases studied in this work
reflects an instability to physical perturbations. While the pseudospectrum
captures the possible migration of modes under any type of perturbation to the
evolution operator $L$, perturbing only the effective potential is sufficient
to trigger the instability. In Fig. 7 we show explicitly that this also
applies to the $\ell=0$ case. We plot an unperturbed QNM spectrum, along with
the spectrum after a perturbation to the effective potential in the form of
the addition of a random vector with components of magnitude $10^{-3}$ to its
discretized version. For other values of $\ell$, such a vector perturbation
has a corresponding energy norm of the same order, $10^{-3}$, suggesting that
this is likely a very small perturbation in any physically reasonable, well-
defined norm. In other words, we treat $\epsilon$ as a free parameter,
extrapolating the likely order of magnitude of the perturbation from other
known cases. We observe that for such a perturbation, overtones are indeed
destabilized, independently of the family they belong to.

###  III.3 Fundamental mode (in)stability

In this final subsection, we revisit another open issue put forward in Ref.
[39] regarding the spectral stability of the fundamental QNM with respect to
(deterministic) perturbations to the effective potential. In what follows, we
will present quantitative results which show that the mode which is initially
the fundamental one always appears to be spectrally stable, i.e. its relative
displacement in the complex plane for perturbations of up to order
$\epsilon=10^{-1}$ (in the energy norm) is less than $\epsilon$, in agreement
with [39]. However, we note that in certain regions of the parameter space,
for which the imaginary part of a PS mode frequency is very close to one of
the purely imaginary modes, and these two are the least damped modes, then an
“overtaking” instability may occur, whereby the mode which is initially the
fundamental one becomes an overtone after the perturbation. Particularly, we
observe this for fine-tuned cases of the parameter space in which the
fundamental mode is initially of the dS family, and gets overtaken by the
first overtone, which is of the PS family.

More precisely, following Refs. [32, 33], we will use
$\Delta\omega^{(\epsilon)}\equiv|\omega^{(0)}_{0}-\omega_{0}^{(\epsilon)}|/|\omega^{(0)}_{0}|$
to quantify the relative migration distance between the unperturbed
fundamental mode $\omega_{0}^{(0)}$ and the perturbed fundamental mode
$\omega_{0}^{(\epsilon)}$. In this definition, the fundamental mode is
understood as the mode with the smallest imaginary part. Ref. [32] identifies
two types of spectral instability for the fundamental mode: (i) a
destabilization via a continuous migration of the initial fundamental mode
itself for which $\Delta\omega^{(\epsilon)}\gg\epsilon$; (ii) a discontinuous
overtaking of the fundamental mode, for which
$\Delta\omega^{(\epsilon)}\gg\epsilon$ due to a migrated overtone becoming the
new fundamental mode. This latter case is the one observed here.

We study deterministic perturbations to the potential of the form (43) and
focus on perturbations with wavenumber $k=20$, which is at the high end of
what is smoothly resolved with our number of discretization points (we use
wavenumber perturbations because these are the ones which destabilize the QNM
spectrum the most, as shown in [24]). We use a normalized perturbation scale
ranging in $\epsilon=10^{-6}-10^{-1}$. For perturbations of a different
variety to the same class of BHs considered here, we refer the reader to Ref.
[33].

Figure 8: Relative value of distance in the complex plane
$\Delta\omega^{(\epsilon)}$ between the unperturbed fundamental QNM and the
perturbed fundamental mode versus the scale $\epsilon$ of the deterministic
perturbation-vector with wavenumber $k=20$ that is added to the effective
potential of various RNdS BHs. The cases depicted correspond to $\ell=1$,
$Q/M=0.5$, and $\Lambda M^{2}=3/121\sim 0.025$ (red squares, overtaking
instability from dS to PS family), $\ell=2$, $Q/M=0.5$ and $\Lambda
M^{2}=0.06$ (green circles, PS dominant family) and, finally, $\ell=2$,
$Q/M=0.5$ and $\Lambda M^{2}=0.006$ (blue triangles, dS dominant family). The
number of grid points used for all cases is $N=150$ which guarantee an
accurate resolution for the deterministic perturbation-vector.

Figure 8 quantifies the migration of the fundamental QNM
$\Delta\omega^{(\epsilon)}$ caused by this effective potential perturbation.
We depict three RNdS cases with fixed charge $Q/M=0.5$ and varying $\Lambda
M^{2}$, as well as $\ell$. We observe that when the BH parameters are such
that the fundamental mode is an oscillatory PS one (green circles in Fig. 8)
or a purely imaginary dS one (blue triangles in Fig. 8), and the first
overtone has a sufficiently distinct imaginary part (as occurs in both of the
aforementioned), then the fundamental mode after the perturbation remains in
the same family. In these cases, we always encounter
$\Delta\omega^{(\epsilon)}\lesssim\epsilon$, thus the fundamental mode remains
spectrally stable.

The third case we analyze is a fine-tuned one for which the original
fundamental mode belongs to the dS family, and the first overtone is of the PS
family, with an imaginary part which is only $\leavevmode\nobreak\ 10^{-3}$
larger (in units of $M$). The parameters which lead to this configuration are
$\ell=1$, $Q/M=0.5$, and $\Lambda M^{2}=3/121\sim 0.025$. We observe that the
fundamental mode is stable until a perturbation of $\epsilon=10^{-1}$, for
which the overtaking occurs (red boxes in Fig. 8). Since the new fundamental
mode is at a distance which is over one order of magnitude larger than the
perturbation scale of $10^{-1}$, an overtaking instability can be said to take
place. The cause for this effect, in this case, seems to be related to the
fact that the PS family has the tendency for larger migration when perturbed,
as shown in [33], and discussed above in regards to the shape of the
pseudospectrum contours.

##  IV Conclusions

We have systematically analyzed several key features of the QNM instability of
(neutral) scalar perturbations on the SdS and RNdS spacetimes, which have
proven to be an important testing ground of the phenomenon due to their richer
spectrum as compared to asymptotically flat BHs [48, 49]. Following Ref. [24],
we have expressed the QNM problem in the hyperboloidal framework, which
imposes the appropriate boundary conditions geometrically, and we have studied
the behavior of the spectrum under perturbations, and calculated the
associated pseudospectrum, in several representative examples. Our analysis
has mainly been focused on three points: (i) revisiting the case of the
peculiar spectral stability of QNM overtones for $\ell=0$ scalar perturbations
found in Ref. [39], (ii) the (in)stability of the fundamental mode for other
angular numbers $\ell\geq 1$, also discussed in [39], and (iii) the
calculation of the pseudospectrum of RNdS scalar perturbations, which had not
been explored in the literature yet.

Regarding the stability of overtones for $\ell=0$, we have found that it is
simply an artefact of an ill-defined norm. Particularly, the $\ell=0$ scalar
QNMs in dS BHs possess a non-trivial zero mode, which turns out to have a
vanishing energy norm. This implies that the norm operator acting on the space
of solutions has a zero eigenvalue, associated with this mode, and it is thus
not invertible. Calculating norms numerically with this operator thus leads to
rapidly diverging results, and attempting to normalize perturbations with such
a norm leads to practically vanishing results. It is thus likely that such an
ill-defined normalization has lead to a nearly vanishing perturbation to be
dubbed “large”, and the overtones within the range of study to be mistaken for
stable. In our analysis, by treating the perturbation scale as a non-
normalized free parameter, we have shown that the $\ell=0$ QNMs are likely to
be spectrally unstable, just like in all other cases studied in the literature
thus far.

We have further explored the energy norm for $\ell=0$, finding that the cutoff
induced by the spatial discretization can lead to qualitatively reasonable
results in the calculation of some quantities, the pseudospectrum in the
vicinity of the overtones being one of them. From it, we also obtain a
qualitative picture of spectral instability, both in SdS and RNdS BHs.

Furthermore, we have examined the pseudospectrum of RNdS for cases with
$\ell\geq 1$, for which the energy norm calculations are formally well-
defined, finding qualitatively the same picture of spectral instability in all
of them. Specifically, we have looked at representative examples in both
weakly-charged and near extremal cases, for which different families of modes
are dominant.

Finally, the stability of the fundamental mode under high-wavenumber
deterministic perturbations to the effective potential (akin to those used in
[24]) has been tested, particularly for cases $\ell\geq 1$, for which the
magnitude of the perturbations can be measured well. We find that when the PS
or dS family is sufficiently dominant with respect to the other families, the
fundamental mode is stable. However, when the fundamental mode is of one
family and the first overtone of another, and they have sufficiently similar
imaginary parts, the migrated overtone can take over as the new fundamental
mode after a sufficiently large perturbation. We have particularly observed
this for a case in which a PS overtone overtakes a dS purely imaginary
fundamental mode, and have dubbed this an “overtaking” instability. We note,
however, that the individual migration of each of the two modes in this case
was within the bounds of spectral stability.

Ultimately, we have taken another step towards confirming the universality of
QNM spectral instabilities of BHs and compact objects in spherical symmetry.
These may become relevant for the BH spectroscopy program when the signal-to-
noise ratio of current and future ground- and space-based detectors is
increased sufficiently, giving us detections of ringdown signals which go
beyond the dominant early-time emission [27, 28].

###### Acknowledgements.

The authors would like to warmly thank Vitor Cardoso and José-Luis Jaramillo
for stimulating discussions and comments. K.D. acknowledges financial support
provided under the European Union’s H2020 ERC, Starting Grant agreement no.
DarkGRA–757480 and the MIUR PRIN and FARE programmes (GW-NEXT,
CUP:B84I20000100001). K.D. also acknowledges hospitality and financial support
from the Center of Astrophysics and Gravitation (CENTRA) where part of this
work was conducted. R.P.M acknowledges support from the Villum Investigator
program supported by the VILLUM Foundation (grant no. VIL37766) and the DNRF
Chair program (grant no. DNRF162) by the Danish National Research Foundation.
This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Sklodowska-Curie grant
agreement No 101131233 V.B. acknowledges support form the Spanish government
projects PID2020-118159GB-C43 and PID-2020-118159GB-C44. This project has
received funding from the European Union’s Horizon 2020 research and
innovation programme under the Marie Sklodowska-Curie grant agreement No
101007855. We thank FCT for financial support through Project No.
UIDB/00099/2020. We acknowledge financial support provided by FCT/Portugal
through grants PTDC/MAT-APL/30043/2017 and PTDC/FIS-AST/7002/2020. This
research project was conducted using the computational resources of “Baltasar
Sete-Sois” cluster at Instituto Superior Técnico.

## References

  * Abbott  _et al._ [2021] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run, (2021), arXiv:2111.03606 [gr-qc] .
  * Kokkotas and Schmidt [1999] K. D. Kokkotas and B. G. Schmidt, Quasinormal modes of stars and black holes, Living Rev. Rel. 2, 2 (1999), arXiv:gr-qc/9909058 .
  * Nollert [1999] H.-P. Nollert, Quasinormal modes: the characteristic ‘sound’ of black holes and neutron stars, Class. Quant. Grav. 16, R159 (1999).
  * Berti  _et al._ [2009] E. Berti, V. Cardoso, and A. O. Starinets, Quasinormal modes of black holes and black branes, Class. Quant. Grav. 26, 163001 (2009), arXiv:0905.2975 [gr-qc] .
  * Konoplya and Zhidenko [2011] R. A. Konoplya and A. Zhidenko, Quasinormal modes of black holes: From astrophysics to string theory, Rev. Mod. Phys. 83, 793 (2011), arXiv:1102.4014 [gr-qc] .
  * Leaver [1986] E. W. Leaver, Spectral decomposition of the perturbation response of the Schwarzschild geometry, Phys. Rev. D 34, 384 (1986).
  * Gundlach  _et al._ [1994] C. Gundlach, R. H. Price, and J. Pullin, Late time behavior of stellar collapse and explosions: 1. Linearized perturbations, Phys. Rev. D 49, 883 (1994), arXiv:gr-qc/9307009 .
  * Dafermos [2005] M. Dafermos, The Interior of charged black holes and the problem of uniqueness in general relativity, Commun. Pure Appl. Math. 58, 0445 (2005), arXiv:gr-qc/0307013 .
  * Dafermos and Rodnianski [2005] M. Dafermos and I. Rodnianski, A Proof of Price’s law for the collapse of a selfgravitating scalar field, Invent. Math. 162, 381 (2005), arXiv:gr-qc/0309115 .
  * Hintz [2022] P. Hintz, A Sharp Version of Price’s Law for Wave Decay on Asymptotically Flat Spacetimes, Commun. Math. Phys. 389, 491 (2022), arXiv:2004.01664 [math.AP] .
  * Dafermos  _et al._ [2021] M. Dafermos, G. Holzegel, I. Rodnianski, and M. Taylor, The non-linear stability of the Schwarzschild family of black holes, (2021), arXiv:2104.08222 [gr-qc] .
  * Carroll [2001] S. M. Carroll, The Cosmological constant, Living Rev. Rel. 4, 1 (2001), arXiv:astro-ph/0004075 .
  * Riess  _et al._ [1998] A. G. Riess _et al._ (Supernova Search Team), Observational evidence from supernovae for an accelerating universe and a cosmological constant, Astron. J. 116, 1009 (1998), arXiv:astro-ph/9805201 .
  * Perlmutter  _et al._ [1999] S. Perlmutter _et al._ (Supernova Cosmology Project), Measurements of $\Omega$ and $\Lambda$ from 42 high redshift supernovae, Astrophys. J. 517, 565 (1999), arXiv:astro-ph/9812133 .
  * Freedman  _et al._ [2001] W. L. Freedman _et al._ (HST), Final results from the Hubble Space Telescope key project to measure the Hubble constant, Astrophys. J. 553, 47 (2001), arXiv:astro-ph/0012376 .
  * Brady  _et al._ [1997] P. R. Brady, C. M. Chambers, W. Krivan, and P. Laguna, Telling tails in the presence of a cosmological constant, Phys. Rev. D 55, 7538 (1997), arXiv:gr-qc/9611056 .
  * Molina  _et al._ [2004] C. Molina, D. Giugno, E. Abdalla, and A. Saa, Field propagation in de Sitter black holes, Phys. Rev. D 69, 104013 (2004), arXiv:gr-qc/0309079 .
  * Dafermos and Rodnianski [2007] M. Dafermos and I. Rodnianski, The Wave equation on Schwarzschild-de Sitter spacetimes, (2007), arXiv:0709.2766 [gr-qc] .
  * Hintz and Vasy [2016] P. Hintz and A. Vasy, The global non-linear stability of the Kerr-de Sitter family of black holes 10.4310/acta.2018.v220.n1.a1 (2016), arXiv:1606.04014 [math.DG] .
  * Hintz [2016] P. Hintz, Non-linear stability of the Kerr-Newman-de Sitter family of charged black holes 10.1007/s40818-018-0047-y (2016), arXiv:1612.04489 [math.AP] .
  * Aguirregabiria and Vishveshwara [1996] J. Aguirregabiria and C. Vishveshwara, Scattering by black holes: a simulated potential approach, Physics Letters A 210, 251 (1996).
  * Nollert [1996] H.-P. Nollert, About the significance of quasinormal modes of black holes, Phys. Rev. D 53, 4397 (1996), arXiv:gr-qc/9602032 .
  * Nollert and Price [1999] H.-P. Nollert and R. H. Price, Quantifying excitations of quasinormal mode systems, J. Math. Phys. 40, 980 (1999), arXiv:gr-qc/9810074 .
  * Jaramillo  _et al._ [2021] J. L. Jaramillo, R. Panosso Macedo, and L. Al Sheikh, Pseudospectrum and Black Hole Quasinormal Mode Instability, Phys. Rev. X 11, 031003 (2021), arXiv:2004.06434 [gr-qc] .
  * Trefethen and Embree [2005] L. Trefethen and M. Embree,  _Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators_ (2005).
  * Gasperin and Jaramillo [2022] E. Gasperin and J. L. Jaramillo, Energy scales and black hole pseudospectra: the structural role of the scalar product, Class. Quant. Grav. 39, 115010 (2022), arXiv:2107.12865 [gr-qc] .
  * Jaramillo  _et al._ [2022] J. L. Jaramillo, R. Panosso Macedo, and L. A. Sheikh, Gravitational Wave Signatures of Black Hole Quasinormal Mode Instability, Phys. Rev. Lett. 128, 211102 (2022), arXiv:2105.03451 [gr-qc] .
  * Berti  _et al._ [2022] E. Berti, V. Cardoso, M. H.-Y. Cheung, F. Di Filippo, F. Duque, P. Martens, and S. Mukohyama, Stability of the fundamental quasinormal mode in time-domain observations against small perturbations, Phys. Rev. D 106, 084011 (2022), arXiv:2205.08547 [gr-qc] .
  * Giesler  _et al._ [2019] M. Giesler, M. Isi, M. A. Scheel, and S. Teukolsky, Black Hole Ringdown: The Importance of Overtones, Phys. Rev. X 9, 041060 (2019), arXiv:1903.08284 [gr-qc] .
  * Franchini and Völkel [2023] N. Franchini and S. H. Völkel, Testing General Relativity with Black Hole Quasi-Normal Modes, (2023), arXiv:2305.01696 [gr-qc] .
  * Destounis and Duque [2023] K. Destounis and F. Duque, Black-hole spectroscopy: quasinormal modes, ringdown stability and the pseudospectrum (2023) arXiv:2308.16227 [gr-qc] .
  * Cheung  _et al._ [2022] 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, 111103 (2022), arXiv:2111.05415 [gr-qc] .
  * Courty  _et al._ [2023] A. Courty, K. Destounis, and P. Pani, Quasinormal mode (in)stability and strong cosmic censorship, (2023), arXiv:2307.11155 [gr-qc] .
  * Barausse  _et al._ [2014] E. Barausse, V. Cardoso, and P. Pani, Can environmental effects spoil precision gravitational-wave astrophysics?, Phys. Rev. D 89, 104059 (2014), arXiv:1404.7149 [gr-qc] .
  * Cardoso  _et al._ [2022a] V. Cardoso, K. Destounis, F. Duque, R. P. Macedo, and A. Maselli, Black holes in galaxies: Environmental impact on gravitational-wave generation and propagation, Phys. Rev. D 105, L061501 (2022a), arXiv:2109.00005 [gr-qc] .
  * Cardoso  _et al._ [2022b] V. Cardoso, K. Destounis, F. Duque, R. Panosso Macedo, and A. Maselli, Gravitational Waves from Extreme-Mass-Ratio Systems in Astrophysical Environments, Phys. Rev. Lett. 129, 241103 (2022b), arXiv:2210.01133 [gr-qc] .
  * Destounis  _et al._ [2023] K. Destounis, A. Kulathingal, K. D. Kokkotas, and G. O. Papadopoulos, Gravitational-wave imprints of compact and galactic-scale environments in extreme-mass-ratio binaries, Phys. Rev. D 107, 084027 (2023), arXiv:2210.09357 [gr-qc] .
  * Destounis  _et al._ [2021] 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, 084091 (2021), arXiv:2107.09673 [gr-qc] .
  * Sarkar  _et al._ [2023] 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, 104002 (2023), arXiv:2304.06829 [gr-qc] .
  * Arean  _et al._ [2023] D. Arean, D. Garcia-Fariña, and K. Landsteiner, Pseudospectra of Holographic Quasinormal Modes, (2023), arXiv:2307.08751 [hep-th] .
  * Cownden  _et al._ [2023] B. Cownden, C. Pantelidou, and M. Zilhão, The pseudospectra of black holes in AdS, (2023), arXiv:2312.08352 [gr-qc] .
  * Boyanov  _et al._ [2023a] V. Boyanov, V. Cardoso, K. Destounis, J. L. Jaramillo, and R. Panosso Macedo, Structural aspects of the anti de-Sitter black hole pseudospectrum, (2023a), to appear.
  * Boyanov  _et al._ [2023b] 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, 064012 (2023b), arXiv:2209.12950 [gr-qc] .
  * Cardoso  _et al._ [2009] V. Cardoso, A. S. Miranda, E. Berti, H. Witek, and V. T. Zanchin, Geodesic stability, Lyapunov exponents and quasinormal modes, Phys. Rev. D 79, 064016 (2009), arXiv:0812.1806 [hep-th] .
  * Konoplya and Stuchlík [2017] R. A. Konoplya and Z. Stuchlík, Are eikonal quasinormal modes linked to the unstable circular null geodesics?, Phys. Lett. B 771, 597 (2017), arXiv:1705.05928 [gr-qc] .
  * Lopez-Ortega [2006] A. Lopez-Ortega, Quasinormal modes of D-dimensional de Sitter spacetime, Gen. Rel. Grav. 38, 1565 (2006), arXiv:gr-qc/0605027 .
  * Jansen [2017] A. Jansen, Overdamped modes in Schwarzschild-de Sitter and a Mathematica package for the numerical computation of quasinormal modes, Eur. Phys. J. Plus 132, 546 (2017), arXiv:1709.09178 [gr-qc] .
  * Cardoso  _et al._ [2018a] V. Cardoso, J. L. Costa, K. Destounis, P. Hintz, and A. Jansen, Quasinormal modes and Strong Cosmic Censorship, Phys. Rev. Lett. 120, 031103 (2018a), arXiv:1711.10502 [gr-qc] .
  * Destounis [2019a] K. Destounis, Superradiant instability of charged scalar fields in higher-dimensional Reissner-Nordström-de Sitter black holes, Phys. Rev. D 100, 044054 (2019a), arXiv:1908.06117 [gr-qc] .
  * Konoplya and Zhidenko [2022] R. A. Konoplya and A. Zhidenko, Nonoscillatory gravitational quasinormal modes and telling tails for Schwarzschild–de Sitter black holes, Phys. Rev. D 106, 124004 (2022), arXiv:2209.12058 [gr-qc] .
  * Richartz and Giugno [2014] M. Richartz and D. Giugno, Quasinormal modes of charged fields around a Reissner-Nordström black hole, Phys. Rev. D 90, 124011 (2014), arXiv:1409.7440 [gr-qc] .
  * Hod [2017] S. Hod, Quasi-bound state resonances of charged massive scalar fields in the near-extremal Reissner–Nordström black-hole spacetime, Eur. Phys. J. C 77, 351 (2017), arXiv:1705.04726 [hep-th] .
  * Cardoso  _et al._ [2018b] V. Cardoso, J. L. Costa, K. Destounis, P. Hintz, and A. Jansen, Strong cosmic censorship in charged black-hole spacetimes: still subtle, Phys. Rev. D 98, 104007 (2018b), arXiv:1808.03631 [gr-qc] .
  * Destounis [2019b] K. Destounis, Charged Fermions and Strong Cosmic Censorship, Phys. Lett. B 795, 211 (2019b), arXiv:1811.10629 [gr-qc] .
  * Destounis  _et al._ [2019] K. Destounis, R. D. B. Fontana, F. C. Mena, and E. Papantonopoulos, Strong Cosmic Censorship in Horndeski Theory, JHEP 10, 280, arXiv:1908.09842 [gr-qc] .
  * Griffiths and Podolsky [2009] J. B. Griffiths and J. Podolsky,  _Exact Space-Times in Einstein’s General Relativity_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, Cambridge, 2009).
  * Panosso Macedo [2023] R. Panosso Macedo, Hyperboloidal approach for static spherically symmetric spacetimes: a didactical introduction and applications in black-hole physics, (2023), arXiv:2307.15735 [gr-qc] .
  * Liu  _et al._ [2019] H. Liu, Z. Tang, K. Destounis, B. Wang, E. Papantonopoulos, and H. Zhang, Strong Cosmic Censorship in higher-dimensional Reissner-Nordström-de Sitter spacetime, JHEP 03, 187, arXiv:1902.01865 [gr-qc] .
  * Mascher  _et al._ [2022] G. Mascher, K. Destounis, and K. D. Kokkotas, Charged black holes in de Sitter space: Superradiant amplification of charged scalar waves and resonant hyperradiation, Phys. Rev. D 105, 084052 (2022), arXiv:2204.05335 [gr-qc] .
  * Zenginoglu [2011] A. Zenginoglu, A Geometric framework for black hole perturbations, Phys. Rev. D 83, 127502 (2011), arXiv:1102.2451 [gr-qc] .
  * Zenginoglu [2008] A. Zenginoglu, Hyperboloidal foliations and scri-fixing, Class. Quant. Grav. 25, 145002 (2008), arXiv:0712.4333 [gr-qc] .
  * Ansorg and Panosso Macedo [2016] M. Ansorg and R. Panosso Macedo, Spectral decomposition of black-hole perturbations on hyperboloidal slices, Phys. Rev. D 93, 124016 (2016), arXiv:1604.02261 [gr-qc] .
  * Panosso Macedo  _et al._ [2018] R. Panosso Macedo, J. L. Jaramillo, and M. Ansorg, Hyperboloidal slicing approach to quasi-normal mode expansions: the Reissner-Nordström case, Phys. Rev. D 98, 124005 (2018), arXiv:1809.02837 [gr-qc] .
  * Panosso Macedo [2020] R. Panosso Macedo, Hyperboloidal framework for the Kerr spacetime, Class. Quant. Grav. 37, 065019 (2020), arXiv:1910.13452 [gr-qc] .
  * Zhu  _et al._ [2014] Z. Zhu, S.-J. Zhang, C. E. Pellicer, B. Wang, and E. Abdalla, Stability of Reissner-Nordström black hole in de Sitter background under charged scalar perturbation, Phys. Rev. D 90, 044042 (2014), [Addendum: Phys.Rev.D 90, 049904 (2014)], arXiv:1405.4931 [hep-th] .