Stable black holes: in vacuum and beyond

Abstract

Black holes are important objects in our understanding of the universe, as they represent the extreme nature of General Relativity. The mathematics behind them has surprising geometric properties, and their dynamics is governed by hyperbolic partial differential equations. A basic question one may ask is whether these solutions to the Einstein equation are stable under small perturbations, which is a typical requirement to be physically meaningful. We illustrate the main conjectures regarding the stability problem of known black hole solutions and present some recent theorems regarding the fully nonlinear evolution of black holes in the case of vacuum and their interaction with matter fields.

1. Physical motivation

Black holes are astrophysical compact objects which are believed to be overwhelmingly present in the universe. Black holes embody the extreme nature of the theory of General Relativity, and the singularities hidden behind them are thought to be the point of contact of Einstein’s theory of gravity with the theory of quantum mechanics, as quantum effects are believed to be relevant at those scales. A resolution for a quantum theory of gravity will most likely involve a better comprehension of black holes and their singularities. Even though they are now considered central objects in the study of astrophysics, it is only very recently that we have direct and indirect proofs of the existence of black holes.

The first indirect observation of black holes was obtained in 2015 by LIGO (the Laser Interferometer Gravitational-wave Observatory) through the first detection of gravitational waves emitted by the merger of binary black hole systems. The first measurement and the ones that have followed opened up the possibilities of a more direct way of observing and studying these objects. The LIGO interferometers detected the first signal emitted by such an event in September 2015 (see 1), and their discovery was awarded the Nobel Prize in Physics in 2017, just a little more than a year after the announcement. The gravitational waves, which were predicted by Albert Einstein precisely 100 years before (in 1915) in his formulation of the theory of General Relativity, came from a collision between two black holes, which took place 1.3 billion light years away. The event can be described as follows: Two black holes rotate around each other. As they rotate, they lose energy through the emission of gravitational waves. In this process they get closer and closer to each other, while moving almost as fast as the speed of light, until they merge into one larger black hole. After the merger, the final black hole settles down to a stationary state, and in doing that it emits energy through gravitational waves.

What LIGO observed are precisely the waves emitted in this process, and they were detected as a wave-like signal by the two interferometers in Washington and Louisiana. Clearly, in order to extract any relevant information about the event (like the mass of the black holes, their rotations or their distance from us) from the raw detected signal, one needs a very precise model to compare it to. One of the crucial components in the LIGO discovery is in fact the statistical process of comparing the signal observed by the instrument to a catalogue of so-called waveforms, which are numerically simulated signals as expected from the prediction of General Relativity. Remarkably, the first code able to simulate the merger of two black holes was only obtained in 2005 by Frans Pretorius 50, just ten years before the actual detection of gravitational waves by LIGO.

The LIGO interferometers have also detected events involving other astrophysical objects, such as neutron stars 2, where the interaction of different matter fields appear, but nevertheless stability properties are present. These kind of events represent a strong motivation for the study of a nonvacuum solution to the Einstein equation.

Another direct proof of the existence of black holes comes from the images produced by the Event Horizon Telescope Collaboration 13 of the black holes M87 in 2019 and Sagittarius A* in 2022. Both images show a ring of light around the black hole, and its shape is due to properties of the so-called photon sphere outside the black hole, which we will see plays a crucial role in the analysis of stability problems.

LIGO has continued to observe mergers of black holes in recent years, and we are starting to have a much better understanding of the population of black holes in the universe, together with their channel formations and parameter properties. The numerical simulations involved in the process—extremely refined and complex on their own—rely on the mathematical understanding of black holes as solutions to the Einstein equation. In particular, the waves emitted in the merger events are the physical fingerprints of the particular dynamics associated to the Einstein equation: the dispersion properties of the wave equation. The mere detection of gravitational waves emitted by those large perturbations of black holes, in addition to probing their existence, gives us a strong argument for the stability of these objects to perturbations, which is a question that has interested physicists and mathematicians alike in the past 50 years.

2. The mathematics of black holes

The mathematics of black holes was formulated as part of the theory of General Relativity by Albert Einstein in 1915. Before then, the Newtonian theory of gravity was the master theory for describing the movement of celestial bodies, even though some observed phenomena, such as the precession of Mercury’s orbit, could not be explained with Newtonian theory.

According to Newton, the spacetime is given by a three-dimensional flat space and a one-dimensional absolute time. In such spacetime, not all observers agree on the value of the spatial coordinates or the velocity of an object (that depends on the choice of coordinates), but they agree on measurements of distances. In other words, the spatial separation

$$\begin{equation*} ds=\sqrt {dx^2+dy^2+dz^2}, \end{equation*}$$

which is the Euclidean distance, is conserved among different observers.

In 1905 Einstein challenged the Newtonian vision of spacetime, and in his Special Relativity he proposed a theory in which the spacetime is thought to be a four-dimensional flat spacetime, with no absolute slices and no privileged family of observers. To each point of the spacetime, there is an associated cone, called the light-cone, along whose edges light can travel, and everything else, massive objects like us, can only travel slower than the speed of light, therefore being forced to move inside the light-cone. Nothing at all is allowed to move faster than the speed of light, and therefore it cannot travel outside the light-cone. He proposed that the spacetime separation

$$\begin{equation*} ds=\sqrt {-c^2 dt^2+dx^2+dy^2+dz^2}, \end{equation*}$$

where $c$ is the speed of light, is the conserved quantity.

This leads to the definition of Minkowski spacetime, or the spacetime of Special Relativity, as the flat metric on $\mathbb{R}^{3+1}$ with signature $(3, 1)$ given by

$$\begin{equation} g_m=- dt^2+dx^2+dy^2+dz^2, {\tag{2.1}} \end{equation}$$

where from now on the speed of light has been set to $c=1$. Minkowski spacetime is the Lorentzian equivalent of Euclidean space in Riemannian geometry, and its light-cones are all uniform and vertical with slope $1$. Minkowski spacetime represents the empty space, where no masses or objects are present.

It took Einstein about ten years to understand what happens in the presence of a massive object (for instance, to describe the spacetime of a star), and this was the core of his theory of General Relativity. In Minkowski spacetime, an object which travels at constant speed would follow a straight line in spacetime. If a star is present instead, it will deform the spacetime so that the object, in following a straight line (i.e., geodesics), would actually bend toward the star. This fact can be summarized in John Wheeler’s words, “Spacetime tells matter how to move; matter tells spacetime how to curve”. The light-cones of the spacetime of a star are therefore not uniform, but they bend toward the star as they get closer to it.

An interesting phenomenon, which even escaped Einstein’s comprehension for a long time, is that the geometry radically changes if the star becomes more and more massive and dense. If this happens, the spacetime gets distorted: the overall geometry of light-cones changes, and a region where not even light can escape forms. The light-cones become tangent to a confined hypersurface, and in particular if a light signal happens to reach that hypersurface, it will in fact never be able to leave the region enclosed by it. The hypersurface is called the event horizon; the region enclosed by it, the black hole.

2.1. The Einstein equation

From a mathematical point of view, black holes are solutions to the main equation governing the theory of General Relativity, called the Einstein equation. According to General Relativity, a spacetime is a four-dimensional manifold equipped with a Lorentzian metric $g$ (i.e., with signature $(3, 1)$) that satisfies the Einstein equation,⁠¹

¹

Here we set the speed of light and the gravitational constant to unity, i.e., $c=G=1$. In addition, we only consider the case of a zero cosmological constant, $\Lambda =0$.

✖

$$\begin{equation} \operatorname {Ric}(g)-\frac{1}{2} \operatorname {R}(g) g=\operatorname {T}, {\tag{2.2}} \end{equation}$$

where

•

$\operatorname {Ric}(g)$ is the Ricci curvature of $g$,

•

$\operatorname {R}(g)$ is the scalar curvature of $g$,

•

$\operatorname {T}$ is called the stress-energy tensor and contains information about the matter fields present in the spacetime.

The Einstein equation 2.2 is the mathematical representation of the fact that “spacetime tells matter how to move; matter tells spacetime how to curve”: the left-hand side of the equation is a particular combination of curvature of the spacetime, which encodes how curved the spacetime is, while the right-hand side describes the behavior of the matter.

At first glance, equation 2.2 may seem like the definition of the right-hand side $\operatorname {T}$, but this is not how one should think about it. The unknown of the equation is the metric $g$, and even though the left-hand side of 2.2 may seem easy to write, it is actually given by a system of second-order partial differential equations (PDEs) in the metric, as in the definition of Ricci curvature. Those second-order PDEs are sourced by the left-hand side of the equation, which is a given expression containing all the information about the matter present in the universe. In its full generality, the Einstein equation is in fact very difficult to solve and study, and in what follows we will concentrate on two important cases.

The simplest scenario one can consider is the case where there are no matter fields present in the spacetime, i.e., $\operatorname {T}\equiv 0$. A vacuum spacetime is then a spacetime satisfying the Einstein vacuum equation

$$\begin{equation} \operatorname {Ric}(g)=0. {\tag{2.3}} \end{equation}$$

The Einstein vacuum equation is considered the benchmark of the study of General Relativity, and rightly so. In fact, solutions to the Einstein vacuum equation can be shown to have most of the geometrical properties of general spacetimes. Nevertheless, analysis of the resulting equations describing the dynamics is simplified in the vacuum case, as gravitational radiation is not sourced by other radiations and does not interact with matter fields. In the case of the interaction of gravitational radiation with other forms of radiation, various different dynamics take place, depending also on the specific matter field involved.

