The central role of spacetime hyperbolas in relativity was recognized already by Minkowski in 1908 [51]. The Milne model from 1933 [52] or Dirac’s point-form of quantum field theory from 1949 are hyperboloidal [53]. In the 1970s, hyperboloidal studies were performed for the analysis of wave equations [54–56] and quantum field theory [57–60]. However, these early studies use a time-dependent formulation in which time freezes at null infinity.

The first hyperboloidal coordinates foliating null infinity are implicit in Penrose’s work on the global causal structure of spacetimes via conformal compactification [61, 62]. Indeed, one can obtain a hyperboloidal surface from any textbook discussing the Penrose diagram simply by looking at the level sets of Penrose time [63]. In the context of numerical relativity, it was recognized that hyperboloidal time functions that asymptotically approach the retarded time should be beneficial for the computation of gravitational waves [64, 65]. Explicit hyperboloidal coordinates in black-hole spacetimes were constructed in the context of the analysis of constant mean curvature foliations [66]. A remarkable but largely ignored paper by Gowdy in 1981 includes many key elements of the hyperboloidal approach used today in black-hole perturbation theory [67], including the height function approach to preserve time-translation symmetry, compactification fixing null infinity (scri-fixing), hyperboloidal solutions to the wave equation, and the structure of time-harmonic solutions relevant for the frequency domain. These ideas were not picked up by the community at the time.

The first systematic study of the hyperboloidal initial value problem for Einstein equations was initiated by Friedrich in 1983 [68]. Friedrich devised a reformulation of the Einstein equations with respect to a conformally rescaled metric that is regular across null infinity. The conformal field equations are well-suited for the analysis of the asymptotic behavior of Einstein’s equations and have led to seminal results such as the nonlinear, semi-global stability of de Sitter-type and Minkowski-type spacetimes [69, 70]. The developments around conformal field equations and attempts to use them numerically are reviewed in [71, 72].

Twenty years after Gowdy’s paper, Moncrief presented the hyperboloidal compactification of Minkowski spacetime using time-shifted hyperboloids in an unpublished talk [73] leading to the first numerical studies using hyperboloidal foliations in Minkowski spacetime [74–77]. Around this time, various suggestions for hyperboloidal coordinates and numerical simulations in black-hole spacetimes were made [78–83].

The construction widely used today in black-hole perturbation theory is based on scri-fixing coordinates with time-shifted hyperboloids presented in 2008 [84]. The idea is to combine the height function technique that preserves the time-symmetry of the underlying spacetime with an explicit radial compactification whose singular Jacobian at infinity is proportional to a prescribed conformal factor. In the following years, this method was used primarily in the time domain for solving wave-propagation problems [85–97].

The translation of the hyperboloidal method to the frequency domain was presented in [27], where it was demonstrated that hyperboloidal time functions regularize the QNM eigenfunctions in the asymptotic domains. Warnick used a related idea in [31] for AdS spacetimes in which spatial slices are naturally hyperbolic (see also [29, 30]). The first detailed analysis of QNMs in asymptotically flat black-hole spacetimes using the hyperboloidal approach was presented in [28]. We summarize the basic ideas of the hyperboloidal approach below.

The construction of globally regular coordinates consists of a time transformation respecting the time symmetry of the background, a suitable spatial compactification, and conformal rescaling [84]. We first introduce the time function τ via [67, 82, 84]. τ = t + h r . (5)

The time transformation implies an exponential scaling in frequency domain [27]. Writing the time-harmonic ansatz in Equation 2 with respect to the new time coordinate in Equation 5, we get u t , r * = e − i ω t R r * = e − i ω τ e i ω h R r * = e − i ω τ R ̄ r * .

The rescaled radial function R ̄ ( r * ) = e i ω h R ( r * ) is regular both near the event horizon and toward null infinity. To see this in an explicit example, consider the height function for the so-called minimal gauge [98–100] h MG r = − r + 2 M ⁡ log r 2 M − 1 − 4 M ⁡ log r 2 M = − r * − 4 M ⁡ log r 2 M .

The minimal gauge height function has the following asymptotic behavior h MG ∼ − r * f o r r → ∞ , h MG ∼ + r * f o r r → 2 M .

