Pseudospectrum and black hole quasi-normal mode (in)stability

We study the stability of quasinormal modes (QNM) in asymptotically flat black hole spacetimes by means of a pseudospectrum analysis. The construction of the Schwarzschild QNM pseudospectrum reveals the following: (i) the stability of the slowest-decaying QNM under perturbations respecting the asymptotic structure, reassessing the instability of the fundamental QNM discussed by Nollert [H. P. Nollert, About the Significance of Quasinormal Modes of Black Holes, Phys. Rev. D 53, 4397 (1996)] as an"infrared"effect; (ii) the instability of all overtones under small-scale ("ultraviolet") perturbations of sufficiently high frequency, which migrate towards universal QNM branches along pseudospectra boundaries, shedding light on Nollert's pioneer work and Nollert and Price's analysis [H. P. Nollert and R. H. Price, Quantifying Excitations of Quasinormal Mode Systems, J. Math. Phys. (N.Y.) 40, 980 (1999)]. Methodologically, a compactified hyperboloidal approach to QNMs is adopted to cast QNMs in terms of the spectral problem of a non-self-adjoint operator. In this setting, spectral (in)stability is naturally addressed through the pseudospectrum notion that we construct numerically via Chebyshev spectral methods and foster in gravitational physics. After illustrating the approach with the P\"oschl-Teller potential, we address the Schwarzschild black hole case, where QNM (in)stabilities are physically relevant in the context of black hole spectroscopy in gravitational-wave physics and, conceivably, as probes into fundamental high-frequency spacetime fluctuations at the Planck scale.

We study the stability of quasi-normal modes (QNM) in asymptotically flat black hole spacetimes by means of a pseudospectrum analysis. The construction of the Schwarzschild QNM pseudospectrum reveals: i) the stability of the slowest decaying QNM under perturbations respecting the asymptotic structure, reassessing the instability of the fundamental QNM discussed by Nollert [1] as an "infrared" effect; ii) the instability of all overtones under small scale ("ultraviolet") perturbations of sufficiently high frequency, that migrate to a universal class of QNM branches along pseudospectra boundaries, shedding light on Nollert & Price's analysis [1,2]. Methodologically, a compactified hyperboloidal approach to QNMs is adopted to cast the QNM calculation as the spectral problem of a non-selfadjoint operator. In such a setting, spectral (in)stability is naturally addressed through the notion of pseudospectrum, that we construct numerically via Chebyshev spectral methods and foster in gravity physics. After illustrating the approach with the Pöschl-Teller potential, we address the Schwarzschild black hole case, where QNM (in)stabilities are physically relevant in the context of black hole spectroscopy in gravitational wave physics and, conceivably, as probes into fundamental high-frequency spacetime fluctuations at the Planck scale.

I. INTRODUCTION: QNMS AND (IN)STABILITY
A. The black hole QNM stability problem: non-selfadjoint spectral (in)stability and the pseudospectrum Structural stability is essential in the modelling and understanding of physical phenomena. In the context of spectral problems pervading physics, the use of 'pseudospectra' [3][4][5][6] is well spread whenever stability issues of non-conservative systems are addressed, from pioneering applications in hydrodynamics [3] to recent advances [7] covering a wide range in physics. Specifically, the notion of pseudospectrum is devised to assess the spectral (in)stability of systems controlled by non-selfadjoint operators. Such is the case of black hole (BH) QNMs, where outgoing boundary conditions define a leaky system whose resonant frequencies can be cast as the eigenvalues of a certain non-selfadjoint operator.
The problem we address here is the spectral robustness of BH QNMs, namely the stability of the resonant frequencies of BH spacetimes under perturbations. Being an intrinsic property of the background, QNM frequencies encode crucial geometric information about BHs and their environment. Thus, they have become a fundamental tool in astrophysics, fundamental gravitational physics, and mathematical relativity in their attempts to probe spacetime geometry through scattering methods. Beyond gravity, QNMs play currently a key role in other areas of physics, as perfectly illustrated by recent breakthroughs in optical nanoresonator QNMs [8,9]. The ultimate significance of QNM frequencies depends directly on the understanding and control of their spectral stability.

B. BH QNM instability: Nollert & Price's pionneer work
Nollert's groundbreaking work in BH QNM spectral stability [1], complemented by the analysis in Nollert & Price's [2], showed evidence of an overall instability of the Schwarzschild QNM spectrum under a class of small scale perturbations. In these works, using a family of step-like approximations to the Schwarzschild curvature potential, it was demonstrated: i) QNM overtones are strongly instable, their instability increasing with damping.
iii) The black hole ringdowns, at late intermediate times, according to the non-perturbed fundamental mode. Only at very late times the ringdown frequency is controlled by the perturbed fundamental QNM mode.
Beyond Nollert & Price's works, research in BH QNM spectral (in)stability has been further pursued in different gravitational physics settings. In astrophysics, the understanding of possible environmental observational signatures in "dirty" BH scenarios has prompted research [10,11] significantly intensified in recent times [12][13][14]. Indeed, in the era of gravitational-wave astronomy, the stability of QNM overtones is paramount for BH spectroscopy [15][16][17][18][19][20][21]. On the other hand, regarding investigations on the fundamental structure of spacetime, the perspective of accessing quantum scales through high-frequency instabilities of QNM overtones has also tantalized the research [22][23][24][25][26][27]. In spite of these efforts, a comprehensive picture of BH QNM (in)stability seems to be lacking. In particular, the stability status of the slowest decaying QNM remains unclear, whereas the elucidation of the lowest overtone subject to high-frequency instability is an open problem. These two points, key in general BH physics, are actually urgent in gravitational wave astrophysics. The implementation of a pseudospectrum analysis (introducing it, up to our knowlegde, in a gravitational setting) permits to address systematically these questions and to provide sound answers to points i) and ii) above. Regarding point iii), some important conclusions can be drawn from the stability properties of the fundamental QNM mode, but the systematic analysis of the detailed relation between QNM frequencies and BH ringdown frequencies lays beyond the scope of the present work. The calculation of BH QNMs has been the subject of systematic study in gravitational physics and there exists a variety of standard approaches to address this problem (cf. e.g. [28][29][30][31]). From a methodological perspective, our discussion lies on two main ingredients at a conceptual level: i) A hyperboloidal approach to QNMs: this geometric approach casts the QNM calculation as a proper eigenvalue problem of a non-selfadjoint operator.
ii) Pseudospectrum: this notion, together with related spectral tools, provides the key tool to study the potential spectral instability of non-selfadjoint operators.
The combination of these two elements permits to develop a systematic treatment of the problem. At an exploratory stage, prior to the full analytical study, the present work addresses pseudospectra in a numerical approach. This sets a challenging problem demanding high accuracy. The third key ingredient in our approach is the use of spectral numerical methods.
The article is structured as follows. Section II reviews the hyperboloidal approach to scattering problems, with a focus on QNMs. Section III introduces the basic elements to study spectral instability of non-selfadjoint operators, in particular the notion of pseudospectrum. Section IV presents the numerical spectral tools to be employed in the present approach. Then, section V illustrates all the previous elements in the toymodel provided by the Pöschl-Teller potential, that also anticipates some of the main results in the BH setting. Section VI contains the main contribution in the present work, namely the construction of the Schwarzschild QNM pseudospectrum and the consequent analysis of BH QNM (in)stability. Conclusions and perspectives are finally presented in section VII.

II. HYPERBOLOIDAL APPROACH TO QNMS
The adopted (compactified) hyperboloidal approach provides a geometric framework to study QNMs, that characterizes resonant frequencies in terms of an eigenvalue problem [32][33][34][35][36][37][38][39][40][41][42]. The scheme geometrically imposes QNM outgoing boundary conditions by adopting a spacetime slicing that intersects future null infinity I + and, in the BH setting, penetrates the horizon. Since light cones point outwards at the boundary of the domain, the boundary conditions are automatically imposed for propagating physical degrees of freedom.

A. Wave equation in the compactified hyperboloidal approach
We focus on the scattering problem of (massless) linear fields on stationary spherically symmetric BH backgrounds. For concreteness, let us first consider a scalar field Φ, satisfying the wave equation We adopt standard Schwarzschild coordinates ds 2 = −f (r)dt 2 + f (r) −1 dr 2 + r 2 dθ 2 + sin 2 θdϕ 2 , and emphasize that t = constant slices correspond to Cauchy hypersurfaces intersecting both the horizon bifurcation sphere and spatial infinity i o . If we consider the rescaling then Eq. (1) rewrites, expanding φ in spherical harmonics with φ m modes and using the tortoise coordinate defined by dr dr * = f (r) with the appropriate integration constant, as where now r * ∈ ] − ∞, ∞[. Remarkably, when considering electromagnetic and (linearized) gravitational fields, the respective geometric wave equations corresponding to Eq. (1) can be cast in the form of Eq. (4) for appropriate effective scalar potentials. Specifically, a scalar field is introduced for each polarization, satisfying Eq. (4) with the suitable potential V . In the gravitational case, the axial polarization is subject to the so-called Regge-Wheeler potential, whereas the polar one is controlled by the Zerilli potential (cf. e.g [28,43,44]). The BH event horizon and (spatial) infinity correspond, respectively, to r * →−∞ and r * →+∞. We extend the domain of r * to [−∞, ∞] and introduce the dimensionless quantities for an appropriate length scale λ to be chosen in each specific setting. More importantly, we consider coordinates (τ, x) that implement the hyperboloidal approach Specifically: i) The height function h(x) implements the hyperboloidal slicing, i.e. τ = constant is a horizon-penetrating hyperboloidal slice Σ τ intersecting future I + .
We note that the compactification is performed only in the spatial direction along the hyperboloidal slice, and not in time, so that the latter can be Fourier transformed. The relevant compactification here is a partial one, and not the total spacetime compactification leading to Carter-Penrose diagrams. The choice of h(x) and g(x) is, as we comment below, subject to certain restrictions. Under transformation (6), Eq. (4) writes where the prime denotes derivative with respect to x. Admittedly, expression (7) appears more intricate than Eq. (4). However, this change encodes a neat geometric structure and, as we shall argue, it plays a crucial role in our construction and discussion of the relevant spectral problem.

B. First-order reduction and spectral problem
The structure in Eq. (7) is made more apparent by performing a first-order reduction in time, by introducing Then, Eq. (7) becomes where the operator L is defined as with and The structure of L 1 is that of a Sturm-Liouville operator. In particular, functions h(x) and g(x) are chosen such that they guarantee the positivity of the weight function w(x), namely w(x) > 0. The operator L 2 has also a neat geometric/analytic structure adapted to the integration by parts, being symmetric in the sense L 2 = 1 w(x) γ(x)∂ x + ∂ x (γ(x)·) .
A key property of coordinate transformation (6) is that it preserves, up to the overall constant λ, the timelike Killing vector t a controlling stationarity In this sense functions t and λτ "tick" at the same pace, namely they are natural parameters of t a , i.e. t a (t) = t a (λτ ) = 1 (the role of the constant λ being just that of keeping proper dimensions). This is crucial for the consistent definition of QNM frequencies by Fourier (or Laplace) transformation from Eqs. (4) and (7), since variables ω respectively conjugate to t and τ then coincide (up to the constant 1/λ).
In other words: the change of time coordinate in (6) does not affect the values of the obtained QNM frequencies.
Performing then the Fourier transform in τ in the first-order form (9) of the wave equation (with standard sign convention for the Fourier modes, u m (τ, x) ∼ u m (x)e iωτ ) we arrive at the spectral problem L u n, m = ω n, m u n, m , (14) or, more explicitly

