Conserved currents for Kerr and orthogonality of quasinormal modes

We introduce a bilinear form for Weyl scalar perturbations of Kerr. The form is symmetric and conserved, and we show that, when combined with a suitable renormalization prescription involving complex r integration contours, quasinormal modes are orthogonal in the bilinear form for different (l, m, n). These properties are apparently not evident consequences of standard properties for the radial and angular solutions to the decoupled Teukolsky relations and rely on the Petrov type D character of Kerr and its t-$\phi$ reflection isometry. We show that quasinormal mode excitation coefficients are given precisely by the projection with respect to our bilinear form. These properties can make our bilinear form useful to set up a framework for nonlinear quasinormal mode coupling in Kerr. We also provide a general discussion on conserved local currents and their associated local symmetry operators for metric and Weyl perturbations, identifying a collection containing an increasing number of derivatives.


I. INTRODUCTION
Quasinormal ringing is the principal gravitational-wave signature of the final black hole after a binary merger.This is described by a spectrum of complex quasinormal frequencies ω lmn , which are uniquely specified in linear perturbation theory by the mass and spin of the Kerr background 1 [2][3][4].Precise measurement of these frequencies therefore characterizes the background [5] and moreover constrains deviations from general relativity (with more than one mode, or when combined with other measurements) [6][7][8][9].Although data today already hint at modes beyond the fundamental [10][11][12][13], future observations with sensitive detectors are sure to enable detailed spectroscopy [7,14,15].
To interpret future observations, however, it will be necessary to understand quasinormal mode interactions.The ringdown follows a highly nonlinear phase (the merger) and although numerical calculations indicate that a sum of modes may be sufficient to represent the gravitationalwave emission [16][17][18], it is not clear that this corresponds to a full nonlinear description.Indeed, nonlinear ringdown effects have been identified in numerical simulations of binary mergers [19,20] as well as in anti-de Sitter black holes [21,22].In other contexts (e.g., perturbations of large anti-de Sitter black holes) quasinormal modes can interact and even become turbulent [23,24].The point of this paper is to introduce some tools that may be helpful when developing a theory of quasinormal mode interactions.
Compared to normal modes, quasinormal modes do not in general form in a straightforward sense a complete "basis" of solutions to the linearized field equations.In fact, black hole perturbations are only described by quasinormal modes for an intermediate time period in their evolution; at early times they are described by a free propagation piece, and at late times by a power law tail [25][26][27].The spatial wavefunction of a decaying quasinormal mode also diverges at the bifurcation surface and at spatial infinity.This makes it hard to write down canonical (conserved) L2 -type inner products based on the usual Cauchy-surfaces of Kerr. 2 .Without an inner product, it is not clear how to project onto quasinormal modes to study nonlinear mode mixing.
The main goal of this paper is to point out an unconventional bilinear form which may take the place, for some purposes, of an inner product on quasinormal modes of Kerr.Before we introduce this notion, we develop a general theory for conserved -under time evolution -bilinear forms for Weyl scalars or metric perturbations.Similar to [32,33], the key idea is to start with a "Klein-Gordon" type current for Weyl scalars or metric perturbations and to apply symmetry operators to the entries of this bilinear expression.As we show, in Kerr spacetimes, such symmetry operators include, besides the obvious ones descending from the Killing symmetry, also an infinite tower of operators built from Carter's Killing tensor.(For the Weyl scalars, the symmetry operator of lowest differential order has two derivatives; for metric perturbations, it has six derivatives).In particular, using a combination of such operators we find an infinite set of new conserved, local, gauge invariant current associated with Carter's constant [34] in Kerr. 3he bilinear form of main interest for this paper is, however, not obtained from such differential symmetry operators but rather the symmetry operator associated with the discrete t-ϕ reflection.We show that gravitational quasinormal modes with different frequencies are orthogonal with respect to this bilinear form.For the reader interested in the main result, the bilinear form is presented explicitly for quasinormal modes in (55).We show furthermore that the quasinormal mode excitation coefficients of a solution are given precisely by the projection of data onto the corresponding modes via the bilinear form.
The plan of this paper is as follows.In section II we recall the standard recipe for constructing conserved bilinear forms for partial differential operators.In section III we introduce symmetry operators (including symmetry operators related to the Killing tensor, see also footnote 3) to construct further conserved bilinear forms, and currents.In section IV we construct the bilinear form ⟨⟨•, •⟩⟩ using the t-ϕ reflection symmetry, which gives orthogonality of quasinormal modes in section V. Finally, in section VI we explain the relation with excitation coefficients.Some technical aspects of this paper are deferred to various appendices.

