# The pseudospectrum and transient of Kaluza-Klein black holes in Einstein-
Gauss-Bonnet gravity

Jia-Ning Chen1,2,3111e-mail address<EMAIL_ADDRESS>Liang-Bi
Wu1,2222e-mail address<EMAIL_ADDRESS>Zong-Kuan Guo3,2,1333e-mail
address<EMAIL_ADDRESS>1School of Fundamental Physics and Mathematical
Sciences, Hangzhou Institute for Advanced Study, University of Chinese Academy
of Sciences, Hangzhou 310024, China  2School of Physical Sciences, University
of Chinese Academy of Sciences, No.19A Yuquan Road, Beijing 100049, China
3CAS Key Laboratory of Theoretical Physics, Institute of Theoretical Physics,
Chinese Academy of Sciences, Beijing 100190, China

(August 28, 2024)

###### Abstract

The spectrum and dynamical instability, as well as the transient effect of the
tensor perturbation for the so-called Maeda-Dadhich black hole, a type of
Kaluza-Klein black hole, in Einstein-Gauss-Bonnet gravity have been
investigated in framework of pseudospectrum. We cast the problem of solving
quasinormal modes (QNMs) in AdS-like spacetime as the linear evolution problem
of the non-normal operator in null slicing by using ingoing Eddington-
Finkelstein coordinates. In terms of spectrum instability, based on the
generalised eigenvalue problem, the QNM spectrum and $\epsilon$-pseudospectrum
has been studied, while the open structure of $\epsilon$-pseudospectrum caused
by the non-normality of operator indicates the spectrum instability. In terms
of dynamical instability, the concept of the distance to dynamical instability
has been firstly introduced in gravity aspects, which plays a crucial role in
bridging the spectrum instability and the dynamical instability. We calculate
the distance, named the complex stability radius, as parameters vary. Finally,
we report the transient effect in the energy norm of the evolution operator,
whose behaviour can be roughly reflected by the three kinds of abscissas in
context of pseudospectrum, and explain the relationship between the waveform
and the transient effect.

##  I Introduction

The black hole spectroscopy program has made significant progress by utilizing
gravitational waves (GWs) generated from the merger of binary systems. These
GW signals provide valuable information about the characteristics of the
merging systems Abbott  _et al._ (2016a, b, 2017a, 2017b, 2023). This
achievement has sparked the development of various space-based initiatives,
including LISA Amaro-Seoane  _et al._ (2017), TianQin Luo  _et al._ (2016);
Gong  _et al._ (2021), and Taiji Hu and Wu (2017). During the ringdown phase,
the remnant black hole emits a distinct signal as it dissipates energy and
settles into its equilibrium state. This signal is characterized by a discrete
collection of damped harmonic oscillations known as quasinormal modes (QNMs)
Kokkotas and Schmidt (1999); Berti  _et al._ (2009); Konoplya and Zhidenko
(2011); Franchini and Völkel (2023).

With improvement of the accuracy of gravitational waves detection, research
has focused on the QNM spectra of black holes in the presence of surrounding
matter, since black holes are not isolated but are always accompanied by an
astrophysical environment Barausse  _et al._ (2014); Cannizzaro  _et al._
(2024). The characteristics of QNMs are closely tied to the concept of
spectrum instability. Initial studies by Nollert and Price Nollert (1996);
Nollert and Price (1999) have demonstrated that even small perturbations can
have a significant impact on the QNM spectrum, leading to what is known as
spectrum instability. This discovery highlights the unexpected sensitivity of
the QNM spectrum to “ultraviole” or small-scale perturbations, challenging the
assumption that a good approximation of the effective potential would yield
minimal deviations in the observed QNMs.

After extensive surveys, we have identified two main categories of methods
used to study spectrum instability. However, it is important to note that
there can be some overlap or connections between these two approaches. The
first approach involves modifying the effective potential based on various
reasons, which can be either physically motivated or artificially constructed.
For instance, studies Daghigh  _et al._ (2020); Qian  _et al._ (2021) have
argued that even minor modifications, such as using a piecewise approximate
potential, can lead to non-perturbative changes in the asymptotic behavior of
the QNM spectrum. Authors in Hirano  _et al._ (2024) utilized the
parameterized black hole quasinormal ringdown formalism and observed larger
deviations from general relativity in the quasinormal frequencies for higher
overtones.

The addition of small wave packets to the existing effective potential is a
commonly encountered method in the study of spectrum instability Cardoso  _et
al._ (2024); Courty  _et al._ (2023); Berti  _et al._ (2022); Cheung  _et al._
(2022); Rosato  _et al._ (2024); Oshita  _et al._ (2024). Notably, studies
Rosato  _et al._ (2024); Oshita  _et al._ (2024) focus on the (in)stability of
black hole greybody factors and find that these factors are stable observables
for detecting gravitational waves, meanwhile Regge Poles are used to construct
stable observables in Torres (2023) considering the unstable nature of the QNM
spectrum. In Cardoso  _et al._ (2024), it is argued that the addition of a
small bump $\epsilon V_{\text{bump}}$ cannot be considered as a small
perturbation in any sense. According to their perspective, the introduction of
any arbitrarily small bump leads to the emergence of a new length scale, which
explains the instability of the fundamental mode. Additionally, the previous
work approaches the instability of spectra from the perspective of equal
perturbation energy Cao  _et al._ (2024). More recently, modified PT effective
potentials have been considered to study spectrum instability Li  _et al._
(2024).

Pseudospectrum analysis Trefethen and Embree (2005), as a second approach, is
utilized to investigate the instability of the QNMs spectrum Jaramillo  _et
al._ (2021); Destounis and Duque (2023); Jaramillo  _et al._ (2022), which
originally comes from the field of hydrodynamics Trefethen  _et al._ (1993).
This method involves analyzing the properties of non-self-adjoint operators in
dissipative systems and provides a visual understanding of the instability in
the spectrum of such operators Trefethen and Embree (2005). In the context of
black holes, pseudospectra have been employed to identify qualitative
characteristics that serve as indicators of spectral instability across a
range of spacetimes. These investigations not only encompass the findings of
Jaramillo  _et al._ (2021), but also extend them to other types of black
holes, including asymptotically flat black holes Destounis  _et al._ (2021);
Cao  _et al._ (2024), asymptotically AdS black holes Areán  _et al._ (2023);
Cownden  _et al._ (2024); Boyanov  _et al._ (2024), asymptotically dS black
holes Sarkar  _et al._ (2023); Destounis  _et al._ (2024), and horizonless
compact objects Boyanov  _et al._ (2023).

Ref.Jaramillo (2022) has recently expanded the applications of pseudospectrum
in the field of gravitational wave physics. This work focuses on dynamic
phenomena, specifically transient growths in linear dynamics and the
occurrence of pseudoresonances, which arise from non-self-adjoint operators.
While the decaying QNMs of the evolution operator govern the later stages of
the linear system’s dynamics, its non-self-adjoint nature can lead to
unexpected transient growth in the early stages. These transient and
pseudoresonant effects are not fully captured by eigenvalues alone. Instead,
they are encoded within the comprehensive structure of the pseudospectrum and
evaluated using dedicated quantities derived from it. Using the hyperboloidal
approach, Ref.Jaramillo (2022) reveals that there is no early-time transient
growth observed in asymptotically flat spacetimes. This result remains
unchanged regardless of the specific potential governing the wave dynamics and
can be mathematically attributed to the escape of the field through the
spacetime boundaries, including the horizon and infinity. The study Boyanov
_et al._ (2023) investigates the pseudospectrum of the linear operator
associated with exotic compact objects (ECOs) featuring reflective surfaces,
aiming to gain insights into transient growth. In a more recent work,
Ref.Carballo and Withers (2024) explores transient dynamics within the
framework of general relativity, including asymptotically dS black holes,
asymptotically AdS black holes, and asymptotically flat black holes. The study
demonstrates that non-normality gives rise to arbitrarily long-lived sums of
short-lived QNMs.

We believe that the study of pseudospectrum and transient effects should not
be limited to the framework of general relativity and the hyperboloidal
coordinates. Instead, it is worthwhile to investigate pseudospectrum and
transient from the perturbation equation within the context of modified
gravity theories. Modified gravity theories offer valuable insights into the
constraints of Einstein’s gravity theory, tackle challenges faced by current
models, and play a crucial role in advancing our comprehension of the
Universe. Einstein-Gauss-Bonnet (EGB) gravity whose Lagrangian contains only
the linear and quadratic terms of spacetime curvature is widely studied among
theoretical physicists. To gain a deeper understanding of the modified gravity
theory, one effective approach is to delve into its black hole solutions. The
Maeda-Dadhich black hole, a type of Kaluza-Klein (KK) black hole, also
constitutes an exact vacuum solution within EGB gravity with a cosmological
constant that holds a distinctive relationship to the Gauss-Bonnet coupling
constant Maeda and Dadhich (2006). Importantly, a generalized master equation
for tensor perturbations of EGB gravity has been derived in our previous two
works Cao and Wu (2021); Cao  _et al._ (2023). This perturbation equation
provides us with an opportunity to study the pseudospectrum and transient
effects under modified gravity theories.

The objective of this study is to investigate the pseudospectrum of the tensor
perturbation equation for Maeda-Dadhich black holes in EGB gravity, with a
focus on understanding spectral and dynamical instability as well as transient
effects. To incorporate the influence of extra dimensions, we introduce the
characteristic tensors of $\mathcal{K}^{n-4}$, allowing us to reformulate the
perturbation equation as a scalar field equation in four-dimensional spacetime
$M^{4}$. This scalar equation differs significantly from the Klein-Gordon
equation, as the coefficients of the second-order covariant derivative are no
longer solely determined by the metric of $M^{4}$. By employing the standard
procedure of separating the angular part, we transform the scalar equation
into a Schrödinger-like equation. We then utilize the ingoing
Eddington–Finkelstein coordinate to obtain a generalized eigenvalue problem,
which differs in structure from the standard eigenvalue problem. By solving
this generalized eigenvalue problem, we obtain pseudospectra that provide
insights into the spectrum and dynamical instability .

We all know that the strong external environment may affect the linear
perturbation, making it unstable in dynamics, that is, resulting in an
exponentially growing mode. So, we can’t help but ask how much external
perturbation can achieve it. With the energy norm, we define the (complex)
stability radius to quantitatively characterize this critical intensity, and
find different parameters correspond to different stability radii. More
importantly, we study the transient effect of linear tensor perturbations of
Maeda-Dadhich black holes in EGB gravity in terms of the norm of the evolution
operator. It is found that the numerical abscissa is larger than zero, which
indicated the transient growth will happen. This transient growth can also
lead to an increase in the maximum value of the waveform, which we illustrate
through numerical simulations.

The remains of this paper are organized as follows. In Sec.II, we give a very
short review of Maeda-Dadhich black hole, meanwhile the Schrödinger-like
equation is derived. Sec.III offers some necessary preparations for numerical
calculation of pseudospectrum, including ingoing Eddington–Finkelstein
coordinates and the associated physical energy norm. In Sec.IV, we show the
QNM frequencies and the pseudospectrum to convey the properties of spectrum
instability. In Sec.V, we clarify the relation among pseudospectrum, stability
radius and transient effect. Sec.VI provides the conclusions and discussions.
In Appendix.A, the definition of the pseudospectrum of a generalized
eigenvalue problem is given. In Appendix.B, we derive the expression of
$\partial_{v}E[\varphi]$. In Appendix.C, the numerical approach of
pseudospectrum is come up. In Appendix.D, the convergence test of QNMs is
shown. In order to obtain the waveform, the numerical method used in this work
is appeared at Appendix.E.

##  II The Maeda-Dadhich black hole and the Schrödinger-like equation

We begin by introducing a Kaluza-Klein black hole proposed by Maeda and
Dadhich Maeda and Dadhich (2006), which is a vacuum solution in Einstein-
Gauss-Bonnet gravity theory. The action for $n\geq 5$ in the $n$-dimensional
spacetime with a metric $g_{MN}$ is given by

$\displaystyle
S=\int\mathrm{d}^{n}x\sqrt{-g}\Big{[}\frac{1}{2\kappa_{n}^{2}}\Big{(}R-2\Lambda+\alpha
L_{\text{GB}}\Big{)}\Big{]}\,,$ (1)

where $\kappa_{n}$ is the coupling constant of gravity which depends on the
dimension of spacetime, $R$ and $\Lambda$ are the $n$-dimensional Ricci scalar
and the cosmological constant, respectively. The Gauss-Bonnet term is given by

$\displaystyle L_{\text{GB}}=R^{2}-4R_{MN}R^{MN}+R_{MNPQ}R^{MNPQ}\,,$ (2)

where the capital letters $\\{M,N,P,Q,\cdots\\}$ are the indices for the
$n$-dimensional spacetime. The symbol $\alpha$ is the coupling constant of the
Gauss-Bonnet term, which is identified with the inverse string tension,
thereby it’s positive definite. The equation of motion is

$\displaystyle G_{MN}+\alpha H_{MN}+\Lambda g_{MN}=0\,,$ (3)

where $G_{MN}$ is the Einstein tensor and $H_{MN}$ is the Gauss-Bonnet tensor.
The so-called Maeda-Dadhich black hole spacetime is locally homeomorphic to
$M^{4}\times\mathcal{K}^{n-4}$ with the metric reading as