Regularity and outgoing boundary conditions
As emphasized at the beginning of this section, a major motivation for the adopted hyperboloidal approach is the geometric imposition of outgoing boundary conditions at future null infinity and at the event horizon: being null hypersurfaces with light cones pointing outwards from the integration domain, physical degrees of freedom (as the scalar polarization fields we consider here) do not admit boundary conditions. How does this translate into the analytic scheme resulting from the change of variables (6)?
The key point is that transformation (6) must be such that p(x) in the Sturm-Liouville operator L 1 in Eq. (11) vanishes at the boundaries of the compactified spatial domain [a, b] This will be illustrated explicitly in the study cases discussed later. Then the elliptic operator L 1 is a 'singular' Sturm-Liouville operator, this impacting directly on the boundary conditions it admits. Specifically, if regularity is enforced on eigenfunctions, then L 1 does not admit boundary conditions. Moreover, such absence of boundary conditions extends to the full operator L in the hyperbolic problem (cf. discussion in [45]). In brief: if sufficient regularity is imposed on the space of functions u m , then wave equations (7), (9) and the spectral problem (14) do not admit boundary conditions, as a consequence of the vanishing of p(x) at the boundaries of [a, b]. This is the analytic counterpart of the geometric structure implemented in the compactified hyperboloidal approach. QNM boundary conditions are in-built, as regularity conditions, in the 'bulk' of the operator L in Eq. (14).
C. Scalar product: QNMs as a non-selfadjoint spectral problem The outgoing boundary conditions in the present setting define a leaky system, with a loss of energy through the boundaries -null infinity and the black hole horizon -so that the system is not conservative. This suggests that the infinitesimal generator of the evolution in Eq. (7), namely the operator L, should not be selfadjoint. This requires the introduction of an appropriate scalar product in the problem. Moreover, such identification of the appropriate Hilbert space for solutions is also key for the regularity conditions evoked above.
Eq. (7) describes the evolution of each mode φ m in a background 1+1-Minkowski spacetime with a scattering potential V . A natural scalar product in this reduced problem, both from the physical and the analytical point of view, is given in terms of the energy associated with such scalar field mode. In the context of the spectral problem (14), we must consider generically a complex scalar field φ m , for which the associated stress-energy tensor writes (dropping ( , m) indices) where η ab denotes the Minkowski metric in arbitrary coordinates and "c.c" indicates "complex − conjugate". In a stationary situation, the "total energy" contained in the spatial slice Σ τ and associated with the mode φ is given [46] by where t a is again the timelike Killing vector associated with stationarity and n a denotes the unit timelike normal to the spacelike slice Σ τ . Writing explicitly the energy in the compactified hyperboloidal coordinates (τ, x) introduced in (6), we get where we identify the functions appearing in the definition of the L 1 operator in (11) and (12). In particular, if g 2 − h 2 > 0 (as we have required above) andV > 0 (this is required for positivity of the norm) then, identifying ∂ τ φ = ψ as in (8), we can write the following norm for the vector u in (8) We refer in the following to this norm as the "energy norm". We notice that γ(x) in Eq. (12), associated with L 2 does not enter in the norm, that is, in the energy. This norm comes indeed from a scalar product. Rewriting, for making its role more apparent, the q (x) function as the rescaled potentialṼ and under the assumption aboveṼ > 0, we can introduce the "energy scalar product" for vector functions u in Eq. (8), as and indeed, it holds ||u|| 2 E = u, u E . This will be the relevant scalar product in our discussion.
The full operator L in (14) is not selfadjoint in the scalar product (22). In fact, the first-order operator L 2 stands for a dissipative term encoding the energy leaking at I + and the BH horizon. One could, at a first look, consider that this is related to the first-order character of L 2 , which makes it antisymmetric when integrating by parts with a L 2 ([a, b], w(x)dx) scalar product on ψ, in contrast with the selfadjoint character of the Sturm-Liouville operator L 1 in L 2 ([a, b], w(x)dx) for φ functions. However, this is misleading and actually would suggest a wrong 'bulk' dissipation mechanism. When calculating the formal adjoint L † of the full operator L with the scalar product (22), one gets where L ∂ is an operator with support only on the boundaries of the interval [a, b], that we can formally write as with L ∂ 2 given by the expression where δ(x) formally denotes a Dirac-delta distribution. This is just a formal expression, that underlines precisely the need of a more careful treatment on the involved functional spaces, but it has the virtue of making apparent that the obstruction to selfadjointness lays at the boundaries, as one expects in our QNM problem, and not in the bulk, as one could naively conclude from the presence of a first-order operator L 2 (cf. discussion above): L ∂ 2 explicitly entails a boundary dissipation mechanism. In particular, we note that L is selfadjoint in the non-dissipative L 2 = 0 case, as expected, but that this has required the introduction of quite a non-trivial scalar product.
As a bottomline, in this section we have cast the QNM problem as the eigenvalue problem of a non-selfadjoint operator. In the following section we discuss the implications of this.

III. SPECTRAL STABILITY AND PSEUDOSPECTRUM
The spectrum of a non-selfadjoint operator is potentially unstable under small perturbations of the operator. Let us consider a linear operator A on a Hilbert space with scalar product ·, · , and denote its adjoint by A † , satisfying A † u, v = u, Av . The operator A is called normal if and only if [A, A † ] = 0. In particular, a selfadjoint operator A † = A is normal. In this setting, the 'spectral theorem' (under the appropriate functional space assumptions) states that a normal operator is characterized as being unitarily diagonalizable. The eigenfunctions of A form an orthonormal basis and, crucially in the present discussion, the eigenvalues are stable under perturbations of A. The lack of such a 'spectral theorem' for nonnormal operators entails a severe loss of control on eigenfunction completeness and the potential instability of the spectrum of the operator A. Here, we focus on this second aspect.
withλ i the complex conjugate of λ i . Let us consider, for > 0, the perturbation of A by a (bounded) operator δA The eigenvalues [48] in the perturbed spectral problem satisfy where the first line generalizes [4,49] the expression employed (for self-adjoint operators, where u i = v i ) in quantum mechanics first-order perturbation theory, the first inequality in the second line is Cauchy-Schwartz inequality and in the second inequality we make explicit use of an operator norm || · || induced from that one in the vector Hilbert space, so that ||δAv|| ≤ ||δA||||v||, and ||δA|| = 1 in (27). Then, defining the condition number κ i associated with the eigenvalue λ i , we can write the bound for the perturbation of the eigenvalue λ i In the normal operator case, u i and v i are proportional (namely, since A and A † commute they can be diagonalized in the same basis). Then, again by Cauchy-Schwartz, κ i = 1 and we encounter spectral stability: a small perturbation of order of the operator A entails a perturbation of the same order in the spectrum. In contrast, in the non-normal case, u i and v i are not necessarily collinear. In the absence of a spectral theorem nothing prevents u i and v i to become close to orthogonality and κ i can become very large: small perturbations of A can produce large deviations in the eigenvalues.
The relative values of κ i control the corresponding instability sensitivity of different λ i 's to an operator perturbation [50].

B. Pseudospectrum
A complementary approach to the study of the spectral (in)stability of the operator A under perturbations consists in considering the following questions: Given the operator A and its spectrum σ(A), which is the set of complex numbers λ ∈ C that are actual eigenvalues of "some" small perturbation A + δA, with ||δA|| < ? Does this set extend in C far from the spectrum of A?
In this setting, if we are dealing with an operator that is spectrally stable, we expect that the spectrum of A + δA will not change strongly with respect to that of A, so that the set of λ ∈ C corresponding to the first question above will not be far from σ(A), staying in its vicinity at a maximum distance of order . On the contrary, if we find a tiny perturbation δA of order ||δA|| < such that the corresponding eigenvalues of A + δA actually reach regions in C at distances far apart from σ(A), namely orders of magnitude above , we will conclude that our operator suffers of an actual spectral instability.

Pseudospectrum and operator perturbations
The previous discussion is formalized in the notion of pseudospectrum. Specifically, given > 0, the -pseudospectrum σ (A) of A can be characterised as [51] σ (A) This notion of -pseudospectrum σ (A) is a crucial one in our study of eigenvalue instability since it implies that points in σ (A) are actual eigenvalues of some -perturbation of A: if σ (A) extends far from the spectrum σ(A) for a small , then a small physical perturbation δA of A can produce large actual deviations in the perturbed physical spectrum. The pseudospectrum becomes a systematic tool to assess spectral (in)stability, as illustrated in the hydrodynamics context [3].
Although the characterization (31) of σ (A) neatly captures the notion of (in)stability of A, from a pragmatic perspective it suffers from the drawback of not providing a constructive approach to build such sets σ (A) for different 's (see however subsection III C below, for a further qualification of this question in terms of random perturbation probes).

Pseudospectrum and operator resolvent
To address the explicit construction of pseudospectra, a second characterization of the same set σ (A) in (31) is very useful. Such second characterization is based on the notion of the resolvent R A (λ) = (λId − A) −1 of the operator A.
An eigenvalue λ of A is a complex number that makes singular the operator (λId − A). More generally, the spectrum σ(A) of A is the set {λ ∈ C} for which the resolvent R A (λ) does not exist as a bounded operator (cf. details and subtleties on this notion in e.g. [6,49]). This spectrum concept is a key notion for normal operators but, due to spectral instabilities discussed above, σ(A) is not necessarily the good object to consider for non-normal operators in our context. The notion of -pseudospectrum enters then in scene. Specifically, an equivalent definition of the set σ (A) in (31) is given by the following characterization [4,6] σ (A) = {λ ∈ C : ||R A (λ)|| = ||(λId − A) −1 || > 1/ }. (32) This characterization captures that, for non-normal operators, the norm of the resolvent R A (λ) can be very large far from the spectrum σ(A). This is in contrast with the normal-operator case, where (in the || · || 2 norm) In the non-normal case, one can only guarantee (e.g. [4]) where κ is also a condition number, different but related to the eigenvalue condition numbers κ i in (30) (κ, associated with the matrix diagonalising A, provides an upper bound to the individual κ i 's; see [4] for details). In the non-normal case, κ can become very large and -pseudospectra sets can extend far from the spectrum of A for small values of . The extension of σ (A) far from σ(A) is therefore a signature of strong nonnormality and indicates a poor analytic behavior of R A (λ).
The important point here is that the characterization (34) of σ (A) provides a practical way of calculating σ (A). If we calculate the norm of the resolvent ||R A (λ)|| as a function of λ = Re(λ) + iIm(λ) ∈ C, this provides a real function of two real variables (Re(λ), Im(λ)): the boundaries of the σ (A) sets are just the 'contour lines' of the plot of this function ||R A (λ)||. In particular, -pseudospectra are nested sets in C around the spectrum σ(A), with decreasing towards the 'interior' of such sets and such that lim

Pseudospectrum and quasimodes
For completeness, we provide a third equivalent characterization of the pseudospectrum in the spirit of characterising λ's in the -pseudospectrum set σ (A) as 'approximate eigenvalues' of A, 'up-to an error' , with corresponding 'approximate (right) eigenvectors' v. Specifically, it can be shown [4,6] that σ (A) can be characterised as σ (A) = {λ ∈ C, ∃v ∈ C n : ||Lv − λv|| < } . (35) This characterisation introduces the notion of " -quasimode" v (referred to as "pseudo-mode" in [4]), a key notion in the semiclassical spectral study of A [6]. On the other hand, this third characterization also clearly indicates the numerical difficulty that may occur when trying to determine the actual eigenvalues of A, since round-off errors are unavoidable. This signals the need of a careful treating, whenever addressing numerically the spectral problem of a non-normal operator A.

Pseudospectrum and choice of the norm
In this subsection we have presented the -pseudospectrum as a notion that may be more adapted to the analysis of nonnormal operators than that of the spectrum. We must emphasize however, that the notion of spectrum σ(A) is intrinsic to the operator A, whereas the -pseudospectrum σ (A) is not, since it also depends on the choice of an operator norm. This is crucial, since it determines what we mean by 'big/small' when referring to the perturbation δA, and therefore critically impacts the assessment of stability: a small operator perturbation δA in a given norm, can be a large one when considering another norm. In the first case, from a large variation δλ in the eigenvalues we would conclude instability, whereas in the second case such variation would be consistent with stability.
In this sense, from a mathematical perspective, the study of spectral (in)stability through pseudospectra amounts, in a good measure, to the identification of the proper scalar product determining the norm, that is, to the identification of the proper Hilbert space in which the operator A acts. However, from a physical perspective we might not have such a freedom to choose a mathematically conveniently rescaled norm, since what we mean by large and small may be fixed by the physics of the problem, e.g. by the size of involved amplitudes, intensities or the energy contained in the perturbations. Then, the choice of an appropriate norm, both from a mathematical and physical perspective, is a fundamental step in the analysis. This is the rationale behind the choice of the energy norm || · || E in (20). Once the norm is chosen, the equivalent characterizations (31), (32) and (35) emphasize complementary aspects of the -pseudospectrum notion and σ (A) sets.

C. Pseudospectrum and random perturbations
When considering the construction of pseudospectra, we have presented the characterization of σ (A) in terms of the resolvent R A (λ), Eq. (32), as better suited than the one in terms of spectra of perturbed operators, Eq. (31). The reason is that the former involves only the unperturbed operator A, whereas the latter demands a study of the spectral problem for any perturbed operator A + δA with small δA: a priori, the difficulty to control such space of possible δA perturbations hinders an approach based on such characterisation (31).
But the very nature of the obstacle suggests a possible solution, namely to consider the systematic study of the perturbed spectral problem under random perturbations δA as an avenue to explore -pseudospectra sets. This heuristic expectation actually withstands a more careful analysis and constitutes the basis of a rigorous approach to the analysis of pseudospectra [6]. From a practical perspective, the systematic study of the spectral problem of A + δA with (bounded) random δA with ||δA|| ≤ , has proven to be an efficient tool to explore the 'migration' of eigenvalues through the complex plane (inside the -pseudospectra) [4]. This is complementary to (and 'technically' independent from) the evaluation of σ (A) from the contour-lines of the norm ||R A (λ)|| of the resolvent. Such complementarity of approaches will prove key later in our analysis of Nollert & Price's high-frequency perturbations of the Schwarzschild's potential and the related QNMs.
Two important by-products of this random perturbation approach to the pseudospectrum are the following: i) Random perturbations help identifying instabilitytriggering perturbations: -pseudospectra and condition numbers κ i are efficient in identifying the instability of the spectrum and/or a particular eigenvalue λ i , respectively. However, they do not inform on the specific kind of perturbation actually triggering the instability. This can be crucial to assess the physical nature of the found instability. The use of families of random operators adapted to specific types of perturbations sheds light on this precise point. We will make critical use of this in our assessment of Schwarzschild's (in)stability.
ii) Random perturbations improve analyticity: a remarkable and apparently counter-intuitive effect of random perturbations is the improvement of the analytic behaviour of R A (λ) in λ ∈ C [6]. In particular, the norm ||R A (λ)|| gets reduced away from σ(A), as for normal operators [cf. (33)], so that the -pseudospectra sets pattern becomes "flattened" (a signature of good analytic behaviour) below the random perturbation scale .
To complement this perspective on the relation between the two given approaches to spectral (in)stability, namely perturbation theory and -pseudospectra, respectively subsections III A and III B, let us connect eigenvalue condition numbers κ(λ i ) with -pseudospectra σ (A). The question we want to address is: how far can the -pseudospectrum σ (A) get away from the spectrum σ(A)? The κ i 's provide the answer. Let us define the 'tubular neighbourhood' ∆ (A) of radius around the spectrum σ(A) which is always contained in the -pseudospectrum σ (A) [4] The key question is about the inclusion in the other direction. Normal operators indeed satisfy [4] where σ 2 (A) indicates the use of a || · || 2 norm. That is, a (||δA|| < ) perturbed eigenvalue of a normal operator can move up to a distance from σ(A). This is what we mean by spectral stability: an operator perturbation of order induces an eigenvalue perturbation also of order . However, in the non-normal case, where κ(λ i ) > 1, it holds (for small ) [4] σ so that σ (A) can extend into a much larger tubular neighbourhood of radius ∼ κ(λ i ) around each eigenvalue, signaling spectral instability if κ(λ i ) 1. This bound is the essential content of the Bauer-Fike theorem relating pseudospectra and eigenvalue perturbations (cf. [4] for a precise formulation).

