# Pseudospectrum of horizonless compact objects: a bootstrap instability
mechanism

Valentin Boyanov1,2, Kyriakos Destounis3, Rodrigo Panosso Macedo4, Vitor
Cardoso2,5, José Luis Jaramillo6 1Departamento de Física Téorica and IPARCOS,
Universidad Complutense de Madrid, 28040 Madrid, Spain  2CENTRA, Departamento
de Física, Instituto Superior Técnico – IST, Universidade de Lisboa – UL,
Avenida Rovisco Pais 1, 1049 Lisboa, Portugal  3Theoretical Astrophysics,
IAAT, University of Tübingen, 72076 Tübingen, Germany  4STAG Research Centre,
University of Southampton, University Road SO17 1BJ, Southampton, UK  5 Niels
Bohr International Academy, Niels Bohr Institute, Blegdamsvej 17, 2100
Copenhagen, Denmark  6Institut de Mathématiques de Bourgogne (IMB), UMR 5584,
CNRS, Université de Bourgogne Franche-Comté, F-21000 Dijon, France

###### Abstract

Recent investigations of the pseudospectrum in black hole spacetimes have
shown that quasinormal mode frequencies suffer from spectral instabilities.
This phenomenon may severely affect gravitational-wave spectroscopy and limit
precision tests of general relativity. We extend the pseudospectrum analysis
to horizonless exotic compact objects which possess a reflective surface
arbitrarily close to the Schwarzschild radius, and find that their quasinormal
modes also suffer from an overall spectral instability. Even though all the
modes themselves decay monotonically, the pseudospectrum contours of equal
resonance magnitude around the fundamental mode and the lowest overtones can
cross the real axis into the unstable regime of the complex plane, unveiling
the existence of nonmodal pseudo-resonances. A pseudospectrum analysis further
predicts that fluctuations to the system may destabilize the object when next
to leading-order effects are considered, as the triggering of pseudo-resonant
growth can break the order-expansion of black-hole perturbation theory.

##  I Introduction

The historic detection of gravitational waves (GWs) [1] paved the way for
precision GW astrophysics and black hole (BH) spectroscopy to develop at an
unprecedented rate. The increase in sensitivity of ground-based
interferometers, as well as the arrival of next generation space-borne
detectors, will further improve our understanding of the gravitational
interaction [2, 3, 4]. BH binaries produce GWs through the coherent motion of
the sources and form ideal testbeds of the gravitational field in extreme
conditions. Together with new radio and deep infrared interferometry [5],
which provide direct images of the dark compact objects lurking at galactic
centers, we are now able to perform novel tests of issues such as the BH no-
hair theorem [6], dark matter and environmental effects [7, 8, 9, 10, 11, 12],
potential quantum effects at the horizon scale [13, 14] as well as the
existence of horizonless exotic compact objects (ECOs) [15, 13, 16, 17].

Our current understanding of stellar evolution suggests that when
astrophysical objects become sufficiently compact, the interactions within
their matter content can no longer counteract their enormous self-gravity. The
natural conclusion to their evolution would therefore be a collapse into a BH.
Nevertheless, the singularity present at the end of this classical collapse
process suggests that this picture is incomplete. Corrections may, on the one
hand, be localized to a small region around the singularity itself, leaving
the exterior as an effectively classical BH. On the other hand, they may
permeate outwards up to the horizon, or even slightly past it, removing the
classical trapped region altogether. This line of thought has inspired the
search for BH mimickers in the form of effectively classical horizonless ECOs,
such as boson stars [18], gravastars [19], or semiclassically sustained
ultracompact objects [20, 21], as well as other modifications with exotic
near-horizon structure such as firewalls [22] and fuzzballs [23, 24, 25].

It is generally accepted that the ringdown signal of a Kerr BH will be
dominated by its quasinormal modes (QNMs) [26, 27], with the QNM spectrum
being characterized by the BH’s mass and spin. This is due to the fact that
the ringdown frequencies coincide (see e.g. [28, 29] and references therein)
with the QNM frequencies (which describe excitations of the photon sphere).
However, it has been shown [13, 30] that the relaxation of a highly-compact
horizonless object is almost identical to that of a BH at early times,
consisting of the usual photon sphere modes [31, 32, 33, 34], even though
these no longer correspond to QNMs. The replacement of the event horizon by a
surface with at least partial reflectivity entraps modes between it and the
photon sphere, which subsequently emerge during the later stages of the
ringdown and manifest as successively damped repetitions of the initial signal
(known as echoes) which describe the proper vibrational spectrum of the ECO
itself. The QNMs of BHs and ECOs can differ dramatically since they are
defined by eigenvalue problems with different boundary conditions.
Particularly, the fundamental mode and first overtones of non-rotating ECOs
with high compactness become extremely long-lived (with very long decay
timescales) [13, 11].

To add to the difficulties in constructing BH mimickers, these long-lived
modes may well turn out to destabilize the object when effects beyond linear-
order perturbations are taken into consideration [35, 36, 37]. This becomes
more clear when rotation is introduced; ultracompact objects with surfaces
within the ergosphere suffer from the ergoregion instability [38, 39] which
turns long-lived modes into unstable ones [40, 41, 35, 42]. Therefore,
constructing stable ECOs which can reproduce astrophysical BH observations
seems not to be as simple as finding ultracompact configurations with the
potential for counteracting their self-gravity in stationary solutions (which
in itself is a highly non-trivial task). In this work we take a novel approach
to studying the possible instability of ECOs by means of a non-selfadjoint
spectral tool recently introduced [43] to the subject area of GW physics: the
pseudospectrum of linear operators [44, 45, 46, 47].

In particular, ref. [43], followed by [48, 49, 50, 51], applied the concept of
pseudospectrum in the gravitational context, and first employed it to provide
a systematic explanation of the lack of robustness of QNMs against small
fluctuations of the underlying BH scattering potential, first described in
refs. [52, 53, 54, 55]. Broadly speaking, the pseudospectrum [44, 45, 46]
provides a characterization of analyticity properties of an operator’s Green
function in terms of a kind of topographic map on the complex plane. Contour
lines of this topographic map are denominated $\epsilon$-pseudospectra and
they delineate the regions in which eigenvalues can potentially “migrate”
under a system perturbation of size $\epsilon$. With this tool, one is able to
reveal whether a non-conservative system exhibits a spectral instability.
Particularly, the QNM spectrum is recovered in the limit $\epsilon\rightarrow
0$, and tightly packed contour lines around these eigenvalues indicate
spectral stability. On the other hand, $\epsilon$-contour lines extending far
from the eigenvalues indicate spectral instability. In practice, this implies
that a minuscule fluctuation in the underlying (QNM) time generator operator
can lead to a disproportionate (with respect to the scale of the fluctuation)
migration of QNMs in the complex plane, an effect already observed in several
seminal works [52, 53, 54, 55]. With the help of the pseudospectrum analysis,
the Schwarzschild and Reissner-Nordström BHs have been shown to suffer from
such an instability in their overtones [43, 49, 50, 48], with possible
implications also on their fundamental modes when the effective potential is
modified at large scales [43, 56].

More recently, ref. [51] has extended the range of applications of the
pseudospectrum in the setting of GW physics. Specifically, whereas refs. [43,
48, 49, 50] focus on the application of the pseudospectrum to the analysis of
QNM spectral instabilities, ref. [51] rather addresses a set of genuinely
dynamical phenomena, namely the possibility of transient growths in linear
dynamics driven by non-selfadjoint operators, on the one hand, and the
possibility of so-called pseudo-resonances, on the other hand. Indeed, while
at late times the linear system is dominated by the decaying QNMs of the
underlying evolution operator, its non-self-adjoint structure may lead to
unexpected early-time transient growth that is neither related to non-linear
effects nor to resonant-like responses of the system to external forcing or to
proper resonant frequencies (namely, so-called ‘pseudo-resonances’), as
observed in hydrodynamics [57, 44, 58]. Such transient and pseudo-resonant
effects are not captured in the eigenvalues, but are encoded in the full
pseudospectrum structure and estimated by a set of specific quantities
constructed out of it. In the gravitational case, ref. [51] shows the absence
of such early-time transient growth for asymptotically flat, spherically
symmetric BH spacetimes. This result is independent of the specific underlying
potential dictating the wave dynamics, and can rather be mathematically
related to the fact that the field exits through the spacetime boundaries at
the horizon and at infinity. Thus, even though the loss of energy through the
BH horizon and at the far wave zone is responsible for the non-self-adjoint
character of the operator, the same property turns out to preclude early-time
transient growth (in spherical symmetry).

Beyond dynamical transients, ref. [51] also emphasizes that the pseudospectrum
framework provides an efficient tool to account for the resonant-like
dynamical growth of perturbations when the non-self-adjoint linear system is
acted upon by external forces. This fact is particularly relevant for the non-
linear gravitational interaction, since the first order perturbation solution
acts as a source for the second order dynamics [59]. In particular, pseudo-
resonances can be triggered by external forces when the (real) frequencies in
the Fourier decomposition of these forces fall close to values where
pseudospectrum contour lines of small $\epsilon$ cross the real axis into the
unstable half of the complex plane. Thus, the pseudospectrum of the non-
selfadjoint second-order perturbation operator (which turns out to have the
same form as the linear-order one) can be used to anticipate potential non-
linear instabilities, in particular when $\epsilon-$pseudospectrum contour
lines for small $\epsilon$ protrude significantly into the part of the complex
plane where unstable modes would reside.

Our goal in this work is to analyze the pseudospectrum of the linear operator
describing the scalar and gravitational perturbations of a spherically-
symmetric ECO with a reflective surface, in search for information regarding
spectral instability, transient growth, and the potential for pseudo-resonant
phenomena. As an introductory example for these concepts, we first consider
the dynamics dictated by the wave equation with a Pöschl-Teller (PT)
potential, together with a reflective boundary condition at some finite
radius. This simplified scenario already captures many features we expect to
encounter in the behavior of linear perturbations on ECOs. Then, to model
actual ECOs we consider a Schwarzschild exterior geometry with a reflective
spherical surface above—but arbitrarily close to—the Schwarzschild radius.

The pseudospectra of the aforementioned configurations reveal, first of all, a
similar overtone spectral instability to the one present in BH geometries [43,
48]. More importantly, we have observed that the pseudospectrum around the
lowest-lying long-lived QNMs of ECOs has an important structural feature:
contour lines of equal resonant magnitude protrude into the unstable region of
the complex plane. We use this feature of the pseudospectrum to study the
possibility of transient growth in the linear homogeneous problem, as well as
the magnitude of the pseudo-resonant amplification of sources with (real)
frequencies close to the long-lived QNMs, that exist due to stable trapping.
Based on the fact that the homogeneous part of the equation governing second-
order perturbations can be written in the same form as that of the first-order
(whose pseudospectrum is identical) in an appropriate gauge [60], we show that
a “bootstrap” instability can be triggered, which incites a major discrepancy
between linear predictions and the full non-linear evolution. Specifically,
the linear-order perturbations excite pseudo-resonances at second-order,
leading to a possible breakdown of modal analysis and to an eventual non-
linear destabilization of the ECO.

This work is structured as follows. Section II introduces the scenarios we
will work with in terms of the equations which describe them. This section
also introduces the hyperboloidal framework used to formulate the systems as
spectral problems of non-self-adjoint operators, as well as the numerical
methods employed in their analysis. Then section III reviews the notion of
pseudospectrum and the relevant quantities introduced in ref. [51] to assess
the early-time transient dynamics. Finally, section IV presents our main
results, before our concluding remarks in section V. In what follows we use
geometrized units, such that $G=c=1$.

##  II Perturbations of spherically-symmetric ECOs

Perturbation theory in spherically-symmetric backgrounds is typically
described by a dynamical field $\Phi$ obeying the wave equation

$-\Phi_{,{t}{t}}+\Phi_{,{x}{x}}-V\Phi=0,$ (1)

where $t$ is a time coordinate parametrizing the evolution, and $x$ is a
radial coordinate defining the spatial domain of integration;
$x\in(-\infty,+\infty)$ for the standard tortoise coordinate. As is usual, a
comma in the subscript denotes partial differentiation. The potential $V$
depends on the spacetime background and the nature of the perturbation field.