$\displaystyle
g_{MN}\mathrm{d}x^{M}\mathrm{d}x^{N}=-f(r)\mathrm{d}t^{2}+\frac{1}{f(r)}\mathrm{d}r^{2}+r^{2}(\mathrm{d}\theta^{2}+\sin^{2}\theta\mathrm{d}\phi^{2})+r_{0}^{2}\gamma_{ij}\mathrm{d}x^{i}\mathrm{d}x^{j}\,,$
(4)

where $r_{0}$ satisfies

$\displaystyle r_{0}^{2}=-2K\alpha(n-4)(n-5)\,,$ (5)

and $\gamma_{ij}$ is unit metric on the constant curvature manifold
$\mathcal{K}^{n-4}$ with sectional curvature $K=-1$. The metric function
$f(r)$ in Eq.(4) is written as

$\displaystyle
f(r)=1+\frac{r^{2}}{2(n-4)\alpha}\Bigg{\\{}1-\Big{[}1-\frac{2n-11}{3(n-5)}+\frac{4(n-4)^{2}\alpha^{3/2}\mu}{r^{3}}-\frac{4(n-4)^{2}\alpha^{2}q}{r^{4}}\Big{]}^{1/2}\Bigg{\\}}\,,$
(6)

in which $\mu$ and $q$ are arbitrary dimensionless constants with $\mu$
referring to the mass of the central object and $q$ being the chargelike
parameter. For simplicity, $\mu\geq 0$ and $q\leq 0$ will be assumed in the
followings to ensure that the metric function $f(r)$ exists at least one zero
point. One can find some other reasons for choosing such ranges in Cao  _et
al._ (2023). Without loss of generality, the units can be fixed, where the
event horizon radius is at $r_{+}=1$. It means that we have $f(1)=0$. For
convenience, the constant $q$ will be represented by other parameters via

$\displaystyle q$ $\displaystyle=$
$\displaystyle\frac{1}{12\alpha^{2}(n-5)(n-4)^{2}}\Big{[}-2n+11+\alpha(-12n^{2}+108n-240)$
(7) $\displaystyle+\alpha^{3/2}(-960\mu+12\mu n^{3}-156\mu n^{2}+672\mu
n)+\alpha^{2}(-12n^{3}+156n^{2}-672n+960)\Big{]}\,.$

In other words, the four parameters $n$, $\alpha$, $\mu$ and $q$ are not
independent due to the requirement $f(1)=0$, and the metric function $f(r)$
changing with different parameters are shown in Fig.1. From now on, the
parameter $\alpha$ should be referred to be dimensionless.

Figure 1: The metric function $f(r)$ varying with the three parameters $n$,
$\alpha$ and $\mu$ are depicted by maintaining the event horizon being
$r_{+}=1$. The varying parameters corresponding to each color line are
indicated in the legends. The benchmark parameters are $n=7$, $\alpha=0.05$,
and $\mu=0.5$, which is the orange line among above three panels.

The master equation of the tensor perturbation based on the metric Eq.(4) have
been given in Cao and Wu (2021), after the separation of the variables Cao
_et al._ (2023); Dotti and Gleiser (2005), the Klein-Gordan-like equation of a
scalar field $\Phi$ on the manifold $M^{4}$ can be obtained

$\displaystyle\left[\frac{4n-22}{(n-4)(n-5)}g^{ab}-4\alpha\cdot{}^{4}G^{ab}\right]D_{a}D_{b}\Phi+\left[\frac{2+\gamma}{(n-4)(n-5)}{}^{4}R+\frac{3(n-6)(2+\gamma)}{\alpha(n-4)^{2}(n-5)^{2}}\right]\Phi=0\,,$
(8)

where ${}^{4}R$, ${}^{4}G^{ab}$ are the Ricci scalar and the Einstein tensor
on the manifold $M^{4}$, respectively. In fact, the scalar $\Phi$ is the
amplitude of the characteristic field $\bar{h}_{ij}$ corresponding to the
characteristic value $\gamma$. One can refer to Cao and Wu (2021); Cao  _et
al._ (2023) for more detailed processes of the above equation. It should be
pointed out here that the constant $\gamma$ appeared in Eq.(8) is the
eigenvalue of the Laplace-Beltrami operator acting on the symmetric,
transverse, and traceless (STT) tensor fields with rank $2$ on the
$(n-4)$-dimensional hyperbolic space $\mathcal{K}^{n-4}$. According to the
reference Camporesi and Higuchi (1994) [see Eq.(2.4) therein], we can obtain

$\displaystyle\gamma=-\Big{[}\zeta^{2}+\Big{(}\frac{n-5}{2}\Big{)}^{2}+2\Big{]}\,,$
(9)

with $\zeta^{2}$ is continuous and positive. Thus from this perspective, the
constant $\gamma$ is limited by $\gamma<-(n-5)^{2}/4-2$.

Separating the variable $\Phi(t,r,\theta,\phi)$ as

$\displaystyle\Phi(t,r,\theta,\phi)=S(r)\varphi(t,r)Y_{lm}(\theta,\phi)\,,$
(10)

where $Y_{lm}(\theta,\phi)$ is the spherical harmonics, we obtain a wave
equation with a potential $V(r)$ as follows

$\displaystyle\Big{[}-\frac{\partial^{2}}{\partial
t^{2}}+\frac{\partial^{2}}{\partial
r_{\star}^{2}}-V(r)\Big{]}\varphi(t,r)=0\,.$ (11)

For Eq.(10) and Eq.(11), the scaling function $S(r)$ and the effective
potential $V(r)$ can be found in Cao  _et al._ (2023), $V(r)$ is a function of
$n$, $l$, $\alpha$, $\mu$, and $\gamma$, and we omit the symbols of these
parameters in the expression $V(r)$. The tortoise coordinate $r_{\star}$ is
defined as $\mathrm{d}r_{\star}=\mathrm{d}r/f(r)$ here. When $r\to+\infty$, we
obtain the asymptotic behavior of the effective potential as follows

$\displaystyle V(r)\sim V_{0}(\alpha,n,\gamma)r^{2}\,,$ (12)

where we only keep the leading order term. The stability requirement demands
that $V(r)$ tends towards positive infinity, i.e., $V_{0}(\alpha,n,\gamma)>0$
Buell and Shadwick (1995), which also limits the range of $\gamma$ Cao  _et
al._ (2023). Combining two limitations of $\gamma$, the ranges of $\gamma$ in
different spacetime dimensions are shown in Tab.1.

Spacetime dimension $n$ | Value range of $\gamma$
---|---
$n=6$ | $\backslash$
$n=7$ | $-4.82843<\gamma<-3$
$n\geq 8$ | $\gamma<-(n-5)^{2}/4-2$
Table 1: The value ranges of $\gamma$ with respect to different spacetime
dimensions $n$, which are solved by two limitations $\gamma<-(n-5)^{2}/4-2$
and $V_{0}(\alpha,n,\gamma)>0$. The symbol $\backslash$ represents the empty
set.

Furthermore, Fig.2 shows the variation of the effective potential $V(r)$ with
different parameters $n$, $l$, $\alpha$, $\mu$, and $\gamma$. We choose these
parameters for four reasons: (i) the ranges of $q$ and $\mu$ are selected as
$q\leq 0$ and $\mu\geq 0$, (ii) the potential $V(r)$ should remain positive
from the event horizon to infinity, (iii) the effective AdS radius Cao  _et
al._ (2023) calculated with these parameters can be compared with the event
horizon $r_{+}$, and it should also ensure that $r_{+}$ is the outermost event
horizon and $r_{+}=1$ if there exists another horizon $r_{-}$, (iv) the
constant $\gamma$ should be chosen as illustrated in Tab.1. The ranges for the
five parameters are fixed as $n=[7,9]$, $l$ starting from $l=2$,
$\alpha=[0.02,0.08]$, $\mu=[1.2,1.8]$, and $\gamma=[-4.8,-3.1]$ in the overall
subsequent calculations to investigate physical phenomena, although they are
not unique ranges for these parameters. There are two sets of parameters that
need to be clarified. One is the set of benchmark parameters, which are fixed
as $n=7$, $l=2$, $\alpha=0.05$, $\mu=1.5$, and $\gamma=-4.5$. In these
parameters, the effective AdS radius is $\ell_{\text{AdS}}=1.01206$, which is
comparable to the event horizon. The other set is when the spacetime dimension
is $n=9$, and currently $\gamma$ is limited to $\gamma<-6$ due to condition
(iv). Therefore, $\gamma=-6.5$ is fixed when $n=9$.

Figure 2: The effective potential $V(r)$ varying with different parameters
$n$, $l$, $\alpha$, $\mu$ and $\gamma$ are depicted by maintaining the event
horizon be $r_{+}=1$. The varying parameters corresponding to each color line
are indicated in the legends. The benchmark parameters are $n=7$, $l=2$,
$\alpha=0.05$, $\mu=1.5$ and $\gamma=-4.5$, which is the orange line among
above five panels.

##  III The set up

###  III.1 The establishment of eigenvalue problems

In this subsection, we utilize the ingoing Eddington-Finkelstein (EF)
coordinates to cast the problem of solving for the QNM frequencies as an
eigenvalue problem Cownden  _et al._ (2024). There are two reasons for
choosing this coordinate system. First, these coordinates exhibit regularity
at the black hole horizon and the AdS boundary is allowed to arrive at. By
utilizing these coordinates, we can impose the necessary boundary conditions
at the extremities of the physical domain, which correspond to the two ends of
$v=\text{constant}$ surfaces. However, it is important to note that this
approach is not feasible in asymptotically flat spacetimes, where surfaces of
constant $v$ lead to past null infinity $\mathcal{J}^{-}$ instead of future
null infinity $\mathcal{J}^{+}$. The motivation for proposing hyperboloidal
coordinates is to address this issue Zenginoglu (2008, 2011); Zenginoğlu
(2024); Jaramillo  _et al._ (2021); Panosso Macedo (2024). However, in our
current analysis, it is found that hyperboloidal coordinates are not necessary
and the ingoing EF coordinates are simpler to work with. The drawback of using
these coordinates is that it leads to a generalized eigenvalue problem, which
has a different structure compared to the standard eigenvalue problem. The
second reason for choosing the ingoing EF coordinates is based on the
discussion regarding the structural aspects of the AdS black hole
pseudospectrum, as presented by Boyanov et al. Boyanov  _et al._ (2024). Their
work suggests that the convergence properties of the pseudospectrum are more
favorable in such case. However, the discussion on the convergence of the
pseudospectrum goes beyond the scope of our work.

Given the two aforementioned reasons, it is reasonable for us to use the
ingoing EF coordinates to handle our problem. For convenience, we can
compactify the radial coordinate $r$ by introducing $z=1/r$. By doing so, the
black hole event horizon is located at $z=1$ and the AdS boundary is located
at $z=0$. Therefore, combining the ingoing EF coordinate and the
compactification coordinate $z$, we have the coordinate transformation, which
is given by

$\displaystyle
v=t+r_{\star}\,,\quad\frac{\mathrm{d}r_{\star}}{\mathrm{d}r}=\frac{1}{f(r)}\,,\quad
z=\frac{1}{r}\,.$ (13)

From the above coordinate transformations, one obtains the metric in ingoing
EF coordinates on the four dimensional manifold $M^{4}$

$\displaystyle\mathrm{d}s^{2}=-f(z)\mathrm{d}v^{2}-\frac{2}{z^{2}}\mathrm{d}v\mathrm{d}z+\frac{1}{z^{2}}(\mathrm{d}\theta^{2}+\sin^{2}\theta\mathrm{d}\phi^{2})\,,$
(14)

where

$\displaystyle f(z)=1+\frac{1}{2(n-4)\alpha
z^{2}}\left\\{1-\Big{[}1-\frac{2n-11}{3(n-5)}+4z^{3}(n-4)^{2}\alpha^{3/2}\mu-4z^{4}(n-4)^{2}\alpha^{2}q\Big{]}^{1/2}\right\\}\,.$
(15)

By using the ingoing EF coordinates, Eq.(11) becomes444For convenience, we
write $V(r(z))$ into $V(z)$.

$\displaystyle\Big{[}-2z^{2}f(z)\partial_{vz}+z^{2}f(z)\partial_{z}(z^{2}f(z)\partial_{z})-V(z)\Big{]}\varphi(v,z)=0\,.$
(16)

In order to regularize the solution of Eq.(16) at the AdS boundary $z=0$, we
proceed to redefine

$\displaystyle\varphi(v,z)=z^{\lambda}\psi(v,z)\,,$ (17)

where the index $\lambda>0$ is confirmed as Cao  _et al._ (2023)

$\displaystyle\lambda=\frac{1}{2}\Bigg{\\{}1+\sqrt{1+\frac{16(n-4)^{2}\alpha^{2}V_{0}}{\Big{[}1-\sqrt{\frac{n-4}{3(n-5)}}\Big{]}^{2}}}\Bigg{\\}}\,.$
(18)

By substituting Eq.(17) into Eq.(16), we have

$\displaystyle z^{\lambda}\Bigg{\\{}-2\lambda
zf(z)\partial_{v}\psi(v,z)-2z^{2}f(z)\partial_{vz}\psi(v,z)+z^{4}f^{2}(z)\partial_{z}^{2}\psi(v,z)$
$\displaystyle+\Big{[}2(\lambda+1)z^{3}f^{2}(z)+z^{4}f(z)f^{\prime}(z)\Big{]}\partial_{z}\psi(v,z)$
$\displaystyle+\Big{[}\lambda(\lambda+1)z^{2}f^{2}(z)+\lambda
z^{3}f(z)f^{\prime}(z)-V(z)\Big{]}\psi(v,z)\Bigg{\\}}=0\,,$ (19)