IV. NUMERICAL APPROACH: CHEBYSHEV'S SPECTRAL METHODS
The present work is meant as a first assessment of BH QNM (in)stability by using pseudospectra. At this exploratory stage, we address the construction of pseudospectra in a numerical approach. As indicated in section III B 3, the study of the spectral stability of non-normal operators is a challenging problem that demands high accuracy. Spectral methods provide well-adapted tools for these calculations [4,52,53].
We discretize the differential operator L in (9)-(14) via Chebyshev differentiation matrices, built on Chebyshev-Lobatto n-point grids, producing L N matrix approximates (we note systematically n = N + 1 in spectral grids, cf. appendix C). Once the operator is discretized, the construction of the pseudo-spectrum requires the evaluation of matrix norms. A standard practical choice [4,52] involves the matrix norm induced from the Euclidean L 2 norm in the vector space C n that, starting from Eq. (32), can be rewritten [4,52] as where σ min (M ) denotes the smallest singular value of M , that is, Although Eq. (40) captures the spectral instability structure of A, the involved L 2 scalar product in C n is neither faithful to the structure of the operator L in Eq. (9), nor to the physics of the BH QNM problem (cf. discussion in section III B 4). Instead, we rather use the natural norm in the problem, specifically the Chebyshev-discretrised version of the 'energy norm' (20), following from the Chebyshev-discretised version of the scalar product (22). Specifically, we write the discretised scalar product in an appropriate basis as (we abuse the notation, since we use ·, · E as in (22), although this is now a scalar product in a finite-dimensional space C n ) where G E ij is the Gram matrix corresponding to (22) (cf. appendix C for its construction) and we note u * =ū t . The adjoint A † of A with respect to ·, · E writes then The vector norm || · || E in C n associated with ·, · E in Eq. (41) induces a matrix norm || · || E in M n (C) (again, we abuse notation by using the same symbol for the norm in C n and in M n (C)). Then, cf. appendix B, the -pseudospectrum σ where s E is the smallest of the "generalized singular values" with M ∈ M n (C) and its adjoint M † given by Eq. (42).

V. A TOY MODEL: PÖSCHL-TELLER POTENTIAL
As presented in the previous sections, in our study of BH QNMs and their (in)stabilities, we exploit the geometrical framework of the hyperboloidal approach to analytically impose the physical boundary conditions at the BH horizon and at the radiation zone (future null infinity). As discussed in section II, a crucial feature of such a strategy is that it allows us to cast the calculation of the QNM spectrum explicitly as the spectral problem of a non-selfadjoint differential operator, which is then the starting point for the tools assessing spectral instabilities presented in section III, namely the construction of the pseudospectrum and the analysis of random perturbations. Finally, spectral methods discussed in section IV are employed to study these spectral issues through a discretisation for the derivative operators. Prior to the study of the BH case, the goal of this section is to illustrate this strategy in a toy model, namely the one given by the Pöschl-Teller potential.

A. Hyperboloidal approach in Pöschl-Teller
The Pöschl-Teller potential [54] V (45) has been widely used as a benchmark for the study of QNMs in the context of BH perturbation theory (e.g. [55][56][57]). Interestingly, QNMs of this potential have been very recently revisited to illustrate, on the one hand, the hyperboloidal approach to QNMs in a discussion much akin to the present one (cf. [42], cast in the setting of de Sitter spacetime) or, on the other hand, functional analysis key issues related to the selfadjointness of the relevant operator [58]. Our interest in Pöschl-Teller stems from the fact that it shares the fundamental behavior regarding QNM (in)stability to be encountered later in the BH context, but in a mathematically much simpler setting. In particular, Pöschl-Teller presents weaker singularities than the Regge-Wheeler and Zerilli potentials in Schwarzschild, that translates in the absence of a continuous part of the spectrum of the relevant operator L (corresponding to the "branch cut" in standard approaches to QNMs). Let us consider the compactified hyperboloids given by Bizoń-Mach coordinates [59,60] or, equivalently In the spirit of the conformal compactification along the hyperboloids described in section II A, we add the two points at (null) infinity (no BH horizon here), namely 1]. Under this transformation the wave equation (4) reads namely the version of Eq. (7) corresponding to the transformation (46). We notice that angular labels ( , m) are not relevant in the one-dimensional Pöschl-Teller problem. If x = 1 we can divide by (1 − x 2 ) and, defining we can write This expression is formally valid for any given potential V (x) (although analyticity issues may appear if the asymptotic decay is not sufficiently fast, as it is indeed the case for Schwarzschild potentials at I + ). If we now insert the Pöschl-Teller expression (45) and notice sech 2 (x) = 1 − x 2 , we get a remarkably simple effective potentialṼ , actually a constant In particular, the Pöschl-Teller wave equation (50) exactly corresponds to Eq. (4) in [42], so that the Pöschl-Teller problem is equivalent to the Klein-Gordon equation in de Sitter spacetime with mass m 2 = V o . In the following, we choose (5), so that we can set Performing now the first-order reduction in time (8)- (9) we get for w(x), p(x), q(x) and γ(x) in Eq. (12) the values and therefore the operators L 1 and L 2 building the operator L in Eq. (10) write, in the Pöschl-Teller case, as As discussed in section II B 1, the function p(x) = 1−x 2 vanishes at the boundaries of the interval [a, b] = [−1, 1], defining a singular Sturm-Liouville operator. This is at the basis of the absence of boundary conditions, if sufficient regularity is enforced on the eigenfunctions of the spectral problem. Regularity therefore encodes the outgoing boundary conditions (see below). Finally, the scalar product (22) writes in this case as B. Pöschl-Teller QNM spectrum

Exact Pöschl-Teller QNM spectrum
Pöschl-Teller QNM spectrum can be obtained by solving the eigenvalue problem (14)-(15) with operators L 1 and L 2 given by Eq. (54). As commented above, no boundary conditions need to be added, if we enforce the appropriate regularity. In this particular case, this eigenvalue problem can be solved exactly. The resolution itself is informative, since it illustrates this regularity issue concerning boundary conditions.
If we substitute the first component of (15) into the second or, simply, if we take the Fourier transform in τ in Eq. (50) (withṼ = 1 from the chosen λ leading to Eq. (52)), we get This equation can be solved in terms of the hypergeometric (see details in appendix D). In particular, for each value of the spectral parameter ω we have a solution that can be written as a linear combination of linearly independent solutions obtained from 2 F 1 (a, b; c; z). Discrete QNMs are obtained only when we enforce the appropriate regularity, that encodes the outgoing boundary conditions. In this case, this is obtained by enforcing the solution to be analytic in x ∈ [−1, 1] (corresponding in z to analyticity in the full closed interval [0, 1]), which amounts to truncate the hypergeometric series to a polynomial. We emphasize that such a need of truncating the infinite series to a polynomial, a familiar requirement encountered in many different physical settings, embodies here the enforcement of outgoing boundary conditions. In sum, this strategy leads to the Pöschl-Teller QNM frequencies (cf. e.g. [56,61]) with QNM eigenfunctions where P (α,β) n are the Jacobi polynomials (see appendix D). Two comments are in order here: i) QNMs are normalizable: QNM eigenfunctions φ ± n (x) are finite and regular when makingx → ±∞, corresponding to x = ±1. This is in contrast with the exponential divergence of QNM eigenfunctions in Cauchy approaches, where the time slices reach spatial infinity i o . This is a direct consequence of the hyperboloidal approach with slices reaching I + . The resulting normalizability of the QNM eigenfunctions can be relevant in e.g. resonant expansions (cf. e.g. discussion in [9]).
ii) QNM regularity and outgoing conditions: In the present case, Pöschl-Teller in Bizoń-Mach coordinates, analyticity (actually polynomial structure) implements the regularity enforcing outgoing boundary conditions. Analyticity is too strong in the general case. But asking smoothness is not enough (see e.g. [35]). In Refs. [39][40][41] this problem is approached in terms of Gevrey classes, that interpolate between analytic and (smooth) C ∞ functions, identifying the space of (σ, 2)-Gevrey functions as the proper regularity notion. The elucidation of the general adequate functional space for QNMs, tantamount of the consistent implementation of outgoing boundary conditions, is crucial for the characterization of QNMs in the hyperboloidal approach.
2. Numerical Pöschl-Teller QNM spectrum Fig. 1 shows the result of the numerical counterpart of the Pöschl-Teller eigenvalue calculation, whose exact discussion has been presented above, by using the discretised operators L, L 1 and L 2 described in section IV and appendix C This indeed recovers numerically the analytical result in Eq. (57) (we drop the "±" label, focusing on one of the branches symmetric with respect the vertical axis).  (57) is far from being a trivial result, as already illustrated in existing systematic numerical studies. This is in particular the case of Ref. [62] (where Pöschl-Teller is referred to as the Eckart barrier potential), where the fundamental mode ω ± 0 in (57) is stable and accurately recovered, whereas all overtones ω ± n≥1 suffer from a strong instability (triggered, according to the discussion in [62,93], by the C 1 -regularity of the approximation modelling the Pöschl-Teller potential) and could not be recovered.
In our setting, a convergence study of the numerical values shows that the relative error between the exact QNM ω n and the corresponding numerical approximation ω (N ) n (obtained at a given truncation N of the differential operator) actually increases with the resolution. This is a first hint of the instabilities to be discussed later. Indeed, the top panel of Fig. 2 displays the error for the fundamental mode n = 0 and the first overtones n = 1, . . . , 4 when the eigenvalue problem for the discretised operator is naively solved with the standard machine roundoff error for floating point operations (typically, ∼ 10 −16 for double precision). It is astonishing how, despite the simplicity of the exact solution, the relative error grows significantly already for the first overtones and, crucially, more strongly as the damping grows with higher overtones. To mitigate such a drawback, one needs to modify the numerical treatment in order to allow for a smaller roundoff error in floating point operations. The bottom panel Fig. 2 shows the error E (N ) n when the calculations are performed with an internal roundoff error according to 5×Machine Precision, i.e. ∼ 10 −5×16 . In this case, the fundamental QNM n = 0 is "exactly" calculated at the numerical level (i.e. the difference between its exact value and the numerical approximation vanishes in the employed precission). The error for the overtones still grows, but in a safe range for all practical purposes. The values displayed in the bottom panel of Fig. 1 were obtained with internal roundoff error set to 10×Machine Precision and we can assure that the errors of all overtones are smaller than 10 −100 .