In asymptotically flat BH spacetimes, $V\rightarrow 0$ as
$x\rightarrow\pm\infty$, and we solve eq. (1) with the left- and right-moving
boundary conditions

$\Phi_{,t}\pm\Phi_{,x}\sim 0,\quad{\rm for}\quad x\rightarrow\pm\infty,$ (2)

which impose that the wave be purely ingoing at the horizon and purely
outgoing at infinity.

When modeling an ECO, however, one is interested in a solution on the domain
$x\in[x_{o},+\infty)$, where a boundary condition involving a combination of
$\Phi(t,x_{o})$ and $\Phi_{,x}(t,x_{o})$ is imposed at the finite point
$x=x_{o}$. Here we make the simple assumption that all incoming waves are
completely reflected at the ECO’s surface through a Dirichlet boundary
condition

$\Phi(t,x_{o})=0.$ (3)

Alternative ways of imposing boundary conditions at the surface have been
explored, for instance, in [61].

The problem can be studied in the frequency domain, where the Fourier
decomposition $\Phi(t,x)\sim\phi(x)e^{i\omega t}$ yields the eigenvalue
problem

${\phi}_{,{x}{x}}+(\omega^{2}-V)\phi=0.$ (4)

QNMs are solutions of eq. (4) which satisfy purely outgoing boundary
conditions at infinity and appropriate boundary conditions at the ECO’s
surface. Particularly, the Fourier space equivalent of the BC (2) at infinity,
and of (3) at the surface, which read

$\displaystyle\phi(x)\sim e^{-i\omega x}\quad$ $\displaystyle{\rm for}$
$\displaystyle\quad x\rightarrow\infty,$ (5)
$\displaystyle\phi(x_{o})=0.\quad$ $\displaystyle{\rm at}$ $\displaystyle\quad
x=x_{o}.$ (6)

###  II.1 Hyperboloidal foliation

In this work, we follow refs. [62, 63, 64, 65, 66, 67, 43] in incorporating
the asymptotic boundary conditions via the hyperboloidal approach to
perturbation theory. By exploiting the freedom in the choice of coordinates,
we deform the time slices and compactify the radial direction via

$t=\lambda\left(\tau-H(\chi)\right),\quad x=\lambda\,g(\chi),$ (7)

with $H(\chi)$ the so-called height function, which ensures that the wave zone
(future null infinity) is approached as $x\rightarrow\infty$, $g(\chi)$ a
function mapping $x\in[x_{o},\infty)$ into the compact interval
$\chi\in[\chi_{{}_{{I}^{+}}},\chi_{o}]$, and $\lambda$ a characteristic
(constant) length scale of the problem.

The wave equation (1) takes the form [43, 48]

$-\Phi_{,\tau\tau}+L_{1}\Phi+L_{2}\Phi_{,\tau}=0,$ (8)

where

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

and the corresponding functions

$\displaystyle
p=\dfrac{1}{|g^{\prime}|},\,\gamma(\chi)=ph^{\prime},\,w=\dfrac{1}{p}\left(1-\gamma^{2}\right),\,\hat{V}=\dfrac{\lambda^{2}}{p}V,$
(10)

where the prime denotes the derivative with respect to $\chi$.

Equation (8) can be written as a matrix evolution problem through a first-
order reduction in time by introducing $\psi=\phi_{,\tau}$ such that

$\dot{u}=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}\phi\\\
\psi\end{array}\right).$ (11)

Given initial data $u(0,\chi)$, eq. (11) has the formal solution

$u(\tau,\chi)=e^{iL\tau}u(0,\chi),$ (12)

in terms of the _evolution operator_ $e^{iL\tau}$. By further performing a
harmonic decomposition $u(\tau,\chi)\sim\hat{u}(\chi)e^{i\omega\tau}$ in eq.
(11) we arrive at the eigenvalue equation

$L\hat{u}_{n}=\omega_{n}\hat{u}_{n},$ (13)

where, imposing the surface boundary condition (3), $\omega_{n}$ becomes an
infinite set of QNMs of the operator $L$, with $n\geq 0$ the overtone number.

Note that the outgoing boundary condition (2) for $x\to\infty$ is
automatically imposed by the hyperboloidal slicing. Indeed, the main
motivation for employing the hyperboloidal approach is the geometric
imposition of outgoing boundary conditions at future null infinity (and at the
event horizon of BH spacetimes). With the light cones pointing outwards from
the integration domain, the causal degrees of freedom must satisfy the
regularity conditions underlying the wave equation. In practical terms, this
property manifests via $p(\chi_{{}_{{I}^{+}}})=0$, i.e. the term proportional
to the second spatial derivative vanishes, and the remaining term imposes the
appropriate condition for the relation between $\phi(\chi_{{}_{{I}^{+}}})$ and
$\phi_{,\chi}(\chi_{{}_{{I}^{+}}})$.

At the ECO boundary $\chi_{o}$, however, the wave equation does not contain
any singularity, and we are free to impose any desired external boundary
condition. In particular, for our simple ECO model we choose to impose eq.
(3), or equivalently eq. (5), which translates into

$u(\chi_{o},\tau)=\hat{u}(\chi_{o})=0.$ (14)

As in ref. [43], we explicitly construct working examples for the cases of
dynamics dictated by either the PT potential or the external Schwarzschild
spacetime potential. While the former provides us with a simple toy model to
exemplify the ideas presented in the work and set a benchmark for the results,
the latter has direct astrophysical relevance when modeling ECOs.

####  II.1.1 Pöschl-Teller Potential

Let us first analyze the simpler case of eq. (1) with the PT potential,

$V(x)=\frac{V_{0}}{\cosh^{2}(x/\lambda)},$ (15)

where $V_{0}$ is a positive constant. We will use units $V_{0}=1$ by setting
the length scale to $\lambda=\sqrt{V_{0}}$ in transformations (7). The
hyperboloidal coordinates adequate for the asymptotic behavior of this
potential are [43]

$t=\lambda\bigg{[}\tau-\frac{1}{2}\ln(1-y^{2})\bigg{]},\quad{x}=\lambda\,\text{arctanh}(y),$
(16)

where $y\in[-1,1]$. Rather than using the whole range of $y$ for our toy ECO
construction, we place a Dirichlet boundary condition at
$y_{o}=1-\mathcal{E}$.

It is, thus, more convenient to introduce a new spatial coordinate
$\chi\in[-1,1]$ via

$y=y_{o}\dfrac{1+\chi}{2}-\dfrac{1-\chi}{2}=\chi-\frac{\mathcal{E}}{2}(1+\chi),$
(17)

which maps the PT-ECO boundary $y_{o}$ to $\chi_{o}=1$, while keeping future
null infinity $y_{{}_{{I}^{+}}}=-1$ at $\chi_{{}_{{I}^{+}}}=-1$. Inserting
(17) into (16), and making use of the definitions in eqs. (10) one obtains the
functions

$\displaystyle p(\chi)$ $\displaystyle=$ $\displaystyle
1-\chi^{2}+\dfrac{\mathcal{E}}{2}(1+\chi)^{2},$ (18)
$\displaystyle\gamma(\chi)$ $\displaystyle=$
$\displaystyle-\chi+\dfrac{\mathcal{E}}{2}(1+\chi),$ (19) $\displaystyle
w(\chi)$ $\displaystyle=$
$\displaystyle\hat{V}(\chi)=1-\dfrac{\mathcal{E}}{2},$ (20)

from which the operators $L_{1}$ and $L_{2}$ follow directly.

As already mentioned, $p(\chi)$ vanishes at $\chi_{{}_{{I}^{+}}}=-1$, but not
at $\chi_{o}=1$ when $\mathcal{E}\neq 0$. Thus, in this case we have the
freedom to specify the boundary condition at the PT-ECO surface. Furthermore,
we recover the expressions from ref. [43] in the limit ${\cal E}\rightarrow
0$, in which $\chi=1$ also corresponds to an asymptotic region.

####  II.1.2 Regge-Wheeler Potential

A similar set of hyperboloidal coordinates is introduced for perturbations
propagating on spherically-symmetric ECO of mass $M$ with a surface located at
a radial position

$r_{0}=(1+\mathcal{E})r_{\rm h},\quad\mathcal{E}\ll 1,$ (21)

where $r_{h}=2M$. The geometry outside the ECO is described by the
Schwarzschild metric. In the standard Schwarzschild coordinates
$(t,r,\theta,\varphi)$, the line element reads

${\rm d}s^{2}=-f(r){\rm d}t^{2}+f(r)^{-1}{\rm d}r^{2}+r^{2}{\rm
d}\Omega^{2}_{2},$ (22)

where $f(r)=1-r_{\rm h}/r$ and ${\rm d}\Omega^{2}_{2}$ is the line element of
the unit $2$-sphere. The radial coordinate $x$ defining the spatial domain in
eq. (1) is the tortoise coordinate

$x=r+r_{\rm h}\log\left(\frac{r}{r_{\rm h}}-1\right),$ (23)

whereas the potential $V$ follows from the Regge-Wheeler-Zerilli formalism
[27]. In this work, as an example, we focus on the Regge-Wheeler potential

$V=\dfrac{f(r)}{r^{2}}\left(\ell(\ell+1)+(1-s^{2})\frac{r_{\rm h}}{r}\right),$
(24)

with $\ell$ the angular number and $s$ the spin parameter characterizing the
scalar ($s=0$), vector ($\pm 1$) and tensorial ($\pm 2$) sectors of the metric
perturbation.

We begin with the hyperboloidal coordinates $(\tau,\sigma,\theta,\varphi)$
introduced in [65, 43]

$t=\lambda\Bigg{[}\tau-\dfrac{r_{\rm
h}}{\lambda}\left(\ln(\sigma)+\ln(1-\sigma)-\sigma^{-1}\right)\Bigg{]},\,r=\dfrac{r_{\rm
h}}{\sigma},$ (25)

which map the BH horizon $r=r_{\rm h}$ to $\sigma=1$ and future null infinity
to $\sigma=0$. Then, according to eq. (21), the ECO surface $r_{o}$ is located
at

$\sigma_{o}=\dfrac{r_{\rm h}}{r_{o}}=\bigg{(}1+\mathcal{E}\bigg{)}^{-1}\approx
1-\mathcal{E},\quad{\cal E}\ll 1.$ (26)

Similar to the PT case, we introduce a final coordinate $\chi\in[0,1]$ via

$\sigma=\sigma_{o}\chi.$ (27)

Thus, we retain the location of future null infinity at
$\chi_{{}_{{I}^{+}}}=0$, while mapping the ECO surface to $\chi_{o}=1$.

The transformations (25) and (27), combined with eq. (23) yield the functions

$\displaystyle p(\chi)$ $\displaystyle=$ $\displaystyle\dfrac{\lambda}{r_{\rm
h}}\dfrac{\chi^{2}\left(1+{\cal E}-\chi\right)}{(1+{\cal E})^{2}},$ (28)
$\displaystyle\gamma(\chi)$ $\displaystyle=$ $\displaystyle
1-\dfrac{2\chi^{2}}{(1+{\cal E})^{2}},$ (29) $\displaystyle w(\chi)$
$\displaystyle=$ $\displaystyle\dfrac{4r_{\rm h}}{\lambda}\dfrac{(1+{\cal
E}+\chi)}{(1+{\cal E})^{2}},$ (30) $\displaystyle\hat{V}(\chi)$
$\displaystyle=$ $\displaystyle\dfrac{\lambda}{r_{\rm
h}}\left[\frac{\ell(\ell+1)}{1+{\cal E}}+\frac{(1-s^{2})\chi}{(1+{\cal
E})^{2}}\right].$ (31)

Once again we emphasize the singular character of the equation at future null
infinity, captured by the property $p(0)=0$. Imposing the regularity of the
perturbation field there, automatically fixes the necessary boundary
condition. At $\chi=1$, the equation does not degenerate, and we are free to
impose the boundary conditions fixed by eq. (14). This freedom is no longer
available for ${\cal E}\rightarrow 0$. In this limit, eqs. (28)-(31) reduces
to the corresponding expression for the Schwarzschild BH spacetime [43], in
which $p(1)$ also vanishes.

