Teukolsky formalism for nonlinear Kerr perturbations

We develop a formalism to treat higher order (nonlinear) metric perturbations of the Kerr spacetime in a Teukolsky framework. We first show that solutions to the linearized Einstein equation with nonvanishing stress tensor can be decomposed into a pure gauge part plus a zero mode (infinitesimal perturbation of the mass and spin) plus a perturbation arising from a certain scalar ("Debye-Hertz") potential, plus a so-called"corrector tensor."The scalar potential is a solution to the spin $-2$ Teukolsky equation with a source. This source, as well as the tetrad components of the corrector tensor, are obtained by solving certain decoupled ordinary differential equations involving the stress tensor. As we show, solving these ordinary differential equations reduces simply to integrations in the coordinate $r$ in outgoing Kerr-Newman coordinates, so in this sense, the problem is reduced to the Teukolsky equation with source, which can be treated by a separation of variables ansatz. Since higher order perturbations are subject to a linearized Einstein equation with a stress tensor obtained from the lower order perturbations, our method also applies iteratively to the higher order metric perturbations, and could thus be used to analyze the nonlinear coupling of perturbations in the near-extremal Kerr spacetime, where weakly turbulent behavior has been conjectured to occur. Our method could also be applied to the study of perturbations generated by a pointlike body traveling on a timelike geodesic in Kerr, which is relevant to the extreme mass ratio inspiral problem.

Metrics describing a perturbed Kerr black hole play an important role in gravitational physics, for instance in order to model the ringdown phase after a black hole merger, or an extreme mass ratio inspiral. Their analysis also presents nontrivial mathematical challenges: one may wish to prove the expected decay properties of the perturbations (see e.g., [1][2][3] and references for recent progress directly relevant to this work), or to understand the complicated but highly special geometrical structure of the perturbation equations on the Kerr background (see e.g., the monograph [4]).
In the case of linear perturbations (i.e., solutions to the linearized Einstein equation around Kerr), one has the Teukolsky formalism [5,6], which in effect reduces the problem to the study of a linear scalar wave equation. Furthermore, this equation can be solved by a separation of variables ansatz. Combined, these two features dramatically simplify the analysis of linear perturbations. As originally derived, Teukolsky's equation applies to the extreme components ψ 0 or ψ 4 of the linearized Weyl tensor, cf. section II. However, it was discovered soon afterwards [7,8] that the equation also applies directly to the metric perturbation itself in a sense. In fact, if one makes an ansatz of the form h ab = Re S † ab Φ, where S † ab is a certain second order differential operator [cf. (24)] constructed using geometrical objects of the background Kerr geometry, then if the complex "Hertz potential" Φ is a solution to the source-free (adjoint) Teukolsky equation, h ab is a solution to the source-free linearized Einstein equation. Furthermore, it is generally accepted-and argued more carefully in section III-that this ansatz in a sense covers all solutions with the exception only of pure gauge ones and "zero modes," by which one means linear perturbations towards a Kerr black hole with different mass and/or spin.
While the Weyl tensor components are sufficient to extract the information about gravitational radiation at null infinity, it is also highly desirable to know directly the metric perturbation itself. For instance, in the self-force approach [9,10] to extreme mass ratio inspiral, the first step is to determine the linearized metric perturbation resulting from a stress tensor of a point particle on a geodesic in the Kerr background. If one is interested in higher order perturbations, one has to solve, for the n-th order perturbed metric, the linearized Einstein equation with a source determined by the metric perturbations up to order (n − 1). Both of these problems involve solving the linearized Einstein equation with a nonvanishing source. But in the presence of a source, the Hertz-potential ansatz h ab = Re S † ab Φ is generically no longer consistent! Furthermore, while the metric perturbation can be extracted from the perturbed Weyl scalars ψ 0 or ψ 4 for perturbations satisfying the source free linearized Einstein equation by means of the Teukolsky-Starobinsky relations, this is no longer true for the Einstein equation with source. Thus, it would seem that the Teukolsky formalism is not useful to study higher order perturbations of Kerr (for example), and the mileage gained from the highly special properties of the background seems lost again.
In this paper, we propose a way around this problem, which takes advantage of the very special geometric properties of Kerr (and in fact, algebraically special solutions). Our starting point is the linearized Einstein equation with source 1 , where E is the linearized Einstein operator on the Kerr background (26). In applications, T ab would be e.g., the stress tensor of a point source on a geodesic in Kerr, or the nonlinear terms in the equation for the n-th order gravitational perturbation. Instead of the Hertz-potential ansatz, which is generically inconsistent with a nontrivial source T ab , we will prove that the metric perturbation can be decomposed in the following manner (see section II for the precise definition of the operators S, O, T ) h ab = (L ξ g) ab +ġ ab + x ab + Re(S † Φ) ab , where: • (L ξ g) ab = ∇ a ξ b + ∇ b ξ a is a pure gauge perturbation (irrelevant). • x ab is called a corrector field. As we will show, its nontrivial tetrad components [cf. (55)] are obtained from T ab by solving three decoupled ordinary differential equations [cf. (107)] along a congruence of outgoing null geodesics aligned with a principal null direction.
• Φ is a complex potential which is a solution to the adjoint Teukolsky equation with a source, The source η is determined in terms of T ab by an equation of the form (T † η) ab = T ab − (Ex) ab . As we will show, this equation can be solved by integrating one ordinary The difficult step in our algorithm is thus to solve the adjoint Teukolsky equation with source O † Φ = η, a scalar wave equation, which can be dealt with by a separation of variables ansatz [6]. By contrast, the determination of the source η and the tetrad components of the corrector field x ab involves the comparatively easy task of solving scalar ordinary differential equations depending on T ab . Thus, our method represents a drastic simplification compared to the nonseparated tensor equation (1), taking advantage to the maximum possible extent of the special geometric features of Kerr. Before presenting the details of our algorithm in the following sections, let us mention three potential applications of our method, which were the motivation of our analysis.
A. Point sources Point sources are described by a distributional stress tensor supported on a timelike geodesic γ, Since our corrector field x ab is obtained by integrating linear ordinary differential equations sourced by appropriate components of T ab along a congruence of outward geodesics from the past horizon to future null infinity, the corrector will only be nonzero for those geodesics intersecting γ. Furthermore, if we impose zero boundary conditions on x ab on H − , it follows that x ab is supported in each spatial slice along a "string" attached to the point source going out to infinity. The locus of this string is where the null geodesics intersecting γ pierce the slice (see figure 1). Thus, our analysis shows that the perturbation h ab is defined in terms of a Hertz potential (as well as possibly trivial zero mode and gauge perturbations), plus a corrector x ab supported on a semi-infinite string emanating at the point particle going to infinity. If instead vanishing boundary conditions on x ab are imposed at I + rather than H − , then the string goes to the horizon rather than to infinity.
The occurrence of such "Dirac-type strings" has been observed also by a number of previous authors [11] who have studied the equation (1) for point sources by a somewhat different approach. In their approach, the exterior of the black hole is divided into two domains M + and M − separated by an interface S containing the particle orbit. In each domain, an ansatz h ± ab = Re(S † Φ ± ) ab is made for the metric, and then the existence of a distributional source is encoded into suitable junction conditions relating Φ + to Φ − on the interface S and in particular at the worldline γ in the form of differentiability and continuity requirements.
The perturbations obtained in this manner can be categorized [12] based on the form of the singularity into the "half-string" classes and the "full-string" class. A "no-string" gauge may be formed by joining together the two regular sides of two opposite half-string perturbations along the interface S through the particle. Such a perturbation, introduced by Friedman and collaborators in [13,14], avoids the stringlike singularities, but has instead a gauge discontinuity (and also delta-function distributions [12]) on the interface S .
Our construction scheme by contrast avoids the consideration of an interface, but the metric perturbation involves the corrector piece x ab , which has a delta-function-like distribution along a string emanating from the particle. Furthermore, our potential Φ will be a solution to the adjoint Teukolsky equation with a source (3) that will be supported along the string. Our construction thereby gives a rather transparent conceptual understanding for the occurrence of such Dirac-type strings that have been observed in other approaches. It is also conceivable that our construction, which avoids the consideration of an interface, might provide a viable alternative scheme for numerical integrations.