If we assume that the gravitational field can interact with electromagnetic radiation (if we want to literally add some color to the equations!), we obtain an electrovacuum spacetime, which satisfies the Einstein–Maxwell equation,

$$\begin{equation} \operatorname {Ric}(g)=2 \operatorname {F} \cdot \operatorname {F} - \frac{1}{2} |\operatorname {F}|^2 g, {\tag{2.4}} \end{equation}$$

where $\operatorname {F}$ is a $2$-form, called the electromagnetic tensor, satisfying the Maxwell equations,

$$\begin{equation} d\operatorname {F}=0, \qquad \operatorname {div}\operatorname {F}=0. {\tag{2.5}} \end{equation}$$

In this case, the Ricci curvature is not identically zero, but rather sourced by a quadratic expression in terms of the electromagnetic field, which itself satisfies the equations of electrodynamics 2.5.

The electromagnetic radiation governed by the Maxwell equations represents one of the interesting matter fields for various reasons. As seen in the next subsection, there are explicit black hole solutions to the Einstein–Maxwell equation with a nontrivial electromagnetic field. Moreover, the Maxwell equations themselves have a similar hyperbolic structure as the Einstein equation, and therefore the interaction of electromagnetic and gravitational radiation is very instructive for subsequent studies of interaction with other matter fields.

2.2. Black hole solutions

We now present some particular solutions to the Einstein vacuum and Einstein–Maxwell equations and their properties.

The simplest solution to the Einstein vacuum equation is the Minkowski spacetime $g_m$, given by the metric 2.1, which can also be written in standard spherical coordinates $(t, r, \theta , \phi )$ as

$$\begin{equation} g_m=-dt^2+dr^2+r^2\left( d\theta ^2+ \sin ^2\theta d\phi ^2\right). {\tag{2.6}} \end{equation}$$

As mentioned above, Minkowski is the spacetime of Special Relativity with no presence of matter.

A more interesting solution was discovered by Karl Schwarzschild in 1916 52, just a few months after Einstein wrote down his equation. Schwarzschild wanted to find an explicit solution to the Einstein vacuum equation 2.3, and he looked for one with a high degree of symmetry. By imposing a spherical symmetry on a solution (the spherical part $r^2\left( d\theta ^2+ \sin ^2\theta d\phi ^2\right)$ is also present in the Minkowski metric 2.6), he reduced the problem to solving an ordinary differential equation (ODE). The ODE can be solved explicitly in terms of the radial coordinate, and in the process of integrating for the solution one obtains a constant of integration $M$. The Schwarzschild solution is then a one-parameter family of solutions, parametrized by $M \in \mathbb{R}$, given by

$$\begin{equation} g_M=-\left( 1-\frac{2M}{r} \right) dt^2+\left( 1-\frac{2M}{r} \right)^{-1} dr^2+r^2\left( d\theta ^2+ \sin ^2\theta d\phi ^2\right). {\tag{2.7}} \end{equation}$$

One can immediately notice that if $M=0$, the above metric reduces to the Minkowski metric 2.6. Even more interestingly, it was clear early on that for large $r$ the Schwarzschild metric $g_M$ gives the gravitational potential according to the Newtonian theory of an isolated body of mass $M$. This means that the Schwarzschild solution gives an explicit formula for the equivalent in General Relativity to the spacetime of a star.

One can easily notice that the metric 2.7 has a singular behavior at $r=2M$, where the coefficient $\left( 1-\frac{2M}{r} \right)$ vanishes. Nevertheless, it took long time before the nature of such singularity became clear. At $r=2M$, there is no geometrical invariant which becomes singular, but it is rather a coordinate singularity. This means that there is a change of coordinates which modifies the metric 2.7 into one which is perfectly well-behaved at $r=2M$. Still, the hypersurface $r=2M$ has an interesting geometrical property: the light-cones of the manifold bend toward the region $r<2M$ and become tangent to $r=2M$. In particular, light rays traveling along the edges of those cones, as well as massive objects traveling in the interior of the light-cones, are not able to leave the region $r\leq 2M$ once they enter it. It is a black hole region, and the hypersurface $r=2M$ is an event horizon.

A region of spacetime where not even light can escape was a very difficult concept to grasp for the physics community at that time, and the Schwarzschild solution was indeed only considered to be valid outside the black region, i.e., for $r > 2M$. It was believed that the presence of the black hole region was an artifact of the symmetry of the solution: as in the real world there is no perfectly spherically symmetric massive object, so there would not be any black hole region. The argument went, if one were to consider a more realistic star, a region where not even light can escape would not form. Einstein himself wrote in 1939 22 that the “Schwarzschild singularities [which is what they called the event horizon] do not exist in physical reality”. The entire physics community was convinced that the black hole region was merely a mathematical property of the very symmetric Schwarzschild solution.

The general expectation about the existence of black hole solutions changed dramatically when Roy Kerr in 1963 37 wrote down another explicit solution to the Einstein equation, parametrized by two real parameters $M$ and $a$, where $a$ is a rotation parameter, with $|a|\leq M$.⁠² The Kerr solution is given by

²

If $|a| >M$, the spacetime is said to contain a naked singularity.

✖

$$\begin{equation} g_{M,a}=-\frac{\Delta }{\rho ^2}\left( dt-a\sin ^2\theta d\phi \right)^2+ \frac{\rho ^2}{\Delta }dr^2+\rho ^2d\theta ^2+\frac{\sin ^2\theta }{\rho ^2}\left( a dt-(r^2+a^2) d\phi \right)^2, {\tag{2.8}} \end{equation}$$

where

$$\begin{equation*} \Delta =r^2-2Mr+a^2, \qquad \rho ^2=r^2+a^2\cos ^2\theta , \end{equation*}$$

and for $a=0$ it reduces to the Schwarzschild metric 2.7. The Kerr solution represents a much more realistic case of a massive object of mass $M$ which rotates around an axis of symmetry with angular momentum $Ma$.

Analyzing its geometrical properties, one can see that the hypersurface $\Delta =0$ is still an event horizon: the light-cones become tangent to it and the region enclosed by it is a black hole region. The derivation of such a metric proved to be a game changer for the concept of a black hole, as its existence proved that the black hole region was clearly not simply an artifact of symmetry, and it did not have to be just a mathematical property of an unrealistic physical object. It could in fact very well be real. The Kerr solution is now considered to be the most fundamental black hole solution, and it is believed to represent the astrophysical black holes present in the universe.

The Einstein–Maxwell equation has two explicit solutions which are directly related to Schwarzschild and Kerr. The Reissner–Nordström solution is parametrized by two parameters $M$ and $Q$, where $Q$ is the electric charge, with $|Q| \leq M$, and is given by

$$\begin{equation} g_{M, Q}=-\left( 1-\frac{2M}{r}+\frac{Q^2}{r^2} \right) dt^2+\left( 1-\frac{2M}{r} +\frac{Q^2}{r^2}\right)^{-1} dr^2+r^2\left( d\theta ^2+ \sin ^2\theta d\phi ^2\right). {\tag{2.9}} \end{equation}$$

It represents a spherically symmetric charged black hole.

Finally, the Kerr–Newman solution 49 is parametrized by three parameters $M$, $a$ and $Q$ with $a^2+Q^2 \leq M^2$, and is given by the same expression as in 2.8, where

$$\begin{equation*} \Delta =r^2-2Mr+a^2+Q^2, \qquad \rho ^2=r^2+a^2\cos ^2\theta . \end{equation*}$$

The Kerr–Newman metric is the most general known explicit black hole solution to the Einstein–Maxwell equation, and it is a three-parameter family which describes the gravitational field around an isolated rotating charged black hole of mass $M$, angular momentum $Ma$, and electric charge $Q$. The Kerr–Newman metric generalizes the Reissner–Nordström solution (for $a=0$), the Kerr solution (for $Q=0$), and the Schwarzschild metric (for $Q=a=0$). As such, the Kerr–Newman spacetime plays a fundamental role in describing the final state of evolution in General Relativity.

2.3. The dynamics of black holes

Kerr–Newman solutions are extremely useful for performing explicit computations, such as calculations of geodesics, curvature invariants, etc. Nevertheless, one can immediately notice that the metric coefficients appearing in 2.8 do not depend on the time variable $t$. We call such solutions stationary. In particular, stationary solutions cannot describe events such as the merger of black hole binaries that LIGO observed. In those events, the black holes rotate around each other and continuously change their configuration, in contrast with a black hole described by a stationary solution, such as 2.8, which sits still forever, rotating around its axis in the same way as time passes by.

In order to describe events that are not stationary, we need a mechanism to study the dynamics of black holes as solutions to the Einstein equation. This tool is precisely given by the Choquet-Bruhat theorem 23 which formulated the Einstein equation as an initial value problem, and effectively started the field of Mathematical General Relativity.

Theorem 2.1 (Choquet-Bruhat 23, 1952).

The Einstein equation in wave coordinates is given by a hyperbolic system of PDEs of the form

$$\begin{equation} \square _g g = N(g, \partial g), {\tag{2.10}} \end{equation}$$

where $\square _g = g^{\mu \nu } D_\mu D_\nu$ is the D’Alembertian operator associated to the metric $g$, $D$ is the covariant derivative associated to $g$, and $N(g, \partial g)$ denotes nonlinear terms in $g$ and its first derivative with initial data is given by $(g|_{\Sigma _0}, k|_{\Sigma _0})$, i.e., the metric and its second fundamental form on a spacelike hypersurface $\Sigma _0$, which satisfy some compatibility conditions, called the constraint equations.