This compactified hyperboloidal approach allows us to cast the relevant
equations in terms of the spectral problem of a non-selfadjoint operator, and
to employ the corresponding tools that recover the results from the theory of
scattering resonances (see refs. in [43, 50]). Section III introduces several
pseudospectrum (or ‘non-modal analysis’ [68]) quantities which we will use in
this context. Hereafter, we solve the underlying equations and evaluate these
quantities numerically. The following section presents the numerical methods
employed.

###  II.2 Numerical Method

We employ a collocation spectral method with the Chebyshev polynomials as a
basis [69, 70, 44, 71]. To that end, we fix a numerical resolution $N$ (number
of spatial discretization points) and introduce the Chebyshev-Lobatto grid

$\chi_{j}=\cos(j\,\pi/N),\quad j=0\cdots N,$ (32)

which allows the discretization of the differential operators from eqs. (9) in
terms of their corresponding spectral approximation [69, 70, 71, 44] (when the
interval $\left[\chi_{{}_{{I}^{+}}},\chi_{o}\right]$ is different form
$[-1,1]$, we employ a linear rescaling). The resulting matrix $L$ in eq. (11)
has dimensions $\mathcal{N}=2(N+1)\times 2(N+1)$.

The reflective boundary condition at the ECO surface is subsequently imposed
by removing the row and column of the discretized matrix $L$ corresponding to
the position of the surface in the discrete spatial grid (see e.g. [69]). Then
it follows straightforwardly from (13) that the spectrum $\sigma(L)$ of this
discretized and dimensionally-reduced (due to the reflective boundary
condition) $L$ consists of the QNMs of the ECO with purely outgoing boundary
conditions at future null infinity and a Dirichlet boundary condition at its
surface.

We note that we explored several different methods of imposing the necessary
boundary conditions, before finding that the above is the most efficient for
numerical computation. For instance, to include future null infinity into the
computational domain, and thus impose the outgoing condition there through
regularity, a spacetime slicing can also be performed along the null outgoing
Eddington-Finkelstein coordinate $u=t-x$ (toward which the hyperboloidal
$\tau$ tends at infinite radius). In BH spacetimes, this slicing method is not
viable because the limit $r\rightarrow r_{h}$ along the surface $u=$const.
leads to the white-hole horizon, where the ingoing boundary condition for QNMs
is not compatible. However, one encounters no problem in the ECO setup because
the domain is cut at the ECO surface $x_{o}$. On the other hand, to impose
boundary at ECO surface, we have explored the alternative methods described in
refs. [69, 72, 44] as well as the possibility of enforcing the boundary
condition at the ECO surface via an auxiliary field $u=(1-\chi)w$.

##  III Pseudospectrum and non-modal analysis: dynamical transients

In this section we introduce the mathematical concepts relevant to our study
of ECOs, closely following the discussion in [51]. In particular, we provide a
brief introduction to the notion of pseudospectrum (for more detail see e.g.
refs. [44, 43]) and we review the concepts of numerical abscissa and Kreiss
constant. As we will explain, these quantities allows us to obtain information
regarding non-modal transient effects, as introduced to the context of GW
physics in [51].

###  III.1 Pseudospectrum

Before we introduce the pseudospectrum of an operator $L$, let us first recall
that the operator’s spectrum $\sigma(L)$ is understood as the set of complex
values $\lambda$ where the so-called resolvent

$R_{L}(\lambda)=(L-\lambda\mathbb{I})^{-1}$ (33)

(with $\mathbb{I}$ the identity operator) is not defined as a bounded
operator. Essentially, the resolvent is the Green function of the eigenvalue
problem $(L-\lambda\mathbb{I})u=0$ . This definition generalises our practical
experience with finite matrices, where the spectrum is obtained by finding
discrete $\lambda$-values such that
$\det\left[L-\lambda\mathbb{I}\right]=\det\left[R_{L}(\lambda)^{-1}\right]=0$.

Building upon this notion, the $\epsilon$-pseudospectrum
$\sigma^{\epsilon}(L)$ is characterised by the set of complex values $\lambda$
such that $\|R_{L}(\lambda)^{-1}\|<\epsilon$. Two aspects are worth
mentioning: (i) the spectrum $\sigma(L)$ is recovered when
$\epsilon\rightarrow 0$, and (ii) in contrast to the spectrum, the
$\epsilon$-pseudospectrum depends on the choice of a scalar product [49]. We
will discuss the latter in more detail in sec. III.3.

In fact, it can be shown that the following definitions of the pseudospectrum
are equivalent [44]

$\displaystyle\sigma^{\epsilon}(L)$ $\displaystyle=$
$\displaystyle\\{\lambda\in\mathbb{C}:\|R_{L}(\lambda)\|>1/\epsilon\\},$ (34)
$\displaystyle=$ $\displaystyle\\{\lambda\in\mathbb{C},\exists\delta
L,\|\delta L\|<\epsilon:\lambda\in\sigma(L+\delta L)\\},$ (35)
$\displaystyle=$ $\displaystyle\\{\lambda\in\mathbb{C},\exists
v\in\mathbb{C}^{n}:\|Lv-\lambda v\|<\epsilon\|v\|\\}.$ (36)

The first is the aforementioned one, where the norm of the resolvent provides
an $\epsilon-$contour level bounding the set $\sigma^{\epsilon}(L)$. This
prescription is particularly useful for numerical computations.

The second definition illustrates how the spectral instability under a
perturbation $\delta L$ is encoded in the pseudospectrum. In other words, when
a perturbation of magnitude $\epsilon$ is added to the operator $L$, the
perturbed spectrum $\sigma(L+\delta L)$ can be displaced only within the
confines of the $\epsilon$-pseudospectral contour lines, and the
$\epsilon$-pseudospectrum can in fact be defined as the set of all possible
displacements of the spectrum under all possible perturbations of equal
magnitude $\epsilon$.

Finally, the third definition introduces the concept of a pseudo-eigenmode $v$
of the operator, which is connected to the notion of _pseudo-resonance_ : a
point in the complex plane (not coincident with the spectral points) where the
resolvent norm has a very large value.

Upon calculating the pseudospectrum, a general criterion for visually
interpreting the result is the following: for each value of $\epsilon>0$, if
$\sigma^{\epsilon}(L)$ is a set of circular regions of radius $\epsilon$
around $\sigma(L)$ in the complex plane, then the operator can be said to be
spectrally stable (see e.g., fig. 4 in [43] and fig. 10 in [48]). Otherwise,
if $\sigma^{\epsilon}(L)$ encompasses regions in the complex plane larger than
these, $L$ is spectrally unstable, as has been found for the QNM operators in
the Schwarzschild and Reissner-Nordström BHs [43, 48, 50].

As discussed in the introduction and in further detail in refs. [57, 44, 68,
51], the pseudospectrum contains more information than just the
characterization of the spectral stability of $L$. In the $\epsilon\to 0$
limit, it determines the late-time evolution of the problem through the QNM
spectrum,111We restrict the discussion to the discrete part of the spectrum
$\sigma(L)$ yielding the QNM ring-down behaviour. As it is well known for
asymptotically flat BH spacetimes, one also encounters a very late time power-
law decay [73]. In the current context, this tail is understood in terms of
the continuous part of the spectrum $\sigma(L)$ [74, 65, 66] but away from
this limit it provides information regarding the early and intermediate-time
evolution of the perturbations, as well as the existence and magnitude of
possible pseudo-resonances, which the eigenvalue analysis fails to identify.

###  III.2 Transient growth of linear perturbations

Following closely the ideas put forward in [51] in the gravitational setting
(and more generally in refs. [57, 44, 68]), we will briefly present the
quantities related to evolution operator $e^{iL\tau}$ which allow us to obtain
information regarding the early-time evolution of perturbations. We recall
that the dynamics of the system is given by the wave equation expressed in the
form (11), and the evolution operator arises when expressing the formal
solution to the system (12). In particular, the operator’s norm
$\|e^{iL\tau}\|$ monitors the maximum growth rate of solutions $u(\tau)$ in
time.222For simplicity, in this section we omit the solution’s dependence on
the spatial coordinate $\chi$.

In the standard modal analysis of the dynamical stability of the wave
equation, one typically focuses on the fundamental QNM $\omega_{0}$ —
particularly, on the sign of its imaginary part. Modal stability is shown when
${\rm Im}(\omega_{0})$ takes on values in the ‘stable’ half of the complex
plane, i.e. ${\rm Im}(\omega_{0})>0$ with the convention used in this work.
The quantity typically used to monitor this late-time stability is the so-
called _spectral abscissa_ 333Note the $-\pi/2$ rotation as compared with the
standard ‘abscissa’, due to the ‘$i$’ factor in the definition of the
infinitesimal time generator, as compared with the notation in e.g. [44],
namely $\alpha(iL)=\sup_{z\in\sigma(iL)}{\rm Re}(z)$. This accounts for the
minus sign in the definition (see [51] for a more complete discussion).

$\alpha(L)=-\inf_{z\in\sigma(L)}{\rm Im}(z),$ (37)

which formally identifies ${\rm Im}(\omega_{0})$ (when a spectral gap exists
between $\omega_{0})$ and the real line).

For self-adjoint operators $\|e^{iL\tau}\|\sim e^{\alpha(L)\tau}$, implying
that the growth rate of solutions $u(\tau)$ is well estimated by the decaying
exponential $e^{-\tau{\rm Im}(\omega_{0})}$. The picture changes when $L$ is
not self-adjoint. In this case, the behavior $\|e^{iL\tau}\|\sim e^{-\tau{\rm
Im}(\omega_{0})}$ manifests only at late times. Throughout the rest of the
evolution, the spectral abscissa only provides the lower bound

$\|e^{iL\tau}\|\geq e^{\alpha(L)\tau}=e^{-\tau{\rm
Im}(\omega_{0})},\quad\forall\tau\geq 0,$ (38)

leaving open the possibility of initial transient growth, even if at late
times the solutions still decays due to ${\rm Im}(\omega_{0})>0$.

The possibility for transient growth has far-reaching implications in physics,
where the linear evolution is often only an approximation, the validity of
which relies on the absence of such growth. For instance, in hydrodynamics
[57] the accumulation of transient growth in perturbations of certain
configurations has been shown to lead to a breakdown of the linear
approximation, from where non-linear instabilities take over — all this in a
system where an eigenvalue analysis would predict stability.

It is therefore important to have a full grasp of the evolution of system
described by (12) at different times. Specifically, one can identify three
distinct phases to the evolution dictated by an operator with stable
eigenmodes (${\rm Im}(\omega_{0})>0$): a possible (i) initial growth, reaching
a (ii) transient maximum, before proceeding to the (iii) late-time decay.
Specific tools are available to monitor each of these stages. While stage
(iii) is described by the spectral abscissa $\alpha(L)={\rm Im}(\omega_{0})$,
we will now introduce the quantities which provide information on stages (i)
and (ii).

####  III.2.1 Numerical abscissa

The _numerical abscissa_ , that is defined in terms of the so-called
‘numerical range’ of $L$ (see [44, 51]), characterizes the initial slope of
the norm of the evolution operator

$\omega(L)=\left.\frac{\text{d}}{\text{d}\tau}\|e^{iL\tau}\|\right|_{\tau=0}.$
(39)

Therefore, a value $\omega(L)>0$ indicates the presence of an initial growth
in the evolution. This quantity can be related to the
$\epsilon-$pseudospectrum in the limit $\epsilon\rightarrow\infty$ [51].
Moreover, this quantity also provides an upper bound on the evolution
operator’s norm at any time $\tau>0$ through the relation

$\|e^{iL\tau}\|\leq e^{\omega(L)\tau},\quad\forall\tau\geq 0.$ (40)

Particularly, one sees that when $\omega(L)=0$ the solution of the evolution
problem cannot grow and, generically, its norm decreases in time (the operator
$L$ in this case is said to be contractive, namely $\|e^{iL\tau}\|\leq
1,\forall\tau>0$), whereas when $\omega(L)>0$ the possibility of transient
growth remains.

For practical purposes, the numerical abscissa can also be characterized and
computed [51] through the eigenvalue the problem