II. BILINEAR FORM -BASIC CONSTRUCTION
Consider a partial differential operator X acting on sections of some vector bundle, E, over a manifold M .We assume that M is equipped with a volume form, ϵ a1...an ; later we will always have a metric g ab , so the volume form is chosen as the one compatible with the metric.Let Ẽ be the dual vector bundle, i.e., each fibre is given by the C-linear maps of the corresponding fiber of E. If ψ is a section of E and ψ is a section of Ẽ, we can pointwise form the scalar ψψ ∈ C. The formal adjoint is the unique differential operator X † defined by the formula where x a [ ψ, ψ] is local, i.e., at any point built from finitely many derivatives of the fields at that point.The divergence operator on the right is defined by our volume form and if it comes from a metric, as we assume from now, it is equal to the usual covariant derivative operator.Said differently, X † ψ is obtained by the usual "partial integration" procedure dropping surface terms as if the above equation were placed under an integral sign.Note that, by contrast to quantum mechanics, † as defined above is C-linear, rather than anti-linear.Now let ( ψ, ψ) be a pair of solutions to X ψ = 0 = X † ψ, and let Σ be a codimension 1 submanifold of M (later to be chosen as a constant t slice of Kerr).Then, by Gauss' theorem, if ψ, ψ have sufficient decay on Σ for the following integral to be suitably convergent (e.g., if they are compactly supported), then the bilinear form is unchanged under local deformations of Σ, and we say that it is "conserved".(Here ⋆ denotes the Hodge dual.)As a simple example, consider X = ∇ a ∇ a −m 2 , the Klein-Gordon operator acting on real-valued functions ψ, so E = Ẽ = R is the trivial bundle.Then X † = X and x a = − ψ∇ a ψ+ψ∇ a ψ is the Klein-Gordon (symplectic) current, which is of course conserved for any pair of solutions.The bilinear form in this case is just the symplectic form for Klein-Gordon theory.It is anti-symmetric under ψ ↔ ψ, but note that in the general case we cannot say that about the bilinear form since the bundles E and Ẽ cannot usually be identified in a natural way.As a second example, let E be the linearized Einstein operator on a Ricci-flat spacetime.It acts on symmetric covariant rank-2 tensors h ab , so E is equal to Sym(T * M ⊗ T * M ) in this case, and the dual bundle Ẽ corresponds to symmetric contravariant rank-2 tensors, Sym(T M ⊗ T M ).The formula is and under the identification of E with Ẽ (by using the metric g ab to raise indices), we have E † = E.As in the Klein-Gordon case, this last relation follows because the linearized Einstein equation arises from an action principle.By explicit calculation, the boundary term where The bilinear form is the symplectic form of General Relativity [40].
Our third, and most important, example concerns the Teukolsky operator(s) for the perturbed Weyl scalars of the Kerr spacetime (M, g ab ), to which we will restrict attention from now on.For this, we shall employ the GHP formalism [33,[41][42][43] in the following, and we now briefly review the essential portions of this formalism which simplifies and also conceptualizes many calculations in the Kerr -or more generally, Petrov type D -geometry.l a and n a are taken to be the repeated principal null directions which are completed to a null tetrad by defining a smooth pair of complex null rays (m a , ma ) that span the remaining dimensions.We choose the normalization l a n a = 1 and m a ma = −1, corresponding to the −2 signature.The metric then takes the form The basic idea is to contract any tensor field on M into the legs of the Newman-Penrose (NP) tetrad (l a , n a , m a , ma ) in all possible ways 4 and to represent the action of the covariant derivative operator ∇ a in terms of these tetrad components, in a way that preserves a natural grading by spin and boost weights.Fields η obtained by contracting with the tetrad are classified according to their spin and boost weights as follows.Under a local rotation that preserves the real null pair, the tetrad transforms as (l a , n a , e iΓ m a , e −iΓ ma ), whereas under a local boost that preserves the directions of the real null pair, it transforms as (Λl a , Λ −1 n a , m a , ma ), where Λ, Γ are smooth real-valued functions.If we combine these functions into the complex function λ 2 = Λe iΓ , then η is said to possess (real) GHP weights (p, q) if under the above combined local rotation and boost of the tetrad, it transforms as We write η ⊜ (p, q) if this is the case.In the GHP formalism, only quantities with the same weight may be added, whereas weights behave additively under multiplication.
From the mathematical viewpoint, the GHP formalism can be understood in terms of principal fibre bundles and their associated vector bundles, as follows.Consider the set of oriented null frames aligned with the given null directions.On each such frame, we may pointwise perform a boost/rotation, which as we described can be combined into a nonzero complex number λ ∈ C × .Thus, we have a multiplicative action of C × on the set of frames which gives this set the structure of a principal G-bundle: A principal bundle is abstractly a bundle P over M such that a group G can act by right multiplication X → X • g in the fibre -in our case X is an NP frame aligned with the principal null direstions and g ↔ λ.Given a principal G-bundle and a representation R of G on some vector space V , there is a canonical construction of an "associated" vector bundle.The sections of this bundle correspond physically to quantities defined on M that "transform in the representation R".More precisely, the elements in this associated bundle are the equivalence classes of pairs (X, v) where X ∈ P and v ∈ V where (X, v) is declared to be equivalent to (X • g, R(g)v).In the present example, R p,q (λ)v = λ p λq v and V = C, which corresponds precisely to the "transformation law" (8).The associated vector bundle is denoted in general by P ⋉ R V and its fibres are isomorphic to V .In our case, we get 1-dimensional complex ("line") bundles L p,q = P ⋉ p,q C over M labelled by the GHP weights (p, q).The number s = 1 2 (p − q) is commonly referred to as the spin.Of course, we could tensor L p,q with the usual tensor bundles T (r,s) M to host objects that have GHP weights and tensor indices at the same time such as l a or R abcd m a m d .
The advantage of the above invariant viewpoint involving associated vector bundles is that we can naturally see what quantities are defined in a frame independent manner, which quantities can naturally be added, etc.This provides not only an extremely useful guiding principle in the -usually very complicated -calculations related to Kerr, but also means that one is always intrinsically dealing with objects that behave in a well-defined manner under a change of frame.To make the formalism really useful, one needs covariant derivative operators on the bundles L p,q .These are given by The Teukolsky operators also feature the "gravitomagnetic potential" which is given by where ρ, τ are related to spin-coefficients [41][42][43]; see appendix F. The Teukolsky operator acts on GHP-scalars of the same weight 5 as the perturbed Weyl scalar ψ 0 , i.e., (p, q) = (4, 0) and is given by with Ψ 2 a background Weyl-scalar.So E = L 4,0 now.Since the dual vector bundle to L p,q is L −p,−q , the adjoint Teukolsky operator O † acts on GHP scalars of weight (−4, 0).It is given by It follows that the boundary term x a [ Υ, Υ] ≡ π a (with Υ ⊜ (4, 0), Υ ⊜ (−4, 0)) is given in the case of the Teukolsky operator by We denote the corresponding bilinear form -formally similar to the Klein-Gordon inner product of a charged scalar field -by 5 For the definition of O and O † for general GHP weights see appendix B.
The Teukolsky equation/operator and the linearized Einstein equation/operator are well-known to be related and this implies that the bilinear forms W and Π as in (6) and ( 14) are related, too.[44] have shown that for Υ a smooth solution to O † Υ = 0 arising from compact support data and h ab a smooth solution to Eh ab = 0, an identity of the following form holds where H ab is a skew symmetric local tensor.Furthermore [45,46] where Z abcd ≡ Z ab Z cd , and , are operators such that the Teukolsky-Wald identity holds: This equation encodes that the action T (h) on a metric perturbation h ab -which equals the perturbed Weyl scalar ψ 0 -gives a solution to Teukolsky's equation Oψ 0 = 0. Conversely, taking an adjoint of (17), i.e., ES † = T † O † , shows that any solution O † Υ = 0 of GHP weight (−4, 0) ("Hertz potential") is such that h ab = Re S † ab Υ is a solution to the linearized Einstein equations.
Ref. [44] did not derive the explicit form for H ab but argued for the above equation (15) to hold on general grounds based on (17).The main use of the above identity (15) is to relate the corresponding bilinear forms W [h, S † Υ] and Π[T h, Υ] for a Cauchy surface Σ of the exterior of Kerr.This identity is obtained by simply integrating the above identity over Σ.If all fields are falling off rapidly at the horizon and spatial infinity, then the boundary term arising from H ab will not contribute; in other cases, H ab will contribute surface terms.Their computation is fairly long and non-trivial and therefore deferred to appendix A. If Σ is a co-dimension one surface with boundary ∂Σ, Υ is a smooth solution to O † Υ = 0 and h ab a smooth solution to Eh ab = 0, then we have where B = ∂Σ H ab dΣ ab .When Σ is a slice of constant t in Boyer-Lindquist coordinates, ∂Σ would correspond to the bifurcation surface at r = r + and the sphere at r = ∞.Using this formula, the reader can readily transfer results on bilinear forms in this paper between the metric perturbation and Teukolsky variables.

III. BILINEAR FORMS FROM INFINITESIMAL SYMMETRY OPERATORS
Consider again a general partial differential operator X acting on sections of some vector bundle, E, over a manifold M .We have the corresponding conserved bilinear form X[ ψ, ψ] defined by (2).Now suppose C is a partial differential operator acting on E mapping solutions to X ψ = 0 to solutions -this is equivalent to the statement that there is a partial differential operator D such that X C = DX .Such an operator is called a "symmetry operator".The symmetry operators form an algebra which is trivial for a generic operator X .If we have a symmetry operator, then X[ ψ, Cψ] is also a conserved bilinear form, i.e., invariant under local changes of the surface Σ in (2), see e.g.[32,35] for a similar observation.
Let us apply this recipe to the linearized Einstein operator E on the Kerr spacetime.The Kerr spacetime has two Killing vector fields, t a , ϕ a corresponding to asymptotic time translations and rotations.The Lie derivatives L t , L ϕ evidently commute with E and thus provide two conserved quadratic forms: They correspond to the canonical energy and canonical angular momentum of the perturbation h ab when Σ is a Cauchy surface stretching between the bifurcation surface and spatial infinity [47].
If we want to repeat a similar construction for the Teukolsky operator O and the corresponding bilinear form Π we face the problem that the Lie-derivative in general is not well-defined on an arbitrary vector bundle (though it is on the usual bundles of tensors over M ).In the GHP formalism, the vector bundles L q,p in question are defined relative to an NP tetrad, and in such a case we can still give a definition of the Lie derivative along a Killing vector field, though not an arbitrary vector field, as we now describe.The point is that if g ab has an isometry φ that preserves the globally defined null directions, then this constitutes an intrinsically defined action on GHP tensors η ⊜ (p, q).More explicitly, if φ preserves the null directions, then it must be the case that it acts on a given null frame as φ * l a = Λl a , φ * n a = Λ −1 n a , and φ * m a = e iΓ m a , for some real functions Λ, Γ on M that depend on the chosen frame and φ.The action of φ on η is then invariantly defined since GHP tensors are functionals of the null tetrads giving rise to the prescribed pair of null directions.In the given null frame, this action amounts to φ GHP * η ≡ λ −p λ−q φ * η, where λ 2 = Λe iΓ and φ * is the standard pushforward on functions (or tensors).In particular, the tetrad vectors are invariant under φ GHP * .Infinitesimally, if φ t is a 1-parameter group of transformations generated by a Killing field χ a with corresponding λ t , then the corresponding "Lie" transport of η ⊜ (p, q) is given by [48] in the given frame.Here, L denotes the standard Lie derivative, and If we introduce the bivector Y ≡ n ∧ l − m ∧ m (for further details on the bivector calculus see, e.g.[45,49]) and use the fact that χ a is a Killing field, so ∇ (a χ b) = 0, then (20) can be manipulated to obtain where L Θ is the standard Lie derivative with ∇ a derivatives replaced by Θ a derivatives.In this notation, the GHP Lie derivative is also defined for GHP-tensors, i.e., sections in a bundle L p,q tensored with T M or T * M .In any case, the GHP Lie derivative defined here is manifestly GHP covariant, and it can be checked that it satisfies the Leibniz rule.The expression for Ł χ in a chosen NP tetrad will depend on that choice.For the Kinnersley tetrad (F5), w = 0, but w can be different from zero for other choices of the frame.
With these definitions, it then follows that Ł χ for χ a either t a or ϕ a commutes with the covariant derivative Θ a and annihlates g ab , n a , l a , m a , ma , B a .Therefore, Ł χ also commutes with the Teukolsky operators, and it thus defines a symmetry operator.The corresponding conserved currents arising from π a (13) using the general construction described above have been discussed by [33].
There exist other symmetry operators in the Kerr (and more generally, Petrov type D-) spacetimes related to the Killing tensor K ab that exists in those spacetimes.The construction of those operators for spin s = 0, 1  2 in the Teukolsky equation goes back to [32,35]; here we present the corresponding symmetry operator for arbitrary GHPweights (p, q).Similar operators have appeared also in [36,37], eq.III.3, for spin s = 1, 2, though not in the GHP covariant form presented here which makes manifest the relationship with the Killing tensor.This tensor is given by where we use the shorthand with ρ one of the spin coefficients in the GHP formalism.The desired symmetry operator K acting on GHP scalars of weights (p, q) is defined as where is proportional to a Killing vector field, and γ = (ζ 2 − ζ2 )/(8ζ).Here and in the following, a prime as in B ′ a means the GHP priming operation n a ↔ l a , m a ↔ ma .In Boyer-Lindquist coordinates and the Kinnersly frame when acting on GHP quantities of weight (4, 0) or (−4, 0), respectively.The proof of this statement is rather nontrivial and deferred to appendix B, where we also prove the commutation property for arbitrary (p, q).It follows from the properties of the GHP Lie derivative that [Ł χ , K] = 0 for any Killing vector field χ a , so we have: Theorem 1. Ł t , Ł ϕ , K generate a commutative, infinitedimensional algebra of symmetry operators for Teukolsky's operator O for any GHP weights (p, q).
Hence, by the general scheme, if we have solutions to O Υ = 0 = O † Υ, and symmetry operators A, B, then the bilinear form Π[A Υ, BΥ], with Π as in (14), is conserved, i.e. unchanged under local deformations of the Cauchy surface Σ.We caution the reader that such bilinear forms can be trivial, i.e., be equivalent to forms that are conserved identically; see appendix C for some discussion.
It is possible to derive symmetry operators also for the linearized Einstein tensor E (and for the Maxwell equations) on Kerr or more generally, a Petrov type D spacetime.Let n = 0, 1, 2, . . .and set as well as where for spin-2 considered here we should take s = 2, and where we use the GHP priming operation.Then where we used twice the Teukolsky-Wald identity (17), the commutation [O † , K] = 0, as well as the intertwining relation When acting on a perturbation h ab , T ′ (h) gives the perturbed Weyl scalar ψ 4 , which is gauge invariant.Therefore, we see that C n (h) = 0 for any gauge perturbation h ab = L ξ g ab .
By the results of appendix B another symmetry operator for E would be with similar proof, see appendix B for the definition of G), and further symmetry operators are obtained by the GHP prime-and overbar operations applied to these C n 's.Finally, by putting s = 1 in the above expressions, and defining T , S so that the analog of the Teukolsky-Wald identity (17) holds for electromagnetic perturbations, where (EA) a = ∇ b ∇ [a A b] , we get similar operators in the electromagnetic case.
As a consequence, in all cases, C n give symmetry operators for E of order 4 + 2n for spin-2 and of order 2 + 2n for spin-1.Regarding our operator C 0 for spin-2, we remark that a very similar looking operator has been considered by [36], Eq.III.14.Regarding our operator C 1 , a similar looking operator has been considered in [36], Eq.III.47 and also in [39], Thm.16.However, closer inspection of the operator6 in [36], Eq.III.47 shows that it is non-local, while our operators are all local and also manifestly GHP covariant.The relation of our operators C 1 to the order 6 symmetry operator asserted in [39] is not completely clear to us and the same goes for our other operators C ′ 1 , etc.For spin-1, symmetry operators of orders 2 and 4 have been discussed in [37,38], and the comparison to ours is qualitatively similar. 7y the general theory, for example (n = 0, 1, 2, . . . ) with W as in ( 6) are conserved for all solutions h ab to the linearized Einstein equations, i.e. unchanged under local deformations of the Cauchy surface Σ.The corresponding conserved currents are with w a as in (4).Note that each j a (n) is local and gauge invariant from the properties of C n .The concrete expressions of j a (n) are very long and contain 2n + 9 derivatives of h ab .For the reason explained below, we call j a (n) the "Carter current(s)".
To gain some insight into the meaning of the conserved quantities χ (n) , we make a WKB (high frequency) analysis similar to [50], see also [36].If the momentum of the sharply collimated WKB wave packet h ab is p a and its amplitudes defined with respect to a suitable basis of polarization tensors are A +,× , the result is where K ab p a p b = Q(p) denotes the Carter constant.See appendix D for more detail on the derivation of this formula and on the precise definitions of the WKB wave functions, polarizations, etc.We can obviously form alternative conserved quantities by other combinations of the various symmetry operators of the linearized Einstein operator described above giving e.g., the GHP primed version of our Carter currents j a′ (n) .We note also that such currents could have alternatively been constructed from π a (13), taking Υ = ζ 4 ψ 4 and Υ = ψ 0 and acting on those with various symmetry operators for the Weyl scalars, as described above.
We finally remark that conserved currents for metric perturbations related to Carter's constant have also been considered in [36], Eqs.IV.14-16.Eq.IV.14 is very similar to our j a (0) but their currents Eqs.IV.15-16 are different from our Carter currents j a (n) or their GHP primes because unlike ours, they are based on non-local currents requiring a mode decomposition of the solutions.