B. Nonlinear perturbations: Decay and turbulent behavior
By iterating our prescription, we can recursively obtain the n-th order gravitational perturbation up to gauge as h ab + x (n) ab + Re(S † Φ (n) ) ab in terms of an n-order potential subject to a sourced adjoint Teukolsky equation ab are determined by transport equations (ODEs) from the nonlinear terms T (n) ab in the n-th order Einstein tensor (involving h (1) ab , . . . , h (n−1) ab ), and where g (n) ab is an n-th order perturbation to another Kerr black hole. While a thorough discussion is outside the scope of the present paper, we here want to ask whether higher order perturbations could be physically relevant in astrophysical situations.
At the linear level, an important aspect of characterizing black hole perturbations is the determination of a quantized set of complex frequencies associated with certain mode solutions to Teukolsky's equation called quasinormal modes [15]. Quasinormal modes are the natural resonances of black hole spacetimes and carry information about the parameters of binary merger remnants. Rapidly spinning Kerr black holes exhibit a special family of weakly damped quasinormal frequencies, which in a sense "corotate" with the horizon. These modes each have the same real part proportional to the horizon frequency and become arbitrarily long-lived in the extremal limit [16,17]. At linear order, the long-lived modes all get coherently excited, and their collective behavior leads to interesting physical effects.
Near the horizon, the superposition of modes results in a transient growth of tidal fields as measured by infalling observers. Away from the horizon, the mode sum yields a transient power-law tail, slowing the asymptotic decay [17,18]. Furthermore, as they all have the same real part of frequency, these modes can be nonlinearly in a state of resonance and can readily transfer energy amongst each other. As such, the nonlinear coupling between multiple overtones is conjectured to lead to turbulent cascades exhibiting a Kolmogorov-type scaling [19]. In this way, such nonlinear features could in principle lead to a unique signature in the ringdown signals of highly spinning black holes, as measured by future gravitational wave detectors.
Generically, however, an arbitrary perturbation to Kerr will extract angular momentum through the mechanism of superradiance, the wave analog of the Penrose process, see e.g., [20]. For nonlinear mode-coupling effects to become important, it is, therefore, necessary that they develop faster than (a) the time scale for linear decay of the quasinormal modes, and (b) the spin-down time scale of the black hole. We now give a heuristic argument, concluding that as long as the perturbation amplitude and surface gravity are appropriately chosen, both of these conditions can be satisfied.
We first estimate the spin-down time scale of the black hole. 2 Suppose we have a mode perturbation to Kerr of amplitude A, frequency ω, and azimuthal number m. The energy and angular momentum fluxes into the black hole scale as [21] where the overdot denotes differentiation with respect to ingoing Kerr time and Ω H is the angular velocity of the outer horizon. For short, with dimensionless frequencies The fractional change in the extremality parameter (surface gravity) κ = 1 − J 2 /M 4 due to absorption at the horizon is given by In the last line, we have used that Ω H = M/(2J) + O(κ). The long-lived "near-horizon" modes of nearly extremal Kerr have ω − mΩ H = O(κ), so this term does not change the order of the final term. Therefore, for near-horizon modes, we have Thus Aω n κ suffices to make the fractional change in κ small over a typical "fast" time scale M . For near-horizon modes, with ω n ∼ κ, the spin-down time scale goes as τ spin-down ∼ M/A 2 .
The time scale for linear decay of a near-horizon mode, meanwhile, is τ decay ∼ M/κ. Thus, for nonlinear mode-coupling to become important, the perturbation amplitude must be such that the mode-coupling time scale τ mode-coupling τ spin-down , τ decay .
The mode-coupling time scale depends on whether the fundamental interaction involves a 3-or 4-mode coupling [22]. Three-mode interactions are stronger, but their presence depends 2 We thank S. Gralla for a discussion in which we made this estimate. on the dispersion relation of the modes and any associated selection rules, analysis of which we do not attempt here. Four-mode interactions are more generic. Three-mode interactions have an interaction time scale that goes as 1/A, whereas four-mode interactions go as 1/A 2 .

3-mode
The overall strength of the interaction depends on the precise overlap integrals, which would have to be studied using the framework developed in [23]. For now, we assume this goes as 1/M . By comparing these time scales to those for dissipation and spin-down, we obtain constraints on the size of the perturbation compared to the spin of the black hole in order for nonlinear mode-coupling effects to become important (see table I).
To summarize this discussion, we conclude that if the leading interaction involves three modes, then provided 1 A κ, nonlinear coupling could become important for nearextremal black holes. If the leading interaction involves four modes and A 2 κ, then the interaction time scale is the same as the spin-down time scale, and such nonlinear effects are possible, although only marginally so. We will analyze this in detail in a future paper [24].

C. Perturbative quantum gravity
Our method should also bring substantial simplifications to perturbative quantum gravity on the Kerr spacetime. For linear perturbations, the corrector tensor x ab is zero, and our results give a decomposition h ab =ġ ab + Re(S † Φ) ab up to gauge. It turns out that, with respect to the canonical symplectic form W for general relativity [cf. (29)], the zero mode perturbationsġ ab (finitely many degrees of freedom) are symplectically orthogonal to the Hertz-potential perturbations 3 (S † Φ) ab , and thus can be "quantized" independently.
Thus, up to the zero mode perturbations (quantization of a finite-dimensional phase space), the quantization of the gravitational field reduces at linearized order to the quantization of the Hertz-potential field. For these, we have the naturally induced symplectic form , which may be used as the basis for the quantization of the field Φ.
For the higher perturbative orders, we should be able to proceed in a recursive manner just as one does when defining the higher orders of the interacting Klein-Gordon quantum field in φ 4 -theory on a curved spacetime, say, see [26,27]. The only difference should be that the sources of the equations for the n-th order perturbation are determined by a certain transport equation as described in section V, and that in addition to the Hertz-potential quantum field Φ, we need to keep track of the corrector x ab at each perturbation order, also as described in section V. We leave the details to a future investigation.
This paper is organized as follows. In section II we review the essential parts of the GHP formalism for the convenience of the reader. In section III we study the decomposition (2) for solutions to the homogeneous linearized Einstein equations. In section IV, we generalize this to solutions to the linearized Einstein equations with source. In section V we point out the fairly obvious application to higher order (nonlinear) perturbations. Some formulae are relegated to the appendix.