$\dfrac{1}{2i}\left(L^{\dagger}-L\right)u_{L}=\lambda u_{L}.$ (41)

In particular, the eigenvalues $\lambda$ of this problem are real numbers (by
construction, since the operator in the left-hand side is selfadjoint), and
$\omega(L)$ is simply the largest one.

As mentioned above, the spectral and numerical abscissa are related to the
$\epsilon\rightarrow 0$ and $\epsilon\rightarrow\infty$ limits of the
$\epsilon-$pseudospectra, respectively. Given the dynamical stages of the
evolution they are each related to, a time scale $\tau_{\epsilon}\sim
1/\epsilon$ was introduced in ref. [51] to provide an estimation of the
relation between different moments of the evolution governed by the operator
(12) and the $\epsilon$-pseudospectral contours.

####  III.2.2 Pseudospectral abscissa and Kreiss constant

To estimate the maximum of the possible dynamical transient, that is
controlled by the maximum growth of the norm of the evolution operator at
intermediate times, we can utilize the so-called _pseudospectral abscissa_ of
$L$, defined in our setting as 444Again, as in eq. (37), the minus sign
follows from the factor “$i$” in the infinitesimal time generator.

$\alpha_{\epsilon}(L)=-\inf_{z\in\sigma^{\epsilon}(L)}{\rm Im}(z).$ (42)

This quantity allows us to put a sharper lower bound on the evolution norm
than eq. (38), giving a lower bound to a possible transient peak. This bound
is obtained from the $\epsilon$-pseudospectrum at intermediate values of
$\epsilon$ through the relation

$\sup_{\tau\geq
0}\|e^{iL\tau}\|\geq\frac{\alpha_{\epsilon}(L)}{\epsilon},\quad\forall\epsilon>0.$
(43)

In other words, if the $\epsilon-$pseudospectral contour lines protrude
significantly (as quantified by $\alpha_{\epsilon}(L)>\epsilon$) into the
‘unstable’ half of the complex plane, then transient growth occurs in the
evolution and the peak corresponds to the intermediate value of $\epsilon$
that maximizes the right-hand side of eq. (43). This naturally leads, upon
maximizing eq. (43) over $\epsilon>0$, to the absolute lower bound

$\sup_{t\geq 0}\|e^{iL\tau}\|\geq{\cal K}(L)$ (44)

in terms of the so-called _Kreiss constant_

$\displaystyle\mathcal{K}(L)$ $\displaystyle=$
$\displaystyle\sup_{\epsilon>0}\frac{\alpha_{\epsilon}(L)}{\epsilon}$
$\displaystyle=$ $\displaystyle\sup_{{\rm Im}(z)<0}\left\\{|{\rm
Im}(z)|\,\cdot\|R_{L}(z)\|\right\\},$

where the characterization in the second line (cf. [44, 51]) provides an
efficient approach to its calculation. If the fundamental QNM were unstable,
i.e. if $\mathrm{Im}(\omega_{0})<0$ and therefore $\alpha(L)=\lim_{\epsilon\to
0}\alpha_{\epsilon}(L)>0$, then the Kreiss constant $\mathcal{K}(L)$ would be
directly divergent (as follows from the first line in (III.2.2) by taking
$\epsilon\to 0$), consistently with the unbounded growth such a mode would
produce. On the other hand, if $\alpha(L)\leq 0$, then ${\cal K}(L)$ can be
either larger than 1, in which case there would be necessarily a transient
growth, or equal to 1. In the latter case, no conclusion can be drawn in
general, since (44) saturates and trivializes for $\tau=0$ and the inequality
is not significant. However, if the numerical abscissa actually vanishes
$\omega(L)=0$ (indicating through (39) the absence of an initial growth), then
the Kreiss constant must indeed be unity ${\cal K}(L)=1$ (cf. discussion [51])
and in this (contractive) case we know a priori that no transient happens
along the evolution.

At a practical level, the calculation of the ratio
$\alpha_{\epsilon}(L)/\epsilon$ comes “for free” once the pseudospectrum is
known. All one needs to do is take the point of an $\epsilon$-pseudospectrum
contour line with minimum imaginary part. Particularly, lines which protrude
into the negative imaginary half-plane are the ones which will provide a
sensible lower bound. The supremum corresponding to the Kreiss constant is
then reached somewhere between the first such protruding $\epsilon$-contour
and the outermost one, i.e. $\epsilon\to\infty$.

Note that all the quantities discussed in this section rely on specifying a
norm for the operators involved. When calculating the pseudospectrum, and
quantities derived therefrom, this norm is used to evaluate $\|R_{L}(z)\|$.
For the explicit calculation of the numerical abscissa via eq. (41), that
assumes a norm associated with a scalar product, such scalar product enters in
the very construction of the adjoint operator $L^{\dagger}$. We now proceed to
define the norm, and the scalar product it is associated with, which we will
employ in this work.

###  III.3 Energy norm

To choose the norm with which to evaluate the above quantities, we follow past
works [43, 48, 50, 49] in appealing to a physical argument for the measurement
of the magnitude of the perturbations described by (1). Particularly, we will
employ the _energy norm_ (see also [75]), which quantifies this magnitude
through its contribution to the system’s energy (for an in-depth discussion of
the role of the norm and, in particular, the energy norm, see [49]).

The energy associated to the state vector $u=\left(\phi,\psi\right)^{T}$
satisfying the wave equation (11) reads

$\displaystyle E(\phi,\psi)$ $\displaystyle=\|u\|^{2}_{{}_{E}}$ (46)
$\displaystyle=\frac{1}{2}\int_{\chi_{{}_{{I}^{+}}}}^{\chi_{o}}\left(\omega(\chi)|\psi|^{2}+p(\chi)|\partial_{\chi}\phi|^{2}+\hat{V}(\chi)|\phi|^{2}\right)d\chi.$

The subscript $E$ denotes the energy norm and the integration limits
correspond to the compactified spatial interval between future null infinity
and the ECO surface in our hyperboloidal scheme. An energy scalar product then
follows as the quadratic form associated to such energy norm, in the form

$\displaystyle\langle
u_{1},\\!u_{2}\rangle_{{}_{E}}=\Big{\langle}\begin{pmatrix}\phi_{1}\\\
\psi_{1}\end{pmatrix},\begin{pmatrix}\phi_{2}\\\
\psi_{2}\end{pmatrix}\Big{\rangle}_{{}_{E}}$ (47)
$\displaystyle=\frac{1}{2}\int_{\chi_{{}_{{I}^{+}}}}^{\chi_{o}}\left(\omega(\chi)\bar{\psi}_{1}\psi_{2}+p(\chi)\partial_{\chi}\bar{\phi}_{1}\partial_{\chi}\phi_{2}+q(\chi)\bar{\phi}_{1}\phi_{2}\right)d\chi,$

such that, by construction, $\|u\|^{2}_{{}_{E}}=\langle u,u\rangle_{{}_{E}}$.
With eq. (47) we can evaluate the pseudospectrum according to definition (34),
and we can also analytically deduce the expression for the adjoint
$L^{\dagger}$ [43, 49].

More specifically, the $\epsilon$-pseudospectrum in the energy norm
$\sigma^{\epsilon}_{{}_{E}}(L)$ can be calculated as the set

$\sigma^{\epsilon}_{{}_{E}}(L)=\\{\omega\in\mathbb{C}:s_{{}_{E}}^{\mathrm{min}}(\omega\mathbb{I}-L)<\epsilon\\},$
(48)

with $s_{{}_{E}}^{\mathrm{min}}$ the minimum of the generalized singular value
decomposition

$s_{{}_{E}}^{\mathrm{min}}(\mathcal{M})=\min\\{\sqrt{\omega}:\omega\in\sigma(\mathcal{M}^{\dagger}\mathcal{M})\\}.$
(49)

Appendix C.3 of ref. [43] details the calculation of the discretized energy
scalar product within the numerical scheme presented in sec. II.2.

By following [43, 49], we also derive an analytic expression for the formal
adjoint $L^{\dagger}$

$\displaystyle\expectationvalue{L^{\dagger}u_{1},u_{2}}$ $\displaystyle=$
$\displaystyle\expectationvalue{u_{1},Lu_{2}}$ (50) $\displaystyle=$
$\displaystyle\expectationvalue{Lu_{1},u_{2}}+\expectationvalue{L^{\partial}u_{1},u_{2}},$
(51)

with $L^{\partial}$ yielding a boundary contribution in the form

$\expectationvalue{L^{\partial}u_{1},u_{2}}=i\left[2\gamma\bar{\psi}_{1}\psi_{2}+p\left(\bar{\phi}_{1}^{\prime}\psi_{2}+\phi_{2}^{\prime}\bar{\psi}_{1}\right)\right]\bigg{|}_{\chi_{{}_{I}^{+}}}^{\chi_{o}}.$
(52)

We emphasize that the term $\expectationvalue{L^{\partial}u_{1},u_{2}}$ is
precisely the actual responsible for the non-self-adjoint character to the
operator $L$.

At future null infinity, the only contribution to $L^{\partial}$ comes from
the term which contains $\gamma(\chi)$, since $p(\chi)$ always vanishes at
$\chi_{{}_{{I}^{+}}}$. At the ECO surface,
$\expectationvalue{L^{\partial}u_{1},u_{2}}$ vanishes altogether as a
consequence of the boundary condition (14). This result makes apparent the
fact that the reflective ECO surface does not contribute to the energy outflow
which makes the operator non-self-adjoint; the only such contribution comes
form waves propagating away towards infinity, and is the same as in the BH
case [43]. We can thus express the boundary operator $L^{\partial}$ in terms
of a Dirac-delta distribution

$L^{\partial}={i}\begin{pmatrix}0&0\\\
0&-2\gamma\delta(\chi-\chi_{{}_{{I}^{+}}})\end{pmatrix}.$ (53)

##  IV Results

We now apply the elements of the pseudospectrum framework introduced in the
previous section to the equations governing the dynamics of perturbations on
ECO spacetimes. As an example and motivation, we first study the toy model of
a wave equation with the (modified) PT potential. This model already captures
much of the results expected for ECO perturbations, but has a cleaner
numerical behavior which allows us to more easily discern physically relevant
results. After this analysis, we proceed to present the results for the
dynamics of perturbations propagating on a Schwarzschild solution with a
boundary condition which represents an ECO.

To summarize briefly, we find that the spectral instability present in the BH
case [43] is also present for ECOs. Furthermore, despite the pseudospectrum
contour lines we find protruding into the unstable half of the imaginary
plane, our numerical results provide strong indications that the evolution
operator $\|e^{iL\tau}\|$ is in fact contractive. This result is further
supported by an analytical calculation of the numerical abscissa. This entails
the absence of a transient growth but, following from the same analysis. we
argue that pseudo-resonances may indeed trigger instabilities when second
order effects are taken into account.

###  IV.1 Toy model: the Pöschl-Teller potential

Figure 1: Left Panel: QNMs (red dots) and pseudospectrum (white curves) of the
$L$ operator for the PT potential, with a Dirichlet boundary condition at
$\mathcal{E}=10^{-3}$ and $N=200$ spatial discretization points. Contour lines
correspond to values between $-4.7$ (uppermost contour) and $-1.1$ (lowest
contour) with steps of $0.3$ in the log scale. The horizontal dashed line
indicates the position where $\text{Im}(\omega)=0$. Right Panel: Lower bound
$\alpha_{\epsilon}(L)/\epsilon$ to the evolution operator norm for different
values of $\epsilon$, computed numerically from the pseudospectrum of the PT
case with Dirichlet condition at $\mathcal{E}=10^{-3}$ and with $N=200$. The
horizontal dashed lines corresponds to $\alpha_{\epsilon}(L)/\epsilon=1$.

The PT potential provides us with a clean setup for the pseudospectrum
analysis. In the case of dissipative boundary conditions at both asymptotic
regions (${\cal E}=0$), the QNMs are known exactly, and they follow a
distribution in the complex plane very similar to the large frequency limit of
BH QNMs. Additionally, the similarities between PT and BH potentials turn out
to be not limited to their spectrum, as both models present a very similar
pseudospectrum as well [43]. The main reason for the PT potential’s cleaner
numerical behavior is its exponential asymptotic decay, which eliminates the
branch cut present in the BH case. As discussed in ref. [43], this branch cut
manifests numerically as non-convergent eigenvalues on the imaginary axis,
which affect the behavior of the pseudospectrum in their vicinity.