where the prime denotes derivative with respect to $z$. Therefore, through
multiplying the factor $\Big{[}z^{\lambda+1}f(z)\Big{]}^{-1}$ into Eq.(III.1),
we obtain the equation of evolution in terms of $\psi(v,z)$, which is given
by555We may call $H=L_{2}^{-1}L_{1}$ is the Hamiltonian of our present system.

$\displaystyle L_{2}\partial_{v}\psi(v,z)=L_{1}\psi(v,z)\,,$ (20)

where the operators $L_{1}$ and $L_{2}$ are satisfied with

$\displaystyle
L_{1}=z^{3}f(z)\partial_{z}^{2}+\Big{[}2(\lambda+1)z^{2}f(z)+z^{3}f^{\prime}(z)\Big{]}\partial_{z}+\Big{[}\lambda(\lambda+1)zf(z)+\lambda
z^{2}f^{\prime}(z)-\frac{V(z)}{zf(z)}\Big{]}\,,$ (21) $\displaystyle
L_{2}=2z\partial_{z}+2\lambda\,.$ (22)

In order to obtain QNM frequencies, we perform a Fourier transform like
$\psi(v,z)\sim e^{i\omega v}\psi(z)$, then we are going to get a generalised
eigenvalue problem

$\displaystyle(L_{1}-i\omega L_{2})\psi(z)=0\,.$ (23)

By using the Chebyshev grid to discretize two operators $L_{1}$ and $L_{2}$ ,
one gets its matrix versions $\mathbf{L}_{1}$ and $\mathbf{L}_{2}$. The
generalised eigenvalues $\omega$ (QNM frequencies) of such matrix pencil are
derived from

$\displaystyle\text{det}(\mathbf{L}_{1}-i\omega\mathbf{L}_{2})=0\,.$ (24)

The pseudospectrum is also based on Eq.(23). One can find more details in
terms of the pseudospectrum of such generalized eigenvalue problems with
matrix version in Appendix.A. So far, we have completed the establishment of
the eigenvalue problem. Finally, to increase readability, we give the reason
for choosing such $\lambda$. Note that $f\sim f_{0}/z^{2}$ and $V(z)\sim
V_{0}/z^{2}$ as $z\to 0$, in which

$\displaystyle
f_{0}=\frac{1}{2(n-4)\alpha}\Big{[}1-\sqrt{\frac{n-4}{3(n-5)}}\Big{]}\,.$ (25)

It is easy to find that the first two terms of $L_{1}$ and the two terms of
$L_{2}$ are all regular in the interval $[0,1]$. The third term in $L_{1}$ is
given by

$\displaystyle\lambda(\lambda+1)zf(z)+\lambda
z^{2}f^{\prime}(z)-\frac{V(z)}{zf(z)}=\frac{1}{z}\Big{[}\lambda(\lambda-1)f_{0}-\frac{V_{0}}{f_{0}}\Big{]}+\mathcal{O}(z)\,,\quad
z\to 0\,.$ (26)

In order to keep the above term regular, the coefficient of $1/z$ has to be
vanished. Therefore, the index $\lambda>0$ should be chosen as Eq.(18).

###  III.2 The energy norm

In order to quantify the size of perturbations, we should define an
appropriate norm $\lVert\cdot\rVert_{E}$, in which the energy norm is more
natural. As is well known or empirically speaking, a wide class of
perturbations around a black hole can be reduced into such form like Eq.(11).
Such equation has $1+1$ covariant form

$\displaystyle\eta^{\alpha\beta}\nabla_{\alpha}\nabla_{\beta}\varphi-V\varphi=0\,,$
(27)

where $\eta_{\alpha\beta}$ is an effective metric with ingoing EF coordinates
($v$, $r_{\star}$) in $2$-dimensional Minkowski spacetime. On the other hand,
it is easy to find that Eq.(27) stems from the action

$\displaystyle
S\sim-\frac{1}{2}\int\mathrm{d}^{2}y\sqrt{-\eta}\Big{(}\eta^{\alpha\beta}\nabla_{\alpha}\bar{\varphi}\nabla_{\beta}\varphi+V\bar{\varphi}\varphi\Big{)}\,.$
(28)

The energy momentum tensor of the scalar field $\varphi$ from the action (28)
reads as Jaramillo  _et al._ (2021); Cownden  _et al._ (2024)

$\displaystyle
T_{\alpha\beta}=\frac{1}{2}\nabla_{\alpha}\bar{\varphi}\nabla_{\beta}\varphi+\frac{1}{2}\nabla_{\alpha}\varphi\nabla_{\beta}\bar{\varphi}-\frac{\eta_{\alpha\beta}}{2}(\nabla_{\mu}\varphi\nabla^{\mu}\bar{\varphi}+V\varphi\bar{\varphi})\,,$
(29)

where the bar means complex conjugation. There are two Noether conserved
quantities associated with this action (28). To be specific, energy is the
invariance quantity with respect to the time translation, and the other
Noether charge $Q$ is related to $U(1)$ gauge transformation. Since we are
focusing on the null slicing, then the energy $E[\varphi]$ on the null slicing
can be written as

$\displaystyle
E[\varphi]=\int_{\Sigma_{v}}T_{\alpha\beta}\xi^{\alpha}\kappa^{\beta}\mathrm{d}\Sigma_{v}=\frac{1}{2}\int_{a}^{b}\mathrm{d}r_{\star}\Big{(}\partial_{r_{\star}}\varphi\partial_{r_{\star}}\bar{\varphi}+V\varphi\bar{\varphi}\Big{)}\,,$
(30)

where $\xi=\partial_{v}$ is the Killing vector, $\kappa=-\partial_{r_{\star}}$
is the normal vector to the null hypersurface $\Sigma_{v}$, and $a$ and $b$
are the lower and upper boundary in tortoise coordinate $r_{\star}$. The
derivative of $E[\varphi]$ with respect to $v$ is given by

$\displaystyle\partial_{v}E[\varphi]=-\frac{1}{4}|\partial_{r_{\star}}\varphi|^{2}\Big{|}^{b}_{a}+\frac{1}{4}V|\varphi|^{2}\Big{|}^{b}_{a}+\int_{a}^{b}\mathrm{d}r_{\star}\Big{(}\frac{1}{2}V\partial_{v}|\varphi|^{2}-\frac{1}{4}\frac{\mathrm{d}V}{\mathrm{d}r_{\star}}|\varphi|^{2}\Big{)}\,.$
(31)

The above equation is deducted from Eqs.(30) by using Eq.(27) with coordinate
$(v,r_{\star})$ and integration by parts, and one can find more details in
Appendix.B. It should be emphasized that the sign of $\partial_{v}E[\varphi]$
can not be identified directly, which is pretty different from the
hyperboloidal case Carballo and Withers (2024). Furthermore, the energy norm
of $\psi$ based on Eq.(30) is given by Cownden  _et al._ (2024); Jaramillo
_et al._ (2021)

$\displaystyle\lVert\psi\rVert_{E}^{2}\equiv\langle\psi,\psi\rangle_{E}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\int_{\Sigma_{v}}\Big{(}\partial_{r_{\star}}\bar{\varphi}\partial_{r_{\star}}\varphi+V\bar{\varphi}\varphi\Big{)}\mathrm{d}r_{\star}$
(32) $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{0}^{1}\Big{[}z^{2\lambda+2}f(z)(\partial_{z}\bar{\psi})(\partial_{z}\psi)+\lambda
z^{2\lambda+1}f(z)\psi\partial_{z}\bar{\psi}+\lambda
z^{2\lambda+1}f(z)\bar{\psi}\partial_{z}\psi$
$\displaystyle+z^{2\lambda-2}\Big{(}\lambda^{2}z^{2}f(z)+\frac{V(z)}{f(z)}\Big{)}\bar{\psi}\psi\Big{]}\mathrm{d}z\,.$

In the second line, we have used Eq.(13) and Eq.(17). The upper and lower
limits of integration come from the range of compactified coordinate $z$.
Since one can proof that $V(z=1)=0$ Cao  _et al._ (2023), the definition of
the inner product

$\displaystyle\langle\psi_{1},\psi_{2}\rangle_{E}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{0}^{1}\Big{[}z^{2\lambda+2}f(z)(\partial_{z}\bar{\psi}_{1})(\partial_{z}\psi_{2})+\lambda
z^{2\lambda+1}f(z)(\psi_{2}\partial_{z}\bar{\psi}_{1}+\bar{\psi}_{1}\partial_{z}\psi_{2})$
(33)
$\displaystyle+z^{2\lambda-2}\Big{(}\lambda^{2}z^{2}f(z)+\frac{V(z)}{f(z)}\Big{)}\bar{\psi}_{1}\psi_{2}\Big{]}\mathrm{d}z\,.$

is regular for $z\in[0,1]$, which is important for numerical calculation. When
$\psi_{1}=\psi_{2}=\psi$, Eq.(33) reduces into Eq.(32). One can see Gasperin
and Jaramillo (2022) for further discussions of the inner product. In the
Appendix.C, we provide a detailed introduction to the numerical technique of
using the energy norm defined above to determine the pseudospectrum, which in
turn will inform us of the instability of the QNM frequencies.

##  IV The results of the spectrum and the pseudospectrum

In following sections, we will analyze the properties of linear perturbation
dynamics in Maeda-Dadhich black hole in view of pseudospectrum framework. The
results of the spectrum and the pseudospectrum are shown in this section. One
can find the definition of $\epsilon$-pseudospectrum in Appendix.A for our
generalised eigenvalue problem. Technically, based on the spectral method, the
non-normal differential operator that describes the linear system will becomes
a $(N+1)\times(N+1)$ discrete matrix through Chebyshev ploynomials, if we
choose that the resolution is $N$. Besides, the invariant subspace method is
exploited to improve the computational efficiency for pseudospectrum, which
has already been used in Cao  _et al._ (2024).

The fundamental mode varying with different parameters is shown in Fig.3, and
the spectrum and pseudospectrum are displayed in Fig.4 with benmark
parameters, where we solve them with resolution $N=100$ and $100$ digits of
precision. The convergence of such eigenvalues in view of the generalised
eigenvalue problem are discussed in detail in Appendix.D.

In Fig.3, we depict the variation of the fundamental mode with different
parameters in the Maeda-Dadhich black hole. The ranges of each parameter are
specified in Sec.II. Each color and shape of the points in the panel represent
a different parameter, as indicated in the legend. For example, the blue
triangle points represent the fundamental mode variation with respect to
$\alpha$, and the direction of the colored arrow indicates the trend of the
fundamental mode as the corresponding parameter increases. In general, a set
of parameters $(n,l,\alpha,\mu,\gamma)$ determines the effective potential,
which in turn affects the fundamental mode. The fundamental mode corresponding
to the benchmark parameters is located at the intersection of various color
and shape points. We amplify the data around this intersection in the subgraph
indicated by the black arrow. As $\alpha$ and $\mu$ increase, the effective
potential is depressed, resulting in a smaller real part of the fundamental
mode. For $\alpha$, the imaginary part of the fundamental mode becomes
smaller, while for $\mu$, it becomes larger. On the other hand, as $l$ and
$\gamma$ increase, the effective potential is lifted, leading to a larger real
part of the fundamental mode. For $l$, the imaginary part of the fundamental
mode becomes smaller, while for $\gamma$, it becomes larger. Especially, it is
found that for spacetime dimensions $n=8$ and $n=9$, the fundamental modes are
both located on the imaginary axis, corresponding to $1.523601i$ and
$1.788444i$, respectively, which are not shown in the Fig.3.

Figure 3: The fundamental mode changes with different parameters
$(l,\alpha,\mu,\gamma)$ with numerical resolution $N=100$ and $100$ digits of
precision. The benchmark parameters are $n=7$, $l=2$, $\alpha=0.05$,
$\mu=1.5$, $\gamma=-4.5$, and the corresponding fundamental mode is
represented by $\star$ in the panel and the amplified subgraph. The ranges of
parameters used are: $l=[2,5]$, $\alpha=[0.02,0.08]$, $\mu=[1.2,1.8]$ and
$\gamma=[-4.8,-3.1]$, for last three parameters we select $11$ discrete data
points. Each color and shape of the points represents each parameter, as
illustrated in the legend.

In Fig.4, we present the spectrum and pseudospectrum with the benchmark
parameters. In the left panel, the green and yellow points represent the QNM
frequencies in the complex plane. The family of yellow points is located on
the imaginary axis, which is meaningful and corresponds to physical modes.
However, in flat spacetime, these points are introduced due to the
discretization process and are considered meaningless. The right panel shows
the $\epsilon$-pseudospectrum around the fundamental mode666The overtone
number is denoted by $\bar{n}$ and it start counting from $1$, the fundamental
mode is defined by the smallest value of imaginary part of overall
eigenvalues, including eigenvalues on the imaginary axis.
$\omega_{1}=4.051108+1.659013i$. The open sets formed by the contour lines of
the $\epsilon$-pseudospectrum indicate spectrum instability. Specifically, it
implies that the QNM frequencies can migrate a further distance in the complex
plane than what is possible under a given perturbation size. This property is
intrinsic to non-normal operators and is different from the case of normal
operators, where the contour lines of the $\epsilon$-pseudospectrum form
closed sets. It is worth noting that this situation is different from the case
of the hyperboloidal coordinates in global AdS4 Cownden  _et al._ (2024),
where the pseudospectrum structure consists of perfectly closed spherical
contour lines, indicating the stability of the spectrum.