Conventions:
We use the (+ − −−) signature convention for the metric and the conventions of [20] otherwise. We also make extensive use of the Geroch-Held-Penrose (GHP)formalism [28] and the associated standard notation.

II. GHP FORMALISM AND TEUKOLSKY EQUATION
Our construction makes heavy use of the special geometric features of the background Kerr geometry. The key simplifying feature behind most arguments is in fact the algebraically special nature of the background geometry and, of course, the vacuum Einstein equation assumed throughout this article. We recall (see [20] for details) that the Weyl tensor C abcd of a 4-dimensional Lorentz metric generically has four distinct "principal null directions" at each point, where a principal null direction is a null vector k a such that Definition 1. A 4-dimensional spacetime is called algebraically special (or Petrov type II) if two principal null directions coincide at every point, k a 1 = k a 2 = l a , and in this case C abc[d l e] l b l c = 0. An algebraically special spacetime is called Petrov type D if there is another pair of coinciding principal null directions k a 3 = k a 4 = n a .
In an algebraically special spacetime, it is true that l a is tangent to a congruence of shear free null-geodesics by the Goldberg-Sachs theorem [29]. Examples of such spacetimes are the Kerr-and Robinson-Trautman families of metrics [30]. The former ones are actually of type D, while in the latter the twist of l a vanishes. In order to take full advantage of the algebraically special property, it is useful to consider null tetrads (l a , n a , m a ,m b ) such that l a is aligned with the principal null direction (and n a with the other one in type D). Our conventions for the tetrad are: l a n a = 1, m am a = −1 and all other inner product vanish. From the null tetrad, we define as usual the Weyl components, 2 C abcd (l a n b l c n d + l a n b m cmd ), Ψ 3 = −C abcd l a n bmc n d , and the spin coefficients ρ, ρ , σ, σ , τ, τ , κ, κ , β, β , , ; see appendix B. The real part of ρ corresponds to the expansion of the null-geodesic congruence tangent to l a , the imaginary part to the twist, and the real and imaginary parts of σ correspond to the shear of l a . κ measures the failure of the congruence to be geodesic in general, so by the Goldberg-Sachs theorem [29] and its corollaries, The corresponding primed quantities correspond to sending l a → n a and m a →m a . By choosing n a appropriately, one can set further quantities to zero, but there are different options: n1) In any type II spacetime, we can perform a null rotation l a → l a , m a → m a + Al a , n a → Am a +Ām a + AĀl a , under which τ transforms as τ → τ + Aρ. Then, if ρ = 0, we can adjust A so that τ = 0, and this choice fixes the null vector n a up to a rescaling.
Making further use of the Einstein field equation, one then sees that τ = σ = 0 in addition to (13); see [31] for details.
n2) In type D, we may alternatively choose n a to be the second principal null direction, which results in the primed version of the Goldberg-Sachs theorem κ = σ = Ψ 3 = Ψ 4 = 0 in addition to (13). In Kerr, an example of such a tetrad is the Kinnersley frame; see appendix B.
Except for the case of Schwarzschild spacetime, one must in general choose between these two possibilities, and we shall make use of either one of them in the following depending on the purpose.
The GHP formalism [28] allows one to handle in a very efficient manner the equations obtained for the spin coefficients and Weyl components. Furthermore, that formulation has a geometrical basis which renders its equations automatically invariant under the remaining permissible rescalings of the tetrad (l a , n a , m a ,m b ). These rescalings consist of (a) a local boost, sending l a → λλl a , n a → λ −1λ−1 n a , or (b) a local rotation m a → λλ −1 m a , where λ is a nonzero complex number depending on the spacetime point. A GHP quantity η is said to have weights (p, q) if transforms as under a combined local boost and rotation.
More mathematically speaking, the construction can be seen as follows. Given two null directions-rather than vectors-aligned with the given l a , n a at each point, one can defined the bundle of all null frames (l a , n a , m a ,m b ) over M such that l a , n a point along the prescribed directions and such that m a ,m a are positively oriented and span the complexified orthogonal complement of the plane determined by l a , n a . This set is seen to define a principal fibre bundle, P, over M with structure (gauge) group the boosts and rotations, i.e., , isomorphic also to the multiplicative group G = C × of nonzero complex numbers; see [32] for the basic definitions related to the notion of a principal fibre bundle.
The gauge group G acts precisely by local boosts and rotations on P. G = C × also acts by multiplication with λ pλq on the 1-dimensional vector space C, and this gives a representation π p,q . The GHP quantities η (p, q) are precisely the sections of the so-called associated complex line bundle L p,q = P πp,q C. We may also consider "mixed" tensorial/GHP objects of GHP weight (p, q) and tensorial rank (r, s), which are sections in the vector bundles The key idea of the GHP formalism is to view the spin coefficients ρ, ρ , κ, κ , τ, τ , σ, σ as GHP quantities of the appropriate weight as specified in appendix B. Likewise, the tetrad components of an ordinary tensor are GHP quantities of an appropriate weight. For instance, if ξ a is a covector field then ξ l = ξ a l a (1, 1), or ξ n = ξ a n a (−1, −1), or ρ (1, 1), etc. The remaining spin coefficients β, β , , by contrast are used to define a covariant derivative operator (i.e., a connection) called Θ a of the vector bundle L p,q ⊗ T r,s M , where Ordinary tensor fields are identified with sections of L 0,0 ⊗ T r,s M , i.e., mixed GHP-tensor fields of weight (0, 0). The GHP-covariant directional derivatives along the tetrad legs are denoted traditionally by These operators shift the GHP weights by the amounts Þ : (1, 1), Þ : (−1, −1), ð : (1, −1), ð : (−1, 1). Apart from mathematical elegance, this geometric viewpoint is useful because we may encounter null tetrads that behave in a singular way at certain points of M , such as the horizon, infinity, the north pole of a sphere, etc. The components of a mixed GHP quantity in such a singular null frame will therefore also be singular, but such a singularity is obviously absent in the geometrical viewpoint where it is simply ascribed to a bad choice of gauge.
By working throughout with invariantly-defined GHP quantities, the formalism therefore automatically takes care of such artificial singularities.
As noted by Wald [33], the essence of the Teukolsky formalism [5,6] and its extension by [7,8] may be succinctly encoded in an operator relation. Let (M , g ab ) be a type II vacuum solution. To state the operator relation, we introduce the linearized Einstein operator (26), as well as partial differential operators S, T , O defined as follows. S, T act on symmetric rank-two covariant tensor fields and give a GHP scalar of weight (4, 0), O, the wave operator appearing in Teukolsky's master equation [5,6], takes a (4, 0) GHP scalar to another (4, 0) GHP scalar and is defined by The operator relation is and can be exploited as follows. The quantity T h = ψ 0 is by construction equal to the linearization of the Weyl component Ψ 0 ; see (12). Then, if h ab is a solution to the linearized field equation Eh ab = 0, it follows that ψ 0 solves the homogeneous Teukolsky equation Oψ 0 = 0. We can also take the formal adjoint of this operator relation and use that [The formal adjoint 4 P † of a partial differential operator P taking sections in L p,q ⊗ T r,s M to sections in L p ,q ⊗ T r ,s M is a partial differential operator taking sections for some vector field w a constructed locally out of η, η and their derivatives of GHP weight zero.] This results in the relation As a consequence, suppose valued solution to the linearized Einstein equation, Eh ab = 0. By taking the real part and using that E is a real operator, we thereby get a real tensor field satisfying the linearized Einstein equation from the "Hertz potential" Φ. The adjoint Teukolsky operator O † is where the second line is the coordinate expression in the Kinnersley tetrad in the Kerr spacetime, see appendix B. The adjoint operators T † , S † are explicitly and As one sees from this expression, the real perturbation h ab = Re(S † Φ) ab automatically is in a gauge where We follow the tradition in which this gauge is referred to as "ingoing radiation gauge" (IRG), although in a general type II vacuum spacetime this terminology has little physical meaning.