The height function regularizes the QNM eigenfunctions in the asymptotic domains. The regularity of the QNM eigenfunctions is directly related to the regularity of the minimal gauge at the asymptotic boundaries near the horizon and near infinity (see Figure 1). The minimal gauge is unique in its simplicity and appears in different setups as a natural construction [101, 102]. Surprisingly, the minimal gauge was implicitly used by Leaver in his papers on QNMs in BH spacetime [103, 104]. Related hyperboloidal regularization procedures have been suggested over the years by various authors without an explicit recognition of the geometric background of their construction [105–111].

In [27], it was shown that the time translation must be combined with a suitable rescaling to arrive at a regular representation of QNMs. The rescaling takes into account the asymptotic fall-off behavior of the QNM eigenfunctions toward the BH and toward infinity. The resulting equations have short-range potentials suitable for compactification of the exterior black region from the radial coordinate r ∈ r h , ∞ — or equivalently, r ∗ ∈ ( − ∞ , + ∞ ) — into a compact domain σ ∈ [ σ h , σ I + ] . This rescaling is related to the conformal compactification of black-hole spacetimes.

The external boundary conditions (Equation 4) are automatically satisfied in terms of a radially compact hyperboloidal coordinates ( τ , σ ) [27, 100], when the underlying function R ̄ ( σ ) is regular at the black-hole horizon ( σ h ) and future null infinity ( σ I + ) . Thus, we no longer need to impose boundary conditions to the wave equation by hand. The boundary condition is replaced by a regularity condition on the underlying solution R ̄ ( σ ) in the entire domain σ ∈ [ σ h , σ I + ] . In practical terms, one derives the regularity condition at the boundary directly from the hyperboloidal differential equation. When formulating the frequency-domain problem in coordinates ( τ , σ ) , the resulting differential equation equivalent to Helmoltz Equation 3 assumes a generic form [28, 99, 100]. α 2 σ d 2 d σ 2 + α 1 σ d d σ + α 0 σ R ̄ σ = 0 , (6) with coefficients α 2 , α 1 and α 0 depending on the particular choice for the hyperboloidal height function. The most important property of the above equation is that it is a singular ordinary differential equation, i. e., its principal part behaves as α 2 ∼ ( σ − σ I + ) 2 ( σ − σ h ) and therefore α 2 ( σ h ) = α 2 ( σ I + ) = 0 . Hence, at the boundaries, Equation 6 provides us directly with the relation between the field and its first σ -derivatives serving as boundary data ensuring a bounded solution. From the spacetime perspective, and the resulting wave equation, the same condition α 2 ( σ h ) = α 2 ( σ I + ) = 0 ensures that, at the boundaries, the light cones point outwards the numerical domain, or equivalently, that the characteristic speeds of incoming modes vanish [27, 100].

The finite behaviour of the function R ̄ ( σ ) is one of the most important aspects in the hyperboloidal approach. As we discuss below, this feature allows us to unveil new properties of the QNM eigenfunctions, develop novel numerical algorithms and attack new problems relevant to black-hole physics and gravitational wave astronomy.

The hyperboloidal framework regularizes solutions to the Helmholtz problem Equation 3. In fact, any bounded solution satisfying the singular ordinary differential Equation 6 automatically fulfills the Sommerfeld conditions (Equation 4). One would naively think that bounded solutions exist only at the QNMs frequencies. However, we saw in Section 2 that the conditions (Equation 4) are necessary, but not sufficient to specify the QNM problem uniquely. In the traditional formulation, the QNM eigenfunctions grow asymptotically. The complex plane spanned by the frequency ω might contain regions with solution satisfying (Equation 4), but contaminated by unwanted solutions that decrease at the boundaries. Removing the asymptotic blow-up allows us to peek directly into these unphysical solutions, which exist in the entire half-plane I m ( ω ) < 0 [28, 31].

The left panel of Figure 2 shows two solutions to Equation 6 which are bounded in the entire exterior BH domain, from σ = 0 representing future null infinity and σ = 1 the BH horizon: the solution in blue is obtained at a given QNM frequency ω QNM , whereas the solution in red corresponds to a given frequency in the half-plane I m ( ω ) < 0 , but with ω ≠ ω QNM . These solutions are obtained with two different numerical approaches: the solid line results from solving the ODE with a Chebyshev collocation point spectral method [112], whereas the dotted points arise when using Leaver’s Taylor expansion, which corresponds to a frequency domain hyperboloidal formulation in the minimal gauge [28, 99, 113]. Both strategies yield the same results. At first glance, there is nothing special in the behavior of the solutions that allows us to distinguish a QNM from a non-QNM eigenfunction. Indeed, it is even possible to specify a hyperboloidal initial data such that the corresponding time evolution has an arbitrary exponentially damped oscillation [28].