IV. BILINEAR FORM FROM t-ϕ REFLECTION
In the previous section, we combined the basic conserved bilinear form (14) with symmetry operators, which arise in particular from the Killing vector fields of Kerr.One naturally expects that a similar construction should be possible for the discrete isometry of Kerr, namely the t-ϕ reflection map J : (t, ϕ) → (−t, −ϕ) where here and in the following we refer to Boyer-Lindquist coordinates.However, just as for Killing vectors, some care has to be taken when defining the action of J on GHP scalars with nontrivial weights (p, q).So we first turn to this issue.
The map J swaps the null directions l a and n a and changes the orientation on the orthogonal complement of these null directions spanned by m a , ma .There must thus be Λ, Γ depending on the null tetrad such that J * l a = −Λn a , J * n a = −Λ −1 l a , and J * m a = e iΓ ma , where we have defined J to act on tensors by the push-forward.By analogy with the previous case of isometries which are continuously deformable to the identity, it is then natural to define for η ⊜ (p, q) a GHP reflection in the given frame.The operator J is evidently a GHP priming operation combined with t → −t, ϕ → −ϕ, and is therefore easily seen to be GHP covariant (i.e.defined intrinsically as a map from sections in L p,q to sections in L −p,−q , irrespective of the chosen frame), but, by contrast to the "pull-back" arising from isometries continuously connected to the identity as considered above, it changes the GHP weights.In this sense it is similar to the CPT operator arising in quantum field theory.It is clear that J 2 = 1 and one can relatively easily show the "anti-commutation" relations Ł t J = −J Ł t , Ł φ J = −J Ł φ with the GHP Lie-derivative defined above.We also note an important intertwining property of the t-ϕ reflection operator J with the Teukolsky operator and its adjoint, namely, where we used basic properties of gravito-magnetic field B a and its GHP prime B ′ a , as well as the relation In the Kinnersley frame and Boyer-Lindquist coordinates (see appendix F), the J operator corresponds to sending t → −t, ϕ → −ϕ and multiplication according to (36) by appropriate powers of λ, λ, where λ is given in this case explicitly by We are now in a position to define the bilinear form.For simplicity, we restrict at first to entries having compact support on the Cauchy surface Σ in order to avoid any convergence problems.
Lemma 2. Under the conditions of the definition, we have for t a the time translation Killing field, and (iv) ⟨⟨Υ 1 , Υ 2 ⟩⟩ is independent of the chosen Cauchy surface Σ.
Before we prove this lemma, we remark that, e.g. by (46), the bilinear form may be viewed as defined on the initial data of the Teukolsky equation on the Cauchy surface Σ.On an initial data set Ł t corresponds to the action of a suitably definined Hamiltonian operator H. Then item (iii) corresponds to the statement that i.e. to the fact that the Hamiltonian operator is symmetric with respect to our bilinear form.We refer the interested reader to appendix E for details on the Hamiltonian formulation of the Teukolsky equation.
We also note that although we defined our bilinear form on s = −2 GHP scalars (i.e., solutions to the adjoint Teukolsky equation), we could also define a bilinear form on s = +2 solutions to the original Teukolsky equation.In this case, we set ⟨⟨ Υ1 , It can be shown that the s = +2 bilinear form satisfies all the same properties as the s = −2 form.Proof.(i) This is obvious from the definition.
(ii) By explicit calculation, we have with π abc = ϵ abcd π d and π a as in (13), using J 2 = 1 and (38).Now integrate over Σ.Since J reverses the orientation of Σ, the claim follows.
(iii) We first remark that, by Cartan's magic formula, we have that on solutions (where π = π abc dx a ∧ dx b ∧ dx c ), if dπ = 0. Integrating over Σ and using Stokes's theorem, as, for compact support data, the contribution on ∂Σ evaluates to zero.In our case, ) is indeed closed, dπ = 0. On the other hand, we have, since background quantities are all GHP-Liederived by t a = M 1/3 ξ a , and since Inserting this into the left hand side of (44) evalu-ated on the solutions Ψ We end this section with an explicit expression of our bilinear form in Boyer-Lindquist coordinates and the Kinnersley frame: where we refer to appendix F for the definitions of Σ, ∆, and Λ.8 V. QUASINORMAL MODE ORTHOGONALITY