III. DECOMPOSITION OF h ab FOR HOMOGENEOUS EQUATION
Although we are generally interested in the inhomogeneous linearized Einstein equation (1), we will in this section first inspect the homogeneous equation As we have recalled, on a type II vacuum background, if Φ is a (IRG) Hertz potential, i.e., a weight (−4, 0) GHP scalar satisfying the (adjoint) homogeneous Teukolsky equation and we would now like to argue that, in a sense, this is the most general solution in the case of the Kerr background, or more precisely, when (M , g ab ) is the exterior region of the black hole in the Kerr spacetime.
We shall give an argument for (27) following that given by Prabhu and Wald [25] in the case of Schwarzschild, with suitable modifications for Kerr. Their argument is of a conceptual nature as it is embedded in the phase-space (Hamiltonian) formulation of general relativity.
One begins by defining the symplectic current for vacuum general relativity. Let h ab , h ab be two symmetric tensor fields, not necessarily solutions to the linearized Einstein equations.
By the definition of the adjoint, since E † = E, we have where w a is in the present context called the "symplectic current." Evidently, it is conserved if h ab , h ab solve the linearized Einstein equation. The explicit form of w a can be found e.g., in [34]. Following that reference, one defines the symplectic form of general relativity to be where Σ is a Cauchy surface with unit future pointing normal ν a , stretching from the bifurcation surface S = H − ∩ H + to spatial infinity. Since w a is conserved, the definition of W is invariant under deformations of Σ leaving S and spatial infinity fixed. If ξ a is a smooth vector field vanishing at S and approaching as r → ∞ either a time translation, spatial translation, or rotation [in the asymptotically Cartesian coordinate system (t, determined by (t, r, θ, φ)], the symplectic inner products with the corresponding pure gauge define the perturbed ADM mass, linear momentum, respectively angular momentum, see [34].
From a solution to the linearized Einstein equation h ab we may determine its initial data (q ab , p ab ), where q ab is a smooth symmetric tensor field on Σ representing the perturbed intrinsic metric, and where p ab is a smooth symmetric tensor density on Σ representing the conjugate momentum (related to the perturbed extrinsic curvature). In terms of these, the symplectic form reads The symplectic form is automatically finite on the space of square integrable pairs (q ab , p ab ), and we denote the corresponding L 2 -type norm by . .
Following [35] one can consider in this L 2 space of (unconstrained) initial data, the symplectic complement V c of the subspace W c of all gauge perturbations (L ξ g) ab generated by smooth vector fields ξ a that become an asymptotic translation or rotation at infinity and whose projection onto S vanishes. Such initial data satisfy the constraints weakly, and havė Based on results of [35], Prabhu and Wald [25] have shown that any (q ab , p ab ) representing a smooth perturbation with compactly supported intial data on Σ can be approximated in for some h ab ∈ V c implies that such a h ab is pure gauge. Furthermore, [25] have given strong arguments-but not a complete proof-that this must in fact be true in the Schwarzschild background. Thus, in the Schwarzschild spacetime, this gives a strong argument that L 2 initial data (q ab , p ab ) havingṀ =Ṗ =J = 0 and satisfying the constraints weakly can be approximated up to a pure gauge perturbation by the intial data of a smooth perturbation Re(S † Φ) ab arising from a Hertz potential. Since perturbations with a given nonzeroṀ ,Ṗ ,J can be obtained by adding a perturbationġ ab towards another, perhaps boosted/rotated, Kerr black hole this gives an argument-and precise sense-in which (27) should hold true.
We now would like to give a corresponding argument for the sub-extremal Kerr metrics, i.e., for 0 < a ≤ M . [25] have shown that in a general type II background, the relation (32) implies T h = 0, so ψ 0 = 0. The argument is then completed by the following theorem, which shows that h ab must be pure gauge (since it cannot be a perturbation to another Kerr black hole in view ofṀ =Ṗ =J = 0).

Theorem 1. Let h ab be a weak solution to the linearized Einstein equation on a Kerr
background with 0 < a ≤ M , with initial data (q ab , p ab ) ∈ L 2 × L 2 , having ψ 0 = 0. Then also ψ 4 = 0 and h ab is up to gauge equal to an infinitesimal perturbation towards another Kerr black hole.
Before we give the proof of this theorem, we point out that Wald [36] has given a proof for smooth mode solutions. Our proof thus generalizes that of Wald in the sense that our solutions need not be smooth nor mode solutions.
Proof. By the Teukolsky-Starobinsky identities which hold in a general type D background (see e.g., [37], [38]), where ξ a = Ψ Because ψ 0 = 0, these relations evidently give strong constraints on ψ 4 , and we shall now argue that, in fact, we learn ψ 4 = 0. From (33a), we first conclude in view of Þρ = ρ 2 that Here and in the following, a degree mark • on a quantity indicates that this quantity is annihilated by Þ, so x • means that Þx • = 0. We now substitute (35) into (34). To take full advantage of Þα • i = 0, we employ a formalism invented by Held [31,40]. His formalism has three main ingredients. The first is a set of new operators, denoted byð,ð andÞ which replace ð, ð and Þ and contrary to the latter have the property thatð x • ,ð x • and Þ x • are quantities annihilated by Þ. The second is a list of identities, generated from the Einstein equations and Bianchi identities, expressing any background quantity in terms of ρ and quantities annihilated by Þ, together with a table of how the new operators act on these. These two features enable one to write any expression appearing in (34) as a Laurent polynomial in ρ with coefficients that are annihilated by Þ. The third ingredient is a simple lemma showing that if GHP quantities a • i satisfy [Proof: Divide by ρ n , then act with Þ n and use repeatedly Þρ = ρ 2 to obtain a • 0 = 0. Substitute this relation and iterate the process.] We begin by recalling the precise definition of the new operators [40] on weight (p, q) GHP quantities, Using the vacuum Einstein equations in GHP form, [40] next shows that on a type D background of the Case II in Kinnersley's classification [41] (which includes Kerr) with a null tetrad satisfying n2), with Ω • (−1, −1). These identities are complemented by identities for the action of the new operators on ρ and on Ω • , ρ • , τ • , Ψ • 2 given in the table II extracted from [40]. Now we substitute (35) into (34), replacing the GHP operators in terms of Held's operators (37), and expanding all quantities using (38). Each of the resulting terms is given by a product of integer powers of ρ andρ times a quantity annihilated by Þ. We then divide byρ 2 in order to get rid of any positive power ofρ, and for the negative powers ofρ we In deriving these identities, we have used that Ψ • 2 and ρ • are real in Kerr. It is already clear that without hidden linear dependencies between these 7 equations, it will be impossible to get a nonzero solution. In fact, this is very easy to see for smooth h ab . First, we note the values of Ω • , τ • , Ψ • 2 in the Kinnersley tetrad in Kerr: Now (46) gives α • 0 = 0 away from the "equatorial plane" θ = π/2-the only place where Ω • vanishes (here we must use 5 a > 0) and the north-and south poles θ = 0, π-the only places where τ • vanishes. Since α • 0 = 0 is assumed to be smooth, it vanishes everywhere. Next, (45) combined with α • 0 = 0 gives α • 1 = 0 by the same argument and similarly (44) and (43) then subsequently give α • 2 = 0 and α • 3 = 0. So, ψ 4 = 0 in view of (35). Thus, h ab is a perturbation towards another type D metric. These have been analyzed by [36] who has shown that h ab is modulo gauge either a perturbation towards another Kerr metric, or a perturbation towards a rotating C-metric or a NUT metric (in total 4 real linear independent perturbations). However, by inspection the latter two do not have initial data that are in L 2 × L 2 in any gauge, so the proof is complete when h ab is smooth.
When h ab is a distributional solution, the argument as given only shows that α • i , hence ψ 4 , are distributions supported on the equatorial plane and the north-and south poles. Since Þ = l a ∇ a − p − q¯ and = 0 in Kerr in the Kinnersley frame, and since l a is asymptotically tangent to a congruence of outgoing null geodesics, a GHP quantity annihilated by Þ may be viewed geometrically as function on I + (or H − ); see figure 2. If we define outgoing Kerr-Newman coordinates (u, r, θ, ϕ * ) as in (B4), then (u, θ, ϕ) are good coordinates on I + (or H − ) away from the poles. On I + , the Held operatorsð,ð andÞ arẽ in the Kinnersley tetrad. In outgoing Kerr-Newman coordinates, r is an affine parameter along the orbits of l a (see figure 2), so we can write α • i near the equatorial plane in the general form for some k < ∞. In fact, we must have k < 2, for otherwise the metric perturbation h ab -being two derivatives up from the curvature component ψ 4 -would have a δ-function singularity on the equatorial plane, contradicting our assumption that the initial data are in L 2 × L 2 . If we substitute the expression into (40), we see using (48) that the derivative operatorð ð will generate terms containing second and third derivatives of the delta-functions inside α • 3 , which cannot be compensated by the other terms in view of k < 2. Therefore, we see that α • 3 = 0. We substitute this result into (41), and then the same type of argumentation yields α • 2 = 0. Repeating this subsequently with equations (42) and (43) finally shows that also α • 1 = α • 0 = 0. A completely analogous argument can be made for the poles. The proof is complete.
Remark: In the case of Schwarzschild spacetime, Ω • = τ • = 0. In that case, only the four equations (40) -(43) are non-trivial, and become These equations can be conveniently analyzed by going to a specific tetrad and decomposing the quantities α • i into (spin-weighted) spherical harmonics. An analysis leads to the solutions already given in eq. 3.41 of [25], which correspond to the algebraically special perturbations of Schwarzschlid found previously by Couch and Newman [42]. Thus, for Schwarzschild, there are non-trivial solutions, but these fail to give metric perturbations in L 2 as argued by [25]. They are hence excluded by the assumptions of the theorem.

IV. DECOMPOSITION OF h ab FOR INHOMOGENEOUS EQUATION
Now we would like to study solutions h ab to the sourced linearized Einstein equation where E is the linearized Einstein operator (26). The first part of our analysis will be entirely local and works in any type II spacetime and with any, possibly distributional, T ab . The final part of the analysis involves global features of Kerr, and to avoid difficult issues of analytical nature distracting from the main points, we will assume that the the source T ab is smooth and of compact support. Of course, for there to be any solutions h ab , we must have Ideally, we would like to be able to write the perturbation up to gauge as h ab = Re(S † Φ) ab in terms of a Hertz potential in order to take advantage of the special features of type II backgrounds. However, it is easy to see that this will not be possible for a generic T ab .
Indeed, such a h ab is automatically in the ingoing radiation gauge h nl = h ll = h ml = h mm = 0, and the ll component of (51) together with (A1) for (Eh) ll then gives T ll = 0. Thus, if T ll = 0 somewhere-as happens, e.g., for a the stress tensor of a point particle (4)-then the perturbation can definitely not be written as h ab = Re(S † Φ) ab in terms of a Hertz potential (up to gauge).
The main idea of this paper is that we may nevertheless write h ab = Re(S † Φ) ab + x ab up to gauge for a relatively simple "corrector tensor" x ab determined from T ab . Furthermore, the Hertz potential Φ will be a solution to the adjoint Teukolsky equation with a certain source also determined from T ab . To begin, it is easy to see that without the use of any Einstein equation, a metric perturbation h ab may be put in a gauge where h ab l b = 0. Indeed, if h ab is not in this gauge to begin with, we can solve for a gauge vector field ξ a satisfying (h ab − L ξ g ab )l a = 0.
Since Þ = l a ∇ a − p − q¯ , the first equation is an ordinary differential equation (ODE) for ξ l along the orbits of l a . Substituting a solution ξ l into the second equation, this similarly gives an ODE for ξ m along the orbits of l a . Finally substituting both solutions ξ l , ξ m into the third equation gives an ODE for ξ n along the orbits of l a . Having determined the null-tetrad components (ξ l , ξ n , ξ m , ξm) of ξ a we can redefine h ab → h ab − L ξ g ab in order to obtain a perturbation such that h ab l b = 0.
The corrector field x ab is now determined in such a way that with the idea to eliminate any l component from T ab . Our ansatz for x ab is where x nm =x nm , so there are 4 real independent components encoded in x mm , x mn , x nn by which we attempt to satisfy the 4 real independent equations (54). We now first transvect (54) with l a and use the ll component of the Einstein operator E, (A1), to obtain Next, we transvect (54) with m a and use the ml component of the Einstein operator E, (A4), to obtain: Finally, we transvect (54) with n a and use the nl component of the Einstein operator E, (A3), to obtain: To simplify, we could set τ = τ = σ = 0 by an appropriate choice of null tetrad; see n1). ODE for x nn along the orbits of l a . We can thereby obtain x ab by solving three scalar ODEs, a comparatively manageable task in practice.
Having determined x ab , we now adjust h ab → h ab − x ab . The adjusted h ab still has h ab l b = 0 and additionally satisfies What x ab has thereby achieved is to remove the l components from the new source S ab , i.e., S ab l b = 0. At this stage, we can apply a result by [43] showing that there exists a gauge vector field η a such that after a gauge transformation h ab → h ab − L η g ab , the perturbation h ab satisfying (59) is in ingoing radiation gauge, h nl = h ll = h ml = h mm = 0.
To summarize, we have obtained from the original perturbation satisfying (51) by a shift h ab → h ab − L ξ g ab − L η g ab − x ab a new perturbation in ingoing radiation gauge satisfying (59). Furthermore, the new source S ab has no l components. Since we have described a constructive procedure to obtain x ab , S ab from the original source T ab , we can from now work with the shifted h ab satisfying (59), remembering to add x ab back at the end.
At this stage, we would like to show that the shifted h ab satisfying (59) It is natural to use the first equation (60) in order to define Φ. Since (60) can be viewed as a second order ODE along the geodesic tangent to l a , integration leaves us with two constants per geodesic. A change of those integration constants entails the change where A • (−4, 0), B • (−1, 3), and, as throughout this article, a degree • marks a quantity annihilated by Þ. We shall have to use this freedom momentarily when trying to satisfy the second equation, (61). In order to understand this relation, we consider the lm component of the Einstein equation (59), (Eh) lm = 0; see (A4). Into this equation, we substitute (60), and we use the following GHP operator identity: which follows from the GHP commutators, Bianchi identity, and background Einstein equation in any type II spacetime 6 . Then, if y is defined as the left hand side of (61) minus the right side, we get This is again a second order ODE along the geodesics tangent to l a , which integrates using where a • (−3, −1), b • (0, 2) are undetermined GHP scalars annihilated by Þ. In order to compensate these, we now use (63), under which y changes as On the right side, we now substitute for the GHP operators, ð, ð and Þ Held's operatorsð,ð , andÞ (37), and we use the table III of operations in type II spacetimes and the following definitions/results (68) taken from [31], taking advantage also of the special tetrad choice described in n1): Then it is found after a calculation that the change (67) will set y as in (66) to zero, provided A • , B • can be chosen so that We shall postpone the discussion of these equations until after theorem 2 and assume that they can be satisfied. Then both (61) which follows from the GHP commutators, Bianchi identity, and background Einstein equation in any type II spacetime 7 . Then, if z is defined as the left hand side of (62) minus the right side, we get: This is a first order ODE along the geodesics tangent to l a , which integrates using Þρ = ρ 2 to: where c • (1, 1) is a real, undetermined GHP scalar annihilated by Þ. So the remaining task is to eliminate this c • by whatever freedom we have left. As noted by [43] there exist gauge vector fields ζ a such that after a gauge transformation h ab → h ab − L ζ g ab , the ingoing radiation gauge perturbation h ab satisfying (59) is still in ingoing radiation gauge, h nl = h ll = h ml = h mm = 0. These so called "residual gauge vector fields" can be characterized as follows [43]. Denoting by ζ l , ζ n , ζ m , ζm the covariant tetrad components of ζ a in a tetrad satisfying n1), we have where ζ • l , ζ • n , ζ • m are GHP quantities annihilated by Þ which must satisfỹ These conditions leave ζ • m undetermined, and we shall now use this freedom to remove c • . First, under a residual gauge transformation h ab → h ab − L ζ g ab , the Hertz potential Φ as defined above changes to whereas h nn changes to (76) 7 It is the ln component of (21). Now, z is defined as the left hand side of (62) minus the right side, related to c • by (72).
From this, it follows after a lengthy calculation using formulas from table III that, provided that we can chose ζ • l , ζ • n , ζ • m so that then performing the gauge transformation h ab → h ab − L ζ g ab with ζ • l , ζ • n , ζ • m subject to (74), will in effect cancel c • . Therefore, (62) will be satisfied. If we combine the gauge vector fields defined so far into X a = η a + ξ a + ζ a , then we can summarize our findings to this point as follows.