The left panel of fig. 1 shows the spectrum (red dots) and pseudospectrum
(color code) for the PT-ECO with a Dirichlet boundary at
$\mathcal{E}=10^{-3}$. We have used $N=200$ grid points in our discretization
method, and have checked the robustness of the results by varying $N$ to
larger values. The internal precision for all calculations is set to $\sim
200$ digits. For comparison, the figure also displays the original PT QNMs (in
blue dots) obtained with dissipative boundary conditions, given by
$\omega_{n}=\pm\sqrt{3}/2+i\left(n+1/2\right)$. As opposed to the original PT
QNMs, one observes that the fundamental mode of the PT-ECO resides very close
to the real axis. Decreasing $\mathcal{E}$, this mode moves even closer to the
real axis. This is an illustration of stable trapping and the appearance of
long-lived modes, ubiquitously present in ECO configurations.

The proximity of QNMs to the real axis can be quantified in terms of the
pseudospectrum: contour lines of $\epsilon\simeq 0.02$ and larger values cross
this axis into the region of unstable modes (i.e. the negative imaginary half-
plane, in our convention for Fourier decomposition). This is suggestive of a
perturbation-triggered instability. According to eq. (34), one may presume
that adding a minuscule $\|\delta L\|$ to the underlying operator could make
the fundamental mode migrate into the unstable region. However, this migration
is not guaranteed to occur in realistic physical scenarios. To understand
this, take for instance the pseudospectrum of any self-adjoint operator, whose
eigenvalues lie on the real axis. The $\epsilon-$contour lines are circles
around these eigenvalues, and thus they necessarily spread out into parts of
the complex plane with both positive and negative imaginary parts. While this
does suggest that a perturbed version of such an operator could have modes
growing boundlessly, it does not speak to the physicality of such a
perturbation (i.e. this perturbation may need to change the problem
structurally in a way which, for physical systems, does not occur naturally).

The most direct method of checking how this instability manifests itself is by
explicitly adding small perturbations $\delta L$ to the original operator and
calculating the modified spectrum. Ref. [43] explored this possibility by
modifying solely the potential $\hat{V}$ within the operator $L_{1}$ in eq.
(9). In particular, ad hoc modifications of this potential in terms of random,
or oscillatory perturbations were explored. Employing the same strategy here,
we reached the same conclusions: the fundamental mode is not affected by
ultraviolet perturbations to the potential, at least ones which preserve the
asymptotic structure. To actually destabilize the fundamental mode and make it
migrate to the unstable half of the complex plane with this type of small-
scale perturbation turns out to require a more generic modification of the
original operator as a whole. Changes in the fundamental mode are also
expected when infrared effects are introduced, as noted originally in ref.
[54] (see also discussions on infrared modifications in [43, 76, 77]) and
recently systematically addressed in ref. [56]. Analyzing whether any such
perturbations can be related to a physically reasonable modification of our
system is an intriguing question, but it falls beyond the scope of this work.

As discussed, the pseudospectrum also brings information about early time
transients in the dynamical evolution under the unperturbed operator $L$. The
notable protrusion of the pseudospectrum lines past the real axis is very
suggestive for a transient period of growth for solutions to the homogeneous
time-evolution problem (12). The pseudospectral abscissa (42) quantifies the
size for which such protrusion may impact the initial dynamics, in particular
via the ratio $\alpha_{\epsilon}/\epsilon$. The right panel of fig. 1 presents
the dependence of the ratio $\alpha_{\epsilon}/\epsilon$ on the parameter
$\epsilon$. We observe a tendency toward a value $\sim 1$, which suggests a
Kreiss constant $\mathcal{K}(L)=1$ achieved through (III.2.2) at
$\epsilon\to\infty$. This result strongly suggests that the evolution operator
$\|e^{iL\tau}\|$ is in fact contractive, and therefore does not bring about
any transient growth.

On the other hand, such protrusion of the pseudospectrum into the lower half-
plane does indeed indicate is the ease with which a pseudo-resonance can be
triggered if the evolution equation had a time dependent external source,
whose Fourier transform had a (real) frequency close to the pseudospectrum
protrusion region, as discussed in [57, 51]. We will elaborate on the
implications this may have in combination with the non-linearity of
gravitational dynamics in the end of this section. First we proceed to present
the analogous of the above results for perturbations on an ECO spacetime with
an external Schwarzschild geometry.

###  IV.2 Pseudospectrum of ECOs

Figure 2:  ECO’s QNMs (red dots) and pseudospectra for several configurations.
The plots also display the respective BH QNMs (cyan squares), as well as non-
convergent eigenvalues of the operator resulting from a discretised the branch
cut spreding from $\omega=0$. _Top Panel_ : Scalar perturbation for $\ell=2$
with ECO compactness $\mathcal{E}=1$, calculated with $N=200$ grid points. The
contour lines range from $-44$ (uppermost contour) to $-3.5$ (lowest contour)
with steps $4.5$ in the $\log$ scale. _Middle Panel_ : Same as top panel but
with ECO compactness $\mathcal{E}=10^{-3}$. The contour lines range from
$-6.5$ (uppermost contour) to $-1.7$ (lowest contour) with steps $0.6$ in the
$\log$ scale. _Bottom Panel_ : Same as middle panel, but for the axial modes
of the gravitational field. For ECO’s surfaces far from the horizon (top
panel), the ECO QNMs lie far from the real axis. For smaller values of
$\mathcal{E}$ (middle and bottom panels), the fundamental mode becomes long
lived with a very small imaginary part. The pseudospectrum indicates and QNM
instability for all cases.

We now present the results for the pseudospectrum, and quantities derived
therefrom, in the case of ECOs with a reflective surface of compactness
dictated by the parameter $\mathcal{E}$. To ensure the validity of our
results, we have first performed convergence tests for the QNM spectrum with
an increasing number of grid points, and we have found that the value $N=200$
is sufficient for accuracy in the region of the complex plane we analyze (as
long as $\mathcal{E}\gtrsim 10^{-3}$). We have also compared the resulting
QNMs with those obtained by a shooting method and found significant agreement.
The internal precision for all calculations is set to $\sim 300$ digits. As
opposed to the PT potential, the algebraic decay of the Regge-Wheeler
potential towards infinity implies the existence of a branch cut, which in our
notation is set to lie along the positive half of the imaginary axis. With our
discretized operator, this branch cut yields non-converging (in the numerical
resolution $N$) modes spreading out from the origin.

####  IV.2.1 QNM instability

Figure 2 shows several pseudospectra which are representative for the ECO
case, with varying compactness and different spin of the perturbed field. The
top and middle panels compare the results for the $\ell=2$ mode of a scalar
field ($s=0$) with two different compactness parameters ${\cal E}$. For less
compact objects, such as ${\cal E}=1$ (top panel), the ECO QNMs (red dots)
have a spectral abscissa (the imaginary part of the fundamental mode) which is
larger than that of a BH with the same mass (blue squares), though the order
of magnitude stays approximately the same. As the object becomes more compact
($\mathcal{E}=10^{-3}$, middle panel), the spectral abscissa can decrease by
orders of magnitude, i.e. the fundamental QNM of ECOs can become extremely
long-lived (for the $\mathcal{E}=10^{-3}$ case the spectral abscissa becomes
of the order $10^{-4}$ in units of the Schwarzschild radius). The exact same
behavior is found for gravitational ($s=2$) perturbations as well (bottom
panel) [13]. The mechanism behind this is attributed to stable trapping of
perturbations between the surface of the ECO and the photon sphere and leads
to the formation quasibound states of [78, 79, 80, 81, 82, 83, 84] that
repeatedly excite the photon sphere modes, giving rise to echoes in the
ringdown signal [30], as well as an overall slower decay characterized by the
actual ECO’s QNMs.

In all these cases, the opening of the $\epsilon$-pseudospectrum contour lines
across the complex plane is indicative of QNM instability, just as observed in
the BH case [43, 48], as well as in the PT-ECO toy model from the previous
section. The main difference with respect to the BH case is precisely the
closeness of the fundamental mode to the real axis. As was the case in the PT-
ECO model, for sufficiently small values of ${\cal E}$,
$\epsilon-$pseudospectral contour lines for small $\epsilon$ cross the real
axis and protrude into the region of the complex plane where unstable QNMs
would reside. In the examples of fig. 2 the crossing happens for
$\epsilon\gtrsim 10^{-3}$. Moreover, the protrusion is quite visibly not only
for the long-lived fundamental QNM, but also for the first overtone, where it
occurs for contour lines $\epsilon\gtrsim 10^{-2}$.

To better understand this instability, we have once again analyzed
perturbations $\delta L$ (with $\|\delta L\|\simeq\epsilon$) which consist of
adding a random or high-frequency sinusoidal function to the effective
potential in $L_{1}$, analogously to the analysis performed in [43]. The
fundamental QNM of ECOs appears to be stable all the way up to perturbations
of energy norm $\epsilon\gtrsim 10^{-1}$ and higher, suggesting that a QNM
migration to the unstable side of the complex plane is not triggered by
ultraviolet effects. Akin to the discussion in the PT case, a possible
destabilization of the fundamental mode into the unstable region requires a
more generic modification on structure of the discretized $L$ operator, and a
general disturbance of any kind to its position requires at least the
inclusion of infrared effects [43, 56].

When the ECO is less compact and the reflective surface reaches the photon
sphere (e.g. top panel of fig. 2, where $\mathcal{E}=1$), then we observe a
different response to perturbations. Particularly, when the ECO’s fundamental
mode acquires a larger imaginary part, the $\epsilon$-pseudospectrum contour
lines no longer protrude into the unstable part of the complex plane for small
values of $\epsilon$, though their opening in the complex plane still
indicates QNM instability. Interestingly, the fundamental mode for these
configurations are susceptible to destabilization by an ultraviolet
perturbation in the effective potential. This appears to be the case when the
spectral abscissa is larger than that of the Schwarzschild BH (of the same
mass). In other words, the stability of the fundamental mode against these
types of operator perturbations (with energy norm up to $10^{-1}$) is present
only when the spectral abscissa is of the same (or smaller) order of a
Schwarzschild BH.

As evidence for this, fig. 3 shows the spectrum of the $L$ operator with a
perturbation added to the effective potential which maximizes the frequency
that can be resolved with the particular resolution used. The perturbation
physically translates to a set of random numbers added to the potential which
simulate very high frequency sinusoidal perturbations and solely depend on the
energy norm of the perturbation [43]. For the first two cases of fig. 2
($\ell=2$ component of a scalar perturbation on $\mathcal{E}=1$ and
$\mathcal{E}=10^{-3}$ ECOs) we perform a QNM calculation of the perturbed
potential and find the following: in the first case, the fundamental mode is
clearly unstable under a rather small perturbation of energy norm $10^{-2}$,
while in the second case, even perturbations of up to $10^{-1}$ are
insufficient to disturb the fundamental mode and even the first overtone. We
also observe that, much like in the BH case [43], the branches of the
perturbed modes approximately follow the pseudospectral contour lines.

Figure 3: Top two panels of fig. 2 with the addition of the spectrum of a
perturbed version of the effective potential in the $L$ operator. The red,
cyan and yellow points designate the ECO QNMs, Schwarzschild BH QNMs and the
QNMs of the ECO with a perturbed effective potential, respectively. Left
Panel: ECO with compactness $\mathcal{E}=1$, and a random perturbation to the
effective potential of energy norm $10^{-2}$. We observe that the fundamental
mode is spectrally unstable. Right Panel: ECO with compactness
$\mathcal{E}=10^{-3}$, and a random perturbation to the effective potential of
energy norm $10^{-1}$ (the largest reasonable perturbation). We observe that
the fundamental mode and the first overtone are spectrally stable.