A. Quasinormal modes
Consider modes of the form with m ∈ Z and ω ∈ C, in the Kinnersley frame.This form leads to separation of the spin-s Teukolsky equation [51], OΥ = 0 (for any integer spin s), into an angular equation, and a radial equation, with H ≡ (r 2 + a 2 )ω − am.Here K is a separation constant.Imposing regularity at the poles θ = 0, π, the angular equation leads to a discrete set of modes s S ℓmω and separation constants s K ℓmω , both of which are indexed by ℓ ∈ Z ≥max(|m|,|s|) .The functions s S ℓmω (θ)e imϕ are known as spin-weighted spheroidal harmonics [51].For ω ∈ R, the angular problem reduces to a Sturm-Liouville eigenvalue problem.Modes with the same s, m, and real ω, but different ℓ are orthogonal, and we normalize them such that Orthogonality can be checked by verifying that the angular operator is symmetric with respect to this product.
To discuss boundary conditions of the radial equation it is convenient to introduce a "tortoise" coordinate dr * = (r 2 + a 2 )/∆dr, see (F3).For fixed s, l, m, ω one considers the solutions R in and R up "defined" by the "boundary conditions" where k ≡ ω − mΩ H , where Ω H is the angular frequency of the outer horizon Ω H = a/(2M r + ), and where the radii of the inner-and outer horizons (roots of ∆) are denoted by r ± , respectively.The conditions (51) correspond physically to the absence of incoming radiation from the past horizon and past null infinity, respectively.As stated (51) do not really pick out uniquely a solution in the case Imω < 0 because we may always add a multiple of the subdominant solution as |r * | → ∞ without affecting the asymptotic behavior.More precisely, mode solutions may be obtained via series expansions [52], involving three-term recurrence relations for the coefficients.Selecting the so-called "minimal solution" [53] of the recurrence relations ensures that the series represenation converges at the horizon (in) or infinity (up). 9Imposing both of these conditions simultaneously 10 gives rise to a discrete set of quasinormal modes ω n ∈ C, where n = 0, 1, 2, . . .are the so-called "overtone" numbers.We restrict to frequencies with Im ω ≤ 0, as modes growing exponentially in time are not in the specturm of Kerr [55].