Condition numbers of QNM frequencies
The growth in the relative error as we move to higher overtones in Fig. 2, suggests an increasing spectral instability in n of eigenvalues ω ± n , triggered by numerical errors related to machine precision, so that this instability can be reduced (but not eliminated) by improving the internal roundoff error.
At the level of the non-perturbed spectral problem (59), and in order to assess more systematically such spectral instability, we can apply the discussion in section III A to the Pöschl-Teller approximates L N . Namely, solving the righteigenvector problem (59) together with left-eigenvector one we can compute the condition numbers κ n || E introduced in Eq. (30). Notice that this is quite a non-trivial calculation, since it involves: first, the construction of the adjoint operator L N † = G E −1 · (L N ) * · G E and, second, the calculation of scalar products ·, · E and (vector) energy norms ||·|| E . These calculations involve the determination of the Gram matrix G E associated with the energy scalar product (55) by implementing expression (C23) in appendix C. These expressions are quite non-trivial and in the following section we provide a strong test to the associated analytical and numerical construction.
The result is shown in the top panel of Fig. 1. The ratio of the condition numbers κ n , relative to the condition number of the fundamental mode κ 0 , grows strongly with n. This indicates a strong and increasing spectral instability consistent with the error convergence displayed in Fig. 2. The rest of this section is devoted to address this spectral stability issue. Precision. Note that missing points for n = 0 correspond to errors that exactly vanish at the employed machine precision.

Motivating the pseudospectrum
As the previous discussion makes apparent, a crucial question that arises after obtaining the QNM spectrum of the operator L in Eq. (10), with L 1 and L 2 in (54) is whether such QNM eigenvalues are stable under small perturbations of L. More specifically for QNM physics, and in the context of the wave equation (4), whether the QNM spectrum is stable under small perturbations of the potential V . The latter is the specific type of perturbation we are assessing in this work.
In the numerical approach we have adopted, perturbations in the spectrum under small pertubations in L may arise either from numerical noise resulting from the chosen discretisation strategy, or they can originate from "real-world sources", namely small physical perturbations of the considered potential V . Ultimately, in the BH setting for which Pöschl-Teller provides a toy model, such physical perturbations could stem  Fig. 1 are superimposed for reference. The color logscale corresponds to log 10 , white lines indicating the boundaries of -pseudospectra sets σ , whose interior extends upwards in ω-C.
from a "dirty" environment surrounding a black hole, and/or emergent fluctuations from quantum-gravity effects. Therefore, the question of whether QNM spectrum instability is a structural feature of the operator L -i.e. not just an artifact of a given numerical algorithm -is paramount for our understanding of the fundamental physics underlying the problem. A pragmatic approach [63] to address this question consists in explicitly introducing families of perturbations, and study their effect on the QNM spectra themselves [10][11][12][13][14]21]. We will make contact with this approach later in section V D, but before that, we apply the pseudospectrum approach described in section III B to the Pöschl-Teller problem. Indeed, one of the main goals of our present work is to bring attention to and emphasise the fact that the unperturbed operator already contains crucial information to assess such (in)stability features. We have already encountered this fact in the evaluation of the condition numbers κ n in Fig. 1, that only depends on the unperturbed operator L, but we develop this theme further with the help of the pseudospectrum notion. Indeed, pseudospectrum analysis provides a framework to identify the (potential) spectral instability, which is oblivious to the particular perturbation employed. Then, in a second stage, actual perturbations of the operators with a particular emphasis on random perturbations along the lines in section III C, can be used to complement and refine such pseudospectrum analysis. Fig. 3 shows the pseudospectrum for the Pöschl-Teller potential in the energy norm of Eq. (20) associated with the scalar product (55). Let us explain the content of such a figure. According to the characterization in Eq. (31) of thepseudospectrum of the operator L, the set σ (L) is the collection of all complex numbers λ ∈ C that are actual eigenvalues for some operator L + δL, where δL is a small perturbation of "size" smaller than a given > 0. Consequently and crucially, introducing a perturbation δL with ||δL|| E < entails an actual ("physical") change in the eigenvalues λ n that can reach up the boundary of the σ (L) set, marked with white lines in Fig. 3. So the key question is to assess if -pseudospectra sets for small can extend in large areas of C or not. This is intimately related to condition numbers κ n controlling the spectral (in)stability of eigenvalues, and explicitly estimated by the Bauer-Fike relation (39) between -pseudospectra sets and 'tubular neighbourhoods' ∆ κ of radii κ n around the spectrum. Let us first discuss a selfadjoint test case and, in a second stage, the actual non-selfadjoint case [64].

Pseudospectrum: selfadjoint case
As discussed in section II C, setting L 2 = 0 in Eq. (10) -while keeping L 1 as in Eq. (54)-leads to a selfadjoint operator L [65]. Therefore the associated spectral problem is, cf. section III A, stable. A typical pseudospectrum in the seldajoint (more generally 'normal') case is illustrated by Fig. 4: a "flat" pseudospectrum with large values of forpseudospectra sets, when moving "slightly" (in the C-plane) from the eigenvalues. Note also, in this case, the horizontal contour lines far from the spectrum, indicating that all eigenvalues share the same stability properties in the energy norm. Let us describe Fig. 4 in more detail. Boundaries ofpseudospectra σ (A) are marked in white lines, with 's corresponding to the values in the color log-scale. Pseudospectra σ (L) are, by construction, "nested sets" around the spectrum (red points in Fig. 4), the latter corresponding to the "innermost set" σ (L) when → 0. In this selfadjoint case, condition numbers in (30) must satisfy κ n = 1, as we have verified and explicitly shown in the top pannel of Fig. 4. Then, and consistently with Eq. (38), the corresponding nested sets σ (L) are actually tubular regions ∆ (L) of "radius " around the spectrum, so that a change δL with a norm of order in the operator L entails a maximum change in the eigenvalues of the same order . Specifically, -pseudospectra sets show concentric circles around the spectra that quickly reach large-epsilon values, i.e. ∼ O(1), when moving away from eigenvalues. As a consequence, one would need perturbations in the op-erator of the same order to dislodge the eigenvalues slightly away from their original values: we say then that L is spectrally stable. Pseudospectra sets with small are then "tightly packed" in "thin throats" around the spectrum, so that green colors are indeed so close to spectrum "red points" that they are not visible in the scale of Fig. 4, giving rise to a typical "flat" pseudospectrum figure of a "single color".
Horizontal boundaries of -pseudospectra, when far from the spectrum, is a feature (in this particular problem) of the use of the energy norm. If another norm is used, e.g. the standard one induced from the L 2 norm in C n , the global "flatness" of the pseudospectrum is still recovered, especially when comparing with the corresponding scales in Fig. 3, indicating already a much more stable situation than the general L 2 = 0 case. But when refining the scale, one would observe that pseudospectra contour lines far from the spectrum are not horizontal but present a slope growing with the frequency. This indicates that, under perturbations of the same size in that L 2 norm, higher frequencies can move further that low frequencies, this being in tension with the equal stability of all the eigenvalues. What is going on is the effect commented in section III B 4 concerning the impact of the norm choice on the notions of "big/small": when using the L 2 norm, we would be marking with the same "small" different perturbations among which there exist δL instances that actually excite strongly the high frequencies, but such a feature is blind to the L 2 norm. If using however a norm sensitive to highfrequency effects, as it is the case of the energy norm that has a H 1 character incorporating derivative terms, those same perturbations δL would have a norm much higher than , the derivative terms in the energy norm indeed weighting more as the frequency grows. What in the L 2 norm was a small perturbation δL, turns out to be a big one in the energy one, so stronger modifications in the eigenvalues are indeed consistent with stability. In practice, in order to construct a given -pseudospectrum set, such "high-energy" perturbations δL need to be renormalized to keep fixed, something that the energy norm does automatically. This is a neat example of how the choice of the norm affects the assessment of spectral stability and, in particular, of the importance of the energy norm in the present work, namely for high-frequency issues. Fig. 4 may appear as a boring figure, but it is actually a tight and constraining test of our construction, both at the analytical and numerical level. First of all, the two panels in Fig. 4 correspond to different calculations: the top panel results from an eigenvalue calculation (actually two, one for L and another for L † ), whereas the "map" in the bottom panel is the result of calculating the energy norm of the resolvent R L (λ) = (λ Id − L) −1 at each point λ ∈ C. Both calculations depend on the construction of the Gram matrix G E , but are indeed different implementations. The κ n = 1 in the top pannel is a most stringent test, since modifications in either the analytical structure of the scalar product (55) or the slightest mistake in the discrete counterpart (C23) spoil the result. As discussed at the end of section V B 3, this provides a strong test both at analytical and the numerical discretization of the differential operator and scalar product. On the other hand, the plain flatness of the pseudospectrum in the bottom pannel is a strong test of the selfadjoint character of L when L 2 = 0 that, given the subtleties of the spectral discretization explained in appendix C, provides a reassuring non-trivial test to the whole numerical discretization scheme.

Non-selfadjoint case: Pöschl-Teller pseudospectrum
In contrast with the selfadjoint case, when considering the actual L 2 = 0 of the Pöschl-Teller case, pseudospectra sets σ (L) with small extend in Fig. 3 into large regions of C (with typical sizes much larger than ) and therefore the operator L is spectrally unstable: very small (physical) perturbations δL, with ||δL|| E < , can produce large variations in the eigenvalues up to the boundary of the now largely extended region σ (L). Such strong variations of the spectrum are not a numerical artifact, related e.g. to machine precission, but they rather correspond to an actual structural property of the non-perturbed operator. Indeed, large values of the condition numbers κ n in the top panel of Fig. 1 entail that the tubular sets ∆ κ (L) in Eq. (36) extend now into large areas in C. This fact on κ n 's is consistent with the large regions in Fig. 3 corresponding to σ (L) sets with very small 's. Such an non-trivial pattern of -pseudospectra is a strong indication of spectral instability, although without a neat identification of the actual nature of perturbations triggering instabilities.

Reading pseudospectra: "topographic maps" of the resolvent
In practice, if one wants to read from pseudospectra -such as those in Fig. 3 or Fig. 4-the possible effect on QNMs of a physical perturbation of (energy) norm of order , one must first determine the "white-line" corresponding to that (using the log-scale). Then, eigenvalues can move potentially in all the region bounded by that line (namely, thepseudospectrum set for the non-perturbed operator L) that, in Fig. 3, corresponds to the region "above" white lines.
Pseudospectra can actually be seen as a "map" of the analytical structure of the resolvent R L (λ) = (λ Id − L) −1 of the operator L, taken as a function of λ. This corresponds to the characterization in Eq. (32), which it is indeed the one used to effectively construct the pseudospectrum (specifically, its realisation (43) in the energy norm; cf. section B 3 for details). In this view, the boundaries of the -pseudospectra (white lines in Figs. 3 and 4) can be seen as "contour lines" of the "height function" ||R L (λ)|| E , namely the norm of the resolvent R L (λ). In quite a literal sense, the pseudospectrum can be read then as topographic map, with stability characterised by very steep throats around eigenvalues fastly reaching flat zones away from the spectrum, whereas instability corresponds to non-trivial "topographic patterns" extending in large regions of the map away from eigenvalues.
In sum, this "topographic perspective" makes apparent the stark contrast between the flat pattern of the selfadjoint case of Fig. 4, corresponding to stability, and the non-trivial pattern of the (non-selfadjoint) Pöschl-Teller pseudospectrum in Fig. 3, in particular indicating a (strong) QNM sensitivity to perturbations that increases as damping grows.

D. Pöschl-Teller perturbed QNM spectra
Pseudospectra inform about the spectral stability and instability of an operator, but do not identify the specific type of perturbation triggering instabilities. Therefore, in a second stage, it is illuminating to complement the pseudospectrum information with the exploration of spectrum instability with "perturbative probes" into the operator, always under the perspective acquired with the pseudospectrum. A link between both pseudospectra and perturbation strategies is provided by the Bauer-Fike theorem [4], as expressed in Eq. (39).