The D’Alembertian operator associated to the Minkowski metric is given by the standard wave operator,

$$\begin{equation*} \square _{g_m}= -\partial _t^2+\partial _x^2+\partial _y^2+\partial _z^2, \end{equation*}$$

and the most important properties satisfied by solutions to the standard wave equation, such as finite speed of propagation and continuous dependence on the initial data, can be proved to be valid for solutions to the general covariant wave operator $\square _g$ for any Lorentzian metric $g$.

In particular, Theorem 2.1 implies local well-posedness and continuous dependence on the initial data for the Einstein equation. More precisely, given a set of compatible initial data for the Einstein equation on some initial time slice $\Sigma _0$, one can uniquely solve locally the Einstein equation in the future of $\Sigma _0$. This gives a mechanism for studying the dynamics of the Einstein equation. Suppose we are given an initial data set that describes two black holes that rotate around each other. Then we can use the Choquet-Bruhat theorem to solve locally for the Einstein equation and describe the behavior of the solution after a very short period of time.

What about the global behavior of the solution, i.e., the behavior of the solution a long time after the initial time slice? This is a much more difficult question that Theorem 2.1 cannot answer. For example, consider a much simpler case than the Einstein equation: a nonlinear scalar wave equation of the form

$$\begin{equation} \square _{g_m}\phi = (\partial _t \phi )^2. {\tag{2.11}} \end{equation}$$

If we consider equation 2.11 with initial data $\phi |_{t=0}=\partial _t \phi |_{t=0}=0$, then, by uniqueness, the solution is given by $\phi (t)=0$ for all $t \geq 0$, and it is therefore global, smooth, and bounded in time. If instead we consider equation 2.11 with initial data $\phi |_{t=0}=\partial _t \phi |_{t=0}=\epsilon$, for some $\epsilon >0$, there exists a finite time $T$ such that $\phi \to \infty$ for $t \to T$. The solution is then not defined globally in time and in fact blows up at finite time.

We say that the trivial solution to equation 2.11 is not stable under small perturbations: if the initial data changes of size $\epsilon$, for any small $\epsilon$, then the behavior of the solution changes dramatically, passing from one which is bounded for all times to one which blows up in finite time.

Since the Einstein equation 2.10, as given by the Choquet-Bruhat theorem, has in principle the same schematic structure as equation 2.11 (albeit more complicated as it involves tensorial quantities and quasi-linear terms), one may worry that perturbations of even the trivial solution (i.e., the Minkowski spacetime) could blow up in finite time. Surprisingly, this does not happen, as proved by Christodoulou and Klainerman in 1993 12.

Theorem 2.2 (Christodoulou and Klainerman 12, 1993).

The Minkowski spacetime $g_m$ is globally nonlinearly stable as a solution to the Einstein vacuum equation.

One of the underlying reasons that the Einstein equation behaves differently than equation 2.11, whose trivial solution is not stable under small perturbations, is the absence of nonlinear terms of the form $(\partial _t \phi )^2$ on the right-hand side. This is called the Klainerman null condition 39. The monumental proof of the Christodoulou–Klainerman theorem does not depend on coordinates, but is rather geometrical as it lays open the structures of the spacetimes which are constructed as perturbations of Minkowski spacetime.

The Christodoulou–Klainerman theorem answers affirmatively the question of stability of the trivial solution to the Einstein equation. What about the global behavior of perturbations of nontrivial solutions to the Einstein equation, such as black holes? If we consider initial data given in a black hole background and a pulse of radiation at a given time, then we know by the Choquet-Bruhat theorem that it will propagate with finite speed as a wave. We can expect some of the radiation to fall inside the black hole and some to disperse toward the far-away region. What about the bulk of radiation in between?

One of the main differences between the Minkowski space and the black hole solutions is the presence of a hypersurface, called the photon sphere, outside the event horizon, where null geodesics (i.e., geodesics which travel along the edges of the light-cones) tend to concentrate. Those are called trapped null geodesics, as they are trapped on this bounded region for a long time. In the case of Schwarzschild spacetime, the photon sphere occurs at $r=3M$, outside the event horizon $r=2M$. In the case of Kerr, the structure of the photon sphere is more involved, as it does not concentrate in one hypersurface in physical space.

Because of the presence of the photon sphere, one may worry that the radiation in perturbations of a black hole will tend to accumulate on the photon sphere region, creating a concentration of energy that will cause the whole spacetime to explode or blow up in finite time.

Surprisingly, once again, this does not happen, as the trapped null geodesics which concentrate around the photon sphere outside the black hole solutions are unstable: in particular, they tend to scatter off and disperse after some time in the region. The particles of light, which start orbiting around the black hole along the photon sphere and then gets scattered off and disperse away, are precisely the ones observed by the Event Horizon Telescope, and they form the luminous ring outside the black hole depicted in their images 13. The instability of the photon sphere is one of the important properties used to prove the stability of the wave equation on black holes.

3. The stability problem for the Einstein equation

The stability problem for black hole solutions concerns the long-time behavior of solutions to the Einstein equation which are perturbations of the known family of black hole solutions. We now present the stability problem as one aspect in the bigger context of the Final State conjecture, which aims to describe the state of the evolution of any reasonable initial data for the Einstein equation.

Conjecture 3.1 (Final State Conjecture 38).

Reasonable initial data for the Einstein equation evolve asymptotically in time to a finite number of Kerr–Newman black holes, moving away from each other.

The Final State Conjecture as stated here is very general, and its proof is out-of-reach in the foreseeable future. Nevertheless, this conjecture implies three sub-conjectures, each of which is interesting in its own right and for which progress has been made.

(1)

Collapse conjecture. The Final State Conjecture predicts that black holes are the state of evolution of any initial data for the Einstein equation. But how do black holes form in the first place? The Collapse conjecture states that there is a mechanism for which large initial data give rise to the formation of a black hole.⁠³ For more detail on the Collapse conjecture, see 44 and 5.

³

Small initial data do not, as implied by the Christodoulou–Klainerman theorem of the stability of Minkowski space.

✖

(2)

Rigidity conjecture. The Final State Conjecture predicts that initial data for the Einstein equation evolve into a very special form of black holes, i.e., the Kerr–Newman family given by the metric 2.8 with $\Delta =r^2-2Mr+a^2+Q^2$. One may ask why one should expect this form of the metric, and not any other one (which we may not have discovered yet!) to be precisely the final state of evolution. The Rigidity conjecture states that the Kerr–Newman family is in fact the only family of stationary solutions to the Einstein equation. The Rigidity conjecture is proved in the case of analytic solutions, and is known as the no-hair theorem. For more detail on the Rigidity conjecture, see 17.

(3)

Stability conjecture. The Final State Conjecture predicts that the evolution of initial data settles down to a member of the Kerr–Newman family, in particular implying that all the extra radiation will not concentrate in the process and instead disperse to infinity. This implies the Stability conjecture, which states that the Kerr–Newman family is stable under small perturbations of the initial data.

We now focus on the Stability conjecture.

Conjecture 3.2 (Stability of the Kerr–Newman family conjecture).

Initial data for the Einstein equation, which are sufficiently close to a Kerr–Newman black hole, evolve asymptotically in time to another member of the Kerr–Newman family.

Observe that the conjecture does not claim that one member of the Kerr–Newman family is stable, but it is instead that the whole family is stable. Physically, this corresponds to the fact that when perturbing a black hole of a certain mass, rotation and charge, one expects the perturbation to modify its mass, rotation, or charge. Nevertheless, the overall structure of its metric is expected to be unchanged.

We now translate the Stability conjecture into the language of PDEs. Let the Einstein equation be represented by the nonlinear operator

$$\begin{equation} \mathcal{P}(\phi )=0, {\tag{3.1}} \end{equation}$$

and let $\phi _{\lambda }$ be a family of stationary solutions (the Kerr–Newman family) parametrized by some parameter $\lambda$ (given by $M, a, Q$), i.e., $\mathcal{P}(\phi _{\lambda })=0$.