B. Bilinear form on quasinormal modes
We would now like to extend our definition of the bilinear form ⟨⟨•, •⟩⟩, originally only for compactly supported solutions/data on the Cauchy surface Σ, to quasinormal modes.The immediate problem is that, according to the boundary conditions on the corresponding solutions to the radial equation, these blow up both at the horizon r = r + and infinity r → ∞.In this subsection, inspired by the work of [56], we show that the Kerr bilinear form can be defined for quasinormal mode data by a suitable 9 This definition is satisfied by a radial solution of the form where with dn coefficients that are a minimal solution to a three-term recursion relation [54] so that the series is uniformly absolutely convergent as r → ∞ 10 The problem is made complicated, however, because ω and K appear in both the angular and radial equations, ω nonlinearly.One must jointly solve both equations to obtain a self-consistent solution of this nonlinear eigenvalue problem.Using Hamiltonian methods (see appendix E) one can recast this as the eigenvalue problem HΥ = iωΥ, i.e., the problem is linear in ω, but the angular and radial problems remain coupled.
deformation of the radial integration into the complex plane. 11onsider the bilinear form acting on two quasinormal modes with quasinormal frequencies ω 1 and ω 2 .The integrand in the bilinear form (55) goes as ∼ e ±i(ω1+ω2)r * as r * → ±∞, and therefore diverges exponentially for Im(ω 1 + ω 2 ) < 0, which is the case for all modes that decay in time.Therefore, we clearly see that the bilinear form as defined for compact support data ( 46) is divergent.
We can obtain a finite bilinear form by analytic continuation in r.The radial mode functions R in/up (r) are analytic with branch points at r = r ± [52], and we take the branch cut as the wiggly line in Fig. 1 going from r + to r − .We take the branch cut for the tortoise coordinate (F3) r * (r) to be identical, so that we can think of both the radial functions R in/up and r * as defined on the same multisheeted covering of the twice cut complex r-plane.The integrand of the bilinear form, given by the 3-form π abc = ϵ abcd π d [see (13)] evaluated on two mode functions as in (40) or equivalently (46), therefore has an analytic continuation on the multi-sheeted complex r-plane.
In (40) or equivalently (46), we now define an integration contour going into this complex r-plane as shown qualitatively in fig. 1.In terms of r * (r), which is a function on the same multi-sheeted complex r-plane, the contour is defined in such a way that 0 < arg((ω 1 +ω 2 )r * ) < π on the right limit, and −π < arg((ω 1 + ω 2 )r * ) < 0 on the left, then as |r * | → ∞, the volume integral will converge exponentially with |r * |.
To achieve this for any Im(ω 1 + ω 2 ) < 0 me may take a snake shaped contour u → r * (u, ϵ) of the radial coordinate in the complex r * plane, with the properties where r * 1 < 0, r * 2 > 0 can in principle be chosen arbitrarily.We give a sketch of this contour, C * , which corresponds to one in terms of r, in the right panel of Fig. 1.The corresponding 3-dimensional submanifold (depending on ϵ > 0 and on t ∈ R) of the analytically continued Kerr manifold M C is denoted by In practice, the angle ϵ > 0 is chosen sufficiently small such that the integral in the following definition of the bilinear form converges, Replacing Σ with the contour Σ C as described in section IV, thanks to the analyticity of the integrand FIG. 1. Left: Sketch of the complex r contour C * defining the bilinear form on quasinormal modes.The contour cannot be pulled back to the real axis because the integrand crosses (an infinite number of) different sheets associated with the branch points r− and r+.Right: Same contour, but in the complex r * plane.Note that this contour cannot be pulled back to the real axis due to the presence of Stokes lines along which the integrand of the bilinear form would diverge.
and its fall off on ∂Σ C , all properties of the bilinear form of of lemma 2 continue to hold on quasinormal modes.In particular, from item (iii) of lemma 2, we get (ω 1 − ω 2 )⟨⟨Υ 1 , Υ 2 ⟩⟩ = 0 for a pair of quasinormal modes with complex frequencies ω 1 , ω 2 .Furthermore, by (iv), the value of the bilinear form is independent of the precise choice of t, details of the complex integration contour such as the asymptotic angle ϵ against the real half-axes and/or r * 1 , r * 2 , as long as the integrand is exponentially decaying.
Corollary 3 (Orthogonality of quasinormal modes).Let Υ 1 and Υ 2 be quasinormal modes for the s = 2 Teukolsky equation with frequencies ω 1 and ω 2 .Then either Our bilinear form takes the following form on quasinormal mode solutions (47).After plugging two s = −2 mode solutions in separated form into (46), we can carry out the ϕ integration to obtain with C * the contour for the r * -integration described above and the Kerr quantities ∆, Σ, Λ as given in appendix F.
The integrands depend on θ and r in a nonfactorizable way, so this expression is the best that can be achieved in general: for Kerr, the orthogonality relation expressed by the previous corollary (vanishing of the above inner product for ω 1 ̸ = ω 2 ) is fundamentally two-dimensional.This has to do with the fact that the orthogonality relation (50) for spin-weighted spheroidal harmonics occurs between modes of different ℓ but the same m and ω; if ω 1 ̸ = ω 2 , then no such relation exists, and one cannot expect to be able to perform the θ integration to obtain a δ ℓ1ℓ2 factor.
In the a → 0 Schwarzschild limit, however, the integral does factorize: the θ dependence of the integrand reduces to the sin θ volume factor on the sphere, the spheroidal harmonics reduce to spherical harmonics (independent of ω), and the θ integral is proportional to δ ℓ1ℓ2 .One is left with a radial integration, which must vanish for ω 1 ̸ = ω 2 .
In Fig. 2, we numerically evaluate the bilinar form on two pairs of Kerr quasinormal modes.We do so along the most convergent contour -along which the integrand is purely damped -for the given pair of frequencies ω 1 , ω 2 : arg r * (u) + arg(ω 1 + ω 2 ) = π/2 for r * → ∞ and arg r * (u) + arg(ω 1 + ω 2 ) = −π/2 for r * → −∞.In other words, the two complex sections of the contour are given by r * (u) = r * ,1,2 + ue −i arg(ω1+ω2)+iπ/2 with u < 0 in the section emanating from r * 1 and with u > 0 in the other (see again Fig. 1).We then integrate u between 0 and a finite u lower, upper , respectively, and study the convergence of the integral as we take u upper, lower → ±∞.As we can see from the figure, the contour integral in the bilinear form converges quite well, which is useful in practice when using it to extract excitation coefficients, as we describe in the next section.Furthermore, since orthogonality is an exact result for quasinormal modes, it can be used potentially as a benchmark check for approximations.For example, we have considered approximations to quasi-