Physical instabilities: perturbations in the potential V
Not all possible perturbations of the L operator are physically meaningful. An instance of this, in the setting of our numerical approach, are machine precision error perturbations δL N to the L N matrix. As discussed in section V B 2, machine precision errors indeed trigger large deviations in the spectrum, consistently with the non-trivial pattern of the pseudospectrum in Fig. 3, but clearly we should not consider such effects as physical. They are a genuine numerical artifact, since the structure of the perturbation δL N does not correspond to any physical or geometrical element in the problem.
The methodology we follow to address this issue is: i) given a grid resolution N , we first set the machine precision to a value sufficiently high so as to guarantee that all nonperturbed eigenvalues are correctly recovered, and ii) we then add a prescribed perturbation with the specific structure corresponding to the physical aspect we aim at studying.
In the present work we focus on a particular kind of perturbation, namely perturbations to the potential V and, more specifically, perturbations δṼ to the rescaled potentialṼ in (21). This is in the spirit of studying the problem in [1]. That is, we consider perturbations δL to the L operator of the form We note that, at the matrix level, the δṼ submatrix is just a diagonal matrix. Therefore, the structure of δL in Eq. (62) is a very particular one. The pseudospectrum in Fig. 3 tells us that L is spectrally unstable, and we know that machine precision perturbations trigger such instabilities, but nothing guarantees that L is actually unstable under a perturbation of the particular form (62). It is a remarkable fact, crucial for our physical discussion, that L is indeed unstable under such perturbations and, therefore, under perturbations of the potential V .
2. Random and high-frequency perturbations in the potential V We have considered two types of generic, but representative, perturbations δL of the form given in Eq. (62): i) Random perturbations δṼ r : we set the perturbation according to a normal Gaussian distribution on the collocation points of the grid. This is, by construction, a high-frequency perturbation. Random perturbations are a standard tool [4] to explore generic properties of spectral instability and there exists indeed a rich interplay between pseudospectra and random perturbations [6].
ii) Deterministic perturbations δṼ d : we have chosen in order to address the specific impact of high and low frequency perturbations in QNM spectral stability, by exploring the effect of changing the wave number k.
Perturbations δṼ are then rescaled so as to guarantee ||δL|| E = . The impact on QNM frequencies resulting from adding these perturbations is shown in Fig. 5. In both random and deterministic cases, the sequence of images in Fig. 5 shows a high-frequency instability of QNM overtones, that "migrate" to new QNM branches. The fundamental (slowest decaying) QNM is however stable under these perturbations.
More generally, such QNM instability is sensitivity with respect to both perturbation's "size"and frequency.
Before we further discuss the details of the QNM instability, namely the nature of the new QNM branches, an important point must be addressed: whether the values obtained correspond to the actual eigenvalues of the new, perturbed operator L + δL, or whether they are an artifact of some numerical noise. As in the non-perturbed case discussed in section V B 2, and as explained above when introducing the employed methodology, results are obtained with a high internal accuracy (10×Machine Precision), so that any numerical noise is below the range of showed values. Proceeding systematically, Fig. 6 presents the convergence tests for a few eigenvalues resulting from the deterministic perturbation (random perturbations do not admit this kind of test) with norm ||δṼ d || E = 10 −8 and frequency k = 20 (bottom right panel of Fig. 5). The relative error is calculated as i.e., in the absence of exact results, we take as reference the values with a high resolution N = 400. As representative QNMs, we have chosen: a) The last "unperturbed" overtone, whose value is actually very close to the (truly) unperturbed QNM ω 4 .
b) The first new QNM on the imaginary axis.
c) Three QNMs along the new branch with values spread in 1 Re(ω n ) 10 and 5 Im(ω n ) 8. Left column: Sequence of QNM spectra for the Pöschl-Teller potential subject to a random perturbation δṼr of increasing "size" (in energy norm). The sequence shows how "switching on" a perturbation makes the QNMs migrate to a new branch (that actually corresponds to a pseudospectrum contour line, compare with Fig. 3), in such a way that the instability starts appearing at highly-damped QNMs and descends in the spectrum as the perturbation grows (unperturbed values, in red, are kept along the sequence for comparison). The top panel corresponds to the non-perturbed potential shown in Fig. 1, the second panel shows how a random perturbation of with (energy) norm ||δṼr|| E = 10 −16 already reaches the 6th QNM overtone, whereas in the third panel a perturbation with ||δṼr|| E = 10 −8 already reaches the 3rd overtone. This confirms the instability already detected in the pseudospectrum, indicating its high-frequency nature. Crucially, to reach the fundamental mode, a perturbation of the same order O(1) as the variation of the eigenvalue is required, this demonstrating the stability of the fundamental QNM in agreement with the pseudospectrum in Fig. 3. Right panel: Sequence of QNM spectra for Pöschl-Teller subject to a deterministic perturbation δṼ d ∼ cos(2π kx). The first panel shows again the unperturbed potential, whereas the second one shows that a "low frequency" (k = 1) perturbation leaves the spectrum unperturbed, in spite of the ||δṼ d || E = 10 −8 norm (compare with the random case with the same norm): this illustrates the harmless character of "low frequency" perturbations. The third panel shows how keeping the norm of the perturbation but increasing its frequency indeed "switches on" the instability, confirming the "high frequency" insight gained from random perturbations. The fourth panel shows how the instability increases with the frequency but less efficiently than with random perturbations of the same norm. One observes a systematic convergence, with the relative error dropping circa 10 orders of magnitudes when the numerical resolution increases [66] from N = 150 to N = 400. This result confirms that the spectrum corresponds indeed to the new, perturbed operator, and is not a numerical artifact. This neatly shows the unstable nature of the QNM spectrum of the unperturbed Pöschl-Teller operator: eigenvalues indeed migrate to new branches under very small perturbations.

Perturbed QNM branches and pseudospectrum
High-frequency perturbations trigger the migration of QNM overtone frequencies to new perturbed QNM branches. Fig. 7 displays the perturbed QNM spectra on the top of the pseudospectra for the unperturbed operator. The remarkable "predictive power" of the pseudospectrum becomes apparent: perturbed QNMs follow the boundaries of pseudospectrum sets. That is, QNM overtones "migrate" to new branches corresponding to the -pseudospectra contour lines. This happens for both random and deterministic high-frequency perturbations. Crucially, no such instability is observed for lowfrequency deterministic perturbations, with small wave number k. Consequently, we shall refer in the following to this effect as an ultraviolet instability of QNM overtones.
Remarkably, such high-frequency QNM instability is not limited to highly damped QNMs but indeed reaches the lowest overtones, the random perturbations being more effective in reaching the slowest decaying overtones for a given norm ||δV || E = . This result is qualitatively consistent with analyses in [2,13] for Dirac-delta potentials (compare e.g., perturbed QNM branches in Fig. 7 here with Fig. 1 in Ref. [13]). These findings advocate the use of pseudospectra to probe QNM instability, demonstrating its capability to capture it al- ready at the level of the non-perturbed operator. At the same time, pseudospectra are oblivious to the nature of the perturbation triggering instabilities. A complementary perturbation analysis, in particular through random perturbations, has been then necessary to identify the high-frequency nature of the instability, confirming its physicality in the sense of being associated with actual perturbations of the potential V .

High-frequency stability of the slowest decaying QNM
The high-frequency instability observed for QNM overtones is absent in the fundamental QNM. The slowest decaying QNM is therefore ultraviolet stable. Such stability is already apparent in the pseudospectrum in Fig. 3, where the order of the 's corresponding to -pseudospectra sets around the fundamental QNM reaches the values in the stable self-adjoint case in Fig. 4. This high-frequency stability is then confirmed in the perturbation analysis. Indeed, Fig. 7 demonstrates the need of large perturbations in the operator in order to reach the fundamental QNM, namely (random) perturbations with a 'size' ||δṼ || E of the same order as the induced variation in ω ± 0 . This behaviour is a tantamount of spectral stability. The contrast between the high stability of ω ± 0 and the instability of overtone resonances ω ± n≥1 has already been evoked in V B 2, when referring to the large condition number ratios κ n /κ 0 , in particular referring to Bindel & Zworski's discussion in [62,93]. This high-frequency stability of the fundamental mode is in tension with the instability found by Nollert in [1] for the slowest decaying mode for Schwarzschild. We will revisit this point in section VI D 3. For the time being, we simply emphasize that the observed stability relies critically on the faithful treatment of the asymptotic structure of the po- FIG. 8. QNMs of Pöschl-Teller "cut" potential. Setting Pöschl-Teller potential to zero outside a certain interval [xmin, xmax] introduces high-frequency perturbations that make QNM overtones migrate to pseudospetrum contour lines, as well as an "infrared" modification that alters the fundamental QNM frequency. Whereas the latter tends to the non-perturbed Pöschl-Teller value as xmin → −∞ and xmax → ∞, QNM overtones remain always strongly perturbed.
tential, that is in-built in the adopted hyperboloidal approach permitting to the capture the long-range structure of the potential up to null infinity I + . It is only when we enforce a modification of the potential at "large distances" that the "low frequency" fundamental QNM is affected. This is illustrated in Fig. 8 (see also [67]), corresponding to a Pöschl-Teller potential set to zero beyond a compact interval [x min , x max ]: such "cut" introduces high frequencies that make migrate the overtones to the new branches and, crucially, alters the asymptotic structure so that the fundamental QNM is also modified. Such "infrared" effect is however compatible with the spectral stability of the fundamental QNM, since such "cut" of the potential does not correspond to a small perturbation in δL.

Regularization effect of random perturbations
Before proceeding to discuss the BH case, let us briefly comment on an apparently paradoxical phenomenon resulting from the interplay between random perturbations and the pseudospectrum. In contrast with what one might expect, the addition of a random perturbation to a spectrally unstable operator L does not worsen the regularity properties of L but, on the contrary, it improves the analytical behaviour of its resolvent R L (λ) [6,[68][69][70]. This is illustrated in Fig. 9, that shows a series of pseudospectra corresponding to random perturbations of the Pöschl-Teller potential with increasing ||δṼ r || E . In addition to the commented migration of QNM overtones to pseudospectra contour lines, we observe two phenomena: i) -pseudospectra sets with > ||δṼ r || E are not affected by the perturbation, whereas ii) the pseudospectrum structure for Im(ω n ) FIG. 9. Pseudospectra of Pöschl-Teller under random perturbations δṼr of increasing norm, demonstrating the "regularizing" effect of random perturbations: pseudospectra sets σ bounded by that "contour line" reached by perturbed QNMs become "flat", a signature of improved analytic behaviour of the resolvent, as illustrated in Fig. 4. Pseudospectra sets not attained by the perturbation remain unchanged. Regularization of R L+δL (λ) increases as ||δṼr|| E grows.
< ||δṼ r || E is smoothed into a "flat pattern". As we have discussed in the setting of Fig. 4, such flat pseudospectra patterns are the signature of spectral stability, a tantamount of regularity of the resolvent R L (λ). The resulting improvement in the spectral stability of L + δL, as compared to L, is indeed consistent with the convergence properties of the respective QNM spectra, as illutrated by the contrast between the corresponding convergence tests, namely Fig. 6 and Fig.2. In brief, random perturbations improve regularity, an intriguing effect with suggestive physical implications in our QNM setting.

VI. SCHWARZSCHILD QNM (IN)STABILITY
We address now the physical BH case, namely the stability of QNMs in Schwarzschild spacetime. Whereas the previous section has been devoted, to a large extent, to discuss some of the technical issues in QNM stability, the spirit in this section is to focus more on the physical implications, in particular in the perspective of assessing the pionnering work in [1,2].

A. Hyperboloidal approach in Schwarzschild
The attempt to implement the QNM stability analysis in the coordinate system employed for Pöschl-Teller, namely the Bizoń-Mach chart (46), is unsuccessful. The reason is the bad analytic behaviour at null infinity of Schwarzschild potential(s) in the corresponding coordinate x. Instead of this, we resort to the 'minimal gauge' slicing [35,36,71], devised to improve regularity in the Schwarzschild(-like) case.
To construct horizon-penetrating coordinates reaching null infinity, one defines a height function h in (6) by first considering an advanced time coordinate built on the rescaled tortoise coordinatex = r * /λ, with r * = r + 2M ln(r/2M − 1), so that the BH horizon is atx → −∞, and then enforcing a deformation of the Cauchy slicing into a hyperboloidal one through the choice of a 'minimal gauge', prescribed under the guideline of preserving a good analytic behavior at I + . In a second stage, the function g in (6) implementing the compactification along hyperboloidal slices is implicitly determined by (note that instead of x in (6), we rather use σ for the spatial coordinate, so as to keep the standard use in [35,36,71]) Choosing λ = 4M in the rescalingx = r * /λ of Eq. (5), the steps above result (see details in [35,36,71]) in the 'minimal gauge' hyperboloidal coordinates for the transformation (6) that, upon addition of the BH horizon and I + points, maps x ∈ [−∞, ∞] to the compact interval σ ∈ [a, b] = [0, 1], with the BH horizon at σ = 1 and future null infinity at σ = 0. Implementing transformation (69) in the first-order reduction in time, cf. Eqs. (8)-(9), we get now for w(σ), p(σ), q (σ) (now explicitly depending on ) and γ(σ) in Eq. (12) leading to the L 1 and L 2 operators building L in Eq. (10) where the rescaled potentialṼ (σ) := q (σ) results, in the respective axial and polar cases, in the explicit expressions Finally, from Eqs. (70) and (22), the energy scalar product is where the weightṼ is fixed by Eq. (71) for each polarization.