Figure 4: The results of spectrum (QNM frequencies) and
$\epsilon$-pseudospectrum with resolution $N=100$ and 100 digits of precision.
The parameters are the mentioned benchmark parameters, i.e. $n=7$, $l=2$,
$\alpha=0.05$, $\mu=1.5$ and $\gamma=-4.5$. The contour lines in both panels
denote the boundary of $\epsilon$-pseudospectrum sets given the perturbation
size $\epsilon$, whose values are $\log_{10}\epsilon$ and are shown in the
legends. The yellow and green points in the left panel are the QNM frequencies
in complex plane. The right panel shows $\epsilon$-pseudospectrum around the
fundamental mode $\omega_{1}=4.051108+1.659013i$. Note that the color in the
right enlarged panel is inconsistent with the left original panel. One should
refer to the actual legend on the right side of each panel.

It is worth mentioning that there are two kinds of instability in view of
pseudospectrum framework, one is the spectrum instability based on the
spectrum method, which is reflected by the open structure of
$\epsilon$-pseudospectrum, the other is the dynamical instability based on the
linear evolution system, which is defined by the eigenvalues entering into the
unstable region from the stable region under perturbations in view of
pseudospectrum. The dynamical instability implies that the eigenmodes of the
linear system will have exponential growth. The spectrum instability has been
shown in this section, the dynamical instability will be discussed in the
following section.

##  V The pseudospectrum and the transient effect

The properties of non-normal matrices and operators encompass various aspects,
notably their ability to generate the transient effect in time-dependent
dynamical systems. The reminded fact is that eigenvalues are unable to
accurately describe all features of the transient effect, a natural question
is that can pseudospectra provide a more accurate physical picture? The answer
is unequivocally affirmative: while pseudospectra may not always furnish a
precise result, they are adept at detecting and measuring the transients that
eigenvalues overlook Trefethen and Embree (2005). In this section, we study
the correlation between the transient dynamics of a non-normal Hamiltonian $H$
which comes from Eq.(20) and its pseudospectrum.

We begin this section by writing Eq.(20) into the following linear evolution
system

$\displaystyle\frac{\mathrm{d}\mathbf{\Psi}}{\mathrm{d}v}=\mathbf{H}\mathbf{\Psi}\,,$
(34)

where $\mathbf{H}=\mathbf{L}_{2}^{-1}\mathbf{L}_{1}$ is a matrix acting on the
function $\mathbf{\Psi}(v)$. Here, for the convenience of subsequent numerical
calculation, we have already discretized $\psi(v,z)$ and the Hamiltonian $H$.
Therefore, the formal solution777If $\mathbf{H}$ is a matrix or bounded
operator, $e^{v\mathbf{H}}$ can be defined by a convergent power series. For
an unbounded operator $\mathbf{H}$, the meaning of $e^{v\mathbf{H}}$ comes
from the mathematical theory of semigroups Trefethen and Embree (2005). of the
problem (34) is

$\displaystyle\mathbf{\Psi}(v)=e^{v\mathbf{H}}\mathbf{\Psi}(0)\,,$ (35)

where $\mathbf{\Psi}(0)$ is the initial condition, and we call
$e^{v\mathbf{H}}$ the evolution operator in null slicing. Within the
gravitational context, the evolution operator, which is adhering to the
concepts proposed in Jaramillo (2022), enables us to gain the evolution
information of field, especially for that in the initial stage.

Specifically, the norm of the operator $e^{v\mathbf{H}}$ tracks the maximum
rate of growth over time for the solutions $\Psi(v)$. This leads to the
definition of $\lVert e^{v\mathbf{H}}\rVert_{E}$, i.e.,

$\displaystyle G(v)\equiv\lVert
e^{v\mathbf{H}}\rVert_{E}=\sup_{\mathbf{\Psi}(0)\neq\mathbf{0}}\frac{\lVert
e^{v\mathbf{H}}\mathbf{\Psi}(0)\rVert_{E}}{\lVert\mathbf{\Psi}(0)\rVert_{E}}\,,$
(36)

where the norm $\lVert\cdot\rVert_{E}$ is defined in Eq.(32). The next
important task is to know estimation on the size of $\lVert
e^{v\mathbf{H}}\rVert_{E}$. It is important to convert all energy norm
$\lVert\cdot\rVert_{E}$ related quantities into $2$-norm related quantities
for the convenience of numerical calculations [see also Eq.(54)]. For the
function $G(v)$, we have

$\displaystyle G(v)=\lVert e^{v\mathbf{H}}\rVert_{E}=\lVert
e^{v\tilde{\mathbf{H}}}\rVert_{2}\,,$ (37)

where $\tilde{\mathbf{H}}=\mathbf{W}\mathbf{H}\mathbf{W}^{-1}$. Note that the
conclusion in Trefethen (1999) has been already utilized in the second
equality, and the matrix $\mathbf{W}$ is derived from the Cholesky
decomposition of Gram matrix $\tilde{\mathbf{G}}^{E}$ (see Appendix.C). Now,
computing $G(v)$ is straightforward since $2$-norm of a matrix is given by its
maximum singular value so that we can do it by Mathematica.

In order to describe the dynamic behavior of the function $G(v)$, three
abscissas are supposed to introduced, i.e., spectral abscissa
$\bar{\alpha}(\mathbf{H})$, $\epsilon$-pseudospectral abscissa
$\bar{\alpha}_{\epsilon}(\mathbf{H})$ and numerical abscissa
$\bar{\omega}(\mathbf{H})$888In this section, for convenience, all the
calculations related to the pseudospectrum are based on $\mathbf{H}$ or
$\tilde{\mathbf{H}}$, which differs from the pseudospectrum calculations of
Sec.IV by the presence of the imaginary unit $i$. We trust that the readers
will understand the connections and differences between these two.. Their
definitions and their relationship with $G(v)$ are as follows Trefethen and
Embree (2005); Jaramillo (2022); Boyanov  _et al._ (2023):

  1. 1.

The spectral abscissa $\bar{\alpha}(\mathbf{H})$ is defined by

$\displaystyle\bar{\alpha}(\mathbf{H})=\sup\text{Re}(\sigma(\mathbf{H}))=\sup\\{\text{Re}(\lambda),\text{with}\,\lambda\in\sigma(\mathbf{H})\\}\,,$
(38)

where $\sigma(\mathbf{H})$ represents the spectra of $\mathbf{H}$, i.e., the
QNM frequencies. In the QNM aspects, $\bar{\alpha}$ is called spectral gap
which is useful in SCCC Cardoso  _et al._ (2018). The lower bound of $G(v)$ is
controlled by the spectral abscissa,

$\displaystyle G(v)\geq e^{v\bar{\alpha}(\mathbf{H})}\,,\quad\forall v\geq
0\,.$ (39)

It should be noted that even if $\bar{\alpha}<0$, i.e., the solution decays at
late times, there exists the possibility of initial transient growth.

  2. 2.

In order to estimate the maximum of the possible dynamical transient, which is
dictated by the maximum growth in the norm of the evolution operator during
intermediate stages of time, we are supposed to use the so-called
pseudospectral abscissa. For each $\epsilon>0$, the $\epsilon$-pseudospectral
abscissa $\bar{\alpha}_{\epsilon}(\mathbf{H})$ is defined by

$\displaystyle\bar{\alpha}_{\epsilon}(\mathbf{H})=\sup\text{Re}(\sigma_{\epsilon}(\mathbf{H}))=\sup\\{\text{Re}(\lambda),\text{with}\,\lambda\in\sigma_{\epsilon}(\mathbf{H})\\}\,,$
(40)

where $\sigma_{\epsilon}(\mathbf{H})$ is the pseudospectrum of $\mathbf{H}$
(see Appendix.A). In terms of the pseudospectral abscissa, the supremum of
$G(v)$ can be estimated, i.e.,

$\displaystyle\sup_{v\geq
0}G(v)\geq\frac{\bar{\alpha}_{\epsilon}(\mathbf{H})}{\epsilon}\,,\quad\forall\epsilon>0\,.$
(41)

As a consequence of maximizing Eq.(41) over $\epsilon>0$, this inherently
results in the definitive lower bound,

$\displaystyle\sup_{v\geq 0}G(v)\geq\mathcal{K}(\mathbf{H})\,.$ (42)

The constant $\mathcal{K}(\mathbf{H})$ is called the Kreiss constant of
$\mathbf{H}$, which is defined by
$\sup_{\epsilon>0}\bar{\alpha}_{\epsilon}/\epsilon$. In Fig.5, we depict its
tendency with benchmark parameters as $\epsilon$ increases, the
$\epsilon$-pseudospectral absicissa $\bar{\alpha}_{\epsilon}$ is calculated by
the criss-cross algorithm Trefethen and Embree (2005); Burke  _et al._ (2003).
Obviously, the Kreiss constant in the set of benchmark parameters equals to
$1$.

Figure 5: The tendency of $\bar{\alpha}_{\epsilon}/\epsilon$ varies with
$\epsilon$ from $0.5$ to $6.5$, and $\delta\epsilon=0.5$. The blue points are
discrete data, which are calculated with the benchmark parameters, and the
horizontal orange line corresponds to $\bar{\alpha}_{\epsilon}/\epsilon=1$.
The resulotion $N=100$ and $100$ digits of precision are used.

  3. 3.

The numerical abscissa $\bar{\omega}(\mathbf{H})$ can be obtained by

$\displaystyle\bar{\omega}(\mathbf{H})=\sup\sigma\Big{(}\frac{1}{2}(\mathbf{H}+\mathbf{H}^{\dagger})\Big{)}\,.$
(43)

Note that the symbol $\dagger$ means the adjoint of an operator (matrix). The
discrete version $\mathbf{H}^{\dagger}$ can be obtained from

$\displaystyle\mathbf{H}^{\dagger}=(\tilde{\mathbf{G}}^{E})^{-1}\cdot\mathbf{H}^{\star}\cdot\tilde{\mathbf{G}}^{E}\,,$
(44)

where the symbol $\star$ stands for the conjugate transpose. The numerical
abscissa in Eq.(43) characterizes the initial slope of the norm of the
evolution operator $e^{v\mathbf{H}}$, which means that

$\displaystyle\bar{\omega}(\mathbf{H})=\frac{\mathrm{d}}{\mathrm{d}v}\lVert
e^{v\mathbf{H}}\rVert_{E}\Big{|}_{v=0}\,.$ (45)

Therefore, if the value $\bar{\omega}(\mathbf{H})>0$, there will exist an
initial growth in the evolution. The upper bound of $G(v)$ is controlled by
the numerical abscissa,

$\displaystyle G(v)\leq e^{v\bar{\omega}(\mathbf{H})}\,,\quad\forall v\geq
0\,.$ (46)

If $\bar{\omega}=0$, $G(v)$ fails to exhibit growth at initial stage, and
diminishes as $v$ increases, which signifies the Kreiss constant must
inherently equal unity, i.e., $\mathcal{K}=1$ Trefethen and Embree (2005). But
the reverse is not necessarily true. The model we are studying is a
counterexample where its numerical abscissa is larger than $0$ but the Kreiss
constant is still $1$. The overall numerical abscissa are shown in Fig.6
varying with different parameters.

Figure 6: The numerical abscissa varies with different parameters $n$,
$\alpha$, $l$, $\mu$ and $\gamma$ with resolution $N=100$ and $100$ digits of
precision. The numerical abscissa corresponding to the benchmark parameters
approximately equals to $20.8$, which has been amplified in the first panel
labeled by the black arrow.

###  V.1 The distance to the dynamical instability

As mentioned earlier, there are two types of instability: spectrum instability
and dynamical instability. Spectrum instability is characterized by the open
structure of the pseudospectrum, as discussed in Section IV. Dynamical
instability, on the other hand, refers to the behavior of a linear system
exhibiting exponential growth. The behavior of a linear system is well-known
to be controlled by the spectral abscissa, as shown in Eq. (39). If the
spectrum abscissa is greater than zero, exponential growth occurs. In other
words, a linear system is considered dynamically stable if the spectrum lies
entirely within the open left half-plane. Conversely, if the spectrum extends
into the right half-plane, the system is considered dynamically unstable.
Actually, these two types of instability are closely related. Spectrum
instability is responsible for the significant shift of the spectrum from its
initial state when subjected to external perturbations. This shift can even
result in the spectrum migrating into the so-called dynamically unstable
region.

Naturally, how much external perturbation would lead to the spectrum migrating
to the dynamically unstable region? Although this is a very physical problem,
we can treat it as a mathematical problem. Therefore, it is natural to
investigate the distance to dynamical instability of
$\tilde{\mathbf{H}}=\mathbf{W}\mathbf{H}\mathbf{W}^{-1}$, which controls the
evolution direction of the entire system. Its definition is given by Trefethen
and Embree (2005)

$\displaystyle\min\Big{\\{}\lVert\mathbf{E}\rVert_{2}:\tilde{\mathbf{H}}+\mathbf{E}\,\text{is
unstable}\Big{\\}}\,.$ (47)