Proving that the family of solutions $\phi _{\lambda }$ is stable under small perturbations as a solution to 3.1 boils down to proving that a solution with initial data close to $\phi _{\lambda }$ converges asymptotically in time to $\phi _{\lambda '}$ with $\lambda '$ close to $\lambda$.

There are various levels of increasing difficulty for the stability problem:

(1)

Consider the linearized equation$$\begin{equation} (d\mathcal{P})|_{\phi _{\lambda }}(\psi )=0, {\tag{3.2}} \end{equation}$$

and prove the following.

•

Separated solutions of 3.2, of the form$$\begin{equation} \psi ( t, r, \theta , \phi )= e^{- i \omega t} e^{i m\phi } R(r) S(\theta ) {\tag{3.3}} \end{equation}$$

for $\omega \in \mathbb{C}$, $m\in \mathbb{Z}$ do not exponentially grow in time. This is called mode stability.

•

All solutions of 3.2 decay in time. This is the full linear stability.

(2)

Prove that all solutions to the fully nonlinear equation 3.1 decay in time. This is the full nonlinear stability.

Why is the Stability portion of the Kerr–Newman family conjecture believed to be true? There have been clues in the past 50 years which have been pointing toward the expectation that the conjecture would hold. The clues concern the following three aspects of the problem.

•

Mode stability. There are no exponentially growing modes for separated mode solutions of the form 3.3 to the Einstein equation in Schwarzschild (Regge and Wheeler 51, Bardeen and Press 6), in Reissner and Nordström (Moncrief 47, Chandrasekhar 11), and in Kerr (Teukolsky 57, Chandrasekhar 11, Whiting 58, Shlapentokh and Rothman 54). No mode stability has been proved for Kerr–Newman spacetime, and we will come back to this anomaly in Section 5.4.

•

Stability for the scalar wave equation on black hole backgrounds. General solutions arising from regular initial data to$$\begin{equation} \square _g \psi =0 {\tag{3.4}} \end{equation}$$

remain bounded and decay in time in Schwarzschild and in Reissner and Nordström (Kay and Wald 36, Blue and Soffer 8, Dafermos and Rodnianski 20), in Kerr and in Kerr and Newman (Dafermos and Rodnianski 20, Tataru and Tohaneanu 56, Andersson and Blue 3, Dafermos, Rodnianski, Shlapentokh, and Rothman 21).

Extensive progress has been made in the last 15 years which allowed us to go beyond mode analysis, to tackling full linear stability (b) for the linear wave equation. A robust geometric interpretation of the redshift effect 18, a physical space analysis of the trapping region and the superradiance 20, and a hierarchy of $r$-weighted decay 19 all contributed to a complete understanding of the boundedness of solutions to the linear wave equation.

•

Nonlinear stability of Minkowski spacetime. Solutions to the fully nonlinear Einstein vacuum equation, which are small perturbations of Minkowski give rise to a complete spacetime which converges to Minkowski space (see Christodoulou and Klainerman 12, Lindblad and Rodnianski 45, and Bieri 7).

3.1. Known results in stability for the Einstein equation

We now collect some of the known results for full linear and nonlinear stability of black hole solutions to the Einstein vacuum and Einstein–Maxwell equation.

The Schwarzschild spacetime has been proved to be linearly stable to gravitational perturbations (Dafermos, Holzegel, and Rodnianski 16, Keller, Hung, and Wang 34, Johnson 35, Hung 33), and nonlinearly stable under the symmetry class of axially symmetric polarized perturbations (Klainerman and Szeftel 42), and up to a codimension-3 submanifold of moduli space (Dafermos, Holzegel, Rodnianski, and Taylor 14).

The Kerr spacetime has been proved to be linearly stable using Newman–Penrose formalism (Andersson, Bäckdahl, Blue, and Ma 4) and using harmonic gauge (Häfner, Hintz, and Vasy 32) in both cases for $|a| \ll M$. Other important results also concern the proof of boundedness and decay for solutions to the Teukolsky equation for $|a| \ll M$ (Ma 46 and Dafermos, Holzegel, and Rodnianski 15) and some progress for $|a| < M$ (Shlapentokh, Rothman, and Teixeira da Costa 55). The nonlinear stability of Kerr spacetime for $|a| \ll M$ has been very recently obtained in a combination of results by Klainerman and Szeftel 40, 41, 43, Shen 53, and Giorgi, Klainerman, and Szeftel 31.

The Reissner–Nordström spacetime has been proved to be linearly stable to electromagnetic-gravitational perturbations for $|Q|<M$ in our series of work, 27, 25, 26, and 24. We will discuss later the partial results on boundedness and decay for solutions to the Teukolsky system in Kerr–Newman spacetime for $|a| \ll M$ as obtained in our recent work, 29 and 28.

We state here the statement of the linear stability of Reissner–Nordström spacetime, as obtained in our 27, in order to illustrate what kind of statement is expected from the stability results quoted above.

Theorem 3.3 (Giorgi 27, 2020).

All solutions to the linearized Einstein–Maxwell equations around Reissner–Nordström spacetime for $|Q|<M$ arising from regular initial data

(1)

remain uniformly bounded on the exterior and

(2)

decay to a linearized Kerr–Newman solution

after adding a pure gauge solution⁠⁴ which can itself be estimated by the initial data.

⁴

Pure gauge solutions are those derived from linearizing the families of metrics that arise from applying the Bondi gauge to Reissner–Nordström families of coordinate transformations which preserve the choice of gauge on the metric (in the case of 27).

✖

We emphasize here two instabilities appearing in the statement of the theorem: the linearized Kerr–Newman and the pure gauge solutions.

We already mentioned that upon perturbing a black hole solution, the parameters of the spacetime can change. In particular, when perturbing a Reissner–Nordström spacetime (i.e., a Kerr–Newman metric with parameter $M$, $Q$ and $a=0$), we should expect the final perturbations to converge to a Kerr–Newman metric with small $a$. In linear theory, this corresponds to solutions to the linearized Einstein–Maxwell equation around Reissner–Nordström which decay to a linearized Kerr–Newman solution.

There is a more serious instability which is proper of the Einstein equation. Being a tensorial equation for the Ricci tensor, the equation is invariant under choice of coordinates or gauge. In particular, given any metric $g$ solution to the Einstein equation, the pullback of the metric through any diffeomorphism is also a solution to the same Einstein equation. For this reason, in order to prove any decay of the solution, we will effectively have to pick a gauge and control the solution in that gauge. In linear theory, this corresponds to solutions to the linearized Einstein equation which decay in time only up to a pure gauge solution which represents the choice of gauge.

In what follows, we will introduce three important ingredients of many stability problems for the Einstein equation: the formalism of Christodoulou and Klainerman, the Teukolsky equation, and the Chandrasekhar transformation. These concepts are all used in the works on nonlinear stability of Kerr solutions and the linear stability of charged black holes.

3.2. The formalism of Christodoulou and Klainerman

Most of the above-mentioned works on the stability of black holes rely on the formalism first introduced by Christodolou and Klainerman 12, which has been very powerful in analyzing geometrical problems in General Relativity. As an alternative to the Choquet-Bruhat approach in wave coordinates, Christodoulou–Klainerman formalism makes use of null frames, i.e., vector fields which are geometrically defined along the manifold and have the following properties.

Let $(M, g)$ be a $3+1$-dimensional Lorentzian manifold solution to the Einstein equation, and let $D$ be the covariant derivative associated to $g$. Suppose that $(M, g)$ can be foliated by spacelike $2$-surfaces $(S, \mathbf{g})$, where $\mathbf{g}$ is the pullback of the metric $g$ to $S$. To each point of $M$, we can associate a null frame $\left\{ e_3, e_4, e_a\right\}$, with $\{ e_a \}_{a=1,2}$ being tangent vectors to $(S,\mathbf{g})$ such that

$$\begin{equation*} g\left(e_3, e_3 \right) = 0, \qquad g\left(e_4, e_4 \right) = 0, \qquad g\left(e_3,e_4\right) = -2, \end{equation*}$$

$$\begin{equation*} g\left(e_3, e_a \right) = 0, \qquad g\left(e_4, e_a\right) = 0, \qquad g \left(e_a, e_b \right) = \mathbf{g}_{ab}. \end{equation*}$$ We denote $\nabla _3=\nabla _{e_3}$ and $\nabla _4=\nabla _{e_4}$ the projection to $S$ of the spacetime covariant derivatives $D_{e_3}$ and $D_{e_4}$, respectively. We then project all geometric quantities along the null frame $\left\{ e_3, e_4\right\}$ and obtain tensors on the sphere $S$. For the Ricci coefficients (or Christoffel symbols) we have for example

$$\begin{equation} \chi _{ab}\coloneq g(D_a e_4, e_b), \qquad \underline{\chi }_{ab}\coloneq g(D_a e_3, e_b), {\tag{3.5}} \end{equation}$$

which can be interpreted as the second fundamental form of the embedding of the sphere with respect to the normal $e_4$ and $e_3$, respectively. The quantities $\chi$ and $\underline{\chi }$ are 2-tensors on the sphere. For the Weyl curvature, we have for example

$$\begin{equation*} \alpha _{ab}=W(e_a, e_4, e_b, e_4), \qquad \underline{\alpha }_{ab}\coloneq W(e_a, e_3, e_b, e_3), \end{equation*}$$

which are the extreme null curvature components of the spacetime, and are 2-tensors on the sphere. For the electromagnetic tensor, the extreme null components are given by

$$\begin{equation*} {}^{(F)}\beta _a\coloneq F(e_a, e_4), \qquad {}^{(F)}\underline{\beta }_{a}\coloneq F(e_a, e_3), \end{equation*}$$

which are 1-forms on the sphere.

In the proof of linear stability of Reissner–Nordström in 27, we make use of the above-described formalism. On the other hand in the case of Kerr or Kerr–Newman, the above formalism, which in particular relies on the foliation of spacetime in spheres, needs to be extended to the case of nonintegrable frame; see Section 4.1 for more detail.

The Christodoulou–Klainerman formalism was introduced in 12 for the proof of nonlinear stability of Minkowski space. In order to study the Einstein equation in this formalism, one makes use of the fact that the Einstein equation is equivalent to the Bianchi identities for the Riemann curvature, and therefore by projecting those equations to the spheres, one obtains a large number of tensorial equations on the spheres.

The proof of stability of Minkowski space is based on a robust vector field method, used to derive quantitative decay based on generalized energy estimates and commutation with an approximate Killing vector field. As mentioned above, for the nonlinear part of the equations, the null condition is used to identify the structure of the nonlinearity which enables the stability mechanism. Finally, the overall setup of the proof is an elaborate bootstrap argument, according to which one makes educated guesses about the behavior of the solution, which is improved by using a long sequence of a priori estimates.

3.3. The Teukolsky equation

In the case of perturbations of black hole solutions, the projected equations for the curvature and the Ricci coefficients are a coupled system of transport, elliptic, and hyperbolic equations with a complicated structure. One may wonder why such an approach, based on null frames, would be more successful or convenient with respect to the one in wave coordinates.

The main reason to use null frames is that the Kerr–Newman family admits a special frame, called a principal null frame, which diagonalizes the Weyl curvature. This means that for

$$\begin{equation} e_3=\frac{r^2+a^2}{\rho ^2}\partial _t-\frac{\Delta }{\rho ^2}\partial _r+\frac{a}{\rho ^2}\partial _\phi , \qquad e_4=\frac{r^2+a^2}{\Delta }\partial _t+\partial _r+\frac{a}{\Delta }\partial _\phi , {\tag{3.6}} \end{equation}$$

the Weyl curvature becomes

$$\begin{equation*} W_{a4b4}=W_{a3b3}=W_{a343}=W_{a434}=0, \qquad W_{abcd}, W_{3434}\neq 0, \end{equation*}$$

with most vanishing components of the curvature. Since it is easier to linearize around zero, it is convenient to write the many equations on the spheres with respect to the principal null frame. In fact, they simplify dramatically and become tractable: more precisely, in vacuum the symmetric 2-tensors on the spheres $\alpha _{ab}$ and $\underline{\alpha }_{ab}$ satisfy a second-order PDE which is wave-like and decouples from all the other components of the linearization.

The symmetric 2-tensor $\alpha$ satisfies the following wave-like equation, called the Teukolsky equation, as discovered by Teukolsky in 1972 57 and given by

$$\begin{equation*} \mathcal{T}(\alpha )\coloneq \square _g \alpha + c_1(r, \theta ) \nabla _{\partial _r} \alpha + c_2(r, \theta )\nabla _{\partial _t} \alpha + c_3(r, \theta )\nabla _{ \partial _{\phi } }\alpha - V(r, \theta ) \alpha =0. \end{equation*}$$

In the physics community, the Teukolsky equation is expressed in terms of a complex scalar of spin $\pm 2$ (for gravitational perturbations) or of spin $\pm 1$ (for electromagnetic perturbations).

The highest-order terms of the Teukolsky equation are given by a D’Alembertian operator $\square _g$, followed by first-order terms in $\partial _r$, $\partial _t$, $\partial _\phi$, and a potential. One may therefore expect that the techniques developed for the standard wave equation 3.4 can be applied to the Teukolsky equation to obtain estimates for its general solutions.

Unfortunately, that is not the case for the following reason. The boundedness of the energy for $\square _{g_m} \psi =0$ in Minkowski spacetime is obtained by multiplying the wave equation by $\partial _t \psi$ and integrating by parts. Schematically, one obtains

$$\begin{align*} 0=\square _{g_m} \psi \cdot \partial _t \psi &= \big (-\partial _t^2 \psi + \triangle \psi \big )\cdot \partial _t \psi \\ &= -\partial _t^2 \psi \cdot \partial _t \psi - \partial _t \nabla \psi \cdot \nabla \psi +\text{boundary terms}\\ &= -\frac{1}{2} \partial _t ( |\partial _t \psi |^2+ |\nabla \psi |^2 )+\text{boundary terms}, \end{align*}$$

which gives conservation of the energy density $|\partial _t \psi |^2+ |\nabla \psi |^2$. Similarly for a wave equation with a positive potential $V$, i.e., $\square _{g_m} \psi - V \psi =0$, one obtains

$$\begin{align*} 0&= \big (\square _{g_m} \psi - V \psi \big ) \cdot \partial _t \psi \\ &= -\frac{1}{2}\partial _t ( |\partial _t \psi |^2+ |\nabla \psi |^2 + V |\psi |^2)+\text{boundary terms}, \end{align*}$$

which gives conservation of the coercive energy density $|\partial _t \psi |^2+ |\nabla \psi |^2 + V |\psi |^2$. A wave equation of the form $\square _g \psi - V \psi =0$ is called a Regge–Wheeler equation.

For a general Teukolsky equation instead, because of the presence of the first-order terms $c_1 \partial _r+ c_2\partial _t + c_3 \partial _\phi$ in the equation, one cannot directly obtain boundedness of the energy as for the standard wave equation.⁠⁵

⁵

Even for the standard wave equation in Kerr spacetime, the boundedness of the energy cannot be obtained simply by using the vector field $\partial _t$ due to superradiance. Such an issue for the scalar wave equation in Kerr can nevertheless be overcome; see 21.

✖

3.4. The Chandrasekhar transformation

This issue appears already in the study of the Teukolsky equation in Schwarzschild spacetime, as one would like to pass from a Teukolsky equation of the form

$$\begin{equation} \square _g \alpha - V(r) \alpha =c_1(r) \partial _r \alpha + c_2(r) \partial _\phi \alpha + c_3(r) \partial _t \alpha {\tag{3.7}} \end{equation}$$

to a Regge–Wheeler equation of the form

$$\begin{equation*} \square _g \mathfrak{q}- V(r)\mathfrak{q}=0 \end{equation*}$$

with a positive potential. Such a transformation was described by Chandrasekhar 11 in the setting of mode stability, and Dafermos, Holzegel, and Rodnianski 16 introduced a physical-space version of the transformation in Schwarzschild, valid for general solutions to the linearized Teukolsky equation.

The Chandrasekhar transformation consists in taking two null derivatives⁠⁶ along the incoming direction $e_3$ of the curvature component $\alpha$. Schematically,

⁶

This surprising transformation, consisting of taking two derivatives of a bad equation and expecting to obtain a good equation, relies on a hidden relation between the equation for curvature perturbations (the Teukolsky equation) and the equation for metric perturbations (the Regge–Wheeler equation) in black hold perturbation theory.

✖

$$\begin{equation*} \mathfrak{q}= f_1(r) \nabla _3 \big ( f_2(r) \nabla _3 \alpha \big ) \end{equation*}$$

for carefully chosen functions $f_1$ and $f_2$.

The generalized form of the Chandrasekhar transformation is a crucial component of the analysis of the nonlinear Teukolsky equation of the slowly rotating Kerr black hole solutions in vacuum as obtained in our joint work with Klainerman and Szeftel 31. It also plays an important role in the proof of the linear stability of the family of charged Kerr–Newman black holes as obtained in 29 and 28.

4. The nonlinear stability of the slowly rotating Kerr black holes

The proof of the nonlinear stability of the slowly rotating Kerr black holes has recently been obtained as a combination of results by Klainerman and Szeftel 40, 41, 43, Giorgi, Klainerman, and Szeftel 31, and Shen 53. The overall setup of the proof is the construction of the evolution as a limit of continuous spacetimes, where a maximal time of existence for the evolution of the perturbation of the initial data for a slowly rotating Kerr black hole is assumed, together with a list of bootstrap assumptions on the finite region of existence. The proof then comprises various parts, all connected with each other, as it is proper of a purely nonlinear problem. The major ingredients of the proof are the following:

(1)

A formalism to derive tensorial versions of the Teukolsky and generalized Regge–Wheeler (gRW) equations in the full nonlinear setting.

(2)

An analytic mechanism to derive combined Morawetz-energy estimates for the gRW equations, based on a far-reaching extension of the Andersson–Blue method 3.

(3)

A dynamical mechanism for finding the right gauge conditions based on GCM (generally covariant modulated) spheres and hypersurfaces, in which convergence to the final state takes place, and which are crucially used to obtain estimates for the gauge-dependent quantities in the evolution.

For more detail on the structure of the proof, see the introductions of 41 and 31.

Here we will focus on the wave equation estimates in ingredient (2) used to improve the bootstrap estimates as obtained in 31. We will first present how the Christodoulou–Klainerman formalism is extended in order to study perturbations of the Kerr solution, and we then describe the generalized Regge–Wheeler equations of Kerr.

4.1. The formalism in Giorgi, Klainerman, and Szeftel

Recall the null frame formalism introduced by Christodoulou and Klainerman, with the foliation of the spacetime in 2-surfaces and null frames associated to each point of the spacetime.

In passing from Schwarzschild to Kerr or to Kerr and Newman, there is a crucial geometrical difference: the principal null frame 3.6 in Kerr(–Newman) is not integrable in Frobenius’s sense. This means that the subspace of the tangent space which is orthogonal to $e_3$ and $e_4$ is not tangent to a surface, but rather a horizontal distribution or 2-plane field.

For this reason, in 30 we extended the Christodoulou-Klainerman formalism described in Section 3.2 to the case of a nonintegrable horizontal distribution. The definitions of Ricci coefficients and curvature components can be carried through in this case, but one has to be careful about some symmetries which do not hold in the case of a nonintegrable frame.

For example, the Ricci coefficient $\chi _{ab}$ defined in 3.5 is not in general symmetric in its indices, but could contain an antisymmetric part, which we denote $\,^{(a)}\text{tr}\chi$ in 30, i.e.,

$$\begin{equation*} \chi _{ab}=\widehat{\chi }_{ab} +\frac{1}{2} \text{tr}\, \chi \delta _{ab} + \frac{1}{2} \,^{(a)}\text{tr}\chi \in _{ab}. \end{equation*}$$

The antisymmetric component $\,^{(a)}\text{tr}\chi$ measures the failure of integrability of the principal null frame. Because of the presence of the antisymmetric terms, it is convenient to complexify all the curvature and Ricci coefficients.

4.2. The generalized Regge–Wheeler equation in perturbations of Kerr

The complexified curvature component $A=\alpha + i * \alpha$ satisfies the Teukolsky equation

$$\begin{equation*} \mathcal{L}[A]=\operatorname {Err}[\mathcal{L}[A]], \end{equation*}$$

where $\mathcal{L}[A]$ is a second-order tensorial wave operator on the perturbed spacetime, and $\operatorname {Err}[\mathcal{L}[A]]$ are nonlinear errors depending on all the linearized Ricci and curvature coefficients. The above Teukolsky equation reduces to the Teukolsky equation 3.7 in the case of linear perturbations of Schwarzschild.

In the case of linear perturbations of Kerr, Dafermos, Holzegel, and Rodnianski 15 and Ma 46 defined the Chandrasekhar transformation, given schematically by

$$\begin{equation*} \mathfrak{q}= f_1(r, \theta ) \nabla _3 \big ( f_2(r, \theta ) \nabla _3 A\big ), \end{equation*}$$

for suitably chosen functions $f_1$ and $f_2$. In linear theory, the quantity $\mathfrak{q}$ satisfies the following gRW equation:

$$\begin{equation*} \begin{split} &\square _g \mathfrak{q}- V(r, \theta ) \mathfrak{q}-i \frac{4a\cos \theta }{\rho ^2} \nabla _{\partial _t} \mathfrak{q}\\ &\qquad =a \big [ d_1(r, \theta ) \nabla _{\partial _\phi } \nabla _3 A +d_2(r, \theta ) \nabla _{\partial _\phi } A + d_3(r, \theta ) \nabla _3 A+d_4(r, \theta ) A \big ], \end{split} \end{equation*}$$

where $V(r, \theta )$ is a positive real potential. Despite the presence of the first-order term $i \frac{4a\cos \theta }{\rho ^2} \nabla _{\partial _t} \mathfrak{q}$, the left-hand side of the gRW equation still has good divergence properties as the standard Regge–Wheeler equation in Schwarzschild, while the right-hand side of the equation can be treated as lower-order terms for $|a| \ll M$.

Just as in linear theory, to be able to control $A$ we need to perform transformations $\mathfrak{q}= \mathfrak{q}[A]$ which take solutions $A$ of the Teukolsky equation into solutions of nonlinear tensorial versions of Regge–Wheeler equations, which we call gRW equations. In nonlinear theory, the definition of the transformation $\mathfrak{q}[A]$ is more subtle to avoid nonacceptable error terms in the process, and it satisfies

$$\begin{equation} \square _{\mathbf{g}} \mathfrak{q}- V \mathfrak{q}-i \frac{4a\cos \theta }{\rho ^2} \nabla _{T} \mathfrak{q}=L_\mathfrak{q}[A]+\operatorname {Err}[\square _{\mathbf{g}}\mathfrak{q}], {\tag{4.1}} \end{equation}$$

where $T$ is a vector field defined on the perturbed spacetime which reduces to $\partial _t$ in the Kerr case, $V$ is a real positive potential, and $L_\mathfrak{q}[A]$ is linear in $A$ and has important specific properties. Finally, the error terms $\operatorname {Err}[\square _{\mathbf{g}}\mathfrak{q}]$ depending on all linearized Ricci and curvature coefficients are acceptable error terms, i.e., they verify important structural properties, reminiscent of the null condition.

Due to the presence of the linear terms in $A$ on the right-hand side of 4.1, one has to view the wave equations in 4.1 as coupled with the defining equations for $\mathfrak{q}$, that is coupled with second-order transport type equations in $A$. The most demanding part in the analysis of gRW equations 4.1 is to derive global Morawetz-energy type estimates for $(\mathfrak{q},A)$.

Even for solutions to the scalar wave equation in Kerr, the complicated nature of the trapping region and the presence of the ergoregion creates additional problems in the analysis. To treat these difficulties in linear theory, 46 and 15 rely on methods first used in the context of the scalar wave equation in Kerr, based on mode decompositions and construction of vector fields adapted to different modes. It is however not clear how to extend this method, without loss of derivatives, to general perturbations of Kerr. In our work 31, we rely instead on a physical space method introduced by Blue and Andersson in 3 in the context of the scalar wave equation in Kerr for small angular momentum.

In 3 Andersson and Blue develop a generalized vector field method which allows for commutations with second-order differential operators, and then apply it to the Carter differential operator $\mathcal{K}$ and its elliptic counterpart, together with the Killing vector fields of the Kerr metric, to derive energy and Morawetz estimates for the solution for $|a| \ll M$. The crucial new idea in 3 is to supplement the two Killing vector fields of the spacetime, i.e., $\partial _t$ and $\partial _\phi$ with a second-order operator

$$\begin{equation*} \mathcal{K}=D_\alpha ( K^{\alpha \beta } D_\beta ), \end{equation*}$$

where $K^{\alpha \beta }$ is a Killing tensor, which commutes with the scalar wave operator. The remarkable thing about the Kerr spacetime, discovered by Carter, is that is has an additional Killing tensor that is not given by a composition of the Killing vector fields.

Relative to the principal null frame, the Carter tensor takes the form

$$\begin{equation*} K=-a^2\cos ^2\theta g + O, \end{equation*}$$

where $O=|q|^2 (e_1 \otimes e_1+ e_2 \otimes e_2)$. The associated second-order operator

$$\begin{equation*} \mathcal{O}=D_\alpha ( O^{\alpha \beta } D_\beta ) \end{equation*}$$

verifies commutation properties with $\square _g$, i.e., $[\mathcal{O}, |q|^2\square _g]=0$. In their work, Andersson and Blue 3 introduce a set of second-order operators coming from the Killing vector fields and the Carter tensor, and they use them as symmetry operators for the wave equation.

In our work, we substantially adjust their method to be applied to the nonlinear gRW equation for symmetric traceless 2-tensors for $|a|\ll M$, by being careful about obtaining acceptable error terms by commuting with the approximate symmetry operators in perturbations of Kerr spacetime.

5. The stability of charged black holes

We now present some recent results about the stability problem for the Einstein equation coupled with matter fields, in particular for charged black hole solutions to the Einstein–Maxwell equations.

5.1. Why study nonvacuum solutions?

We start by giving some motivation on why one should study nonvacuum solutions to the Einstein equation.

In August 2017, LIGO made another breakthrough discovery when it detected the first merger of two neutron stars 2. Neutron stars are the collapsed core of massive supergiant stars, mostly composed by neutrons and supported against their collapse by neutron degeneracy pressure. Just two seconds after the gravitational wave signal was detected by LIGO, a flash of gamma-rays was detected by the FERMI satellite in space, coming from the same tiny corner of the cosmos. This observation started an entirely new field of gravitational astronomy, called multi-messenger astronomy, where the same cosmic event is observed through two different signals: one in gravitational waves (detected by LIGO) and one in electromagnetic rays (detected by the FERMI satellite).

The merger of two neutron stars is expected to collapse into a rotating and charged black hole 48. In fact, even though the magnetic field potentially surrounding a black hole is expected to disperse very fast, in the case of a rotating collapse, the momentum of the magnetic field produces a current of electric charge which is expected to remain after the magnetic field has dispersed.

In addition to that, there have been recent works 9, 10 pointing out that even the merger of two black holes could be charged! In fact, the first detection in 2015 (which is still the strongest one observed by LIGO) is compatible with a charge-to-mass ration as high as 0.3, which makes the presence of the charge in analysis of these scenarios not negligible.

5.2. The equations

We now look at how the interaction of black holes with matter fields change the analysis of the equations governing the perturbations. As described in Section 2.1, the interaction between gravitational and electromagnetic fields in spacetime is governed by the Einstein–Maxwell equation,

$$\begin{equation} \operatorname {Ric}(g)=2 \operatorname {F} \cdot \operatorname {F} - \frac{1}{2} |\operatorname {F}|^2 g. \tag{5.1} \end{equation}$$

More precisely, the left-hand side of the equation, involving the Ricci curvature, encodes information about gravitational radiation, while the right-hand side of the equation, involving the electromagnetic tensor $F$, encodes information about electromagnetic radiation.

Gravitational radiation is transported by the Weyl curvature $W_{\mu \nu \gamma \lambda }$ which is a 4-tensor whose extreme null component $\alpha _{ab}=W(e_a, e_4, e_b, e_4)$ is a 2-tensor on the horizontal structure. Electromagnetic radiation is transported by the electromagnetic tensor $F_{\mu \nu }$, which is a 2-form whose extreme null component ${}^{(F)}\beta _a=F(e_a, e_4)$ is a $1$-tensor on the horizontal structure.

The gravitational and electromagnetic fields, when taken independently, satisfy the Teukolsky equation 57 of spin $s$, for $s=\pm 2,\pm 1$, respectively,

$$\begin{align*} \mathcal{T}^{[2]} ( \alpha ) &= 0, \\ \mathcal{T}^{[1]} ( {}^{(F)}\beta ) &= 0, \end{align*}$$

where the operators $\mathcal{T}$ are schematically given, as above, by

$$\begin{equation*} \mathcal{T}= \square _g+ c_1(r, \theta ) \nabla _{\partial _r} + c_2(r, \theta )\nabla _{\partial _t} + c_3(r, \theta )\nabla _{ \partial _{\phi } }- V(r, \theta ) =0. \end{equation*}$$

On the other hand, in coupled electromagnetic-gravitational perturbations of a black hole, there is interaction between the gravitational and the electromagnetic radiation. Instead of having two independent Teukolsky equations, one would expect to have equations which couple the 2-tensor $\alpha _{ab}$ with the 1-tensor ${}^{(F)}\beta _a$, such as for example⁠⁷

⁷

In fact, the operators $\nabla \widehat{\otimes }$ and $\mathrm{div}$, respectively, transform a 2-tensor into a 1-tensor and vice versa.

✖

$$\begin{align*} \mathcal{T}^{[2]} ( \alpha ) &= Q \cdot \nabla \widehat{\otimes }{}^{(F)}\beta , \\ \mathcal{T}^{[1]} ({}^{(F)}\beta ) &= Q \cdot \mathrm{div} \alpha , \end{align*}$$

where $Q$ is the charge of the black hole and

•

$( \nabla \widehat{\otimes }\xi )_{ab}=\nabla _a \xi _b + \nabla _b \xi _a-\delta _{ab} ( \mathrm{div} \xi )$ for a 1-tensor $\xi$,

•

$( \mathrm{div} \theta )_a=\nabla ^b \theta _{ab}$ for a symmetric traceless 2-tensor $\chi$.

But things are not quite as simple as one may expect. In fact, in a charged black hole with $Q \neq 0$, the 1-tensor ${}^{(F)}\beta$ is not gauge-invariant. Gauge-invariance is related to the instability due to the choice of gauge described above, and can be expressed in terms of null frames in the following way.

A general linear frame transformation of the null frame $\{ e_3, e_4, e_a\}$ is of the form

$$\begin{equation} \begin{aligned} e_4'&=\lambda ( e_4 + \mu _a e_a), \\ e_3'&= \lambda ^{-1}(e_3+\underline{\mu }_a e_a ), \\ e_a'&= e_a + \frac{1}{2} \underline{\mu }_a e_4 +\frac{1}{2} \mu _a e_3, \end{aligned} {\tag{5.2}} \end{equation}$$

with $\lambda =1+O(\epsilon )$, $\mu , \underline{\mu }=O(\epsilon )$, where $\epsilon$ is the size of the perturbation. In particular, neglecting terms of size $\epsilon ^2$, the transformed frame $\{ e_3', e_4', e_a'\}$ defined by 5.2 is also a valid null frame. We say that a tensor $\Psi$ on the horizontal structure is (quadratically) gauge-invariant if it only changes quadratically with a general linear frame transformation 5.2, i.e.,

$$\begin{equation*} \Psi '= \Psi + O(\epsilon ^2), \end{equation*}$$

where $\Psi '$ is the tensor evaluated in the transformed frame $\{ e_3', e_4', e_a'\}$. We are interested in gauge-invariant quantities as they are believed to represent physical quantities, such as gravitational and electromagnetic waves, which should not depend on the chosen coordinates or gauge.

Even though $\alpha$ is still gauge-invariant for charged black holes, it is not as useful in view of the equation it satisfies. We define instead (see 29) new gravitational and electromagnetic radiation quantities $\mathfrak{f}$ and $\mathfrak{b}$, a 2-tensor and a 1-tensor, respectively, which are gauge-invariant and are related to the Weyl curvature $\alpha$.

As a consequence of the Einstein–Maxwell equation, the $1$-tensor $\mathfrak{b}$ satisfies a Teukolsky-type equation coupled with $\mathfrak{f}$, and the $2$-tensor $\mathfrak{f}$ satisfies a Teukolsky-type equation coupled with $\mathfrak{b}$. We therefore have a system of the schematic form 29

$$\begin{align*} \mathcal{T}^{[1]} ( {}^{(F)}\beta )&=Q \cdot c_4(r, \theta ) (\mathrm{div} \mathfrak{f}), \\ \mathcal{T}^{[2]} (\mathfrak{f})&= Q \cdot d_4(r, \theta )( \nabla \widehat{\otimes }\mathfrak{b}). \end{align*}$$

It is tempting to ask ourselves, How do we identify the gravitational and electromagnetic radiation from a physical point of view? The Einstein–Maxwell equation, which governs the interaction between the two radiation types, does not clearly distinguish between them though, as the whole perturbation is governed by the coupled electromagnetic-gravitational perturbations. Historically, we know that gravitational radiation is a spin-2 quantity while electromagnetic radiation is spin-1. In the case of coupled radiation, the spin-2 quantity $\mathfrak{f}$ is also defined in terms of the electromagnetic tensor, and, vice versa, the spin-1 quantity $\mathfrak{b}$ is also defined in terms of the curvature. Both quantities encode part of the perturbations, and they satisfy the master equations describing the coupled perturbation.

We found that there exists a Chandrasekhar transformation in the case of charged black holes too. The Chandrasekhar-transformed quantities of $\mathfrak{b}$ and $\mathfrak{f}$, given schematically by

$$\begin{align*} \mathfrak{p}&= f_1(r, \theta ) \nabla _3 \big ( f_2(r, \theta )\mathfrak{b} \big ), \\ \mathfrak{q}&= g_1(r, \theta ) \nabla _3 \big ( g_2(r, \theta )\mathfrak{f} \big ) \end{align*}$$

for suitably chosen functions $f_1$, $f_2$, $g_1$, $g_2$, satisfy a symmetric system of

•

Regge–Wheeler equations in Reissner–Nordström spacetime,

•

generalized Regge–Wheeler equations in Kerr–Newman spacetime.

Observe that in this case the transformation only involves one derivative in the ingoing null direction $e_3$.

5.3. The case of Reissner–Nordström spacetime

In Reissner–Nordström, we obtain the following symmetric system of Regge–Wheeler equations.

Theorem 5.1 (Giorgi 27).

The gauge-invariant quantities $\mathfrak{p}$ and $\mathfrak{q}$ representing electromagnetic and gravitational radiations for linear perturbations of Reissner–Nordström satisfy the following coupled system of wave equations:

$$\begin{align*} \square _g\mathfrak{p}-V_1(r) \ \mathfrak{p}&=\frac{Q}{r^2} \ \mathrm{div} \ \mathfrak{q},\\ \square _g \mathfrak{q}-V_2(r) \ \mathfrak{q}& =\frac{Q}{r^2} \ \nabla \hat{\otimes } \mathfrak{p}, \end{align*}$$

where $V_1(r)= \frac{1}{r^2} \left(1-\frac{2M}{r}+\frac{6Q^2}{r^2} \right)$, $V_2(r)= \frac{4}{r^2} \left(1-\frac{2M}{r}+\frac{3Q^2}{2r^2} \right)$ are positive potentials.

The system is symmetric, as the respective right-hand sides of the equations are adjoint operators on the sphere, i.e., for a $1$-tensor $\xi$ and a $2$-tensor $\theta$,

$$\begin{equation*} \int _{S} \xi \cdot (\mathrm{div} \ \theta ) = \int _{S} (\nabla \widehat{\otimes }\xi ) \cdot \theta . \end{equation*}$$

The symmetric structure on the right-hand side of the system implies that it is possible to derive energy estimates for the system by summing the estimates for the two equations. Upon integration on the sphere, the coupling terms cancel out reducing to boundary terms.

In physical-space terms, this is equivalent to the fact that we can define a combined energy-momentum tensor for the system $\mathcal{Q}_{\mu \nu }[\mathfrak{p}, \mathfrak{q}]$ as

$$\begin{equation*} \mathcal{Q}_{\mu \nu }[\mathfrak{p}, \mathfrak{q}]\coloneq \mathcal{Q}_{\mu \nu }[\mathfrak{p}]+\mathcal{Q}_{\mu \nu }[\mathfrak{q}]-\frac{Q}{r^2} \left(\nabla \widehat{\otimes }\mathfrak{p}\cdot \mathfrak{q}\right) g_{\mu \nu }, \end{equation*}$$

where

$$\begin{equation*} \mathcal{Q}_{\mu \nu }[\psi ]=\nabla _\mu \psi \cdot \nabla _\nu \psi -\frac{1}{2} g_{\mu \nu } \left(\nabla _\lambda \psi \cdot \nabla ^\lambda \psi + V|\psi |^2 \right) \end{equation*}$$

is the energy-momentum tensor associated to the Regge–Wheeler equation $\square _g \psi - V \psi =0$. By applying the vector field method to the current $\mathcal{Q}_{\mu \nu }[\mathfrak{p}, \mathfrak{q}] X^\nu$ for a vector field $X$, we obtain 27

•

energy estimates, as above, for $X=\partial _t$,

•

Morawetz estimates for $X=f(r) \partial _r$, with $f(r)$ vanishing at the photon sphere.

In this case, one obtains in the bulk a quadratic form which can be proved to be positive definite in the exterior region for $|Q|<M$. This is done by separating the exterior region in subregions where the function $f(r)$ is defined differently.

5.4. The mode stability of Kerr–Newman

Recall that one of the clues to the validity of the Stability conjecture for the Kerr–Newman family is the proof of mode stability of Schwarzschild, Reissner–Nordström, and Kerr, as obtained by the physics community.

Quite strikingly, the Kerr–Newman solution stands up as genuinely different from the similar cases of Kerr or Reissner–Nordström in the simplest possible form of stability. As stated by Chandrasekhar in 11, Section 111, “the methods that have proved to be so successful in treating the gravitational perturbations of the Kerr spacetime do not seem to be applicable (nor susceptible to easy generalizations) for treating the coupled electromagnetic-gravitational perturbations of the Kerr–Newman spacetime.”

The techniques applied in those early works, which relied on decomposition in frequency modes of perturbations of the solutions, failed to be extended to the case of Kerr–Newman spacetime. The issue in the analysis of a coupled system comes from the decomposition in modes. In mode stability analysis, one study solutions of the Teukolsky equation $\mathcal{T}^{[s]}(\psi ^{[s]})=0$ of the separated form

$$\begin{equation*} \psi ^{[s]}(t, r, \theta , \phi ) =e^{-i \omega t} e^{i m \phi } R^{[s]}(r)S^{[s]}_{m\ell } ( a \omega , \cos \theta ). \end{equation*}$$

From the Teukolsky equation, one derives an angular ODE for $S$ which defines the spin $s$-weighted spheroidal harmonics $S^{[s]}_{m\ell } ( a \omega , \cos \theta )$, as eigenfunctions of the spin $s$-weighted Laplacian

$$\begin{equation*} \triangle ^{[s]}= \frac{1}{\sin \theta } \partial _\theta (\sin \theta \partial _\theta ) - \frac{m^2+2ms \cos \theta +s^2}{\sin ^2\theta }+a^2\omega ^2\cos ^2\theta -2a\omega s \cos \theta . \end{equation*}$$

For $a=0$, they reduce to the spherical harmonics $S^{[s]}_{m\ell } ( 0, \cos \theta ) = Y^{[s]}_{m\ell } (\cos \theta )$. Spin-weighted spherical harmonics of different spins are simply related through the angular operators appearing on the right-hand side of the coupled equations, and have the same eigenvalues. Schematically,

$$\begin{equation*} \nabla \widehat{\otimes }Y^{[+1]}_{m\ell } =-\lambda Y^{[+2]}_{m\ell }, \qquad \mathrm{div} \ Y^{[+2]}_{m\ell } =\lambda Y^{[+1]}_{m\ell }. \end{equation*}$$

On the other hand, in the general axisymmetric case, the spin-weighted spheroidal harmonics of different spins are not simply related through those angular operators.

This in fact explains the “apparent indissolubility of the coupling between the spin-1 and spin-2 fields” 11 for electromagnetic-gravitational perturbations of Kerr–Newman, in contrast with Reissner–Nordström or Kerr. In treating the coupled electromagnetic-gravitational perturbations of Kerr–Newman spacetime, the decomposition in modes of the equations, which had the objective of simplifying the analysis of the perturbations, actually makes them unsolvable.

5.5. The physical-space analysis in Kerr–Newman

Our approach to solving this issue is to abandon the decomposition in modes and perform a physical space analysis of the equations. Following the road map that mathematicians have taken in the last decade in interpreting in physical space the mode analysis done by the physics community, the Kerr—Newman solution may be the case where a physical space approach could succeed where the mode analysis in physics failed. Observe that our proof of boundedness of a general solution through a physical space analysis will in particular imply the absence of exponentially growing modes, therefore proving mode stability.

We derived 29 the following equations governing the linear stability of Kerr–Newman spacetime to coupled electromagnetic-gravitational perturbations.

Theorem 5.2 (Giorgi 29).

The gauge-invariant quantities $\mathfrak{p}$ and $\mathfrak{q}$ representing electromagnetic and gravitational radiations for linear perturbations of Kerr and Newman satisfy the following coupled system of wave equations:

$$\begin{align*} \square _g \mathfrak{p}- V_1(r, \theta ) \mathfrak{p}-i \frac{2a\cos \theta }{\rho ^2} \nabla _{\partial _t} \mathfrak{p}&=\frac{Q\overline{\Gamma }^3}{\rho ^5} ( \mathrm{div} \ \mathfrak{q})+ a \nabla ^{\leq 1}( \mathfrak{b}, \mathfrak{f}), \\ \square _g \mathfrak{q}- V_2(r, \theta ) \mathfrak{q}-i \frac{4a\cos \theta }{\rho ^2} \nabla _{\partial _t} \mathfrak{q}&= \frac{Q\Gamma ^3}{\rho ^5}\big ( \nabla \widehat{\otimes }\mathfrak{p}-\frac{3}{2} (\eta -\underline{\eta }) \widehat{\otimes }\mathfrak{p}\big )+ a \nabla ^{\leq 1}( \mathfrak{b}, \mathfrak{f}), \end{align*}$$

where $V_1$, $V_2$ are positive real potentials and $\Gamma =r+i a \cos \theta$, $\overline{\Gamma }=r-i a \cos \theta$.

The system is symmetric, as the respective right-hand sides of the equations are adjoint operators with respect to the spacetime integral, i.e., for a $1$-tensor $\xi$ and a $2$-tensor $\theta$,

$$\begin{equation} \xi \cdot (\mathrm{div} \ \overline{\theta })= ( \nabla \widehat{\otimes }\xi ) \cdot \overline{\theta } +( (\eta +\underline{\eta } ) \widehat{\otimes }\xi )\cdot \overline{\theta } -D_\mu (\xi \cdot \overline{\theta })^\mu . {\tag{5.3}} \end{equation}$$

In the above $\eta$ and $\underline{\eta }$ are Ricci coefficients which satisfy $\nabla \rho ^2=(\eta + \underline{\eta }) \rho ^2$. Because of this property, in deriving energy estimates for the system by summing the estimates of the two equations and integrating by parts in $\nabla$, the interaction between the angular derivative of the function $\frac{\overline{\Gamma }^3}{\rho ^5}$ and the lower-order terms $-\frac{3}{2} (\eta -\underline{\eta }) \widehat{\otimes }\mathfrak{p}$ in the system precisely cancels out 28, by giving rise to spacetime boundary terms, according to relation 5.3.

As mentioned above, most of the results for scalar, electromagnetic, and gravitational perturbations of Kerr or Kerr–Newman spacetimes rely on the separability in modes and the frequency-decomposition of the solution. Even though these methods are very effective (and they are at the present moment the only ones that allow for the analysis in the subextremal range for general solutions), they are nevertheless not well suited for the analysis of coupled electromagnetic-gravitational perturbations of Kerr–Newman spacetime, as separability in modes cannot be obtained in that case. The notable exception among the above-mentioned methods is the physical-space analysis for the wave equation in slowly rotating Kerr by Andersson and Blue 3, which makes crucial use of the Carter tensor in Kerr and the fact that the differential operator associated to the Carter tensor commutes with the D’Alembertian operator in Ricci-flat spacetimes.

Nevertheless, the main obstruction to the application of the Andersson–Blue method 3 to the case of Kerr–Newman spacetime is that, even though the metric admits a Killing tensor, its associated differential operator does not in general commute with the wave equation.

In 28, we prove that even though the Kerr–Newman spacetime is not Ricci-flat, the Carter differential operator $\mathcal{K}$ associated to the Carter tensor still commutes with the D’Alembertian operator of Kerr–Newman. Interestingly enough, the commutation property is not a direct consequence of the Einstein–Maxwell equations, but rather of peculiar properties of the curvature and electromagnetic components in Kerr–Newman.

Such properties then allow us to extend the physical-space analysis of Andersson and Blue 3 to Kerr–Newman spacetime 28. More precisely, in 28 we prove that boundedness of the energy for the system of gRW equations in Kerr–Newman can be obtained, provided that Morawetz estimate bounds are derived. Once those are obtained, boundedness and decay for the solutions to the Teukolsky system can be deduced and, subsequently, used to obtain bounds on the gauge-dependent quantities. These latter estimates are the object of future work.

6. Conclusions

One of the fundamental problems in General Relativity is to understand the final state of evolution of initial data for the Einstein equation. Through gravitational collapse and dispersion of gravitational waves, the geometry to which solutions to the Einstein equation are expected to relax outside the event horizon of a black hole is the one given by the known stationary and axisymmetric explicit solutions: the Kerr and the Kerr–Newman black hole, in the vacuum and electrovacuum cases, respectively.

Rigorous analysis of the wave equation in a black hole background by the mathematics community has allowed us to go beyond the simplest form of stability— mode stability—to have a better picture of the behavior of linear and nonlinear perturbations of black hole solutions. In the case of the Einstein vacuum equation, the Teukolsky equation satisfied by a curvature component cannot be used directly to obtain boundedness of general solutions. The Chandrasekhar transformation allows us to transform it into a well-behaved equation—the (generalized) Regge–Wheeler equation—to which standard techniques for the wave equation can be applied.

In both the cases of nonlinear stability of the Kerr spacetime and in coupled electromagnetic-gravitational perturbations of Kerr–Newman, the physical space analysis of the gRW is crucial to deriving the main estimates.

In the case of nonlinear stability of Kerr spacetime, the Andersson–Blue physical space method 3 can be adapted to the case of the nonlinear gRW equation for $|a| \ll M$, which is a tensorial wave equation coupled with transport equations. The derivation of such an equation and its analysis necessitates the introduction of a new formalism 30 which allows for a nonintegrable horizontal structure.

In the case of the Einstein–Maxwell equation, some nontrivial coupling of gravitational and electromagnetic radiations takes place, and it is not a priori clear what are the relevant quantities that one should consider. It turns out that new definitions of gauge-invariant quantities are needed to obtain a favorable system of generalized Regge–Wheeler equations, where the coupling terms appear in a symmetric fashion, allowing for derivation of spacetime estimates.

This is particularly important for the most general case of perturbations of the Kerr–Newman family, where in order to overcome the issue of “indissolubility of the coupling between spin-1 and spin-2 fields” 11, we necessarily need nonseparated equations, which should be analyzed in physical space.

The equations governing perturbations of Kerr–Newman are a system of two coupled gRW equations, obtained from the Teukolsky equations through the Chandrasekhar transformation, which crucially have

•

real potentials and first-order terms of the form $i \partial _t$ (as in Kerr),

•

adjoint coupling terms (as in Reissner–Nordström).

The decomposition in modes prevented us from seeing this good structure of the equations. On the other hand, a physical space analysis of the above system allows us to derive boundedness of the energy and spacetime energy decay estimates, giving a stronger result than mode stability, and therefore solving the long-standing problem of stability of the most general black hole solution.