B. Schwarzschild QNM spectrum
As discussed in section II B 1, outgoing boundary conditions have been translated into regularity conditions on eigenfunctions. Specifically, as we have seen in the Pöschl-Teller case, the operator L 1 in (70)  But there is a key difference between the Pöschl-Teller and the BH case: whereas in Pöschl-Teller the function p(x) = (1 − x)(1 + x) vanishes linearly at the boundaries, and therefore x = ±1 are regular singular points, in Schwarzschild this is true for σ = 1 (BH horizon) but not for σ = 0 (I + ), due to the quadratic σ 2 term. Null infinity is then an irregular singular point. This is the counterpart, in our compactified hyperboloidal formulation, of the power-law decay of Schwarzschild potentials responsible for the branch cut in the Green function of Eq. (4), with its associated "tails" in late decays of scattered fields. In the context of our spectral problem for the operator L, this translates in the appearance of a ("branch cut") continuous part in the spectrum. This has an important impact on the numerical approach, since the continuous branch cut is realized in terms of actual eigenvalues of the discretised approximates L N . Such eigenvalues are not QNMs and can indeed be unambiguously identified, but their presence has to be taken into account when performing the spectral stability analysis, that becomes a more delicate problem than in Pöschl-Teller. In this context, the latter becomes a crucial benchmark to guide the analysis in the BH case.
The Schwarzschild (gravitational) QNM spectrum (for = 2) is shown in Fig. 10, that presents the result of the numerical calculation of the spectrum of the L operator defined by (70). This is obtained either for the Regge-Wheeler or the Zerilli rescaled potentials in (71), corresponding respectively to potentials (65) and (66). This provides a crucial internal consistency check for the analytical and numerical construction, since both potentials are known to be QNM-isospectral (see below in section VI D 2). The branch cut structure is apparent in the eigenvalues along the upper imaginary axis. Such "branch cut" points can be easily distinguished from the algebraicly special QNM corresponding to ω n=8 , also in the imaginary axis, simply by changing the resolution: branch points move "randomly" along the vertical axis, whereas ω n=8 stays at the same frequency (see later VI D 1 for a more systematic approach to establish the "non-branch" nature of ω n=8 , when we will consider high-frequency perturbations to QNMs).
Due to the lack of an exact expression for the Schwarzschild QNMs, one must compare the obtained values against those available in the literature via alternative approaches -see, for instance [28][29][30][31][74][75][76]. An estimative for the errors when QNMs are calculated with the methods from this work is found in Ref. [71]. From the practical perspective, and regardless of the numerical methods, it is well known that the difficulty to accurately calculate numerically a given QNM overtone ω ± n increases significantly with n. For instance, convergence and machine precision issues similar to the ones commented above are reported in Refs. [77][78][79], a control of the internal roundoff accuracy being required. Alternatively, iterative algorithms such as the Leaver's continued fraction method [80] require an initial seed relatively near a given QNM, which must be carefully adapted when dealing with the overtones [81]. The bottomline is that the calculation of BH QNM overtones is a challenging and very delicate issue.
In our understanding, the latter challenge is not a numerical hindrance but the consequence of a structural feature of the underlying analytical problem, namely the spectral instability of the Schwarzschild QNM problem. This is manifested already at the present stage of analysis, namely the calculation of QNM frequencies of non-perturbed Schwarzschild, in the eigenvalue condition numbers κ n 's shown in the top panel of Fig. 10: we encounter again the pattern found in the Pöschl-Teller case, cf. Fig. 1, with a growth of the spectral instability as the damping increases, with the notably anomaly of an enhanced stability for the algebraicly special QNM frequency, with n = 8. We devote the rest of the section to explore this spectral instability with the tools employed for Pöschl-Teller.

C. Schwarzschild pseudospectrum
The pseudospectrum of Schwarzschild is presented in Fig.11. As illustrated in Pöschl-Teller, the pseudospectrum provides a systematic and global tool to address QNM spectral instability, already at the level of the unperturbed potential. A "topographic map" of the analytic structure of the resolvent, where regions associated with small -pseudospectra (light green) correspond to strong spectral instability, whereas regions with large (namely O( ) ∼ 1, dark blue) indicate spectral stability. The superposition of the QNM spectrum shows the respective spectral stability of QNM frequencies.
We can draw the following conclusions from Fig. 11: i) The Schwarzschild pseudospectrum indicates a strong instability of QNM overtones, an instability that grows fast with the damping. White-line boundaries corresponding to -pseudospectra with very small 's extend in large regions of the complex plane. This is compatible with the results in [1], providing a rationalealready at the level of the unperturbed potential-for the QNM overtone instability discovered by Nollert.
ii) The slowest decaying QNM is spectrally stable. Fig. 11 tells us that changing the fundamental QNM frequency  Fig. 3), though presenting an enhanced spectral instability indicated by the smaller values of -pseudospectra contour lines (cf. range in color log-scale for log 10 in Fig. 3).
requires perturbations in the operator of order ||δL|| E ∼ 1. This corresponds to spectral stability and it is tension with the results in [1], where the fundamental QNM is found to be unstable. We will address this point below.
iii) Schwarzschild and Pöschl-Teller potentials show qualitatively the same pseudospectrum pattern, with large "green regions" producing patterns in stark contrast with the flat selfadjoint case. On the one hand, this reinforces the usage of Pöschl-Teller as a convenient guideline for understanding the stability structure of BH QNMs and, on the other hand, it points towards an instability mechanism independent, at least in a certain measure, on some of the details of the potential.
We can conclude that Fig. 11 demonstrates -at the level of the unperturbed operator-the stability structure of the BH QNM spectrum, namely the QNM overtone instability and the fundamental-QNM stability. The pseudospectrum does not inform us about the particular type of the perturbations that trigger the instabilities. This is addressed in the following subsection.

D. Perturbations of Schwarzschild potential
Once the Schwarzschild pseudospectrum, together with the condition numbers κ n , have presented evidence of QNM spectral instability at the level of the unperturbed operator, in this section we address the question about the actual physical character of perturbations triggering such instabilities.

Ultraviolet instability of BH QNM overtones
The qualitative agreement between Pöschl-Teller and Schwarzschild pseudospectra, cf. Figs. 3 and 11, together with the experience gained in the study of Pöschl-Teller perturbations regarding the high-frequency instability of all QNM overtones and the stability of the fundamental QNM, guide our steps in the analysis of the BH setting.
a. Random perturbations: spoils from the "branch cut". The presence of a "branch cut" in the Schwarzschild spectrum, discussed in section VI B, translates into a methodological subtlety when considering random perturbations in the BH case, as compared with the Pöschl-Teller one. The difficulty stems from the fact that not only the QNM eigenvalues, but also the eigenvalues associated with the discretized version of the branch cut, are sensitive to random perturbations δṼ r of the potential. As a consequence, the possible contamination from eigenvalues from the branch cut complicates the analysis of the impact of random perturbations on QNM frequencies. This is an artifact of our particular numerical approach, and not a problem of the differential operator itself, but it limits our capability to assess the triggering by random perturbation of the QNM migration to new branches, that was observed in the Pöschl-Teller case (cf. left column of Fig. 5). Other tools, either numerical refinements and/or analytical methodologies, are required to address this specific issue in Schwarzschild.
This does not mean that random perturbations have no use in our BH discussion. An illustrative example is the study of the stability of the algebraicly special Schwarzschild QNM ω n=8 . Whereas random perturbations move "branch cut" eigenvalues away from the imaginary axis, the algebraicly special QNM stays stable. This methodology provides a powerful and efficient tool to probe the "physicality" of specific eigenvalues in very general settings (cf. e.g. Fig. 4 in [82]).
b. Deterministic perturbations. Given the limitations for random δṼ s 's, in the present study we have focused on the class of deterministic perturbations to the potential δṼ d provided by Eq. (63). Crucially, such perturbations do not perturb the "branch eigenvalues" as (much as) random δṼ r do, by-passing then the associated spectral instability contamination. Despite their simplicity, they provide a good toymodel to explore the effects of astrophysically motivated perturbations (assessment of "long range/low frequency" versus "small scale/high frequency" perturbations), as well as those arising from generic approaches to quantum gravity ("small scale/high frequency" effective fluctuations). They are, therefore, conveniently suited to address these instability issues.
The left column in Fig. 12 depicts (with ||δṼ d || E ∼ 10 −8 ) the stability of the first overtones against low frequency perturbations (k = 1, top-left panel) in contrast with the instability resulting from high-frequency perturbations (k = 20, bottom-left panel). Pushing along this line, the right column in Fig. 12 zooms in to study the very first overtones, which are paramount for the incipient field of black-hole spectroscopy. Assessing the (in)stability of the very first overtones is therefore crucial for current research programs in gravitational astronomy. It becomes apparent that the first overtones, this including the very first overtone, are indeed affected without any extraordinary or fine tuned perturbations δṼ d . In particular, and taking the left column as a reference, the first overtone is reached: i) either by considering a "slightly" more intense perturbation (||δṼ d || E ∼ 10 −4 , k = 20), or ii) perturbations with sufficiently high frequency (||δṼ d || E ∼ 10 −8 , k = 60). From this perturbation analysis of the BH potential we conclude: i) all QNM overtones are ultraviolet unstable, as in Pöschl-Teller, the instability reaching the first overtone for sufficiently high frequency; ii) QNMs are stable under low frequency perturbations, this illustrating that spectral instability does not mean instability under "any" perturbation, in particular long-wave perturbations not affecting the QNM spectrum; iii) the slowest decaying QNM is ultraviolet stable, a result in tension with the instability of the fundamental QNM found in [1]. We revisit this point in section VI D 3 below.