We complete the discussion by comparing the ECOs QNMs and pseudospectra for
different angular $\ell$-modes. The key motivation behind this comparison is
the argument that ECO structural instability is triggered in regime of high
angular modes $\ell\gg 1$ [36, 35]. For numerical convergence reasons, we move
up to the case of a slightly less compact ECO, with ${\cal E}=10^{-2}$. The
left panel of fig. 4 shows results for the scalar field with $\ell=2$, where
the result is much the same as for the $\mathcal{E}=10^{-3}$ case, though the
spectral abscissa is slightly larger. However, this is in strong contrast to
the case $\ell=10$ (right panel of fig. 4), where modes become substantially
longer-lived, leading to a plethora of protrusions of pseudospectrum contours
into the unstable regime. In particular, we observe this for the first four
modes (the fundamental mode and first three overtones) in the figure.

Figure 4:  Left Panel: $\ell=2$ scalar QNMs (red dots) and pseudospectra
(white curves) of an ECO with $\mathcal{E}=10^{-2}$, calculated with $N=200$
grid points. The contour lines range from $-21.7$ (uppermost contour) to
$-1.7$ (lowest contour) with steps $2$ in the $\log$ scale. For completeness,
we include the $\ell=2$ scalar QNMs of a Schwarzschild BH (cyan squares).
Right Panel: $\ell=10$ scalar QNMs (red dots) and pseudospectra (white curves)
of an ECO with $\mathcal{E}=10^{-2}$, calculated with $N=200$ grid points. The
contour lines range from $-6.5$ (uppermost contour) to $-1.7$ (lowest contour)
with steps $0.6$ in the $\log$ scale. For completeness, we include the
$\ell=10$ scalar QNMs of a Schwarzschild BH (cyan squares).

To better quantify how these protrusions into the unstable part of the complex
plane affects the field’s dynamics, we now once again proceed to apply the
tools introduced in sec. III.2.

####  IV.2.2 Transient effects

As discussed above and in further detail in [51], the multifaceted nature of
pseudospectrum analysis allows one to extend the study beyond the mere
spectral instability effects, and obtain information on the early (and
intermediate) time evolution governed by (11). Repeating the calculation
performed for the PT example in sec. IV.1, we calculate the ratio
$\alpha_{\epsilon}/\epsilon$ appearing on the right-hand side of eq. (43),
which quantitatively characterizes the effects of the
$\epsilon-$pseudospectral protrusion to the the unstable half of the complex
plane.

As representative examples, fig. 5 shows the results obtained from the cases
with largest protrusion of pseudospectral lines for the scalar field: with
parameters ${\cal E}=10^{-3}$, ${\ell=2}$ (top panel) and ${\cal E}=10^{-2}$,
${\ell=10}$ (bottom panel); the calculation in other cases leads to the same
qualitative result. The plots provide strong indications that the ratio
$\alpha_{\epsilon}/{\epsilon}$ tends asymptotically to unity from below as
$\epsilon\rightarrow\infty$. This result agrees with the one observed in the
PT-ECO toy model, and indicates that the Kreiss constant appears to attain the
value ${\cal K}(L)=1$ through this ratio at $\epsilon\to\infty$, which is
consistent with the time evolution operator $e^{iL\tau}$ being contractive.
This is indeed confirmed by the analytical calculation of the numerical
abscissa presented below. The results for gravitational perturbations are
qualitatively the same.

Figure 5: Lower bound $\alpha_{\epsilon}(L)/\epsilon$ to the evolution
operator norm for different values of $\epsilon$. _Top Panel_ : scalar field
with the $\ell=2$ with ECO compactness $\mathcal{E}=10^{-3}$. _Bottom Panel_ :
scalar field with the $\ell=2$ with ECO compactness $\mathcal{E}=10^{-2}$.
Both cases indicate that $\alpha_{\epsilon}(L)/\epsilon\rightarrow 1$
(horizontal dashed line showing the limit), implying a Kreiss constant
$\mathcal{K}=1$ reached at $\epsilon\to\infty$, which occurs when the
evolution operator is $e^{iL\tau}$ contractive.

It is also worth noting that, as discussed in [51], the pseudospectrum at
large values of $\epsilon$ generally contains information about the early-time
evolution of the system, as can be seen by the fact that the
$\epsilon\to\infty$ limit constrains initial growth. The similarities between
the pseudospectra of BHs and ECOs in this limit (which can be seen by the
analogous behavior for the numerical and pseudospectral abscissa between the
two cases) may well be what encodes the similarity between the early-stage
ringdown of the two objects, which has thus far been a mystery form the point
of view of spectral analyses. This phenomenon is further discussed in point
3.5. of ref. [51].

###  IV.3 Contractive operator and pseudo-resonances

Our numerical results strongly suggests that the ECO evolution operator
$e^{iL\tau}$ is contractive, i.e. $\|e^{iL\tau}\|\leq 1,\forall\tau>0$) (more
specifically, the found unity Kreiss constant ${\cal K}(L)=1$ is a necessarily
condition for contractibility, actually strongly supporting the latter),
regardless the potential dictating the dynamics (PT or Regge-Wheeler with
different values of $\ell$). Here, we prove this result analytically by
evaluating the numerical abscissa and confirming its vanishing $\omega(L)=0$,
that actually implies the contractive character of $L$ [44, 51]. We also argue
that the $\epsilon-$pseudospectral protrusions into the unstable half of the
complex plane are indicative of pseudo-resonances on second order
perturbations, possibly triggered by first-order perturbations and which can
themselves trigger resonant-like instabilities breaking down the perturbative
analysis.

####  IV.3.1 Numerical abscissa: an analytical result

Following the procedure laid out in [51] we use the formal expression (53) of
the boundary operator $L^{\partial}$, which determines the adjoint
$L^{\dagger}$, to analytically calculate the numerical abscissa $\omega(L)$.
In accordance with eq. (41), we obtain $\omega(L)$ from the largest eigenvalue
of the operator

$\displaystyle\dfrac{1}{2i}\left(L^{\dagger}-L\right)$ $\displaystyle=$
$\displaystyle\frac{L^{\partial}}{2i}$ (54) $\displaystyle=$
$\displaystyle\begin{pmatrix}0&0\\\
0&-\gamma\delta(\chi-\chi_{{}_{{I}^{+}}})\end{pmatrix}.$

From eqs. (19) and (29), we obtain $\gamma(\chi_{{}_{{I}^{+}}})=1$ in both the
scenarios modelled by the PT and Regge-Wheeler potentials. In fact, we expect
this result to be valid for any potential, since $\gamma(\chi_{{}_{{I}^{+}}})$
is determined directly from its asymptotic decay, i.e. the asymptotic flatness
of the spacetime in the ECO case.

The eigenvalues of this diagonal operator (54) are given by its diagonal
entries, which either vanish or are negative. Since the numerical abscissa
corresponds to the maximum eigenvalue, we then conclude

$\omega(L)=0.$ (55)

This result also implies ${\cal K}(L)=1$, as commented above, and as confirmed
numerically. It allows us to conclude that the evolution operator is indeed
contractive and exhibits no transient growth, in complete analogy to the BH
case [51].

We emphasize that this result is a direct consequence of the fully reflective
boundary condition imposed at the ECO surface via eq. (14), which eliminates
part of the boundary terms which would appear in (53). Whether this is still
the case for more generic conditions at the ECO surface is an open question.

####  IV.3.2 Pseudo-resonances and non-linear bootstrap instability

The pseudospectrum is a tool which describes the evolution of perturbations on
BH and ECO spacetimes in more detail than a simple modal analysis can, as a
consequence of the non-selfadjoint character of the time evolution generator
$L$. However, the evolution it pertains to is still only at the linear level.
Even though in hydrodynamics the presence of linear-level transient effects
have been used to infer the triggering of non-linear instabilities [57], in
the ECO systems we analyzed here no such transient effects have been observed.
This analysis therefore appears insufficient to capture the expected non-
linear dynamical instabilities discussed in refs. [35, 36].

That said, the pseudospectrum of these systems turns out to be useful even
beyond linear order. In particular, the transitive protrusion of
pseudospectral lines through the real axis, along with the notion of pseudo-
resonance ultimately associated with the third pseudospectrum definition (36),
can be combined with GR perturbation theory at second order [85, 60, 86] to
give rise to a potential instability mechanism. Such a mechanism is described
in appendix A of ref. [51]. Although that mechanism does not prove useful in
the context of [51], it seems perfectly suited to the current ultra-compact
object setting. We follow appendix A in [51] in the discussion below.

Due to the non-linearity of GR, the dynamics of a perturbation at second order
$u^{(2)}$ is always sourced by the first order field $u^{(1)}$ and, in
particular gauges, (the homogeneous part of) the second order evolution
operator is the same as the first order one. Thus, we can extend the linear
system eq. (11) to second order as [51]

$\displaystyle\left(\partial_{\tau}-iL\right)u^{(1)}$ $\displaystyle=$
$\displaystyle 0$ (56) $\displaystyle\left(\partial_{\tau}-iL\right)u^{(2)}$
$\displaystyle=$ $\displaystyle S[u^{(1)}].$ (57)

Because the field $u^{(1)}$ does not undergo any transient growth in its
evolution, we follow [57, 68, 51] to model the source as $S\sim
e^{i\omega\tau}s(\chi)$. Then, a formal solution to eq. (57) follows by acting
with the Green function on the source term via

$\displaystyle u^{(2)}(\tau,\chi)$ $\displaystyle=$ $\displaystyle
ie^{i\omega\tau}\left(L-\omega{\mathbb{I}}\right)^{-1}s(\chi)$ (58)
$\displaystyle=$ $\displaystyle ie^{i\omega\tau}R_{L}(\omega)s(\chi).$

Note that in the second line we have used the fact that the resolvent
$R_{L}(\omega)$ defined in eq. (33) is precisely the Green function of the
system. Thus, an upper bound estimate for the response in the solution
$u^{(2)}$ to the oscillating source term reads

$\displaystyle\|u^{(2)}\|\lesssim e^{-{\rm
Im}(\omega)\tau}\|R_{L}(\omega)\|\|s_{n}\|.$ (59)

Remarkably, as a consequence of the perturbation theory structure in GR (57),
the resolvent norm $\|R_{L}(\omega_{n})\|$ which appears here is precisely the
one used to define the pseudospectrum of the linear operator at first order,
precisely the one we have studied in detail in previous sections.

For standard resonances, the field’s maximal response will be at the QNM
spectral frequencies $\omega_{n}$, since
$\|R_{L}(\omega_{n})\|\rightarrow\infty$. However, in the non-selfadjoint case
and due to the spread of $\epsilon$-pseudospectra sets in the complex plane,
we may also encounter resonant-like behaviour, namely pseudo-resonances, for
which the resolvent’s norm is still large, though not strictly divergent. Such
pseudo-resonances are present in every pseudospectrum of a non-selfadjoint
operator. But what makes the pseudospectrum of ECOs stand out is the fact that
for very compact objects ($\mathcal{E}\ll 1$) there exist strong protrusions
in the unstable part of the complex plane crossing the real line. As a
consequence, the norm of the resolvent $R_{L}(\omega)$ becomes very large for
certain values along the real line. This is the key feature for the pseudo-
resonant mechanism. Specifically, considering a time-depending external source
to $u^{(2)}$ (here provided by $u^{(1)}$), through its Fourier transform, it
can be seen as a superposition of real frequency harmonic contributions. As a
consequence, if the Fourier-transformed external source is peaked at the
(real) frequencies where the pseudospectrum protrudes into the unstable
complex half-plane, a pseudo-resonance occurs. Quantitatively, the maximal
response growth-factor $R_{\mathrm{max}}$ that can be attained [51] is given
by

$\begin{split}R_{\mathrm{max}}=\sup_{\omega\in\mathbb{R}}R_{\mathrm{max}}(\omega)&=\sup_{\omega\in\mathbb{R}}\|(\omega-L)^{-1}\|\\\
&=\sup_{\omega\in\mathbb{R}}\|R_{L}(\omega)\|.\end{split}$ (60)

This pseudo-resonance mechanism may well lead to a growth of $u^{(2)}$,
sourced by $u^{(1)}$, which breaks down the perturbative expansion and leads
to the non-linear dynamical instabilities referred to in [35, 36, 37]. This is
in analogy to what occurs in hydrodynamics [57], though without the transient
growth: whereas in the hydrodynamical case it is the non-selfadjoint linear
transient growth that triggers the non-linear response in a bootstrap
mechanism, in the gravitational case the bootstrap aspect is provided by the
GR perturbative structure that source gravitational second-order perturbations
with gravitational first-order ones.