The interpretation of Eq.(47) in terms of pseudospectra is straightforward:
what is the smallest $\epsilon$ for which the boundary of
$\sigma_{\epsilon}(\tilde{\mathbf{H}})$ touches the imaginary axis? In Fig.7,
the shaded area represents all possible ranges of the spectra after being
disturbed by an external intensity of $\epsilon$. As $\epsilon$ increases, we
always see that there exists such $\epsilon$ that the boundary of
$\sigma_{\epsilon}(\tilde{\mathbf{H}})$ is just tangent to the imaginary axis.
At this time, we call such $\epsilon$ is the (complex) stability radius of
$\tilde{\mathbf{H}}$. Given a linear system whose Hamiltonian is given by
$\tilde{\mathbf{H}}$, then the theorem Trefethen and Embree (2005) below
provides us with a method for calculating the radius of (complex) stability.

Figure 7: The instability threshold for $\tilde{\mathbf{H}}$ corresponds to
the minimal value of $\epsilon$ at which the boundary of set
$\sigma_{\epsilon}(\tilde{\mathbf{H}})$ intersects the imaginary axis.

###### Theorem 1

The $2$-norm complex stability radius of a matrix $\tilde{\mathbf{H}}$ is
equal to

$\displaystyle\Bigg{(}\sup_{\text{Re}z=0}s_{1}\Big{(}\mathbf{R}(z)\Big{)}\Bigg{)}^{-1}\,,$
(48)

where $s_{1}$ denotes the largest singular value and
$\mathbf{R}(z)=(z-\tilde{\mathbf{H}})^{-1}$.

In this section, we use Eq.(48) to calculate the complex stability radius and
present some typical results varying with different parameters in Fig.8. By
examining the five panels, we can observe certain patterns. In the second to
fourth panels, we find that the stability radius decreases with increasing
$l$, $\alpha$, $\mu$. However, in contrast, the stability radius increases
with increasing $\gamma$. In the first panel, the stability radius for two
different spacetime dimensions $n=7$ and $n=8$ can be compared due to the same
$\gamma$, and the former is more dynamically stable than the latter. For the
case of $n=9$, since we choose $\gamma=-6.5$ as illustrated before, it can not
be compared with other cases. The complex stability radius provides a
quantitative measure of dynamical instability of the linear system. Based on
the results, we can infer that increasing $n$, $l$, $\alpha$, and $\mu$ will
stabilize the dynamical system, while increasing $\gamma$ will destabilize it.

Figure 8: The complex stability radius $\epsilon$ changes with different
parameters $n$, $l$, $\alpha$, $\mu$ and $\gamma$ with resolution $N=100$ and
$100$ digits of precision. The each parameter has eleven points in the second
panel from the fifth panel while the range of spacetime dimension $n$ is
limited into $n=7,8,9$ in the first panel.

In addition, the stability radius can also be determined by the pseudospectral
abscissa $\bar{\alpha}_{\epsilon}(\tilde{\mathbf{H}})=0$ by definition. The
criss-cross algorithm is a commonly used approach to calculate the
pseudospectral abscissa. More details about the criss-cross algorithm can be
found in references Burke  _et al._ (2003); Trefethen and Embree (2005). In
general, we can employ the bisection method to iteratively approach the
desired (complex) stability radius. On the other hand, this also can serve as
a method to verify the reliability of the stability radius results obtained
from Eq.(48). To evaluate its accuracy, we can calculate the pseudospectral
abscissas where $\epsilon$ is given by the stability radius, thereby examining
the errors. The pseudospectral abscissas with parameters that are same with
calculations in stability radius are shown in Fig.9. From this figure, we
observe that the pseudospectral abscissas are controlled within the order of
magnitude of $10^{-6}$, indicating relatively small errors, that provides
further validation of the stability radius results.

Figure 9: The pseudospectral abscissa $\bar{\alpha}_{\epsilon}$ changes with
different parameters $n$, $l$, $\alpha$, $\mu$ and $\gamma$ with resolution
$N=100$ and $100$ digits of precision, in which $\epsilon$ is selected as the
complex stability radius.

###  V.2 The transient effect and the waveform

In this subsection, we will investigate the transient effect and the waveform
of Maeda-Dadhich black hole in EGB gravity theory. The transient effect is
characterized by the growth of the norm of the evolution operator
$e^{v\mathbf{H}}$, where the norm is given by the energy norm [see Eq.(32)].
Therefore, it is crucial to have a comprehensive understanding of the
evolution of the system described in Eq.(35) across different stages. The
evolution governed by the operator $e^{v\mathbf{H}}$ can be divided into three
distinct phases in view of its norm. Firstly, the spectral abscissa
$\bar{\alpha}(\mathbf{H})$ dominates the asymptotic behavior of the linear
system in the limit of late time $v\to\infty$. Secondly, the numerical
abscissa $\bar{\omega}(\mathbf{H})$ controls the initial growth as $v\to 0$
[see Eq.(46)]. The $\epsilon$-pseudospectrum concept plays a key role in
bridging the late-time dynamics governed by $\bar{\alpha}(\mathbf{H})$ and the
early-time behavior dictated by $\bar{\omega}(\mathbf{H})$. It provides a
means to control the growth at intermediate times using the pseudospectral
abscissa $\bar{\alpha}_{\epsilon}(\mathbf{H})$ [see Eq.(41)].

On the other hand, the three abscissas have the relationship with perturbation
size $\epsilon$. In the pioneering work Jaramillo (2022), the author study the
BBH merger waveform as a transient phenomenon under the (strong) hypothesis of
effective BBH non-normal linear dynamics, and illustrate that the (late)
inspiral phase corresponding to the limit $\epsilon\to\infty$, the final
ringdown characterised by $\epsilon\to 0$ and the actual merger by
intermediate values of $\epsilon$. However, within GWs context in view of
pseudospectrum, we consider that the limit $\epsilon\to\infty$ corresponds the
phase after the merger is completed, based on the following reasons. Firstly,
pseudospectrum describes the properties of operator $H$ in Eq.(34) in linear
evolution system, which origins from the master equation (11). And the master
equation describes the properties of a single black hole distorted by some
perturbations, which doesn’t involve how the binary black holes evolve
specifically. Thus, the limitation $\epsilon\to\infty$ should correspond to
the phase after the merger is completed, where a single event horizon has been
formed, at that moment the black hole is at intense distortion. Then
naturally, the late time dominated by ringdown phase corresponds to the
limitation $\epsilon\to 0$, whose evolution behaviour is characterised by
QNMs, i.e. spectral abscissa. We will use a schematic diagram to illustrate
this process, see Fig.10.

Figure 10:  The sketch map that describe the relationship between
$v_{\epsilon}$ and $\epsilon$.

However, the three types of abscissas provide rough estimates of $G(v)$. In
order to obtain a more accurate grasp of the behavior of $G(v)$, we can
utilize equation (37) to establish the relationship between $G(v)$ and time.
It is important to clarify that there is another formulation presented in an
other work Carballo and Withers (2024). In that work, they define $G(\tau)$
using the norm-squared, resulting in $G^{2}(v)$ in our work. However, despite
the difference in notation, both formulations capture the same physical
meaning as they utilize the energy norm. It is also worth noting that $G(v)$
represents the maximum possible energy growth of the field at time $v$.

Figure 11: The relationship between the norm $G(v)$ of the evolution operator
$e^{v\mathbf{H}}$ and the time $v$ within the benchmark parameters.

As a typical example, we depict the variation of $G(v)$ with respect to $v$
using the benchmark parameters, as shown in Fig.11. The blue line represents
the function $G(v)$. In the lower left amplified subgraph, the orange line
represents the tangent line of $G(v)$ at $v=0$, and its slope corresponds to
the numerical abscissa (see Fig.6). In the initial stage, there is a
noticeable growth of $G(v)$ controlled by the numerical abscissa, with values
larger than $1$, as indicated in the two amplified subgraphs. This behavior is
fundamentally different from the hyperboloidal case presented in Carballo and
Withers (2024), where $G(\tau)$ is less than or equal to $1$. The reason for
this difference lies in the fact that in the hyperboloidal case,
$\partial_{\tau}E(\tau)$ can be proven to be less than or equal to zero.
However, in our case, it is difficult to determine the positivity or
negativity of $\partial_{v}E[\varphi]$ [see Eq.(31)]. As a result, there are
no such limitations on the behavior of $G(v)$. In the intermediate stage, we
observe a platform in the $G(v)$ curve. These two properties of $G(v)$
strongly suggest the existence of the transient effect. The reason behind
these properties is that although the presence of an event horizon ensures
energy overflow, the eigenmodes (i.e., QNMs) are not completely orthogonal due
to the non-normality of the operator $H$. This lack of orthogonality allows
for energy exchange or transfer between non-orthogonal modes, resulting in the
transient effect, even though the system remains at a linear level. In the
late stage, starting at $v_{\star}\approx 4.7$, we observe a decline in
$G(v)$, which is dominated by the spectral abscissa.

In traditional gravitational wave observations, there is a focus on how the
waveform $|\Psi|$ changes over time. The waveform typically refers to the
variation of the field or its absolute value with respect to time. Motivated
by this, we depict the waveform as it varies with different initial positions,
initial widths of the wave packet, and observation positions, while keeping
the amplitude unchanged. This is shown in Fig.12. The initial condition
$\psi(0,z)$ is set as a Gaussian wave packet, given by:

$\displaystyle\psi(0,z)=h\exp\Big{[}-\frac{(z-z_{0})^{2}}{2\lambda^{2}}\Big{]}\,,$
(49)

where $z_{0}$ is the initial position, $h$ is the amplitude and $\lambda$ is
the width, respectively. After spatial discretization, we solve the linear
evolution Eq.(34) by using the time-symmetric numerical method Markakis  _et
al._ (2019); O’Boyle  _et al._ (2022); Markakis  _et al._ (2023) instead of
the traditional Runge-Kutta method, see more details in Appendix.E.

In Fig.12, the first column panels depict the same waveform, where the initial
conditions are given by $z_{0}=0.95$, $h=1$, and $\lambda=\sqrt{10}/100$. The
waveform is observed at null infinity, $z_{\text{obs}}=0$. It is observed that
after an amplification in amplitude, which begins at approximately $v\approx
2.2$, the wave packet starts to exponentially decay. The first row panels show
the results when varying the initial positions of the Gaussian wave packet.
The three panels from left to right correspond to $z_{0}=0.95$, $z_{0}=0.5$,
and $z_{0}=0.05$. The qualitative characteristics of the waveform in the
second and third figures are consistent with those in the first figure, but
there are two notable features. Firstly, the closer the initial position of
the wave packet is to the observation position, the earlier the wave packet is
observed. Secondly, when the initial wave packet is located at $z=0$, the
maximum amplitude of the wave packet is observed at null infinity. The second
row panels show the results when varying the initial width of the Gaussian
wave packet. The three panels from left to right correspond to
$\lambda=\sqrt{10}/100$, $\lambda=\sqrt{5}/100$, and $\lambda=\sqrt{1}/100$.
These results indicate that the width of the Gaussian wave packet affects the
duration of the wave packet at the observation position, which is consistent
with Baibhav  _et al._ (2023). The third row panels show the results when
varying the observation positions. The three panels from left to right
correspond to $z_{\text{obs}}=0$, $z_{\text{obs}}=0.5$, and
$z_{\text{obs}}=0.96$. From these panels, we observe that the closer the
observation position is to null infinity, the greater the amplitude of the
waveform obtained.

Figure 12: The nine panels are the waveform calculated with Eq.(34) by using
time-symmetric numerical method (see Appendix.E). The initial condition is set
as a Gaussian wave packet with amplitude $h=1$. The first and second row
panels are the results of waveform observed at null infinity
$z_{\text{obs}}=0$ varying with the initial positions and widths of wave
packet, where parameters are $z_{0}=0.95,0,0.05$ and
$\lambda=\sqrt{10}/100,\sqrt{5}/100,\sqrt{1}/100$ from left to right,
respectively. In third row panels, we change the observation positions,
corresponding to $z_{\text{obs}}=0,0.5,0.96$ from left to right, respectively.
Therefore, the first column panels are the same waveform observed at null
infinity $z_{\text{obs}}=0$, where parameters of initial Gaussian wave packet
are $\lambda=\sqrt{10}/100$ and $z_{0}=0.95$.

It is important to note that the moment when the exponential decay begins for
$G(v)$ does not necessarily correspond to the moment when the exponential
decay begins for the waveform, which is basic question in BBH dynamics that
when is ringdown excited. The waveform depends on the spatial observation
position, and further discussions about the excitation of the ringdown in
terms of the magnitude $|\Psi|$ at fixed observation positions can be found in
Baibhav  _et al._ (2023), while $G(v)$ does not since it has integrated out
the spatial part. However, the existence of transient growth in $G(v)$ does
provide some insights about the waveform observed at a fixed spatial position.
It can at least be concluded that the existence of transient growth allows for
the amplitude of the wave received at infinity to be larger than the initially
given amplitude. The initial width and position of the wave packet do not
change the fact that the amplitude is amplified at null infinity.

##  VI conclusions and discussion

In this paper, we investigate the pseudospectrum and the transient effect of
the tensor perturbation for Maeda-Dadhich black holes in EGB gravity theory.
The topology of the so-called Maeda-Dadhich black holes with dimension $n$ is
the product of a usual four-dimensional spacetime with a negative constant
curvature space with dimension $n-4$. Upon examining the asymptotic expansion
of the metric function, it is observed that the behavior of the metric
function resembles the one of the Reissner-Nordström anti–de Sitter spacetime,
albeit in the absence of the Maxwell field.