FIG. 2. Numerical check of the orthogonality between two
Kerr quasinormal modes with the same l = m = 2 and different n = 0, 1 (upper panel) and modes with the same n = 0, m = 2 and different l = 2, 3 (lower panel).We show the result of the numerical evaluation of the bilinear form (55) along the most convergent contour (black points), integrating the complex sections of the contour up to a finite u lower, upper .Because of the difficulty in handling the branch cut in Mathematica, the lower integration limit u lower sets the overall accuracy of this numerical evaluation of the bilinear form.We also show an exponential fit converging to zero as u lower → −∞ (red line).We use the mode solutions (normalized at the horizon) provided by the Black Hole Perturbation Toolkit [61].In this example, we set M = 1, a = 0.7, r * 1 = −8, and r * 2 = 5.
normal modes based on a matched asymptotic expansion for near-extremal black holes, and have found that the orthogonality relation is typically satisfied to a very high accuracy.
We remark that our "norm" on quasinormal modes has some similarities with the "norm" of resonant state wave functions in quantum mechanics defined by [62].Rather than taking the integral of |ψ| 2 , the "norm" used by [62] also involves ψ 2 , whereas our bilinear is complex linear in both arguments as opposed to an inner product (antilinear in the first argument, complex linear in the second).Our regularization procedure differs from that proposed by [62] but was rather inspired by the investigations [56] in the context of leaky optical one-dimensional cavities, and on Schwarzschild black holes in [63].In [64,65] it was recognized that phase space was the natural setting for the bilinear form. 12In fact, in several ways our work was inspired by some of these papers: we work within the Teukolsky formalism, we arrive at the bilinear form starting from the symplectic form, and we recognize the fundamental importance of the t-ϕ reflection symmetry.

VI. EXCITATION COEFFICIENTS
If ⟨⟨•, •⟩⟩ were an honest to God scalar product in a Hilbert space and { s Υ ℓmn } an orthonormal basis, then an arbitrary wave function Υ s could evidently be expanded as where the excitation coefficients are Here ℓmn denotes In the present context, ⟨⟨•, •⟩⟩ is of course only a symmetric bilinear form on solutions to the spin s Teukolsky equation (for the case of interest in this paper, s = −2).It is neither positive definite, nor is the set of quasi-normal modes, while being orthogonal, in any obvious mathematical sense a complete basis for a reasonable function space in as far as we can see.
Inspired by [56], we will nevertheless show in this section that for solutions Υ s to the adjoint Teukolsky equation with compact support on a Cauchy surface Σ, the above expansion can formally be "derived" in the Laplace transform formalism [3,26] for the retarded propagator if we deform the frequency integration contours into the complex plane and collect only contributions from the quasinormal mode frequencies.Thus, (57), while not an exact equality, is expected to capture the transient behavior of the solution Υ s .

A. Laplace transform
The Laplace transform f (ω) = Lf (t) of a function f (t) is given by where Im ω > 0. The Laplace transform is related to the Fourier transform Ff = ∞ −∞ e iωt f (t)dt by sending f (t) → f (t)θ(t), where θ(t) is the Heaviside distribution.
Sufficient conditions for the existence of f (ω) are that the function f (t) be Riemann integrable (continuous except on sets of measure zero) on every closed sub-interval of the path of integration and that it be of exponential order; i.e. at any t one can find constants a and N such that |e −at f (t)| < N .If the Laplace integral exists for some value of ω = ω 0 , then it also exists for all ω with Im ω > Im ω 0 .The lowermost Im ω 0 where convergence occurs is called the abscissa of convergence and the region above this line called the convergence region.The function f (ω) is analytic in the convergence region.
The Laplace transform formalism is naturally adapted to the study of causal dynamics of linear second-order systems, as it incorporates the initial data into a source by taking time derivatives into field values at the initial time The Laplace transform Υs of the spin-s master function Υ s is given by This decomposed into modes in the usual way The inverse transform is given by where c > 0 is chosen such that the integral contour lies within the convergence region.
To formulate the initial data problem within the mode decomposition, we take the Laplace transform of the Teukolsky master equation and substitute (62).We then collect the terms in the master equation with transformed time derivatives to the right-hand side and project onto the angular mode function.This yields a sourced equation for the radial function Here, L is given in ( 49) and the source s I ℓmω = s I ℓmω (r) is comprised of (ℓ, m)-projected initial data [68]: which we take to be of compact support.Imposing the outgoing boundary conditions (51) fixes the freedom of homogeneous solutions to (64) and therefore defines a radial Green's function s g ℓmω (r, r ′ ) = s R in ℓmω (r < ) s R up ℓmω (r > )/W where r < (r > ) is the lesser (greater) of r and r ′ .Here, for any two solutions of the radial equation at fixed s, m, ℓ, ω, the ("∆-scaled") Wronskian has been defined which is independent of r.If R 1 and R 2 are linearly dependent, then the Wronskian vanishes.Thus, if we take R 1 → R in , R 2 → R up , the Wronskian vanishes when ω attains a quasinormal frequency.The quasinormal mode contribution to Υ s can be found by closing the contour of the Laplace integral in the lowerhalf complex ω plane, and can be expressed as a discrete sum over the residues of the radial Green's function arising at the points in the complex frequency plane where the Wronskian vanishes: where s I ℓmn = s I ℓmω | ω=ωn .In considering only the poles when closing the contour, we are effectively ignoring the early-time "direct" contribution from the large-ω arc and the late-time "tail" contribution resulting from the branch point at zero frequency.Thus, the = sign in the above equation is not actually justified and should be understood as meaning this approximation.On a quasinormal mode, R in is a constant multiple of R up , and either may be moved outside the radial integral to write the field as where we have isolated the familiar form of the excitation coefficient B. Equivalence between (69) and (57) To begin, for a given radial function R, and ℓ, m, ω, we can define a s = −2 GHP scalar Υ ℓmω in separated form (47) by appending a spin-weighted spheroidal harmonic and e −iωt time-dependence.For the time being, we do not require R to satisfy any equation.We have the following lemma relating the Wronskian to t • π integrated over the 2-sphere.Lemma 4. Let Υ 1 , Υ 2 ⊜ (−4, 0) be two GHP scalars in separated form (47), with the same m, ℓ, ω, where S 1 , S 2 are normalized spin-weighted spheroidal harmonics solving the angular equation, but where R 1 , R 2 are not necessarily solutions to the radial equation.Then where S 2 (t, r) is a sphere of constant t and r, π abc = ϵ abcd π d and π a as in (13).
Proof.Consider the Cauchy surface Σ(t) = {t = const.}in Boyer-Lindquist coordinates.The future directed normal to Σ and induced area element on S 2 (t, r) are given by, respectively From the first relation on can read off the lapse function N of t a from ν a = (t a − N a )/N .The action of the reflection reverses ν a and from this fact and the formula for π a , see (13), one can deduce that where r a is the normal to S 2 (r, t) inside Σ(t).An explicit calculation shows that in the Kinnersley frame, Using this, as well as expressions for N , dA, and J [using (39)], we obtain Finally, perfoming the integration and using the normalization (50) for the angular functions, we obtain the result.
Next, we take R 1 and R 2 to be solutions ingoing at the horizon and outgoing at infinity.Considered as a function of ω, the Wronskian vanishes at quasinormal frequencies ω n , because at these frequencies the two solutions become linearly dependent.The first derivative with respect to ω, however, is proportional to the "norm" of the quasinormal mode.
Lemma 5. Let R in ω , R up ω be solutions to the radial equation for fixed s = −2, ℓ, m, and allowing ω to vary, that are ingoing at the horizon and outgoing at infinity, respectively, as in (51).Construct Υ in ω , Υ up ω ⊜ (−4, 0) as mode solutions based on the radial functions.Then the derivative of the Wronskian at a quasinormal frequency ω n can be written Proof.Consider the current π evaluated at a generic frequency ω and a quasinormal frequency ω n .By Cartan's magic formula and the fact that π is closed on solutions, where in contrast to the previous lemma the right side does not vanish on account of the different frequencies.
We integrate over a partial Cauchy surface S, and apply Stokes's theorem, Next, we differentiate this equation with respect to ω and take the limit ω → ω n .On the right side, we trivially get The left side can be expressed as three terms, namely d dω ω=ωn l.h.s. of (78) By lemma 4 we can write the second of these terms as the derivative of the Wronskian, Summarizing our results so far, we have shown that As S → Σ C , the boundary integrals vanish exponentially, so the right hand side of (82) reduces to −i⟨⟨Υ in ωn , Υ up ωn ⟩⟩.
The desired equivalence between ( 69) and ( 57) can be seen immediately by substituting (65) for s I ℓmn , comparing with (46), and applying lemma 5.
We checked numerically that the excitation factors obtained via the bilinear form, Eq. ( 76), are in good agreement (for some a, n within a few %, which we believe arises from inaccuracies of our mode functions close to the horizon) with the ones computed via direct differentiation of the Wronskian in Refs.[69,70].