##  V Discussion

The pseudospectrum analysis offers a profound insight into the evolution
dictated by linear operators and it has been widely applied in physical
systems [57, 75, 44, 68]. Its ability to capture not just the poles, but the
whole ‘topographical structure’ of the Green’s function provides a powerful
tool to assess underlying nonmodal effects, such as spectral instabilities,
isospectrality properties, early-time transient phenomena, turbulence and
pseudo-resonances — features that an eigenvalue analysis of the system’s
operator does not account for. Only recently has the pseudospectrum analysis
been introduced to the BH literature [43, 50, 48, 49, 51]. The pseudospectrum
tool was first used to assess BH QNM spectral instabilities in Schwarzschild
[43, 50] and Reissner-Nordström [48] BHs. Such instabilities may offer new
challenges to GW data analysis [50, 49]. Going beyond the study of spectral
instabilities, ref. [51] has shown the wide variety of applications of the
pseudospectrum analysis in gravitational settings, such as the study of binary
BH systems by means of an effective (non-selfadjoint) linear description.

This work studies horizonless ECOs with the tools provided by pseudospectrum
analysis. Our goal was twofold. First, expanding on the results from ref. [43,
50, 48, 49] and searching for signatures of QNMs instabilities. Second,
applying the formalism laid out in ref. [51] to assess possible initial
transients growths that could trigger the non-linear dynamical instability
identified in [35, 36].

To this purpose, we have followed ref. [43] and we initially considered a toy-
model built from the PT potential. Despite its simplicity, this model actually
shares a lot of the qualitative features of an ECO case, while also providing
a useful testing ground for our numerical method.

We found that ECOs with any compactness suffer from the same ultraviolet
spectral instability in the overtones as BHs do. Additionally, we observed new
aspects concerning the stability of the fundamental mode under ultraviolet
modifications in the potential.

A well-known phenomenon is that the ECO QNMs $\omega_{n}^{\rm ECO}$ differ
significantly from those of a BH with the same mass $\omega_{n}^{\rm BH}$. For
instance, the fundamental mode $\omega^{\rm ECO}_{0}$ for larger, less compact
ECOs may acquire values with an imaginary part (spectral abscissa) larger than
in the BH case. As the ECO becomes more compact, $\omega^{\rm ECO}_{0}$, as
well as the overtones, get closer to the real axis, becoming long-lived as a
consequence of the trapping which occurs between the ECO surface and the light
ring potential barrier [13, 30]. We observe that the ECO fundamental mode is
unstable under ultraviolet modifications of the potential only when ${\rm
Im}(\omega^{\rm ECO}_{0})>{\rm Im}(\omega^{\rm BH}_{0})$. As $\omega^{\rm
ECO}_{0}$ goes down towards the real axis, it becomes ultraviolet stable
approximately when ${\rm Im}(\omega^{\rm ECO}_{0})\simeq{\rm Im}(\omega^{\rm
BH}_{0})$. From there on, destabilizing the fundamental mode requires
modifications in the whole operator, or the introduction of infrared effects
in the potential [43, 56].

This result is of particular importance for ultracompact ECOs, whose
fundamental mode lies very close to the real axis. For these objects, the
$\epsilon$-pseudospectrum contour-lines cross the real axis into the unstable
region of the complex plane for relatively small values of $\epsilon$. In
principle, the fundamental mode of a slightly perturbed version of the
operator could then cross into the unstable half of the complex plane.
However, this is not observed for any perturbation to the effective potential.
The perturbing operator which could cause this migration would therefore
likely have to substantially change the structure of the underlying problem
(e.g. introducing new derivative terms, albeit with small multiplicative
coefficients), which is not directly relatable to physical sources of
perturbation.

Even though the fundamental mode is not straightforwardly destabilized, the
$\epsilon$-pseudospectral protrusion into the unstable region of the complex
plane might contribute to an initial transient growth, as observed in in
hydrodynamics [57]. However, calculating the Kreiss constant and numerical
abscissa [51] we find an absence of such initial growth. More specifically, we
conclude that the evolution operator $e^{iL\tau}$ is contractive, i.e. that
the solutions to (12) can only have a constant or decreasing norm. This
result, in particular the analytical calculation of the numerical abscissa
$\omega(L)=0$, turns out to strongly depend on the field boundary conditions
we have chosen for the surface of the ECO (14). An interesting open question
is what occurs with different boundary conditions.

However, an initial transient growth is not the only mechanism able to trigger
a non-linear dynamical instability. The $\epsilon$-pseudospectrum protrusion
into the unstable region of the complex plane shows the existence of pseudo-
resonances at and below the real axis which, when appropriately excited, can
lead to a dynamical instability. In particular, due to the fact that the Green
function for the second-order perturbation field is precisely the resolvent
$R_{L}(\lambda)$ of the first-order operator (from which the pseudospectrum is
obtained) [51], pseudo-resonant excitations may be observed at second order,
with the first-order solution acting as a source. The second-order
contribution to perturbations has been shown to be non-negligible (and
potentially measurable) even in the case of a BH [87], where long-lived modes
and pseudo-resonances are generally absent. Here, we can estimate the amount
by which second-order effects are amplified just through the pseudospectrum we
have calculated. Our conclusion is that for ultracompact horizonless objects,
for which the fundamental modes and first overtones lie close to the real
axis, second-order perturbations will be amplified significantly—potentially
enough to lead to the breakdown of linear perturbation theory, and
consequently to the destabilization of the astrophysical object itself.

Overall, the pseudospectrum framework and non-modal analysis provide a large
range of tools to study wave phenomena at linear level for non-conservative
systems driven by non-selfadjoint time evolution operators, beyond standard
(self-adjoint) modal analysis [51]. Potential future works in this direction
include introducing this methodology to a wider range of BH and ECO setups,
with different asymptotic structures and symmetries (see e.g. [32, 33, 82, 88,
89, 90, 91]).

###### Acknowledgements.

The authors thank Emanuele Berti for helpful discussions. V.B. is funded by
the Spanish Government fellowship FPU17/04471. R.P.M acknowledges financial
support provided by COST Action CA16104 via the Short Term Scientific Mission
grant and from STFC via grant number ST/V000551/1. V.C. is a Villum
Investigator and a DNRF Chair, supported by VILLUM FONDEN (grant no. 37766)
and by the Danish Research Foundation. V. B. and V.C. acknowledge financial
support provided under the European Union’s H2020 ERC Advanced Grant “Black
holes: gravitational engines of discovery” grant agreement no.
Gravitas–101052587. 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 acknowledge financial
support provided by FCT/Portugal through grants 2022.01324.PTDC, PTDC/FIS-
AST/7002/2020, UIDB/00099/2020 and UIDB/04459/2020. This research project was
conducted using the computational resources of “Baltasar Sete-Sois” cluster at
Instituto Superior Técnico. This work was supported by the French
“Investissements d’Avenir” program through project ISITE-BFC (ANR-15-IDEX-03),
the ANR “Quantum Fields interacting with Geometry” (QFG) project
(ANR-20-CE40-0018-02), the EIPHI Graduate School (ANR-17-EURE-0002) and the
Spanish FIS2017-86497-C2-1 project (with FEDER contribution).