Isospectrality loss: axial versus polar spectral instability
Regge-Wheeler and Zerilli potentials for axial and polar perturbations are known to be isospectral in the QNM spectrum (cf. [43,[83][84][85]; see also [44]). In particular, Chandrasekhar identified (cf. point 28 in [43]) a necessary condition for two (one-dimensional) potentials V 1 (x) and V 2 (x), withx ∈]−∞, ∞[, as the rescaled tortoise coordinate, to have the same transmission amplitude and present the same QNM spectrum. Specifically, both potentials must render the same values when evaluating an infinite hierarchy of integrals with These quantities turn out to be the conserved quantities of the Korteweg-de Vries equation and connect the Schwarzschild QNM isospectrality problem to integrability theory through the inverse scattering transform of Gelfand-Levitan-Marchenko (GLM) theory (cf. [86]; see [85] for an alternative approach in terms of Darboux transformations).
The key point for our spectral stability analysis of L is that axial and polar QNM isospectrality is the consequence of a subtle and "delicate" integrability property of stationary BH solutions, so we do not expect it to be robust under generic perturbations of V . In particular, given the non-linear dependence in V of the conserved quantities C n in (73), we would expect either random δṼ r or deterministic δṼ d perturbations to render different values of C n , therefore resulting in a loss of QNM isospectrality. Fig. 13 confirms this expectation: whereas the fundamental QNM mode remains stable under high-frequency perturbations, isospectrality is broken for the overtones with a slight, but systematic, enhanced damping in the axial case. This provides an interesting probe into QNM instability, with potential observable consequences, and will be the subject of a specifically devoted study elsewhere.

"Infrared instability" of the fundamental QNM
Both the pseudospectrum and the explicit perturbations of the potential indicate a strong spectral stability of the slowest decaying Schwarzschild QNM. This is tension with the results in [1,2], where the instability affects the whole QNM spectrum, this including the slowest decaying QNM. This is a fundamental point to establish, since it directly impacts the dominating frequency in the late BH ringdown signal.
In our understanding, and as it was the case of the Pöschl-Teller potential discussed in section V D 4, the instability of the fundamental QNM frequency found by [1] is an artifact of the implemented perturbations, namely step-like approximations to the Schwarzschild potential (in particular Regge-Wheeler, but the same applies for Zerilli) that modify the potential at large distances. Specifically, V is set to zero beyond [x min , x max ], fundamentally altering the long-range nature of Schwarzschild potential that becomes of compact support. What we observe in Fig. 12 is that keeping a faithful treatment of the asymptotic structure at infinity through the compactified hyperboloidal approach keeps spectral stability.
To test this idea (see also the recent [67]), and as we did in Pöschl-Teller, we have implemented a "cut Schwarzschild" potential in our hyperboloidal approach, setting the potential to zero from a given distance (both towards null infinity and the BH horizon). The result is shown in Fig.14, showing a similar qualitative behaviour to Pöschl-Teller in Fig. 8. Overtones are strongly perturbed into the QNM branches already observed in Fig. 12, consistently with the high-frequencies introduced by the Heaviside cut. But, crucially, now the fundamental QNM is indeed also modified, in contrast with its stable behaviour under high-frequency perturbations. This reinforces the understanding of this effect as a consequence of the "suppression" of the large-scale asymptotics of the potential [87]. However, the observed modification of the fundamental QNM frequency is not as dramatic as the one found in [1]. We do not have a good explanation for this, but it may relate to the fact that the analysis in [1,2] deals directly with Eq. (4), in particular in the setting of a Cauchy slicing getting to spatial infinity i o . Such asymptotic framework may be more sensitive to the modification of the potential that the hyperboloidal one, related to null infinity I + . In this setting, and lacking a better expression, we refer to this effect as an "infrared instability" of the fundamental QNM.
Enforcing the compact support nature of V is naturally mo- The sequence of figures shows a zoom into the perturbation of lowest = 2 axial and polar QNM overtones (the branch cut has been removed), with δṼ d fixed to a value reaching the first overtone, and then increasing the frequency. The breaking of axial and polar isospectrality is demonstrated, with perturbed axial overtones slightly more damped than polar perturbed counterparts, though both laying over the same perturbed QNM branches (actually following the pseudospectra contour lines, cf. Fig. 15 below). The fundamental QNM remains unchanged, consistently with its stability, so the dominating ringdown frequency remains "isospectral".
tivated in physical contexts such as optical cavities, and will be studied systematically in such settings [88]. In gravitation FIG. 14. "Infrared" modification of the Schwarzchild fundamental QNM. As in the Pöschl-Teller case, cutting the Schwarzschild potential ( = 2, either Regge-Wheeler or Zerilli) outside a compact interval [xmin, xmax] modifies the fundamental QNM, this accounting for its "instability" found in [1]. All QNM overtones are strongly perturbed due to the high-frequencies in the Heaviside cut, whereas (only) the fundamental QNM is recovered as xmin, xmax → ∓∞.
the physicality of such an effect is more difficult to assess, since gravity is a long-range interaction that, in contrast to the electromagnetic one, is not screened. In any case, insofar as a pertinent gravitational scenario may be envisaged for a such "cut potential", then the "infrared instability" shown for the first time in [1] would constitute a physical effect.

E. Nollert-Price BH QNM branches: instability and universality
We revisit the results in [1,2] (see also [14,67]), under the light of the elements introduced for the study of QNM spectral stability. Fig. 2 in [1] presents the migration of Schwarzschild QNMs to new branches, as the result of perturbing the (Regge-Wheeler) Schwarzschild potential with a step-like approximation with an increasing number "N st " of steps (cf. Fig. 1 in [1]). A salient feature of Nollert's Fig. 2, further analysed with Price in [2], is that the new QNM branches distribute in a perfectly structured family of lines in the complex plane, unbounded in the real part of the frequency, that "move down" in the complex plane as N st (i.e. the frequency in the perturbation) increases [89]. A comparison with Schwarzschild's pseudospectrum in our Fig. 11 shows two remarkable features: i) the pattern of new the branches found and studied by Nollert and Price is consistent with the contour lines of -pseudospectra, ii) the effect of increasing the frequency perturbation indeed corresponds to an increment in the of the corresponding contour line (namely the "energy size" of the pertubation that, as a H 1 norm, includes the frequency). In other words, Nollert and Price's BH QNM branches seem in-deed to correspond to -pseudospectrum contour lines.
In order to test this picture, we bring our perturbation analysis in section VI D into scene. Fig. 15 presents the superposition of perturbed QNM spectra in Fig. 12 onto the Schwarzschild pseudospectrum in Fig. 11. As in the Pöschl-Teller case, perturbed QNMs distribute along -pseudospectra lines, demonstrating the insight gained above on Nollert's QNM instability by using the pseudospectrum: Nollert-Price QNM branches are identified as actual probes into the analytical structure of the non-perturbed wave operator. Moreover, the correlation of -contour lines with the "size/frequency" of the perturbations, endows the pseudospectrum not only with a explicative but also with a predictive power, as a tool to calibrate the relation between spacetime perturbations and QNM frequency changes. The conceptual frame encoded in Fig. 15 is, in our understanding, the main contribution in this work.

QNM structural stability, universality and asymptotic analysis
Building on Nollert and Price's work, our analysis strongly suggests that BH QNM overtones are indeed structurally unstable under high-frequency perturbations: BH QNM branches migrate to a qualitatively different class of QNM branches. Noticeably and in contrast with this, the pseudospectrum analysis combined with the perturbation tools also suggests that the new class of "Nollert-Price BH QNM branches" presents structural stability features pointing to a kind of 'universality' in the QNM overtone migration pattern.
a. "Universality" in the high-frequency perturbations. The QNM migration pattern seems independent of the detailed nature of the high-frequency perturbation in the Schwarzschild potential. First, such universality is manifested by the similar QNM perturbation pattern produced by very different perturbations: step-like perturbations in [1], the sinusoidal deterministic ones showed in Fig. 15 and also random perturbations (not presented it here due to "blurring", consequence of the "branch cut" contamination). Second, the new branches are consistent with the pseudospectra contour lines, a key point in this universality discussion, since it is completely prior to and independent of perturbations.
b. "Universality" in the potential. Perhaps more importantly, universality seems to go beyond the insensitivity to the nature of the perturbation: it seems to be shared by a whole class of potentials. First, the same pattern of perturbed branches is found in Pöschl-Teller, cf. Fig. 7. More dramatically, Nollert and Price's analysis in [2] is particularly illuminating in this respect. They considered a toy model capturing the effect of a (Dirac-delta) high-frequency perturbation on a BH-like potential, referred to as "truncated dipole potential" (TDP), that contains only two QNMs. Adding the singular (high-frequency) "spike" creates an infinite number of QNMs, again following the QNM branch pattern compatible with our pseudospectra contour lines (cf. Fig. 5 in [2] and see below).
But more noteworthy, and again noticed by Nollert [1], beyond the BH setting the new BH QNM branches are strikingly similar to (curvature) w-modes in neutron-star QNMs (cf. e.g. Fig. 3 in [28] and, more systematically, Ref. [90]). This is remarkable, suggesting that exact but unstable BH QNMs migrate to perturbed but stable QNMs branches whose qualitative pattern may be shared by generic compact objects [91] c. Asymptotic analysis and universality. How to address systematically a possible universality in the qualitative pattern of the perturbed QNM branches? Asymptotic analysis provides a sound approach. The study of the spiked TDP QNMs by Nollert and Price provides an excellent illustration, with the identification of the large-n asymptotic form of perturbed QNM branches, according to the logarithm dependence It is suggestive that this makes direct contact with the possible universality of perturbed BH QNMs and (non-perturbed) QNMs of compact objects evoked above. Indeed, as shown in Ref. [90], w-modes of (a class of) neutron stars present exactly this logarithm pattern [92]. Even more, this makes (an unexpected) contact with Pöschl-Teller, where the spectral instability discussed in section V B 2 is explained [62,[93][94][95] in terms of so-called broad "Regge resonances" (not to confuse with "Regge poles"), precisely described by such a logarithmic dependence [96] and explained in terms of an underlying reduced C 1 regularity. It is tempting to refer to the perturbed BH QNM branches as Nollert-Price-(Regge) QNMs, but this requires an elucidation of the role of the reduced C 1 regularity in the generic perturbations we have studied, in particular including regular (high-frequency) sinusoidal deterministic perturbations (63). In sum, the asymptotic pattern (75) provides a starting point to probe, in gravitational wave signals, the physical properties (e.g. energy, frequency) of small perturbations.
Beyond specific models, this kind of universal behaviour, independent of the high-frequency perturbation detailed nature and for a large class of potentials, invites for systematic semi-classical analyses of highly-damped scattering resonances, in terms of the wave operator principal part [97] and crucially including boundaries. In the spirit adopted in this work, we expect asymptotic tools in the semiclassical analy-sis of the pseudospectrum to provide a systematic approach to assess the universality of perturbed BH QNM branches [98].

Overall perspective on Schwarzschild QNM instability
The main result of this article is summarized in Fig. 15. Specifically, it combines Figs. 10, 11 and 12 to demonstrate QNM spectral (in)stability through the respective three distinct calculations: i) the calculation of the eigenfunctions of the exact spectral problem to calculate condition numbers κ n 's, ii) the evaluation of operator matrix norms to generate the pseudospectrum, and iii) the calculation of eigenvalues of the perturbed spectral problem. Calculations i) and ii) work at the level of the unperturbed problem, whereas iii) deals with the perturbed problem. The three calculations fit consistently through the Bauer-Fike theorem that constrains through Eq. (39) the relation between the pseudospectrum and the tubular regions around the spectrum. They lead to these main results: i) QNM overtones: i.1) QNM overtones are ultraviolet unstable, including the lowest overtones. The pseudospectrum provides a systematic explanatory and predictive framework for QNM spectral instability, confirming the result by Nollert and Price [1,2]. Such instability is indeed realised by physical high-frequency perturbations in the effective potential V , reaching the first overtone for sufficiently high frequencies and/or amplitudes in the perturbation. ii) Slowest decaying (fundamental) QNM: ii.1) The slowest decaying QNM is ultraviolet stable. This feature critically relies on keeping a faithful description of the asymptotic structure at infinity through the compactified hyperboloidal approach. This result is in contrast with conclusions in [1,2], but no contradiction appears since the latter implement a step-potential approximation fundamentally modifying V at large distances, resulting rather in an "infrared probe" into QNMs. ii.
2) The slowest decaying QNM is stable under low and intermediate frequency perturbations in the potential. This property is shared by the whole QNM spectrum. ii.
3) The slowest decaying QNM is "infrared unstable". The instability of the fundamental QNM observed in [1,2] is physical inasmuch as fundamental modifications of the large-distance structure of the potential are allowed.
iii.1) 'Nollert-Price BH QNM branches' as pseudospectrum contour lines. The QNM BH spectrum is ultraviolet structurally unstable, migrating to perturbed branches following -contour lines of pseudospectra. Such migration pattern is largely independent of the detailed nature of the high-frequency perturbations and the potential. Once on such 'Nollert-Price branches', QNMs are spectrally stable. These structural stability properties result in the universality of perturbed QNM branches.
iii.1) QNM isospectrality ultraviolet loss. High-frequency perturbations spoil the integrability of Regge-Wheeler and Zerilli potentials, resulting in a slightly enhanced damping of axial modes with respect to polar ones.

A. Conclusions
We have demonstrated: i) the fundamental BH QNM is stable under high-frequency (ultraviolet) perturbations, while unstable under (infrared) modifications of the asymptotics, the latter consistent with [1]; ii) (all) BH QNM overtones are unstable under high-frequency (ultraviolet) perturbations, quantifiable in terms of the energy content (norm) of the perturbation, extending results in [1,2] to show isospectrality loss; and iii) pseudospectrum contour lines provide the rationale underlying the structurally stable pattern of perturbed 'Nollert-Price QNM BH branches'. Pseudospectra, together with tools from the analysis of non-selfadjoint operators, have revealed the analytic structure underlying such (in)stability properties of BH QNMs, offering an integrating and systematic approach to a priori disparate phenomena. The soundness of the results relies on the use of a compactified hyperboloidal approach to QNMs, combined with accurate spectral numerical methods.

Caveats in the current approach to QNM (in)stability
Beyond the soundness of the results, key questions remain: i) How much does the instability depend on the hyperboloidal approach? In other words, is the instability a property of the equation or rather of the employed scheme to cast it? This is a legitimate and crucial question, requiring specific investigation. In spite of this, we are confident in the soundness of our conclusions: as discussed in detail, the same qualitative behaviour is found systematically by other studies not lying on the hyperboloidal approach, in particular Nollert and Price's pioneer work. Details may change from scheme to scheme, but the (in)stability properties seem robust.
ii) A numerical demonstration is not a proof. Moreover, numerical discretizations introduce their own difficulties and limitations. In particular, spectral issues in the passage from matrix approximations to the actual differential operator is a most delicate question. Again, we are confident in our results, as a consequence of mutual consistency of existing results and non-trivial tests like the ones described in the text. Definitely, proofs will require the use of other methods and techniques.
iii) Could the observed QNM spectral instability be an effect of regularity loss, namely a C 1 effect? It may be, but it is difficult to conclude at this stage. C 1 regularity provides indeed a sufficient condition for logarithmic branches (75) that can be traced to works by Regge [96], Berry [99] or Zworski [93] and manifests in our setting in Nollert & Price's analysis of BH QNM instability [2], broad "Regge resonances" in Pöschl-Teller QNM instability [62,94,95] or in neutron star w-modes [90] (cf. also [100] in related Regge poles). But we also attest the same instability phenomenon for regular sinusoidal perturbations of sufficiently high-frequency. Moreover, the pseudospectrum already informs of the instability (cf. contour lines) at the unperturbed "regular" stage. If high-frequency is actually the basic mechanism, then C 1 would provide a sufficient, but not necessary condition for QNM instability. This point must be addressed.