With regard to technical details, we have devised the pseudospectrum
computation by using ingoing Eddington-Finkelstein coordinates, which is also
used in Cownden  _et al._ (2024); Boyanov  _et al._ (2024). The equation
governing the perturbations undergoes a reduction to a first-order temporal
form, with the time derivative appearing in a mixed derivative term [see
Eq.(20)]. This transition naturally prompts a reformulation of the
pseudospectrum as a generalized eigenvalue problem. The adoption of these
coordinates, coupled with an appropriate rescaling of the field, not only
drastically simplifies the problem and facilitates the imposition of boundary
conditions, but also remarkably decreases the dimension of the related
matrices (Within the same resolution $N$, the matrix dimension in the
hyperboloidal setting is twice that in the null setting.), thereby lessening
the computational demands of the numerical implementation. Based on the
analysis of the convergence of the hyperbolidal slicing and the null slicing
for the pseudospectra in Boyanov  _et al._ (2024), we are convinced that the
calculation of pseudospectra under the null slicing is appropriate.

The spectrum and pseudospectrum of the model we studied is given in Sec.IV. We
have depicted the behaviour of the fundamental mode with parameters changing
in Fig.3. The decay time of fundamental mode is decreasing as $\mu$ and
$\gamma$ increase, but increasing as $\alpha$ and $l$ increase, while the
fundamental mode is purely imaginary when $n=8$ and $n=9$. In the Fig.4, the
“inverted bridge” structure of pseudospectrum is displayed. Such type of
structure reflects the spectral instability.

In order to further study the pseudospectrum and the transient effect of our
model, where the transient effect is described by $G(v)$. Moreover, three
kinds of abscissas, i.e., the spectral abscissa $\bar{\alpha}$, the
$\epsilon$-pseudospectral abscissa $\bar{\alpha}_{\epsilon}$ and the numerical
abscissa $\bar{\omega}$, are used to describe the dynamics of $G(v)$. First,
the pseudospectral abscissa is used to define the distance to the dynamical
instability, which is also called the (complex) stability radius. The larger
stability radius represents that making the spectral abscissa into dynamical
unstable region requires larger external perturbation energy. We calculate
such stability radii within different parameters, increasing $n$, $l$,
$\alpha$ and $\mu$ reduces the dynamical stability, but increasing $\gamma$
enhances the dynamical stability. It is notable that there is no guaranteed
relationship between the decay time and the stability radius.

More importantly, we firstly report the transient effect in null slicing with
ingoing EF coordinates. In addition, the numerical abscissas are larger than
$0$ that confirms the the behaviour of $G(v)$ at initial stage, which is
consistent with the analysis of $E[\varphi]$, since there is a relationship
between them that the numerical abscissa equals to the time derivative of
$G(v)$ at $v=0$ [see Eq.(45)]. The case is essentially different with the
hyperboloidal case Carballo and Withers (2024). Then we explain the
relationship between the transient effect and waveform, which is solved by
Eq.(34) with a initial Gaussian wave through the numerical method. The
transient effect enables the amplitude of the wave received at null infinity
to exceed the initially given amplitude. Although the wide of Guassian wave
packet will suppress the duration of waveform observed, the initial conditions
of the wave packet do not change the qualitative fact that the amplitude is
amplified at null infinity.

There are some things that can be further explored. Here, we list some
investigations that can be done in the future in terms of pseudospectrum for
the gravitational physics.

  1. 1.

In the general relativity, the positive definiteness of the effective
potential outside the event horizon guarantees the positive definiteness of
the Gram matrix $\mathbf{G}$, so that the Cholesky decomposition can always be
done. However, in some modified gravity theories, this is not always the case.
The discussion of pseudospectra in cases with the large potential well
occurring will be done in the future.

  2. 2.

In the definition of (complex) pseudospectrum, the elements of the
perturbation matrix $\mathbf{E}$ are generally complex. However, considering
the physical reality, the elements of matrix $\mathbf{E}$ should be real,
which imposes an important restriction on pseudospectrum. In this case, the
pseudospectrum is referred to as the real pseudospectrum, and its structure
differs significantly from that of the complex pseudospectrum mentioned
earlier. Then, extending to a more general case, we have the following
definition of structured $\epsilon$-pseudospectrum Trefethen and Embree
(2005):

$\displaystyle\bigcup_{\mathbf{E}\,\text{structured}\,,\lVert\mathbf{E}\rVert<\epsilon}\sigma(\mathbf{H}+\mathbf{E})\,.$
(50)

More work regarding structured pseudospectra can be referenced in Graillat
(2006); Noschese and Reichel (2016). In gravitational physics, exploring the
differences between complex pseudospectra and structured pseudospectra is
worthwhile, and such work will be completed in the future.

## Acknowledgement

We are grateful to Li-Ming Cao, Yu-Sen Zhou, Long-Yue Li and Xia-Yuan Liu for
helpful discussions. This work is supported in part by the National Key
Research and Development Program of China Grant No.2020YFC2201501, in part by
the National Natural Science Foundation of China under Grant No.12075297 and
No.12235019.

##  Appendix A The pseudospectrum of the generalized eigenvalue problems

In this appendix, we will follow the monograph Trefethen and Embree (2005) and
give some definitions of the pseudospectrum of the generalized eigenvalue
problems. Here, the generalized eigenvalue problems is given by999For our
convention, the matrix $\mathbf{A}$ is chosen as $\mathbf{L}_{1}/i$, and the
matrix $\mathbf{B}$ is chosen as $\mathbf{L}_{2}$, where $\mathbf{L}_{1}$ and
$\mathbf{L}_{2}$ are the matrix version of the operator $L_{1}$ and $L_{2}$,
respectively.

$\displaystyle(\mathbf{A}-\omega\mathbf{B})\mathbf{v}=0\,.$ (51)

The parameter-dependent matrix $\mathbf{A}-\omega\mathbf{B}$ is known as a
matrix pencil. The matrix $\mathbf{B}$ is called the “mass matrix” and the
matrix $\mathbf{A}$ is called the “stiffness matrix”. Not like the standard
eigenvalue problem, there are many definitions of pseudospectra denoted as
$\sigma_{\epsilon}(\mathbf{A},\mathbf{B})$ for a matrix pencil. We select two
primary definitions as follows.

Definition 1\. If $\mathbf{B}$ is nonsingular, one can define define
$\sigma_{\epsilon}(\mathbf{A},\mathbf{B})=\sigma_{\epsilon}(\mathbf{B}^{-1}\mathbf{A})$.
The boundary of $\sigma_{\epsilon}(\mathbf{A},\mathbf{B})$ is the
$\epsilon^{-1}$ level curve of
$\lVert(z-\mathbf{B}^{-1}\mathbf{A})^{-1}\rVert$. Therefore,
$\sigma_{\epsilon}(\mathbf{A},\mathbf{B})$ is defined as

$\displaystyle\sigma_{\epsilon}(\mathbf{A},\mathbf{B})=\\{z\in\mathbb{C}:\lVert(z-\mathbf{B}^{-1}\mathbf{A})^{-1}\rVert>1/\epsilon\\}\,.$
(52)

Definition 2\. The boundary of $\sigma_{\epsilon}(\mathbf{A},\mathbf{B})$ is
the $\epsilon^{-1}$ level curve of
$\lVert(z\mathbf{B}-\mathbf{A})^{-1}\rVert$. Therefore,
$\sigma_{\epsilon}(\mathbf{A},\mathbf{B})$ is defined as

$\displaystyle\sigma_{\epsilon}(\mathbf{A},\mathbf{B})=\\{z\in\mathbb{C}:\lVert(z\mathbf{B}-\mathbf{A})^{-1}\rVert>1/\epsilon\\}\,.$
(53)

Note that when the matrix $\mathbf{B}=\mathbf{I}$, these two definitions are
equivalent. Similar to the ordinary eigenvalue problem, we still need to
convert the general norm (Here, it is the energy norm $\lVert\cdot\rVert_{E}$
in our case.) to a $2$-norm for calculating the pseudospectrum. If we arrive
at the Gram matrix $\mathbf{G}$, which is Hermitian and positive definite,
then the norm is given by
$\lVert\mathbf{u}\rVert=\lVert\mathbf{W}\mathbf{u}\rVert_{2}$, where the
nonsingular matrix $\mathbf{W}$ is obtained from the Cholesky decomposition of
$\mathbf{G}$, i.e., $\mathbf{G}=\mathbf{W}^{\star}\mathbf{W}$. For the first
definition of the pseudospectrum, we can compute

$\displaystyle\lVert(z-\mathbf{B}^{-1}\mathbf{A})^{-1}\rVert$ $\displaystyle=$
$\displaystyle\sup_{\mathbf{u}\neq\mathbf{0}}\frac{\lVert\mathbf{W}(z-\mathbf{B}^{-1}\mathbf{A})^{-1}\mathbf{u}\rVert_{2}}{\lVert\mathbf{W}\mathbf{u}\rVert_{2}}=\sup_{\mathbf{u}\neq\mathbf{0}}\frac{\lVert\mathbf{W}(z-\mathbf{B}^{-1}\mathbf{A})^{-1}\mathbf{W}^{-1}\mathbf{u}\rVert_{2}}{\lVert\mathbf{W}\mathbf{W}^{-1}\mathbf{u}\rVert_{2}}$
(54) $\displaystyle=$
$\displaystyle\lVert\mathbf{W}(z-\mathbf{B}^{-1}\mathbf{A})^{-1}\mathbf{W}^{-1}\rVert_{2}=\lVert(z-\mathbf{W}\mathbf{B}^{-1}\mathbf{A}\mathbf{W}^{-1})^{-1}\rVert_{2}\,.$

For the second definition of the pseudospectrum, we compute in the similar way

$\displaystyle\lVert(z\mathbf{B}-\mathbf{A})^{-1}\rVert=\lVert\mathbf{W}(z\mathbf{B}-\mathbf{A})^{-1}\mathbf{W}^{-1}\rVert_{2}=\lVert(z\mathbf{W}\mathbf{B}\mathbf{W}^{-1}-\mathbf{W}\mathbf{A}\mathbf{W}^{-1})^{-1}\rVert_{2}\,.$
(55)

Although the pseudospectrum have different definitions for generalized
eigenvalue problems, they share the same spirit. Hence, in this paper, we
chose definition $1$ to calculate the pseudospectrum.

##  Appendix B The deduction of $\partial_{v}E[\varphi]$

In this appendix, we focus on the deduction of $\partial_{v}E[\varphi]$. It
directly can be written as

$\displaystyle\partial_{v}E[\varphi]$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{a}^{b}\mathrm{d}r_{{\star}}\Big{[}\partial_{vr_{\star}}\varphi\partial_{r_{\star}}\bar{\varphi}+\partial_{vr_{\star}}\bar{\varphi}\partial_{r_{\star}}\varphi+V\partial_{v}(\varphi\bar{\varphi})\Big{]}\,,$
(56)

where $\partial_{vr_{\star}}$ represents $\partial_{v}\partial_{r_{\star}}$.
On the other hand, the equations of motion for $\varphi$ and $\bar{\varphi}$
can be written explicitly within $(v,r_{\star})$,

$\displaystyle(\partial^{2}_{r_{\star}}+2\partial_{vr_{\star}}-V)\varphi=0\,,\quad\text{and}\quad(\partial^{2}_{r_{\star}}+2\partial_{vr_{\star}}-V)\bar{\varphi}=0\,.$
(57)

By substituting
$\partial_{vr_{\star}}\varphi=(V-\partial_{r_{\star}}^{2})\varphi/2$ and
$\partial_{vr_{\star}}\varphi=(V-\partial_{r_{\star}}^{2})\bar{\varphi}/2$
into Eq.(56), we are going to get Eq.(31) as follows

$\displaystyle\partial_{v}E[\varphi]$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{a}^{b}\mathrm{d}r_{{\star}}\Big{[}\partial_{vr_{\star}}\varphi\partial_{r_{\star}}\bar{\varphi}+\partial_{vr_{\star}}\bar{\varphi}\partial_{r_{\star}}\varphi+V\partial_{v}(\varphi\bar{\varphi})\Big{]}$
(58) $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{a}^{b}\mathrm{d}r_{{\star}}\Bigg{\\{}\Big{[}\frac{1}{2}(V-\partial_{r_{\star}}^{2})\varphi\Big{]}\partial_{r_{\star}}\bar{\varphi}+\Big{[}\frac{1}{2}(V-\partial_{r_{\star}}^{2})\bar{\varphi}\Big{]}\partial_{r_{\star}}\varphi+V\partial_{v}(\varphi\bar{\varphi})\Bigg{\\}}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\int_{a}^{b}\mathrm{d}r_{\star}\Big{[}-\frac{1}{2}\partial_{r_{\star}}(\partial_{r_{\star}}\varphi\partial_{r_{\star}}\bar{\varphi})+\frac{1}{2}V\partial_{r_{\star}}(\varphi\bar{\varphi})+V\partial_{v}(\varphi\bar{\varphi})\Big{]}$
$\displaystyle=$
$\displaystyle-\frac{1}{4}|\partial_{r_{\star}}\varphi|^{2}\Big{|}^{b}_{a}+\frac{1}{4}V|\varphi|^{2}\Big{|}^{b}_{a}+\int_{a}^{b}\mathrm{d}r_{\star}\Big{(}\frac{1}{2}V\partial_{v}|\varphi|^{2}-\frac{1}{4}\frac{\mathrm{d}V}{\mathrm{d}r_{\star}}|\varphi|^{2}\Big{)}\,,$