## References

  * Abbott  _et al._ [2021] R. Abbott _et al._ (LIGO Scientific, Virgo), GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run, Phys. Rev. X 11, 021053 (2021), arXiv:2010.14527 [gr-qc] .
  * Barack  _et al._ [2019] L. Barack _et al._ , Black holes, gravitational waves and fundamental physics: a roadmap, Class. Quant. Grav. 36, 143001 (2019), arXiv:1806.05195 [gr-qc] .
  * Barausse  _et al._ [2020] E. Barausse _et al._ , Prospects for Fundamental Physics with LISA, Gen. Rel. Grav. 52, 81 (2020), arXiv:2001.09793 [gr-qc] .
  * Amaro-Seoane  _et al._ [2022] P. Amaro-Seoane _et al._ , Astrophysics with the Laser Interferometer Space Antenna, (2022), arXiv:2203.06016 [gr-qc] .
  * Doeleman  _et al._ [2008] S. Doeleman _et al._ , Event-horizon-scale structure in the supermassive black hole candidate at the Galactic Centre, Nature 455, 78 (2008), arXiv:0809.2442 [astro-ph] .
  * Cardoso and Gualtieri [2016] V. Cardoso and L. Gualtieri, Testing the black hole ‘no-hair’ hypothesis, Class. Quant. Grav. 33, 174001 (2016), arXiv:1607.03133 [gr-qc] .
  * Macedo  _et al._ [2013a] C. F. B. Macedo, P. Pani, V. Cardoso, and L. C. B. Crispino, Into the lair: gravitational-wave signatures of dark matter, Astrophys. J. 774, 48 (2013a), arXiv:1302.2646 [gr-qc] .
  * Annulli  _et al._ [2020a] L. Annulli, V. Cardoso, and R. Vicente, Stirred and shaken: Dynamical behavior of boson stars and dark matter cores, Phys. Lett. B 811, 135944 (2020a), arXiv:2007.03700 [astro-ph.HE] .
  * Annulli  _et al._ [2020b] L. Annulli, V. Cardoso, and R. Vicente, Response of ultralight dark matter to supermassive black holes and binaries, Phys. Rev. D 102, 063022 (2020b), arXiv:2009.00012 [gr-qc] .
  * Vicente and Cardoso [2022] R. Vicente and V. Cardoso, Dynamical friction of black holes in ultralight dark matter, Phys. Rev. D 105, 083008 (2022), arXiv:2201.08854 [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._ [2022] 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 (2022), arXiv:2109.00005 [gr-qc] .
  * Cardoso  _et al._ [2016a] V. Cardoso, E. Franzin, and P. Pani, Is the gravitational-wave ringdown a probe of the event horizon?, Phys. Rev. Lett. 116, 171101 (2016a), [Erratum: Phys.Rev.Lett. 117, 089902 (2016)], arXiv:1602.07309 [gr-qc] .
  * Conklin and Afshordi [2021] R. S. Conklin and N. Afshordi, Boltzmann Meets Lorentz: A Surrogate Model for Black Hole Echoes, (2021), arXiv:2201.00027 [gr-qc] .
  * Macedo  _et al._ [2013b] C. F. B. Macedo, P. Pani, V. Cardoso, and L. C. B. Crispino, Astrophysical signatures of boson stars: quasinormal modes and inspiral resonances, Phys. Rev. D 88, 064046 (2013b), arXiv:1307.4812 [gr-qc] .
  * Giudice  _et al._ [2016] G. F. Giudice, M. McCullough, and A. Urbano, Hunting for Dark Particles with Gravitational Waves, JCAP 10, 001, arXiv:1605.01209 [hep-ph] .
  * Cardoso and Pani [2019] V. Cardoso and P. Pani, Testing the nature of dark compact objects: a status report, Living Rev. Rel. 22, 4 (2019), arXiv:1904.05363 [gr-qc] .
  * Liebling and Palenzuela [2012] S. L. Liebling and C. Palenzuela, Dynamical Boson Stars, Living Rev. Rel. 15, 6 (2012), arXiv:1202.5809 [gr-qc] .
  * Mazur and Mottola [2001] P. O. Mazur and E. Mottola, Gravitational condensate stars: An alternative to black holes, (2001), arXiv:gr-qc/0109035 .
  * Carballo-Rubio [2018] R. Carballo-Rubio, Stellar equilibrium in semiclassical gravity, Phys. Rev. Lett. 120, 061102 (2018), arXiv:1706.05379 [gr-qc] .
  * Arrechea  _et al._ [2021] J. Arrechea, C. Barceló, R. Carballo-Rubio, and L. J. Garay, Semiclassical constant-density spheres in a regularized Polyakov approximation, Phys. Rev. D 104, 084071 (2021), arXiv:2105.11261 [gr-qc] .
  * Kaplan and Rajendran [2019] D. E. Kaplan and S. Rajendran, Firewalls in General Relativity, Phys. Rev. D 99, 044033 (2019), arXiv:1812.00536 [hep-th] .
  * Mathur [2005] S. D. Mathur, The Fuzzball proposal for black holes: An Elementary review, Fortsch. Phys. 53, 793 (2005), arXiv:hep-th/0502050 .
  * Mathur [2008] S. D. Mathur, Fuzzballs and the information paradox: A Summary and conjectures, (2008), arXiv:0810.4525 [hep-th] .
  * Ikeda  _et al._ [2021] T. Ikeda, M. Bianchi, D. Consoli, A. Grillo, J. F. Morales, P. Pani, and G. Raposo, Black-hole microstate spectroscopy: Ringdown, quasinormal modes, and echoes, Phys. Rev. D 104, 066021 (2021), arXiv:2103.10960 [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 .
  * 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] .
  * Zworski [2017] M. Zworski, Mathematical study of scattering resonances, Bulletin of Mathematical Sciences 7, 1 (2017).
  * Dyatlov and Zworski [2019] S. Dyatlov and M. Zworski,  _Mathematical Theory of Scattering Resonances_, Graduate Studies in Mathematics (American Mathematical Society, 2019).
  * Cardoso  _et al._ [2016b] V. Cardoso, S. Hopper, C. F. B. Macedo, C. Palenzuela, and P. Pani, Gravitational-wave signatures of exotic compact objects and of quantum corrections at the horizon scale, Phys. Rev. D 94, 084031 (2016b), arXiv:1608.08637 [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] .
  * Cardoso  _et al._ [2018a] V. Cardoso, J. a. 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] .
  * 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 [2019a] K. Destounis, Charged Fermions and Strong Cosmic Censorship, Phys. Lett. B 795, 211 (2019a), arXiv:1811.10629 [gr-qc] .
  * Cardoso  _et al._ [2014] V. Cardoso, L. C. B. Crispino, C. F. B. Macedo, H. Okawa, and P. Pani, Light rings as observational evidence for event horizons: long-lived modes, ergoregions and nonlinear instabilities of ultracompact objects, Phys. Rev. D 90, 044069 (2014), arXiv:1406.5510 [gr-qc] .
  * Keir [2016] J. Keir, Slowly decaying waves on spherically symmetric spacetimes and ultracompact neutron stars, Class. Quant. Grav. 33, 135009 (2016), arXiv:1404.7036 [gr-qc] .
  * Cunha  _et al._ [2022] P. V. P. Cunha, C. Herdeiro, E. Radu, and N. Sanchis-Gual, The fate of the light-ring instability, (2022), arXiv:2207.13713 [gr-qc] .
  * Friedman [1978] J. L. Friedman, Ergoregion instability, Communications in Mathematical Physics v.63, 243–255 (1978) 63 (1978).
  * Schutz and Comins [1978] B. Schutz and N. Comins, On the ergoregion instability, Proceedings of the Royal Society of London A, v.364, 211-266 (1978) 364 (1978).
  * Cardoso  _et al._ [2008] V. Cardoso, P. Pani, M. Cadoni, and M. Cavaglia, Ergoregion instability of ultracompact astrophysical objects, Phys. Rev. D 77, 124044 (2008), arXiv:0709.0532 [gr-qc] .
  * Chirenti and Rezzolla [2008] C. B. M. H. Chirenti and L. Rezzolla, On the ergoregion instability in rotating gravastars, Phys. Rev. D 78, 084011 (2008), arXiv:0808.4080 [gr-qc] .
  * Moschidis [2018] G. Moschidis, A Proof of Friedman’s Ergosphere Instability for Scalar Waves, Commun. Math. Phys. 358, 437 (2018), arXiv:1608.02035 [math.AP] .
  * 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).
  * Davies [2007] E. Davies,  _Linear Operators and their Spectra_, Cambridge Studies in Advanced Mathematics (Cambridge University Press, 2007).
  * Sjöstrand [2019] J. Sjöstrand,  _Non-Self-Adjoint Differential Operators, Spectral Asymptotics and Random Perturbations_, Pseudo-Differential Operators (Springer International Publishing, 2019).
  * Ashida  _et al._ [2021] Y. Ashida, Z. Gong, and M. Ueda, Non-Hermitian physics, Adv. Phys. 69, 249 (2021), arXiv:2006.01837 [cond-mat.mes-hall] .
  * 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] .
  * 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] .
  * Jaramillo [2022] J. L. Jaramillo, Pseudospectrum and binary black hole merger transients, (2022), arXiv:2206.08025 [gr-qc] .
  * Aguirregabiria and Vishveshwara [1996] J. M. Aguirregabiria and C. V. Vishveshwara, Scattering by black holes: A Simulated potential approach, Phys. Lett. A 210, 251 (1996).
  * Vishveshwara [1996] C. V. Vishveshwara, On the black hole trail …, in  _18th Conference of the Indian Association for General Relativity and Gravitation_ (Institute of Mathematical Science Report, by Madras Univ. Inst. Math. Sci., Madras, India, 1996) pp. 11–22.
  * 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 .
  * Cheung  _et al._ [2022a] 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 (2022a), arXiv:2111.05415 [gr-qc] .
  * Trefethen  _et al._ [1993] L. N. Trefethen, A. E. Trefethen, S. C. Reddy, and T. A. Driscoll, Hydrodynamic stability without eigenvalues, Science 261, 578 (1993).
  * Schmid [2007a] P. J. Schmid, Nonmodal stability theory, Annual Review of Fluid Mechanics 39, 129 (2007a), https://doi.org/10.1146/annurev.fluid.38.050304.092139 .
  * Cunningham  _et al._ [1980a] C. T. Cunningham, R. H. Price, and V. Moncrief, Radiation from collapsing relativistic stars. III - Second order perturbations of collapse with rotation, Astrophysical Journal 236, 674 (1980a).
  * Gleiser  _et al._ [1996] R. J. Gleiser, C. O. Nicasio, R. H. Price, and J. Pullin, Colliding black holes: How far can the close approximation go?, Phys. Rev. Lett. 77, 4483 (1996), arXiv:gr-qc/9609022 .
  * Maggio  _et al._ [2020] E. Maggio, L. Buoninfante, A. Mazumdar, and P. Pani, How does a dark compact object ringdown?, Phys. Rev. D 102, 064053 (2020), arXiv:2006.14628 [gr-qc] .
  * Zenginoglu [2008] A. Zenginoglu, Hyperboloidal foliations and scri-fixing, Class. Quant. Grav. 25, 145002 (2008), arXiv:0712.4333 [gr-qc] .
  * Zenginoglu [2011] A. Zenginoglu, A Geometric framework for black hole perturbations, Phys. Rev. D83, 127502 (2011), arXiv:1102.2451 [gr-qc] .
  * Zenginoğlu [2021] A. Zenginoğlu, A null infinity layer for wave scattering, (2021), arXiv:2111.14217 [math.NA] .
  * 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).
  * Schmid [2007b] P. J. Schmid, Nonmodal Stability Theory, Annual Review of Fluid Mechanics 39, 129 (2007b).
  * Trefethen [2000] L. N. Trefethen,  _Spectral Methods in Matlab_ (Society for Industrial and Applied Mathematics, Philadelphia, PA, 2000).
  * Boyd [2001] J. P. Boyd,  _Chebyshev and Fourier spectral methods_ (Courier Corporation, 2001).
  * Canuto  _et al._ [2007] C. Canuto, M. Hussaini, A. Quarteroni, and T. Zang,  _Spectral Methods: Fundamentals in Single Domains_, Scientific Computation (Springer Berlin Heidelberg, 2007).
  * Jackiewicz and Renaut [2002] Z. Jackiewicz and R. A. Renaut, A note on stability of pseudospectral methods for wave propagation, J. Comput. Appl. Math. 143, 127–139 (2002).
  * 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 .
  * Leaver [1986] E. W. Leaver, Spectral decomposition of the perturbation response of the Schwarzschild geometry, Phys. Rev. D 34, 384 (1986).
  * Driscoll and Trefethen [1996] T. A. Driscoll and L. N. Trefethen, Pseudospectra for the wave equation with an absorbing boundary, Journal of Computational and Applied Mathematics 69, 125 (1996).
  * Qian  _et al._ [2021] W.-L. Qian, K. Lin, C.-Y. Shao, B. Wang, and R.-H. Yue, Asymptotical quasinormal mode spectrum for piecewise approximate effective potential, Phys. Rev. D 103, 024019 (2021), arXiv:2009.11627 [gr-qc] .
  * Liu  _et al._ [2021] H. Liu, W.-L. Qian, Y. Liu, J.-P. Wu, B. Wang, and R.-H. Yue, Alternative mechanism for black hole echoes, Phys. Rev. D 104, 044012 (2021), arXiv:2104.11912 [gr-qc] .
  * Dolan [2007] S. R. Dolan, Instability of the massive Klein-Gordon field on the Kerr spacetime, Phys. Rev. D 76, 084001 (2007), arXiv:0705.2880 [gr-qc] .
  * Rosa and Dolan [2012] J. G. Rosa and S. R. Dolan, Massive vector fields on the Schwarzschild spacetime: quasi-normal modes and bound states, Phys. Rev. D 85, 044043 (2012), arXiv:1110.4494 [hep-th] .
  * Cardoso  _et al._ [2011] V. Cardoso, S. Chakrabarti, P. Pani, E. Berti, and L. Gualtieri, Floating and sinking: The Imprint of massive scalars around rotating black holes, Phys. Rev. Lett. 107, 241101 (2011), arXiv:1109.6021 [gr-qc] .
  * Pani  _et al._ [2012] P. Pani, V. Cardoso, L. Gualtieri, E. Berti, and A. Ishibashi, Perturbations of slowly rotating black holes: massive vector fields in the Kerr metric, Phys. Rev. D 86, 104017 (2012), arXiv:1209.0773 [gr-qc] .
  * Destounis [2019b] K. Destounis, Superradiant instability of charged scalar fields in higher-dimensional Reissner-Nordström-de Sitter black holes, Phys. Rev. D 100, 044054 (2019b), arXiv:1908.06117 [gr-qc] .
  * Vieira and Kokkotas [2021] H. S. Vieira and K. D. Kokkotas, Quasibound states of Schwarzschild acoustic black holes, Phys. Rev. D 104, 024035 (2021), arXiv:2104.03938 [gr-qc] .
  * Vieira  _et al._ [2022] H. S. Vieira, K. Destounis, and K. D. Kokkotas, Slowly-rotating curved acoustic black holes: Quasinormal modes, Hawking-Unruh radiation, and quasibound states, Phys. Rev. D 105, 045015 (2022), arXiv:2112.08711 [gr-qc] .
  * Cunningham  _et al._ [1980b] C. T. Cunningham, R. H. Price, and V. Moncrief, Radiation from collapsing relativistic stars. III - Second order perturbations of collapse with rotation, Astrophys. J. 236, 674 (1980b).
  * Miller  _et al._ [2016] J. Miller, B. Wardell, and A. Pound, Second-order perturbation theory: the problem of infinite mode coupling, Phys. Rev. D 94, 104018 (2016), arXiv:1608.06783 [gr-qc] .
  * Cheung  _et al._ [2022b] M. H.-Y. Cheung _et al._ , Nonlinear effects in black hole ringdown, (2022b), arXiv:2208.07374 [gr-qc] .
  * Destounis  _et al._ [2020a] K. Destounis, R. D. B. Fontana, and F. C. Mena, Accelerating black holes: quasinormal modes and late-time tails, Phys. Rev. D 102, 044005 (2020a), arXiv:2005.03028 [gr-qc] .
  * Destounis  _et al._ [2020b] K. Destounis, R. D. B. Fontana, and F. C. Mena, Stability of the Cauchy horizon in accelerating black-hole spacetimes, Phys. Rev. D 102, 104037 (2020b), arXiv:2006.01152 [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] .
  * Destounis  _et al._ [2022] K. Destounis, G. Mascher, and K. D. Kokkotas, Dynamical behavior of the $C$-metric: Charged scalar fields, quasinormal modes, and superradiance, Phys. Rev. D 105, 124058 (2022), arXiv:2206.07794 [gr-qc] .