B. Perspectives
While the pseudospectrum framework is already employed in physics (cf. e.g. [3][4][5][6][7]), there seems to be (up to our knowledge) no systematic application in gravitational physics. The introduction of pseudospectra in gravitational physics opens an avenue to interbreed the study of (in)stability and transients with other domains in physics, by using pseudospectrum analysis as a common methodological frame. In the following we mention some possible lines of exploration in different gravitational settings, from astrophysics to mathematical relativity.

Astrophysics and cosmology
The astrophysical status of the ultraviolet QNM overtone instability, that reaches the lowest overtones for generic perturbations of sufficiently high frequency and energy, requires to assess whether actual astrophysical (and/or fundamental spacetime) perturbations are capable of triggering it. Some problems in which this question is relevant are the following: a) BH spectroscopy. If such instability is actually present, this should be taken into account in current approaches to BH spectroscopy. The stability of the slowest decaying QNM guarantees that the dominating ringdown frequency is unaltered. But regarding QNM overtones, note that we have not referred at all to late time ringing frequencies, but to QNM frequencies: since such two sets of frequencies can actually decouple [1,2,14,[101][102][103] and, as already noticed by Nollert [1], the scattered field itself is not much affected by high-frequency perturbations, finding the signature of perturbed QNMs in the gravitational wave signal may pose a very challenging problem. Awareness of this potential effect may however lead to specifically tailored data analysis tools.
b) BH environment. The arrangement of perturbed QNM branches along (a priori known) -contour lines of pseudospectra opens the possibility of probing, in an 'inverse scattering' spirit, environmental BH perturbations. One can envisage to read the "size" of the physical perturbations by comparing observational QNM data with the "a priori" calibrated pseudospectrum. This may help to assess "dry" versus "wet" BH mergers, a point of cosmological relevance in LISA science.
c) Universality of compact object QNMs. The combination of the "universality" of the perturbed "Nollert-Price QNM BH branches" with Nollert's remark on their similarity to neutron star "w-modes", together with the demonstrated loss of BH QNM isospectrality between polarizations, poses a natural question: do QNM spectra of all generic compact objects share a same pattern?
Schemes such as [21] may provide a systematic frame for the analysis of the astrophysical implications. d) BH QNM (in)stability in generic BHs. A natural and necessary extension of the present work is the study of QNM (in)stability in the full BH Kerr-Newman family, in particular understanding how it intertwines with superradiance instability and the approach to extremality.

Fundamental gravitational physics
We note some possible prospects at the fundamental level: a) (Sub)Planckian-scale physics. Planck scale spacetime fluctuations seem a robust prediction of different models of quantum gravity. They represent "irreductible" ultraviolet perturbations potentially providing a probe into Planck scale physics that, given the universality of BH QNM overtone instability, may be 'agnostic' to an underlying theory of quantum gravity. Such a search of quantum gravity signatures in BH gravitational wave physics is akin to [104]. Actually, it would suffice that a Planck scale "cut-off" induces an effective C 1 regularity in the otherwise smooth low-energy description, to trigger the instability phenomenon. BH QNM instability might then provide a particular probe into 'discreteness' of spacetime (e.g. [105] are references therein).
b) QNMs and (strong) cosmic censorship. In the setting of cosmological BHs, the assessment of the extendibility through the Cauchy horizon in Reissner-Nordström de Sitter is controlled by the parameter β = α/κ − , where α is the spectral gap (the imaginary part of the fundamental QNM in our setting) and κ − is the surface gravity of the Cauchy horizon [106,107]. Therefore, a good understanding in this setting of the (in)stability properties of the slowest decaying QNM, and more generally of the QNM spectrum, may be enlightening to assess the thresholds for Cauchy horizon stability. c) Random perturbations and spacetime semiclassical limit. The "regularization effect" of random perturbations [6,[68][69][70] in the scattering Green's function is an intriguing phenomenon that may play a role in the transition to a semiclassical smooth effective description of fundamental gravitational degrees of freedom. Again, the universality of the phenomenon may play a key role.

Mathematical relativity
The presented numerical evidences need to be transformed into actual proofs. Some mathematical issues to address are: a) Regularity conditions and QNM characterization: the mathematical study of QNMs entails subtle functional analysis issues. In the present hyperboloidal approach this involves, in particular, the choice of appropriate regularity conditions and the associated functional space. This connects our pseudospectrum study with the identification in [35] of the full upper-complex plane as the actual QNM spectrum, if general C ∞ eigenfunctions are allowed. More regularity must therefore be enforced. An analysis along the lines in [39][40][41], where Gevrey classes are identified as the proper functional spaces to define QNMs, is therefore required. Likewise, a systematic comparison with QNM stability in the framework of [38,106] is needed (cf. also [94,95]).
b) Semiclassical analysis and QNM (in)stability. The interest of asymptotic tools, in the study of QNM stability, is twofold. On the one hand, an "asymptotic reasoning" [108] built on the semiclassical analysis of QNMs (a subject taken to full maturity in Sjöstrand's works [109][110][111][112][113]) with a small parameter defined in terms of highly-damped QNM frequencies, can help to assess universality patterns of perturbed Nollert-Price BH QMN branches. On the other hand, asymptotic analysis provides powerful tools to prove rigurously spectral instability and non-trivial pseudospectra (cf. e.g. [114]). In particular, the recent work [115] provides an explicit example of scattering resonance (or QNM) instability, sharing much of the spirit of the discussion in this work.
the whole u = (φ, ψ) vector, this norm is an L 2 -norm coming from the energy scalar product ·, · E (for q(x) > 0) that coincides with (22) upon identification q(x) =Ṽ . Note that γ(x) plays no role in the energy scalar product ·, · E .

Adjoint operator L †
A very important object in our discussion of QNM spectral instability and the pseudospectrum construction is the adjoint L † of the operator L. The definition of L † depends on the choice of scalar product and we shall adopt here the energy scalar product (A7). The full construction of the adjoint L † requires a discussion of its domain of dependence. This is a delicate question intimately linked with the boundary and regularity conditions determining the functional space on which L and L † are defined. This functional analysis issues will be addressed elsewhere, and here we focus on the construction of the so-called "formal adjoint", formally satisfying the relation for all u 1 = (φ 1 , ψ 1 ) and u 2 = (φ 2 , ψ 2 ). Taking into account the definition in Eq. (10) of the operator L, this writes where we have used the expressions for L 1 and L 2 in Eq. (11). Using the energy scalar product (A7) and integrating by parts where we have used p(a) = p(b) = 0, the real character of w(x), p(x), q(x) and γ(x) and the Dirac-delta δ(x) distribution to formally evaluate the boundary terms. This allows us to rewrite so that, introducing the operator L ∂ 2 as in Eq. (25) we can write the formal adjoint in Eqs. (23) and (24) (A13) In general γ(x) does not vanish at the boundaries, so L is not even symmetric and therefore cannot be selfadjoint. Eq.
(A13) identifies neatly the loss of selfadjointness with such non-vanishing γ(x), specifically linking spectral instability with a boundary phenomenon, formally cast through the presence of the Dirac-delta terms. This form also explains the (formal) selfadjoint case L 2 = 0 discussed in section V C 2.
More generally, evaluation of adjoints play a key role in all aspects of our discussion of spectral instability: i) calculation of conditions numbers κ n 's, involving the spectral problem of the adjoint L † , cf. Eq. (26); ii) evaluation of the pseudospectrum, involving the calculation of (generalized) sin-gular values of R L (λ) and therefore the spectral problem of R † L (λ)R L (λ), cf. Eqs. (43) and (B21); and iii) the prescription of the norm ||δṼ || E to in the exploration on perturbed spectral QNM problems, again involving the spectral problem of the operator δṼ † δṼ . Details of the calculation of adjoints in our discretised approach are given in appendices B and C.
Appendix B: Pseudospectrum in the energy norm We derive here the relevant expressions for the construction of pseudospectra in the discretised version of the energy norm.

Scalar product and adjoint
Let us consider a general hermitian-scalar product in C n as with G a positive-definite Hermitian matrix where * denotes conjugate-transpose, i.e. u * =ū t and G * = G t (we notice that in the problem studied in this work, the Hermitian positive-definite matrix G is actually a real symmetric positive-definite matrix G t = G, but we keep the discussion in full generality). Using (B1) and (B2) in the relation characterising the adjoint A † of A with respect to the scalar product (B1), we immediately get 2. Induced matrix norm from a scalar product norm The (vector) norm || · || G in C n associated with the scalar product ·, · G in (B1), namely Indeed, we can write The matrix A † A is unitarily diagonalisable and non-negative definite (that is, x, A † Ax G ≥ 0, ∀x ∈ C), so that we can find an orthonormal basis of eigenvectors {e i } A † Ae i = λ i e i , e i , e j G = δ ij , with real non-negative eigenvalues λ i that we order as Expanding x = i x i e i for an arbitrary x ∈ C n , we write that we can recast as Inserting this in Eq. (B8), we conclude To prove that the inequality is actually saturated, it suffices to show that there exits a vector x, ||x|| G = 1, that realizes the equality, i.e. ||Ax, Ax|| 2 G = ρ(A † A). If we consider x = e n ||Ae n || 2 G = Ae n , Ae n G = A † Ae n , e n G = λ n = ρ(A † A) , and we can finally conclude

Characterization of the pseudospectrum
Given an invertible matrix A ∈ M n (C) and a nonvanishing eigenvalue λ, then 1/λ is an eigenvalue of A −1 and where in the passage from the first line to the second we have used (B16) and the definition (B7) of the spectral radius, whereas in the last equality we have used that a matrix AB has the same eigenvalues as the matrix BA.
We consider now the -pseudospectrum characterisation (32) applied to the discretised energy norm || · || G When choosing the energy scalar product in section IV, that is with G = G E (see explicit expression in appendix C), we recover expression (43) for σ E (A). When using the canonical L 2 product we recover the standard σ 2 (A) in (40) with i ∈ {0, 1, . . . , N }. In the construction of our differential operator L, the interpolant of the product of two functions f and g is obtained then by multiplication on grid points, that is In addition to that, we need an expression for the interpolant of the derivative f N (x) = df dx N . This is determined by with These expressions define the Gram matrix G E for the discretised version of the energy scalar product (22), in the basis determined from the Chebyshev-Lobatto spectral grid.

a. Grid interpolation
An important aspect to observe when performing the numerical integration is that Eq. (C4) is exact whenever the original function f (x) is a polynomial of order ≤ N . With this in mind, and assuming that f (x) and g(x) are polynomials, Eq. (C13) is exact only for the case where product (f g)(x) yields polynomials of order ≤ N . In practical terms, the procedure described above hampers the accuracy of the scalar product's numerical integration whenever the order gets > N .
As an illustrative example, take f (x) = P (x) and g(x) = P (x), with P (x) the Legendre polynomials. Then the integral (C12) -with µ(x) = 1 omitted of the expressionyields I(f, g) = 2δ , /(2 + 1). If we now consider the discrete version I N (f, g) given by Eq. (C13), one observes that the exact result is obtained only for the cases + ≤ N , even though each individual function f (x) and g(x) is exactly represented for ≤ N and ≤ N , respectively.
To mitigate this issue, we modify the integration matrix C N µ -or equivalently the Gram matrix G E -by incorporating the following interpolation strategy.
Given an interpolant vector f N (x i ) associated with a Chebyshev-Lobatto grid {x i } N i=0 , one can obtain a second interpolant vector fN (x i ) associated with another Chebyshev-Lobatto grid {x i }N i=0 with a resolution N =N via Components Iī i of the interpolation matrix I are obtained by evaluating Eq. (C3) at the grid {x i }N i=0 , with the coefficients {c i } N i=0 expressed in terms of f N (x i ) via Eq. (C7). Then Note that the interpolation matrix I has sizeN × N , which reduces to a square matrix only ifN = N . In this case, Eq. (C19) is actually the identity matrix as expected.
Then, for a fixed N , we consider the discrete integration (C13) in terms of a higher resolutionN = 2N and interpolate the expression back to the original resolution N . In other words, defining I N µ (f, g) := IN µ (f, g), we can consider the grid-interpolated new discrete integration where C N µ = I t · CN µ · I or, in terms of its components Going back to the illustrative example where f (x) = P (x) and g(x) = P (x), we now obtain I N (f, g) = 2δ , /(2 +1) exactly whenever , ≤ N . In the same way, we grid-interpolate the Gram matrices that allows to perform the scalar product (C16) via