where the last line we have use integration by parts, and $|\varphi|^{2}$
represents $\varphi\bar{\varphi}$.

##  Appendix C The numerical approach of pseudospectrum

In this appendix, we give details associated with the discretization method in
order to calculate pseudospectrum. First, we employ the Chebyshev spectral
method to discretize the differential operator $L_{1}$ and $L_{2}$, resulting
in their matrix representations $\mathbf{L}_{1}$ and $\mathbf{L}_{2}$. In
order to better apply the formulas related to the Chebyshev polynomial which
is naturally defined on the interval $[-1,1]$, it is convenient to map the
radial coordinate $z\in[0,1]$ into $x\in[-1,1]$ via a linear transformation

$\displaystyle x=2z-1\,,\quad z=\frac{x+1}{2}\,.$ (59)

The Chebyshev-Lobatto grid is used in discretization, given resolution $N$, it
contains $N+1$ points $\\{x_{j}\\}_{j=0}^{N}$ with

$\displaystyle x_{j}=\cos\Big{(}\frac{j\pi}{N}\Big{)}\,,\quad
j=0,1,\cdots,N\,.$ (60)

The variable $\psi$ is discretized as $\mathbf{\Psi}$ given by

$\displaystyle\mathbf{\Psi}=(\psi(x_{0}),\psi(x_{1}),\cdots,\psi(x_{N-1}),\psi(x_{N}))^{T}\,.$
(61)

Furthermore, in order to calculate the pseudospectrum, we should turn the
inner product (33) into the discretized version, i.e.,

$\displaystyle\langle\psi_{1},\psi_{2}\rangle_{E}=\mathbf{\Psi}_{1}^{\star}\cdot\tilde{\mathbf{G}}^{E}\cdot\mathbf{\Psi}_{2}\,,$
(62)

where the symbol $\star$ means the conjugate transpose of a vector. The matrix
$\tilde{\mathbf{G}}^{E}$ is known as the Gram matrix corresponding to the
inner product. Finally, we need to talk about an important technical note,
i.e., the integral quadrature of product of two functions. Consider that
functions $\psi_{1}$ and $\psi_{2}$ are approximated by $N+1$ Chebyshev
polynomials, a grid of $2N+2$ points is required to achieve numerical
accuracy. On a practical level, it can be achieved by evaluating at a double
resolution $\bar{N}=2N+1$, and then interpolating down to the original
resolution $N$. Essentially, it means that the matrix $\tilde{\mathbf{G}}^{E}$
is given by

$\displaystyle\tilde{\mathbf{G}}^{E}=\mathbf{P}^{T}\cdot\mathbf{G}^{E}\cdot\mathbf{P}\,.$
(63)

The matrix $\mathbf{P}$ with size $(\bar{N}+1)\times(N+1)$ is called the
interpolation matrix, which can be found as follows

$\displaystyle\mathbf{P}_{\bar{i}i}=\frac{1}{c_{i}N}\Big{[}1+\sum_{j=1}^{N}(2-\delta_{j,N})T_{j}(x_{\bar{i}})T_{j}(x_{i})\Big{]}\,,\quad\bar{i}=0,\cdots,\bar{N}\quad\text{and}\quad
i=0,\cdots,N\,,$ (64)

where $T_{j}$ is the Chebyshev polynomial of order $j$, and

$\displaystyle c_{i}=\begin{cases}2\,,&i=0\quad\text{or}\quad N\,,\\\
1\,,&\text{otherwise}.\end{cases}$ (65)

Moreover, the matrix $\mathbf{G}^{E}$ with size $(\bar{N}+1)\times(\bar{N}+1)$
in Eq.(63) is expressed as

$\mathbf{G}^{E}=\mathbf{D}^{T}_{\bar{N}}\cdot\mathbf{C}_{2}\cdot\mathbf{D}_{\bar{N}}+\frac{1}{2}\Big{(}\mathbf{C}_{1}\cdot\mathbf{D}_{\bar{N}}+\mathbf{D}^{T}_{\bar{N}}\cdot\mathbf{C}_{1}\Big{)}+\frac{1}{4}\mathbf{C}_{V}\,,$
(66)

where the matrices $\mathbf{C}_{1}$, $\mathbf{C}_{2}$ and $\mathbf{C}_{V}$ are
all diagonal. Their diagonal elements are

$\displaystyle(\mathbf{C}_{l})_{\bar{i}\bar{i}}=\frac{2\mu_{l}(x_{\bar{i}})}{c_{\bar{i}}\bar{N}}\Big{[}1-\sum_{k=1}^{\lfloor\frac{\bar{N}}{2}\rfloor}T_{2k}(x_{\bar{i}})\frac{2-\delta_{2k,\bar{N}}}{4k^{2}-1}\Big{]}\,,\quad\bar{i}=0,\cdots,\bar{N}\quad\text{and}\quad
l=1\,,2\,,V\,,$ (67)

where $\lfloor a\rfloor$ is the floor function, i.e., the largest integer that
is less than or equal to $a$. Three weight functions $\mu_{l}$ are

$\displaystyle\mu_{2}(x)$ $\displaystyle=$
$\displaystyle\Big{(}\frac{1+x}{2}\Big{)}^{2\lambda+2}f\Big{(}\frac{1+x}{2}\Big{)}\,,$
$\displaystyle\mu_{1}(x)$ $\displaystyle=$
$\displaystyle\lambda\Big{(}\frac{1+x}{2}\Big{)}^{2\lambda+1}f\Big{(}\frac{1+x}{2}\Big{)}\,,$
$\displaystyle\mu_{V}(x)$ $\displaystyle=$
$\displaystyle\Big{(}\frac{1+x}{2}\Big{)}^{2\lambda-2}\Bigg{\\{}\lambda^{2}\Big{(}\frac{1+x}{2}\Big{)}^{2}f\Big{(}\frac{1+x}{2}\Big{)}+\frac{V[(1+x)/2)]}{f[(1+x)/2]}\Bigg{\\}}\,,$
(68)

respectively. Combined Appendix.A, we have obtained the necessities for the
pseudospectrum of a generalised eigenvalue problem so far. More specifically,
the matrix $\mathbf{G}$ appearing at Appendix.A is nothing but the matrix
$\tilde{\mathbf{G}}^{E}$.

##  Appendix D The convergence test of QNM spectra

A simple way to observe the convergence of eigenvalues is to plot the
eigenvalues under different resolutions on the same graph. If the eigenvalues
are indeed convergent, then these eigenvalues will be highly overlapped under
different resolutions. In Fig.13, we compare the eigenvalues with the
benchmark parameters but from two different resolutions with $N_{1}=90$ and
$N_{2}=100$. The convergence of QNM spectra for the first few eigenvalues can
be demonstrated from Fig.13.

However, in this appendix, we will consider the convergence of eigenvalues
from a more sophisticated perspective to avoid misjudgment caused by eyes.
According to the Eigenvalue Rule-of-Thumb, numerous eigenvalues of the
discretized matrix often bear no meaningful relation to the eigenvalues of the
fundamental differential eigen problem, amounting to mere numerical artifacts
Boyd  _et al._ (2000). Motivated by previous statement, we give a convergence
test of the QNM spectra in this appendix in the wake of Boyd  _et al._ (2000);
Boyd (1996).

How can we distinguish between the “good” and “bad” eigenvalues? Considering
the feature of the QNM problem, the spectra are symmetric about the imaginary
axis. Therefore, we only need to care about the convergence of the eigenvalues
on one side of the imaginary axis, including the eigenvalues on the imaginary
axis. Only those eigenvalues that exhibit minimal differences or “resolution-
dependent drift” can be considered reliable. For two different resolutions
$N_{1}$ and $N_{2}$ with $N_{1}<N_{2}$, we use $\bar{n}_{1}$ and $\bar{n}_{2}$
to represent the number of eigenvalues on the right side of the imaginary axis
at different resolutions, respectively, and we count the numbers of
eigenvalues according to the size of the value of their imaginary parts. For
lower resolution $N_{1}$, it is recommend to define the two kinds of drifts
which are defined as

$\displaystyle
r_{\bar{n},\text{ordinal}}=\frac{\min(|\omega_{\bar{n}}|,\sigma_{\bar{n}})}{\delta_{\bar{n},\text{ordinal}}}\,,\quad\bar{n}=1,2,\cdots,\bar{n}_{1}-1,\bar{n}_{1}\,,$
(69) $\displaystyle
r_{\bar{n},\text{nearest}}=\frac{\min(|\omega_{\bar{n}}|,\sigma_{\bar{n}})}{\delta_{\bar{n},\text{nearest}}}\,,\quad\bar{n}=1,2,\cdots,\bar{n}_{1}-1,\bar{n}_{1}\,,$
(70)

in which the weight function $\sigma_{\bar{n}}$ is defined as

$\displaystyle\sigma_{\bar{n}}=\begin{cases}|\omega_{1}-\omega_{2}|\,,&\bar{n}=1\,,\\\
|\omega_{\bar{n}_{1}-1}-\omega_{\bar{n}_{1}}|\,,&\bar{n}=\bar{n}_{1}\,,\\\
\frac{1}{2}\Big{(}|\omega_{\bar{n}-1}-\omega_{\bar{n}}|+|\omega_{\bar{n}+1}-\omega_{\bar{n}}|\Big{)}\,,&\bar{n}=2,\cdots,\bar{n}_{1}-1\,.\end{cases}$
(71)

The difference between the two definitions origins from the denominator,
depending on what pairs of eigenvalues in different resolutions should be
compared, the first definition is named as “ordinal drift”

$\displaystyle\delta_{\bar{n},\text{ordinal}}=|\omega_{\bar{n}}^{(N_{1})}-\omega_{\bar{n}}^{(N_{2})}|\,,$
(72)

and the second definition is named as “nearest neighbor drift”

$\displaystyle\delta_{\bar{n},\text{nearest}}=\min_{k=1,\cdots,\bar{n}_{2}}|\omega_{\bar{n}}^{(N_{1})}-\omega_{k}^{(N_{2})}|\,,$
(73)

where $\omega_{j}^{(N)}$ is the $j$-th eigenvalue extracted from $\bar{n}_{1}$
or $\bar{n}_{2}$ eigenvalues.

When the eigenvalues are accurate enough, the two ratios would be very large.
In Fig.14, we show the two ratios (69) and (70) with two resolutions
$N_{1}=90$ and $N_{2}=100$. Note that the scale of the vertical axis has been
set to the logarithmic value of the drift. The advantage of the drift ratio
measure lies in its resilience against spurious modes stemming from numerical
errors, as such errors are unlikely to persist near meaningful eigenvalues as
the resolution varies.

Figure 13: The discrete eigenvalues calculated with benchmark parameters in
the right side of the imaginary axis with different resolutions $N_{1}=90$
(orange) and $N_{2}=100$ (blue).

Figure 14:  The eigenvalues are computed with benchmark parameters. The left
image represents the result of the nearest neighbor drift, and the right image
represents the result of the ordinal drift. We have taken the natural
logarithm of $\delta_{\bar{n},\text{ordinal}}$ and
$\delta_{\bar{n},\text{nearest}}$.

##  Appendix E The time-symmetric numerical methods

In this appendix, following the Refs.Markakis  _et al._ (2019); O’Boyle  _et
al._ (2022); Markakis  _et al._ (2023), we are going to solve the coupled ODEs
(34) from the methods of lines for PDE (20). We start from the following
linear differential equations

$\displaystyle\frac{\mathrm{d}\mathbf{\Psi}}{\mathrm{d}v}=\mathbf{H}\mathbf{\Psi}\,,$
(74)

where vector $\mathbf{\Psi}$ is the discrete approximation of the field $\psi$
and $\mathbf{H}$ is time-independent. Using the fundamental theorem of
calculus, Eq.(74) tells us that $\mathbf{\Psi}(v)$ at time $v_{n}$ and
$\mathbf{\Psi}(v)$ at time $v_{n+1}$ have the following relation

$\displaystyle\mathbf{\Psi}(v_{n+1})=\mathbf{\Psi}(v_{n})+\mathbf{H}\int_{v_{n}}^{v_{n+1}}\mathbf{\Psi}(v)\mathrm{d}v\,.$
(75)

The classical Runge-Kutta method for ODEs is based on $1$-point Taylor
expansion. The so-called time-symmetric numerical methods is based on
$2$-point Taylor expansion. In fact, the integral in Eq.(75) is computed by
using the Lanczos-Dyche formula:

$\displaystyle\int_{v_{n}}^{v_{n+1}}\mathbf{\Psi}(v)\mathrm{d}v=\sum_{m=1}^{l}c_{lm}(\Delta
v)^{m}\Big{[}\mathbf{\Psi}^{(m-1)}(v_{n})+(-1)^{m-1}\mathbf{\Psi}^{(m-1)}(v_{n+1})\Big{]}+\mathbf{R}_{l}(\mathbf{\Psi})$
(76)

with coefficients given by

$\displaystyle c_{lm}=\frac{l!(2l-m)!}{m!(2l)!(l-m)!}$ (77)

and the remainder term $\mathbf{R}_{l}(\mathbf{\Psi})$ given by

$\displaystyle\mathbf{R}_{l}(\mathbf{\Psi})=(-1)^{n}\frac{(l!)^{2}}{(2n+1)!(n-l)!}(\Delta
v)^{2n+1}\mathbf{\Psi}^{(2n)}(v)\,,\quad v\in[v_{n},v_{n+1}]\,.$ (78)