VII. CONCLUDING REMARKS
We end this paper with some potential applications and alternatives to our formalism.
The main motivation of this work is to provide some tools needed to study the black hole ringdown beyond linear order in perturbation theory.Higher orders are already needed to interpret high-precision numerical relativity simulations of binary mergers [19,20], and could be needed to analyse gravitational wave observations by future detectors.Roughly speaking, we wish to make an ansatz for the solution of the non-linear system as a linear combination of quasinormal modes with time dependent excitation coefficients similar to (56).The idea is that the bilinear form will help us writing down a dynamical system for these coefficients by analogy with wave equations on compact spaces, where the normal modes would be used instead to compute the overlap integrals required for terms in this dynamical system that are non-linear in the modes.
As extremality is approached, it is well known that a family of quasinormal modes becomes arbitrarily longlived, with Re ω ≈ mΩ H [54,[71][72][73][74][75]. With a commensurate frequency spectrum and arbitrarily slow decay, these modes have been conjectured to become turbulent as a → M [76].By considering the nonlinear excitation coefficients near extremality using the approach detailed above, one may hope to establish (or rule out) the emergence of turbulent behavior, or to get a new perspective on the Aretakis instability [77].It would also be interesting to see if our bilinear form can be used directly at extremality.This would involve understanding better the behavior of the complex contour integrals that we employ to regulate the bilinear form at extremality. 13pplications along similar lines could include mode mixing in clouds of ultralight scalar fields that could form outside Kerr black holes [80].Here, once these clouds grow via the superradiant instability, nonlinear interactions between the modes have been conjectured to give rise to a coherent emission of gravitational waves or a bosenova [81], see [82,83] for recent proposals based on heuristic methods.It would be interesting to see whether our methods could be used to conceptualize or shed more light on the theoretical basis of such proposals.
In [31], a different approach is taken to quasinormal modes.Their essential idea is to consider, instead of a time t Cauchy surface intersecting the bifurcation surface and spatial infinity, a "hyperboloidal" slice intersecting the future event horizon and future null infinity.On such a slice, they define a certain space of "almost analytic" functions (a "Gevrey space") encoding somehow the "boundary conditions" (51).Their space is in fact a genuine Hilbert space, and the time evolution is represented on this space by a semigroup whose generator is essentially the Hamiltonian, H (see appendix E).Their inner product is non-canonical -and is not conserved -and the generator of the semi-group is correspondingly not symmetric, as is also not physically expected due to the "dissipative" nature of quasinormal modes.Nevertheless, their analysis shows that quasinormal modes are genuine eigenfunctions HΥ = iωΥ in this space -crucially, by contrast to their restriction to a constant t surface, they do not blow up on the hyperboloidal slice as the horizon or scri are approached.While the quasinormal modes are not orthogonal in their inner product, the definition of our bilinear form with a hyperboloidal slice is also clearly possible and it would be interesting to see whether quasinormal modes, as defined in the setting of [31] (see also [28][29][30]) are still orthogonal in this bilinear form, as we conjecture.This would provide an alternative to our regularization procedure involving complex contours.
Finally, it will be interesting to explore the relation between our bilinear form and the adjoint-spheroidal functions introduced in Ref. [84].
Of course, H ab is defined by the above equation only up to a total local divergence ∇ c C [abc] , or equivalently, ⋆H is only defined up to an exact local 1-form.Therefore, the expression for H ab given below only represents one possible solution.Unfortunatlely, the proof given in [85] is not really useful to actually find an H ab in practice, so we simply try to "peel off" the total divergence from the right side of (A8) by hand.
For this, it is convenient to write out the divergence operator on H ab in terms of GHP operators as We now substitute the above expressions for S † Φ, Eh, O † Φ and T h respectively into w a , σ a , t a and π a on the right side of (A8).Next, by comparing the result with (A9) we should in principle be able to get the components of the 2-form H ab [Φ, h] up to a total divergence.This calculation is extremely tedious and was done using the Xtensor package of Mathematica, along the following lines.First, we neglect the terms containing less than three derivatives on the right side of (A8) and try to "peel off" a total divergence by making suitable choices of the NP components of H ab appearing in (A9).In fact, we cannot peel off a divergence exactly, but only at the cost of terms containing less than three derivatives.Having determined the highest derivative parts of the NP components of H ab , we collect any terms with less than three derivatives that were left over when peeling of the total divergence, and combine them with all terms with less than three derivatives on the right side of (A8) that were neglected so far.Of these remaining terms, we now discard all terms with less than two derivatives and repeat the process until no terms are left.We are ensured on general grounds that this must happen, though it is in practice rather non-trivial and time-consuming to determine all terms.The corresponding Mathematica notebooks are provided at [86].All in all, it is found that a skew symmetric tensor H ab which satisfies (A8) is given by the following expression. where When integrating H ab over a 2-surface we need the Hodge-dual, which is given by We remind the reader that ⋆H is unique only up to an exact local form.Such a contribution could be introduced e.g. to simplify the form of ⋆H.We have not attempted to do this.Finally, we remark that for spin-1, the tensor H ab has already been constructed in [87].It seems to have a much simpler form than in the spin-2 case.
Appendix B: Proof of (29) In this appendix, we will prove (29).Although that relation is used in the body of the paper only for gravitational perturbations, i.e., spin s = (p − q)/2 = ±2, most of the calculations work for arbitrary s, and in fact p and q.We therefore state the results in their most general form in view of their potential applicability to various values of the spin and for general Petrov type D spacetimes, of which Kerr is an example.
We begin by defining where B a was given above in (10), ζ was given in (26), ξ a was given in (28), and where here and in the following we use the GHP priming operation.Strictly speaking, P denotes not one operator, but sevaral operators depending on the chosen (p, q), and the same is understood for other operators below.
The following lemma is checked by direct computation.
Lemma 6.The operators P and K (see (27)) have the following properties: (i) They are self-dual, i.e., K † = K and P † = P, (ii) They are real, i.e., Kη = Kη and Pη = P η for every smooth weighted scalar η, where we mean the GHP overbar operation, where s = (p − q)/2 is the spin of the theory, i.e., s = 1 for electromagnetic-and s = 2 for gravitational perturbations.
Next, we construct the following operators G, R and L by combining K and P in a certain way.
Here K = g ab K ab is the trace of the Killing tensor K ab (25).See (27) for the definition of the "Carter" operator K for general (p, q).We will make use of the following lemma 7 observed by [88].
Lemma 7. Suppose we have partial differential operators A, B (defined on suitably compatible vector bundles) such that for all f , where α, β are certain functions, i.e., multiplication operators.Then This lemma is used to prove the following result.
Theorem 8. Let η ⊜ {p, q}.Then in any type D spacetime Through a long tedious calculation using heavily the type D property of the background, one can show that Consequently, using the lemma 7, we conclude that Further, In addition, we can also conclude that This demonstrates the claimed relations in the theorem.

Appendix C: Trivial conservation laws
An question when deriving conservation laws is whether the current is nontrivially conserved.A conservation law is said to be trivial if either (a) the conserved quantity itself vanishes on solutions, or (b) the quantity is conserved even on non-solutions [89].A generic trivial conserved quantity will be a combination of (a) and (b), e.g., if by adding a term involving the equation of motion to a conserved current, it becomes conserved even on non-solutions, then the conservation law is trivial.
As an example, consider a real Klein-Gordon field ψ, and a second order symmetry operator L 2 ξ for a Killing vector ξ a .The associated current x a (ψ, L 2 ξ ψ) is conserved as follows: (C1) Thus if we add the quantity ξ a [(X ψ)(L ξ ψ) − (L ξ X ψ)ψ]-which vanishes on solutions X ψ = 0-to x a (ψ, L 2 ξ ψ), then the resulting current is identically conserved.Hence the current associated to L 2 ξ is trivial.It is easy to see that in this context symmetry operators L n ξ for odd n (i.e., skew-adjoint operators) generate nontrivial conserved quantities, whereas those with even n generate trivial ones.
In general, for a symmetry to generate a nontrivial conservation law, it must be a symmetry not just of the equation of motion, but also the Lagrangian.A complete analysis of triviality for the conserved currents discussed in the present manuscript is left for future work.

Appendix D: WKB analysis of conserved currents
We make a WKB (high frequency) ansatz for the metric perturbation of the usual form h ab ∼ Re A ab e iωS (D1) with real frequency ω ≫ 1.For h ab to be a solution to the linearized Einstein equation, S should as usual satisfy the eikonal equation g ab ∇ a S∇ b S = 0, so that k a = ∇ a S is tangent to a future directed -by assumption -congruence of affine null geodesics.We may take the (complex) amplitude A ab to be varying slowly on time scales of order one, and supported in a very small tube around such a geodesic.Imposing the usual transverse-traceless gauge conditions, A ab is found to obey a set of transport equations which are solved order by order in ω −1 , see e.g.[50].In addition, k a A ab = g ab A ab = 0, i.e.A ab is a polarization tensor.To be precise, when the expansion for A ab is carried out up to order ω −N , we only get a solution to the linearized Einstein equation up that order, which must be supplemented with a non-perturbative (in ω −1 ) correction of order O(ω −N ).This correction can be chosen such that h ab is still supported near the world-tube of the geodesic in the vicinity of a given Cauchy surface Σ, see [50] for details of this non-trivial construction.
When a derivative ∇ a hits A ab e iωS , the oscillations dominate, corresponding effectively to the substitution ∇ a → ip a ≡ iωk a (D2) familiar from quantum mechanics.Under this substitution (also called the "principal symbol" in the mathematics literature), the operators C n in (30) become form an orthonormal basis of (real) polarization tensors unless p a ∝ n a which we assume for simplicity is not the case.
In particular, we can write the complex amplitude of our WKB approximation as up to a "gauge-transformation"14 , which does not matter since j a is gauge invariant.Making use of such relations in (34) and ( 4), we find that in the high frequency limit, the conserved quantity χ -the flux of Carter current through Σ -on a WKB solution (D1) is given to leading order in ω by The integrand is by construction sharply localized at the point where the geodesic pierces the Cauchy surface Σ, so the integral is basically the integrand at that point.Thus, χ (n) [h] is essentially a power of the Carter constant for a WKB solution [50] in the high frequency limit.An analogous result [50], holds for the canonical energy E[h] (19): In that case we get the energy of a particle in the WKB limit.
Due to the presence of Im(A + Ā× ) coupling the different polarizations, the Carter currents j a (n) are of "zilch"-type in the terminology of [36].

C n ab cd → − 1 2 ζ 4 [ 1 √ 2 1 √ 2
−Q(p)] n e ab (p)e cd (p) ′ , (D3) where Q(p) = K ab p a p b is Carter's constant which is conserved along the null geodesic on account of ∇ (a K bc) = 0, and where e ab (p) and e ab (p) ′ are the polarization tensors e ab (p) = Z ac Z bd p c p d , e ab (p) ′ = Z ′ ac Z ′ bd p c p d (D4) in which the polarization tensor properties follow from the definitions Z = l ∧ m and Z ′ = n ∧ m.Thus, in the high frequency limit, C n is basically the n-th power of Carter's constant, dressed by the tensor structure e ab (p)e cd (p) ′ involving polarization tensors related to the principal null direction l a respectively n a .For the GHP primed symmetry operator C ′ n , the roles of the polarization tensors should be reversed.We now wish to evaluate χ (n) [h] (35) on our WKB solution h ab (D1).Using the definition of Z ab and of Z ab′ , we can see that ēab (p)e ab (p) = (l a p a ) 4 , e ab (p)e ab (p) = 0, ēab (p) ′ e ab (p) ′ = (n a p a ) 4 , e ab (p) ′ e ab (p) ′ = 0, (D5) in addition to the usual properites e ab (p)p a = g ab e ab (p) = 0, and similarly for the primed polarization tensor e ′ ab (p).This means that ϵ + ab (p) = Re e ab (p) ′ /(n c p c ) 2 ϵ × ab (p) = Im e ab (p) ′ /(n c p c ) 2 (D6)