What distinguishes a QNM from a non-QNM solution is their regularity class. By studying the convergence rate of their discrete numerical representation, one can infer that these functions belong to different regularity classes. The middle panel of Figure 2 shows the Chebyshev coefficients from the Chebyshev collocation point spectral method. These coefficients decay exponentially for C Ω analytic functions or algebraically for C k singular functions. We observe an intermediary decay, suggesting the regularity class of these functions is between C k and C Ω . The Chebyshev coefficients for QNM eigenfunctions decay faster than for non-QNM eigenfunctions, indicating QNM eigenfunctions belong to a better regularity class. A similar conclusion arises from the Taylor expansion coefficients (right panel 2). For QNM eigenfunctions, | a k | decays asymptotically. For non-QNM functions, | a k | grows asymptotically. Even though the series does not converge absolutely, it converges conditionally due to oscillations in I m ( a k ) (left panel of Figure 2). These conclusions are formalized by interpreting QNM as a formal eigenvalue problem of the generator of time translations for a null foliation, acting on an appropriate Hilbert space [29–33, 114–116], where QNM eigenfunctions belong to the Gevrey-2 regularity class.

Even though challenge (i) mainly concerns the GW detection’s signal-to-noise ratio, it heavily relies on accurate predictions for the expected QNM excitations [117]. The excitation of each QNM depends on the particular initial perturbation triggering the dynamics. This perturbation also excites the late-time power-law tail decay. Determining these excitation factors has always been challenging due to the blow-up of the underlying modes at the bifurcation sphere and spatial infinity [118]. A common approach to avoid the infinities at the bifurcation sphere when calculating integrals along the physical coordinate r ∗ ∈ ( − ∞ , ∞ ) is to deform the integration path into the complex plane [119–121].

The hyperboloidal formalism offers an alternative strategy to determine such excitation factors due to the globally regular behavior of the QNM eigenfunctions. The direct identification of Leaver’s continued fraction strategy with spacetime solutions defined on hyperboloidal hypersurfaces allows the further development of the Leaver method to calculate QNMs (and tail decay) excitation factors for problems formulated on hyperboloidal slices [28, 99]. While Leaver’s method relies on a Taylor expansion around the horizon for the underlying hyperboloidal functions, the strategy can be adapted to directly solve a linear partial differential equation having the QNM excitation amplitude as an unknown parameter in the equation [122], or alternatively via the so-called Keldysh scheme [123]. The hyperboloidal formalism is also essential for recent advances in the understanding of the role played by the tail decay in BH spectroscopy [124–126].

Since GR is a nonlinear theory, challenge (ii) emphasizes that BH spectroscopy must also account for second-order, quadratic perturbations [21, 127–136]. The quadratic coupling of first-order solutions dictates the dynamics at second order in perturbation theory [137–139]. When formulated in the standard t slices, the blow-up of QNM eigenfunctions at the bifurcation sphere and spatial infinity imposes severe restrictions for second-order studies, both at theoretical and numerical levels. The hyperboloidal framework for black-hole perturbations beyond the linear order becomes indispensable for regular evolutions [23, 132, 140], as well as for studies in the frequency domain [134, 141].

Apart from his groundbreaking work in QNM [34], Vishveshwara also highlighted that the QNM spectra is very sensitive to small modifications in the black-hole potential [142, 143]. At the same time, the QNM spectra destabilisation was also observed by Nollter and Price [144, 145], but the phenomenon’s impact in the BH spectroscopy programme has been largely overlooked over the past decades. Only recently has the challenge (iii) gained a greater attention [112, 146–155].

Small modifications in oscillatory frequencies for wave equations result in minor spectral responses only if the wave operators are self-adjoint. However, the flow of GWs into the BH and out into the wave zone places BH perturbation theory within the framework of non-self-adjoint operators. The successful application of non-self-adjoint operator theory to gravitational systems was only made possible by the hyperboloidal approach to black-hole perturbations [112] (see also Ref. [156] for an alternative approach akin to “complex scaling”). In this approach, one can use the mathematical formalism of pseudospectra [157] to study the QNM spectral instability [112] and perform a non-modal analysis [151] that a traditional mode analysis might overlook. Since the breakthrough offered by the hyperboloidal framework, the analysis of QNM pseudospectra has been performed in several different contexts, from astrophysically relevant scenarios to applications in the gauge-gravity duality [158–167].