For convenience, $\mathbf{\Psi}(v_{n})$ is set to be $\mathbf{\Psi}^{n}$.
Furthermore, by ignoring the remainder $\mathbf{R}_{l}(\mathbf{\Psi})$ of
Eq.(76), one will grasp at the generalized Hermite rule.

For $l=2$ used in this study, Eq.(76) yields the Hermite rule, which is
accurate to fourth order. Then, a fourth scheme for sloving Eq.(74) is found
which means that

$\displaystyle\mathbf{\Psi}^{n+1}=\Big{[}\mathbf{I}-\frac{\Delta
v}{2}\mathbf{H}+\frac{(\Delta
v)^{2}}{12}\mathbf{H}^{2}\Big{]}^{-1}\Big{[}\mathbf{I}+\frac{\Delta
v}{2}\mathbf{H}+\frac{(\Delta
v)^{2}}{12}\mathbf{H}^{2}\Big{]}\mathbf{\Psi}^{n}\,.$ (79)

Here, the time derivatives $\dot{\mathbf{\Psi}}$ has been replaced by
$\mathbf{H}\mathbf{\Psi}$. In order to reduce round-off errors at each time
step, the above scheme is transformed into a more suitable form of numerical
calculation

$\displaystyle\mathbf{\Psi}^{n+1}=\mathbf{\Psi}^{n}+\Big{[}\mathbf{I}-\frac{\Delta
v}{2}\mathbf{H}\Big{(}\mathbf{I}-\frac{\Delta
v}{6}\mathbf{H}\Big{)}\Big{]}^{-1}(\Delta v\mathbf{H})\mathbf{\Psi}^{n}\,.$
(80)

For the choice of $l\geq 3$, the numerical scheme can be found in Markakis
_et al._ (2019); O’Boyle  _et al._ (2022); Markakis  _et al._ (2023).

## References

  * Abbott  _et al._ (2016a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016a), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ (2016b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 241103 (2016b), arXiv:1606.04855 [gr-qc] .
  * Abbott  _et al._ (2017a) B. P. Abbott _et al._ (LIGO Scientific, VIRGO), Phys. Rev. Lett. 118, 221101 (2017a), [Erratum: Phys.Rev.Lett. 121, 129901 (2018)], arXiv:1706.01812 [gr-qc] .
  * Abbott  _et al._ (2017b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017b), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ (2023) R. Abbott _et al._ (KAGRA, VIRGO, LIGO Scientific), Phys. Rev. X 13, 041039 (2023), arXiv:2111.03606 [gr-qc] .
  * Amaro-Seoane  _et al._ (2017) P. Amaro-Seoane _et al._ (LISA), (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Luo  _et al._ (2016) J. Luo _et al._ (TianQin), Class. Quant. Grav. 33, 035010 (2016), arXiv:1512.02076 [astro-ph.IM] .
  * Gong  _et al._ (2021) Y. Gong, J. Luo, and B. Wang, Nature Astron. 5, 881 (2021), arXiv:2109.07442 [astro-ph.IM] .
  * Hu and Wu (2017) W.-R. Hu and Y.-L. Wu, Natl. Sci. Rev. 4, 685 (2017).
  * Kokkotas and Schmidt (1999) K. D. Kokkotas and B. G. Schmidt, Living Rev. Rel. 2, 2 (1999), arXiv:gr-qc/9909058 .
  * Berti  _et al._ (2009) E. Berti, V. Cardoso, and A. O. Starinets, Class. Quant. Grav. 26, 163001 (2009), arXiv:0905.2975 [gr-qc] .
  * Konoplya and Zhidenko (2011) R. A. Konoplya and A. Zhidenko, Rev. Mod. Phys. 83, 793 (2011), arXiv:1102.4014 [gr-qc] .
  * Franchini and Völkel (2023) N. Franchini and S. H. Völkel, (2023), arXiv:2305.01696 [gr-qc] .
  * Barausse  _et al._ (2014) E. Barausse, V. Cardoso, and P. Pani, Phys. Rev. D 89, 104059 (2014), arXiv:1404.7149 [gr-qc] .
  * Cannizzaro  _et al._ (2024) E. Cannizzaro, T. F. M. Spieksma, V. Cardoso, and T. Ikeda, (2024), arXiv:2405.05315 [gr-qc] .
  * Nollert (1996) H.-P. Nollert, Phys. Rev. D 53, 4397 (1996), arXiv:gr-qc/9602032 .
  * Nollert and Price (1999) H.-P. Nollert and R. H. Price, J. Math. Phys. 40, 980 (1999), arXiv:gr-qc/9810074 .
  * Daghigh  _et al._ (2020) R. G. Daghigh, M. D. Green, and J. C. Morey, Phys. Rev. D 101, 104009 (2020), arXiv:2002.07251 [gr-qc] .
  * Qian  _et al._ (2021) W.-L. Qian, K. Lin, C.-Y. Shao, B. Wang, and R.-H. Yue, Phys. Rev. D 103, 024019 (2021), arXiv:2009.11627 [gr-qc] .
  * Hirano  _et al._ (2024) S. Hirano, M. Kimura, M. Yamaguchi, and J. Zhang, (2024), arXiv:2404.09672 [gr-qc] .
  * Cardoso  _et al._ (2024) V. Cardoso, S. Kastha, and R. Panosso Macedo, (2024), arXiv:2404.01374 [gr-qc] .
  * Courty  _et al._ (2023) A. Courty, K. Destounis, and P. Pani, Phys. Rev. D 108, 104027 (2023), arXiv:2307.11155 [gr-qc] .
  * Berti  _et al._ (2022) E. Berti, V. Cardoso, M. H.-Y. Cheung, F. Di Filippo, F. Duque, P. Martens, and S. Mukohyama, Phys. Rev. D 106, 084011 (2022), arXiv:2205.08547 [gr-qc] .
  * Cheung  _et al._ (2022) M. H.-Y. Cheung, K. Destounis, R. P. Macedo, E. Berti, and V. Cardoso, Phys. Rev. Lett. 128, 111103 (2022), arXiv:2111.05415 [gr-qc] .
  * Rosato  _et al._ (2024) R. F. Rosato, K. Destounis, and P. Pani, (2024), arXiv:2406.01692 [gr-qc] .
  * Oshita  _et al._ (2024) N. Oshita, K. Takahashi, and S. Mukohyama, (2024), arXiv:2406.04525 [gr-qc] .
  * Torres (2023) T. Torres, Phys. Rev. Lett. 131, 111401 (2023), arXiv:2304.10252 [gr-qc] .
  * Cao  _et al._ (2024) L.-M. Cao, J.-N. Chen, L.-B. Wu, L. Xie, and Y.-S. Zhou, (2024), arXiv:2401.09907 [gr-qc] .
  * Li  _et al._ (2024) G.-R. Li, W.-L. Qian, and R. G. Daghigh, (2024), arXiv:2406.10782 [gr-qc] .
  * Trefethen and Embree (2005) L. Trefethen and M. Embree,  _Spectra and Pseudospectra: The Behavior of Nonnormal Matrices and Operators_ (Princeton university press, 2005).
  * Jaramillo  _et al._ (2021) J. L. Jaramillo, R. Panosso Macedo, and L. Al Sheikh, Phys. Rev. X 11, 031003 (2021), arXiv:2004.06434 [gr-qc] .
  * Destounis and Duque (2023) K. Destounis and F. Duque (2023) arXiv:2308.16227 [gr-qc] .
  * Jaramillo  _et al._ (2022) J. L. Jaramillo, R. Panosso Macedo, and L. A. Sheikh, Phys. Rev. Lett. 128, 211102 (2022), arXiv:2105.03451 [gr-qc] .
  * Trefethen  _et al._ (1993) L. N. Trefethen, A. E. Trefethen, S. C. Reddy, and T. A. Driscoll, Science 261, 578 (1993).
  * Destounis  _et al._ (2021) K. Destounis, R. P. Macedo, E. Berti, V. Cardoso, and J. L. Jaramillo, Phys. Rev. D 104, 084091 (2021), arXiv:2107.09673 [gr-qc] .
  * Areán  _et al._ (2023) D. Areán, D. G. Fariña, and K. Landsteiner, JHEP 12, 187 (2023), arXiv:2307.08751 [hep-th] .
  * Cownden  _et al._ (2024) B. Cownden, C. Pantelidou, and M. Zilhão, JHEP 05, 202 (2024), arXiv:2312.08352 [gr-qc] .
  * Boyanov  _et al._ (2024) V. Boyanov, V. Cardoso, K. Destounis, J. L. Jaramillo, and R. Panosso Macedo, Phys. Rev. D 109, 064068 (2024), arXiv:2312.11998 [gr-qc] .
  * Sarkar  _et al._ (2023) S. Sarkar, M. Rahman, and S. Chakraborty, Phys. Rev. D 108, 104002 (2023), arXiv:2304.06829 [gr-qc] .
  * Destounis  _et al._ (2024) K. Destounis, V. Boyanov, and R. Panosso Macedo, Phys. Rev. D 109, 044023 (2024), arXiv:2312.11630 [gr-qc] .
  * Boyanov  _et al._ (2023) V. Boyanov, K. Destounis, R. Panosso Macedo, V. Cardoso, and J. L. Jaramillo, Phys. Rev. D 107, 064012 (2023), arXiv:2209.12950 [gr-qc] .
  * Jaramillo (2022) J. L. Jaramillo, Class. Quant. Grav. 39, 217002 (2022), arXiv:2206.08025 [gr-qc] .
  * Carballo and Withers (2024) J. Carballo and B. Withers, (2024), arXiv:2406.06685 [hep-th] .
  * Maeda and Dadhich (2006) H. Maeda and N. Dadhich, Phys. Rev. D 74, 021501 (2006), arXiv:hep-th/0605031 .
  * Cao and Wu (2021) L.-M. Cao and L.-B. Wu, Phys. Rev. D 103, 064054 (2021), arXiv:2101.02461 [gr-qc] .
  * Cao  _et al._ (2023) L.-M. Cao, L.-B. Wu, Y. Zhao, and Y.-S. Zhou, Phys. Rev. D 108, 124023 (2023).
  * Dotti and Gleiser (2005) G. Dotti and R. J. Gleiser, Class. Quant. Grav. 22, L1 (2005), arXiv:gr-qc/0409005 .
  * Camporesi and Higuchi (1994) R. Camporesi and A. Higuchi, J. Math. Phys. 35, 4217 (1994).
  * Buell and Shadwick (1995) W. F. Buell and B. A. Shadwick, American Journal of Physics 63, 256 (1995).
  * Zenginoglu (2008) A. Zenginoglu, Class. Quant. Grav. 25, 145002 (2008), arXiv:0712.4333 [gr-qc] .
  * Zenginoglu (2011) A. Zenginoglu, Phys. Rev. D 83, 127502 (2011), arXiv:1102.2451 [gr-qc] .
  * Zenginoğlu (2024) A. Zenginoğlu, (2024), arXiv:2404.01528 [gr-qc] .
  * Panosso Macedo (2024) R. Panosso Macedo, Phil. Trans. Roy. Soc. Lond. A 382, 20230046 (2024), arXiv:2307.15735 [gr-qc] .
  * Gasperin and Jaramillo (2022) E. Gasperin and J. L. Jaramillo, Class. Quant. Grav. 39, 115010 (2022), arXiv:2107.12865 [gr-qc] .
  * Trefethen (1999) L. N. Trefethen, Acta Numerica 8, 247–295 (1999).
  * Cardoso  _et al._ (2018) V. Cardoso, J. a. L. Costa, K. Destounis, P. Hintz, and A. Jansen, Phys. Rev. Lett. 120, 031103 (2018), arXiv:1711.10502 [gr-qc] .
  * Burke  _et al._ (2003) J. V. Burke, A. S. Lewis, and M. L. Overton, IMA Journal of Numerical Analysis 23, 359 (2003).
  * Markakis  _et al._ (2019) C. M. Markakis, M. F. O’Boyle, D. Glennon, K. Tran, P. Brubeck, R. Haas, H.-Y. Schive, and K. Uryū, (2019), arXiv:1901.09967 [math.NA] .
  * O’Boyle  _et al._ (2022) M. F. O’Boyle, C. Markakis, L. J. G. Da Silva, R. Panosso Macedo, and J. A. V. Kroon, (2022), arXiv:2210.02550 [gr-qc] .
  * Markakis  _et al._ (2023) C. Markakis, S. Bray, and A. Zenginoğlu, (2023), arXiv:2303.08153 [gr-qc] .
  * Baibhav  _et al._ (2023) V. Baibhav, M. H.-Y. Cheung, E. Berti, V. Cardoso, G. Carullo, R. Cotesta, W. Del Pozzo, and F. Duque, Phys. Rev. D 108, 104020 (2023), arXiv:2302.03050 [gr-qc] .
  * Graillat (2006) S. Graillat, Journal of Computational and Applied Mathematics 191, 68 (2006).
  * Noschese and Reichel (2016) S. Noschese and L. Reichel, (2016), arXiv:1611.04415 [math.NA] .
  * Boyd  _et al._ (2000) J. P. Boyd, T. Marilyn, and P. T. S. Eliot,  _Chebyshev and Fourier Spectral Methods_ (2000).
  * Boyd (1996) J. P. Boyd, Journal of Computational Physics 126, 11 (1996).