Theorem 2. Suppose h ab satisfies the inhomogeneous linearized Einstein equation (Eh) ab =
T ab on a Petrov type II background, and suppose the symmetric tensor field x ab is defined as (55), with tetrad components obtained by successively solving the three ODEs (56), (57), (58). Then there exists a gauge vector field X a and a GHP scalar Φ (−4, 0) such that provided that we can solve for given a • , b • , c • the equations (69), (77) subject to (74).
Remark: We note that the transport operators appearing on the left side of equations (56), (57), (58) may be rewritten in any type II background as . of the Kerr black hole with 0 ≤ a ≤ M . We also assume that T ab has compact support and (for simplicity) that it is smooth. Furthermore, we assume that h ab is the retarded solution. Then the support of h ab satisfies supp(h ab ) ⊂ J + (supp(T ab )) in some gauge (e.g., the Geroch-Xanthopoulos gauge [44]), where J + is the causal future of a set, see figure 3.
These operators are equal to the "spin-weighted" invariant operators on S 2 described from a geometric point of view in detail in [45]. More precisely, a cross section S 2 of I + gives rise to bundles E s → S 2 , where E s is the line bundle isomorphic to an appropriate restriction of L p,q to this cross section with "spin" s = 1 2 (p − q), and the operators are invariantly defined elliptic operatorsð,ð mapping sections of E s to sections of E s+1 , E s−1 , respectively. From the relation betweenð,ð and the spin-weighted spherical harmonics s Y ,µ [46], the eigenvalues of this operator are explicitly found to be 1 2 L(L + 1) where L = ( + 1), = 0, 1, 2, . . . , and the eigenfunctions are the 0 Y ,µ -harmonics. Thus, the "heat equation" can be solved by decomposing into such harmonics. Once we have ζ • l , the remaining components ζ • m , ζ • n are obtainable from (74). is orthogonal to the s Y ,µ spin-weighted spherical harmonics with s = −1, l = 1. Now, the components of the perturbed Bondi angular momentumJ of h ab is obtained by taking r times the integral of −1 Y 1,µ h mn over a sphere S 2 at constant u and letting r → ∞, see e.g., [47], and we know thatJ = 0. From this, it can be seen going through the definitions that a • is orthogonal to the −1 Y 1,µ , µ = −1, 0, 1 spin-weighted spherical harmonics. On the other hand, one can see that if we choose Φ as a solution to (60) with zero initial conditions on H − , say, then b • must be zero, because such a term in y as in (66) is inconsistent with the behavior [44] of h ab at I + in the Geroch-Xanthopoulos gauge. Thus, it is consistent to take In Minkowski, M = 0 (77) simplifies to: The operator on the right side is invertible on the subspace orthogonal to the zero mode = 0. The perturbed Bondi-mass of h ab is obtained by taking r times the integral of h nn over a sphere S 2 at constant u and letting r → ∞, see e.g., [47]. Going through the definitions, it follows that if c • is not orthogonal to the constant function, then the perturbation h ab would have a nonzero perturbed ADM (and Bondi-) massṀ , which cannot be. The equation (69) is solved as in Schwarzschild.
In Kerr, one can also solve the equations, but the argument is more involved. Thus, we shall now present an alternative argument which in effect bypasses equations (69) and (77).
In the argument leading to the preceeding theorem, we go back to the definitions of the gauge vector field ξ a and the tensor x ab . We integrate the ODEs required to solve for these quantities, (53), (56), (57), (58) outwards from H − to I + along the geodesics tangent to which integrates to [43] h mm =h • ρ ρ for undetermined GHP scalars h • , k • annihilated by Þ. However, h mm = 0 near H − . By Held's technique we can write (83) as a polynomial equation in ρ as in (36) with coefficients a • i annihilated by Þ. Then these coefficients vanish near H − , and therefore all coefficients a • i = 0 globally since these quantities are annihilated by Þ, which is a directional derivative along the orbits of l a going from H − to I + ; see figure 2. One thereby easily finds h • = k • = 0, i.e., h mm = 0 globally and the metric is automatically in ingoing radiation gauge. The next step in the argument was to integrate the ODE (60) for Φ along the geodesics tangent to l a going from H − to I + . Again, we do this by choosing vanishing initial conditions for Φ.
Then y as in (66) vanishes in a neighborhood of H − . By (39), we can express y asρ 2 ρ −3 Remarks: 1) As stated, our result is not symmetric under time-reflection because we have been working with a tetrad satisfying n1) where l a is a principal null direction while n a is not (except in Schwarzschild). However, one can derive the same result also for a tetrad satisfying n2) in Kerr, where l a and n a are both principal null directions (e.g., the Kinnersley tetrad, see appendix B). The derivation proceeds along the same lines, but some formulas, especially those related to the residual gauge vector field ζ a , given in [43], are more complicated in this case.
2) Solutions to the inhomogeneous linearized Einstein equation (Eh) ab = T ab with other initial conditions differ from the retarded solution by a solution to the homogeneous equation.
We have demonstrated in the previous section that a corresponding decomposition holds to arbitrary precision in a suitable sense, where there can now be an additional pieceġ ab representing an infinitesimal perturbation towards another Kerr solution.
3) The derivation of our result goes through also for distributional T ab of compact support.
In that case, X a , Φ are distributional as well.
We now turn to the equation satisfied by Φ. To this end, we act with the linearized Einstein operator E (26) on (84), and we use the operator identity (21). This results in Thus, if Φ solves the sourced adjoint Teukolsky equation with η such that then h ab as in (84) is a solution to the sourced linearized Einstein equation (Eh) ab = T ab , and as discussed, every solution arises in this way. We view (87) as the defining equation for η and we know that this equation must have solutions. In components, (87) is by (23) The left side is in GHP form, using the linearized Einstein operator in GHP form given in section A.
x nn , x mn , x mm , which are obtained by successively solving in this order the three ODEs  (36) with coefficients a • i annihliated by Þ. Then these polynomials vanish near I + , and therefore all coefficients a • i = 0 globally since these quantities are annihilated by Þ, which is a directional derivative along the orbits of l a going from H − to I + . As a consequence, we find writing out the explicit expressions for a • i in terms of C • , D • , Ω • and their derivatives byð,ð that (89), (90) are satisfied if and only ifðC • =ðD • = 0, so we may pick C • = D • = 0. Therefore, the source η is given by integrating (88) outwards from H − to I + with vanishing initial conditions.
We arrive at the main result of this paper: with radial function −2 R ω m (r) determined by the radial spin −2 Teukolsky equation [5,6] with a source −2 η ω m (r) obtained from η by a similar decomposition, and with −2 S ω m (θ) the spin-weighted spheroidal harmonics. 8 For retarded initial conditions as in the statement of the theorem, the radial function satisfies the conditions −2 R ω m ∼ r 3 e iωr * , r * → ∞, where ∆ has the conventional meaning for the Kerr metric (see appendix B) and Here Ω H is the angular frequency of the outer horizon Ω H = a/(2M r + ).
since T ll = 0 on H + and ϑ l = σ ab l = 0 on H + in the background,θ l vanishes on the entire future horizon H + . As explained in detail in [35], it follows that the perturbation is automatically in a gauge in which the position of H ± remain in their fixed (background) coordinate location to first order. 6) All our equations and constructions remain valid if we work with n a , the second repeated principal null direction, instead of l a . The integrations necessary to solve the analogs of the ODEs (56), (57), (58), (88) are now along n a rather than l a , going from I − to H + , rather than from H − to I + (see figure 2). Thus, we should prescribe trivial initial conditions at I − for these ODEs rather than at H − . The metric obtained in this way will be in the so-called "outgoing radiation gauge." Let us finally clarify the relationship between the equation O † Φ = η in theorem 4 and the inhomogeneous Teukolsky equation for the perturbed Weyl scalar ψ 0 already given in Teukolsky's original papers [5,6]. First, we recall how this is derived: applying the operator relation (19) to a symmetric tensor h ab such that Eh ab = T ab gives the relation Oψ 0 = T 0 , where ψ 0 = T (h) is the perturbed 0-Weyl-scalar (12) and where T 0 = S(T ) [see (17a)] is the source in Teukolsky's equation for ψ 0 [5,6].
Now, we substitute for h ab the decomposition h ab = x ab + Re(S † Φ) ab + L X g ab +ġ ab of theorem 4. 9 Since the perturbed Weyl scalar ψ 0 is gauge invariant on any background in 9 Note that in theorem 4, we are assuming retarded initial conditions. For other initial conditions, there which the corresponding background scalar Ψ 0 = 0, it follows that T (L X g) = 0. Furthermore, T (x) = 0, since the corrector field x ab has vanishing lm, ll, mm components [see (17b),(55)].
Theorem 5. Let ψ 0 be the perturbed Weyl scalar associated with a metric perturbation h ab = x ab + Re(S † Φ) ab + L X g ab as in theorem 4, with Φ satisfying O † Φ = η with source η as described in theorem 4. Then where Remark: According to this theorem, alternatively, having obtained ψ 0 by solving the it remains necessary to compute the corrector tensor x ab .

V. NONLINEAR PERTURBATIONS OF KERR
So far we have limited the discussion to linearized perturbations. However, even in the absence of matter sources, the linearized metric perturbation acts as source for the second order perturbation and so on. To count orders let us write would be another pieceġ ab representing a perturbation towards another Kerr black hole, see remark 2) following thm. 3, which we can include here without problems.
ab is the background Kerr metric. The equation for h (1) ab is of course just the source free linearized Einstein equation in the absence of matter sources, whereas the higher order corrections for n > 1 satisfy an equation of the form In the perturbative setting, gauge transformations are formal diffeomorphisms f = Exp(ξ a ) generated by a formal vector field and correspond to g ab → f * g ab = exp(L ξ )g ab , giving, for instance etc. As described in detail in section III, at zeroth order, we can apply the decomposition ab + (L ξ (1) g (0) ) ab + Re(S † Φ (1) ) ab , where Φ (1) is a Hertz potential, i.e., a (−4, 0) solution to O † Φ (1) = 0, and where in this section, is an n-th order perturbation to another Kerr metric. Now, if we could inductively apply at each order n > 1 the decomposition for h (n) ab described in theorem 4 and the following remarks (see section IV), with an n-th order gauge vector field ξ (n)a and x (n) ab , η (n) determined from T (n) ab , then we could write the metric g ab as the formal series with O † Φ (n) = η (n) , up to a pull back by the formal diffeomorphism generated by ξ a = ∞ n=1 ξ (n)a . Since the correctors x (n) ab and the GHP-sources η (n) are comparatively easy to obtain by solving ordinary transport equations along the orbits of l a , the analysis of nonlinear perturbations has thereby effectively been reduced to solving a Teukolsky equation with source at each order. In view of the powerful methods/results already available in this setting-such as separation of variables method or the recent decay results by [1][2][3]-this is obviously a substantial simplification of matters.
Unfortunately, as stated, theorem 4 only holds for stress tensors of compact support, while the n-th order stress tensors of the gravitational field will clearly not have this property. In order to make our arguments rigorous, we would have to make some approximation of T (n) ab by stress tensors of compact support and then take a limit. In order to control the limit, we must understand the asymptotic behavior of the perturbations h (in terms of certain energy norms on the initial data) has recently been obtained by [1][2][3].
We expect that their results, suitably generalized to the Teukolsky equation with source, can be used to show rigorously the decomposition (106) to all orders, thereby also establishing the asymptotic behavior of higher order gravitational perturbations on Kerr. However, such an analysis goes beyond the present work and is therefore postponed to another paper.

VI. SUMMARY AND OUTLOOK
For the convenience of the reader, we summarize our integration scheme for the sourced linearized Einstein equation (Eh) ab = T ab , which consists of the following steps: Step 1: For given T ab , we integrate, in this order, the ODEs (56), (57), (58) to obtain x mm , x nm , x nn . In the frame (B5) in outgoing Kerr-Newman coordinates (B4) (u, r, θ, ϕ * ), where Þ = ∂/∂r, they have the form [using equations (79)]: x mn = r.h.s. of (57) with right hand sides involving T ll , T ml , T nl , respectively and ρ = −(r − ia cos θ) −1 . For x mm , x nm , x nn , we choose trivial initial conditions at the (past) event horizon r = r + .
Step 2: We integrate the ODE (88) in order to obtain η. In the frame (B5) in outgoing Kerr-Newman coordinates (B4), this equation has the form: with right hand side involving T mm , x mm , x nm , respectively. For η, we choose trivial initial conditions at r = r + .
Step 3: Solve the adjoint Teukolsky equation O † Φ = η with the desired initial conditions.
One may use to this end a separation of variables ansatz with radial function −2 R ω m (r) determined by the radial spin −2 Teukolsky equation [5,6], and with a source −2 η ω m (r) obtained from η by a similar decomposition. For instance, for a retarded perturbation h ab , Φ should be a retarded solution to the adjoint Teukolsky equation, which is equivalent to the boundary conditions (95).
Step 4: The desired solution of the sourced linearized Einstein equation (Eh) ab = T ab is obtained as h ab = x ab + Re(S † Φ) ab , where x ab is constructed from x mm , x nm , x nn as in (55), and where (S † Φ) ab is as in (24).
By iterating this procedure, one can, for example, construct recursively the n-th order gravitational perturbation h (n) ab in terms of an n-th order potential Φ (n) solving a sourced Teukolsky equation and an n-th order corrector tensor x (n) ab . These quantities are constructed from the nonlinear terms in the n-th order Einstein tensor, which act as an effective stress tensor T (n) ab in the n-th iteration step. A different path for constructing higher order (in fact second order) perturbations in a formalism based on the Teukolsky method was proposed some time ago by Campanelli and Lousto [49]. Their starting point is the wave equation for the Weyl tensor, which has the schematic form ∇ e ∇ e C abcd = terms quadratic in C abcd (110) in any Ricci flat (R ab = 0) spacetime. By transvecting this tensor equation in all possible ways with the legs (l a , n a , m a ,m a ) of a Newman-Penrose tetrad, one obtains a system of equations involving the Weyl scalars (12) and the rotation coefficients (see appendix B).
By developing the metric g ab around Kerr, one obtains a Teukolksy type equation for the second order perturbed Weyl scalars sourced by terms involving the first order perturbed quantities. To get a closed system of equations, one must solve for the first order perturbed rotation coefficients, i.e., effectively the first order perturbed metric, in terms of the first order perturbed (extreme) Weyl scalars. This can be done, e.g., by well-known inversion formulas such as obtained by [50] at first order. However, the method breaks down as formulated at higher orders, since the inversion formulae no longer apply. It seems that our method, which works directly with the metric perturbation rather than Weyl scalars, is superior in this sense.
In the future, we would like to use the formalism developed in this paper in combination with the mode projection method outlined in [23] to study the interaction of quasinormal modes in the near horizon region of a near-extremal black hole [24]. It would also be an interesting project in our view to use our formalism in combination with results obtained by [1][2][3] in order to obtain decay properties of higher order gravitational perturbations in Kerr.
What seems to be required here is primarily a detailed analysis of the decay properties for solutions to the Teukolsky equation with a source (step 3), since the transport equations (steps 1,2) would then be fairly easy to treat. The method developed in this paper can also be applied to obtain the gravitational field of a point particle in the Kerr metric to linear order, and it would be interesting to compare this to the current alternative treatments of this problem. Finally, it might be possible to use our decompositions in order to obtain simplifications for perturbative quantum gravity off a Kerr background.
M. Wald for discussions.

Appendix A: Linearized Einstein operator in GHP form
Here we quote from [43] the tetrad components of the linearized Einstein operator (Eh) ab (26) in GHP form. 10 It is assumed that the first null leg l a of the null tetrad (l a , n a , m a ,m a ) is aligned with a principal null direction of a type II spacetime (so that κ = σ = Ψ 0 = Ψ 1 = 0 in view of the Goldberg-Sachs theorem). In a type D spacetime we have additionally the option n2) described in section II to set further GHP scalars to zero, and in type II we can alternatively make the simplifications described in n1) by an appropriate choice of n a .  10 We thank B. Wardell for pointing out several typos in the expressions in [43].
We choose the normalization n a l a = 1 and m am a = −1, corresponding to the "−2" signature.
The Kerr metric may be defined by the Kinnersley tetrad, which is convenient in explicit computations. The covariant Boyer-Lindquist coordinate components in the order (t, r, θ, φ)