Explicit Triangular Decoupling of the Separated Lichnerowicz Tensor Wave Equation on Schwarzschild into Scalar Regge-Wheeler Equations

We consider the vector and the Lichnerowicz wave equations on the Schwarzschild spacetime, which correspond to the Maxwell and linearized Einstein equations in harmonic gauges (or, respectively, in Lorenz and de Donder gauges). After a complete separation of variables, the radial mode equations form complicated systems of coupled linear ODEs. We outline a precise abstract strategy to decouple these systems into sparse triangular form, where the diagonal blocks consist of spin-$s$ scalar Regge-Wheeler equations (for spins $s=0,1,2$). Building on the example of the vector wave equation, which we have treated previously, we complete a successful implementation of our strategy for the Lichnerowicz wave equation. Our results go a step further than previous more ad-hoc attempts in the literature by presenting a full and maximally simplified final triangular form. These results have important applications to the quantum field theory of and the classical stability analysis of electromagnetic and gravitational perturbations of the Schwarzschild black hole in harmonic gauges.


Introduction
It is textbook knowledge that a standard self-adjoint Sturm-Liouville spectral problem (pφ ) − qφ+ω 2 wφ = 0, with p, q, w > 0, has (i) a real and positive ω 2 -spectrum (or just real ω-spectrum) on a natural function space, with (ii) the generalized eigenfunctions providing a resolution of the identity. In some applications (which we describe below) it is crucial to provably know that properties (i) and (ii) hold (or do not hold) for some equations that do not quite fit into this standard class. Unfortunately, in those cases, the standard Sturm-Liouville theory is no longer applicable and each case has to be studied individually. This work is dedicated to a deep study of a specific ordinary differential equation (ODE) system that naturally appears in the context of linearized gravity on the spherically symmetric Schwarzschild black hole, which is non-standard in the above sense: the coefficients p, q, w are now matrices, either no longer selfadjoint or no longer positive definite, and the spectral parameter ω also appears linearly in the q coefficient. Despite these complications and its superficially very unstructured form, this system turns out to be highly special. Our main result is an explicit transformation of this complicated ODE system into a much simpler (sparse, upper triangular) form, with standard self-adjoint Sturm-Liouville operators on the diagonal. From this simplified form, it becomes obvious that properties (i) and (ii) can be proven to hold by standard theory, while as a byproduct a number of interesting geometric features of the equations of linearized gravity Schwarzschild (under the harmonic gauge condition) are discovered, including an alternative variational formulation, non-trivial generalized symmetries and the existence of Debye potentials (all of these results are concisely collected in Section 2). Our methods are rather special to this ODE system, but having synthesized previous ad-hoc results from the literature into an abstract strategy based on ideas from homological algebra, the geometry and algebra of differential equations, and the theory of rational ODEs, they may also be applicable in other cases. Last but not least, we know of no other tools that can produce any of the same results for the given ODE system.
Since the above topics do not often appear together, let alone applied to problems motivated by linearized gravity, we have aimed for the presentation to be self-contained, though as concise as possible. That balance and the intrinsic complexity of the ODE system under study bear responsibility for the length of this work.
Physical motivation. The study of linear metric perturbations around a Schwarzschild black hole was initiated in [40]. Since then, a rich literature has developed on this topic, including extensions to other linear fields and more general kinds of black holes [10,14,23], with important applications, including the modeling of gravitational wave forms [42], rigorous linear stability of black holes [18,46] (geared towards the nonlinear stability problem), and the study of quantum effects around black holes [23].
Due to diffeomorphism invariance, the linearized Einstein equations require gauge fixing to get well-posed initial value or inhomogeneous problems. Different gauge choices have different merits. The mathematical literature on linear metric perturbations of Schwarzschild has been dominated by the question of classical stability of the classical Cauchy problem. Rather decisive, positive results have already been obtained in a gauge-invariant formalism [21], as well as directly for the metric in a special double null-foliation gauge [17], neither approach using mode decompositions. So why bother studying a different gauge, using mode decomposition, as we do below? The answer is that there is a wider class of interesting questions that we can ask about black hole perturbations, where such information is useful. For instance, applications in quantum field theory (QFT) require an explicit construction of the linear Green function, which is difficult to obtain other than by a mode decomposition. Also, QFT favors specific choices of gauge, with mode stability in those gauges being a necessary condition for the existence of a reasonable vacuum state, while stability in one gauge (favored by the classical Cauchy problem), does not imply stability in a different gauge (favored by QFT). See more precise remarks at the end of Section 2.
In general theoretical treatments of linear metric perturbations, a common choice is the socalled harmonic 1 gauge (also known as de Donder gauge, wave (map/coordinate) gauge, or by analogy with electrodynamics as Lorenz gauge) ∇ ν p µν = 0, with p µν = p µν − 1 2 p λ λ g µν the trace reversed metric perturbation, preferred for its local regularity properties [6] and in applications to quantum field theory [8,12,22]. In this gauge, the linearized Einstein equations take the form of a tensor wave equation, the Lichnerowicz wave equation p µν − 2 R (µ λκ ν) p λκ = 0, which is both covariant and hyperbolic, with both properties of significant theoretical importance. However, due to its algebraic complexity, until recently, relatively little work has been done in this gauge for black hole perturbations, compared to Regge-Wheeler gauge and its variations. Notable early uses of harmonic gauge include [5,9,25]. Following [5], this gauge (there mostly called Lorenz ) has found important applications in the gravitational self-force literature [7]. More recently, the work [26] has used harmonic gauge in the proof of global non-linear stability of Kerr-de Sitter black holes, while in [29,30] 2 and [27,28] the vector field method was used, directly in real-space, to study stability and decay of perturbations of Schwarzschild, with [31] an addendum that implicitly considers an upper triangular decoupling of the Lichnerowicz wave equation under the transverse-traceless conditions, ∇ ν p µν and p λ λ = 0, without referring to or comparing with the fully diagonal decoupling of the same system by [9].
Mathematical problem. To date, despite the strong motivations listed above and the known complete separability of the Lichnerowicz d'Alembertian, there still does not exist an explicit mode-level construction of a harmonic gauge Green function for metric perturbations on Schwarzschild. A necessary step in such a construction would be a determination of the spectrum of the separated radial mode equation and a proof of completeness of its generalized eigenfunctions. In this work, we take a significant step in that direction. Namely, we perform a highly non-trivial simplification that allows this spectral problem to be studied by standard Sturm-Liouville theory.
The main obstacle is that the separated radial mode equation on Schwarzschild is a rather complicated system of ODEs. While it can be put into matrix Sturm-Liouville form, it presents a highly non-standard spectral problem. For one, it is naturally self-adjoint only with respect to an indefinite functional inner product. So either it becomes a self-adjoint spectral problem on a Krein space (analogous to a Hilbert space, but with an indefinite inner product), or a non-selfadjoint problem on a Hilbert space (where we artificially positivise the natural inner product). Further, it presents a quadratic eigenvalue problem in the frequency ω, rather than a regular eigenvalue problem with respect to ω 2 . In either case, the standard spectral theory of self-adjoint operators on Hilbert space becomes totally inapplicable, leaving us without any obvious way to prove (or disprove) the reality of the spectrum (absence of modes exponentially unstable in time) or even to check the completeness of (generalized) eigenfunctions. The indefiniteness of the functional inner product is ultimately due to the Lorentzian signature of the metric, meaning that this complication does not occur in the analogous problem in Riemannian geometry [25].
Methodology. In this work, building on the strategy outlined and successfully implemented in [32,34] for the simpler example of the vector wave equation v µ = 0, we prove that the harmonic gauge radial mode equation can be decoupled into a triangular system, where the diagonal blocks are the spin-s Regge-Wheeler equations, with s = 0, 1 or 2. These Regge-Wheeler equations are then of standard scalar Sturm-Liouville form, with very well-understood spectral properties [14,20]. Such a decoupling essentially reduces the non-standard spectral problem of the radial mode equation to the standard and well-understood spectral problem of the Regge-Wheeler equation. The details of such a spectral analysis are left for to future work.
To our knowledge, prior to [32], which was our warm-up for this work, such a triangular decoupling has never been explicitly discussed in the literature on linearized gravity. Indirect hints of it have previously appeared only in [9] (Remark 5.10) and independently in [41] (though only for the vector wave equation in the latter). Very recently, [27,28] has used some of the resulting formulas from [9]. Unfortunately, though pioneering, the original works [9,41] were based on rather ad-hoc, extensive, explicit computations and the full details of how the original radial mode equations transform into the decoupled form and back are not easy to understand from these references. We have strived to distill their overall strategy into conceptually clear terms 3 and to understand why certain key ad-hoc steps were successful, as recorded and implemented in [32,34] for the vector wave equation. In the end, our results also go a step further than [9,41], by reducing the resulting triangular form as much as possible and proving that no further reduction is possible (in the context of rational ODEs).
Contents overview. In Section 3, we start by reviewing a precise notion of morphism between differential equations, as well as of equivalence up to homotopy between morphisms or between equations (or isomorphism), where we synthesize some elementary notions from homological algebra [49], D-modules [43,Section 10.5] and the categorical approach to differential equations [35,Section VII.5]. Essentially, a morphism is a differential operator mapping solutions to solutions, but also equipped with extra data preserving the structure of the equations. 4 Then, we review how a hierarchical separation a differential equations into gauge modes, gauge invariant modes and constraint violating modes can lead to an equivalence with an equation in block upper triangular form. Our general decoupling strategy is to use this step recursively until a full triangular form is reached. The second step in the strategy is to simplify this initial triangular form further by transforming as many as possible off-diagonal components to zero or to some non-zero canonical form by the methods to be presented in Section 4. It is remarkable that all these equations and transformations between them involve only differential operators with rational coefficients. This review is a condensed version of the more detailed discussion given in [32].
In Section 4, we review some tools from the study of rational solutions of rational ODEs (ODEs with rational coefficients) and apply them to the spin-s Regge-Wheeler equations. Namely, we reduce the problem of deciding when a triangularly coupled Regge-Wheeler system can be made diagonal to a finite dimensional linear problem that can be easily solved with computer algebra. Incidentally, for spins s = 0, 1, 2 that are of interest to us, we prove for the first time that the only non-trivial rational morphisms (in the sense of Section 3) between spin-s Regge-Wheeler equations are the identity morphisms for equal s. (In other words, there do not exist rational spin raising and lowering operators [44] at least between radial modes of spins s = 0, 1, 2, and also likely any other spins.) These methods were first applied in [34]. In this work, they are further improved and refined by new results on the characterization of the cokernel of a rational ODE. Finally, in Section 4.1, we also discuss the Chandrasekhar transformation as an equivalence between the s = 2 Regge-Wheeler and Zerilli equations, with special attention to its dependence on the frequency ω.
In Section 5, we first review the complete separation of variables of the scalar, vector and Lichnerowicz wave equations into spherical and time harmonic modes. Our notation and conventions are compared with the literature in Appendix A. Then we present the resulting radial mode equations, related differential operators and identities between them that are needed to implement the decoupling strategy from Section 3. In Section 5.2 we recall from [32] and slightly improve the explicit full triangular decoupling for the radial vector wave equation. In Section 5.3 we give for the first time a explicit full triangular decoupling for the radial Lichnerowicz wave equation, which follows from the decoupling strategy of Section 3 and builds on the results of Section 5.2. Since some of the explicit formulas needed for the even sector of the Lichnerowicz equation are quite lengthy, they are relegated to Appendix B. In both the vector and Lichnerowicz wave equation cases, we improve on the original results of [9] by giving a more complete and simpler triangular form.
In Section 2 we formulate our main results as a concise theorem. Then we state several immediate qualitative corollaries, which indicate important applications of the explicit knowledge of our decoupling formulas. Finally, these applications and other potential further developments are discussed in Section 6.

Results and applications
Here we concisely state our main result (Theorem 2.1) and then discuss several important applications in the form of corollaries that follow straightforwardly from the main theorem. The commentary given along the way can be seen as a guide to reading the rest of the paper. There are of course further potential applications that require less straightforward consequences of Theorem 2.1. They will be pursued in future works.
We must start by introducing the minimum of the notations and definitions needed to state our results. We may be informal here, but will give pointers to the formal details in the body of the paper.
The central objects that enter the hypothesis of our main theorem are specific linear ordinary second order differential operators (equivalently, differential equations). They have matrix coefficients (of given size) whose entries are rational functions depending on parameters, which are M (mass), l (angular momentum quantum number ) and ω (frequency). They are 0 (1 × 1), and 2 e (7 × 7), where notation reflects that they are related to the mode separation of the tensor wave operator = ∇ µ ∇ ν on the static spherically symmetric Schwarzschild black hole (full details in Section 5) for different spins (s = 0, 1, 2) and belong either to the even or odd sectors under antipodal reflection of the spherical orbits of symmetry. Each operator is defined on the interval r ∈ (2M, ∞), with a regular singularity at r = 2M and an irregular one at r = ∞. The dependence on all the parameters is polynomial, with ω considered the spectral parameter. Since the explicit form of some of these operators is rather complicated, we write them only in the sections where they are derived, namely Sections 5.1, 5.2 and 5.3. However, each of them has the following generic form (see Section 3.1 for our conventions regarding differential operators): where P (r), Q(r), iA(r) and W (r) are self-adjoint matrices (supposing that M and have real values) of appropriate size, meaning that E * = E is formally self-adjoint (Section 3.3). As discussed in the Introduction, the matrix W (r) is of indefinite signature (except for the 1 × 1 cases), meaning that it cannot be used to define a Hilbert space (only a Krein space) on which E could define a self-adjoint unbounded operator. Even if that were possible, the dependence of E on both ω and ω 2 implies that we would have needed to consider it as a quadratic eigenvalue problem, where verifying the reality of the ω-spectrum would have already been non-trivial.
In our main result, each of these complicated differential operators will be transformed to simplified upper diagonal form, which involves the family of spin-s Regge-Wheeler operators (Section 4) where for convenience The following convenient combinations of the parameters often accompany the Regge-Wheeler operators: B l := l(l + 1), When s = 2, the alternative Zerilli operator is often used In Section 4.1 we study in detail the equivalence of D 2 and D + 2 , and explain why we prefer not to use the Zerilli operator. The Regge-Wheeler operators are in fact of standard Sturm-Liouville type [20].

Main result
To state Theorem 2.1 below, besides the notations we have just introduced above, one also needs the notion of equivalence between two differential equations. Informally, for our purposes, two differential equations are equivalent when there exist differential operators that map solutions to solutions, in both directions, and are moreover mutually inverse on the solutions. More precisely, our notion of equivalence involves the existence of a number of auxiliary differential operators satisfying some compositional identities, which are conveniently summarized in an equivalence diagram such as (2.1) below. The full details of the definition and the rationale for it are given in Section 3.1.
To prove that such an equivalence exists, it is of course sufficient to write down all the differential operators and verify the required identities between them. That is precisely our proof strategy, with full details presented in Section 5, with Sections 3 and 4 building up preliminary material for it.
The main content of the theorem is that each of our differential equations of interest is equivalent to a much simpler upper triangular form. To appreciate the magnitude of the simplification, it is sufficient to compare the upper triangular forms in (2.2) with the original forms of 0 (5.1c), 1 o (5.3d), 1 e (5.3c), 2 o (5.6d) and 2 e (5.6c).
Theorem 2.1. Let E be one of the 0 , 1 o , 1 e , 2 o or 2 e operators, representing the systems of radial mode equations for the scalar (Section 5.1), vector (Section 5.2) or Lichnerowicz (Section 5.3) wave operators on the Schwarzschild black hole of mass M , mode separated as described in Section 5, in either the odd or even sector. Then the following is true: (i) With respect to the Schwarzschild radial coordinate r, E is a rational ODE, with singular points only at r = 0, 2M, ∞. E also depends polynomially on the frequency ω and angular momentum B l = l(l + 1) spectral parameters, as well as the mass parameter M .
(ii) There exists a rational ODE systemẼ in upper triangular form and an equivalence diagram (Section 3) such that the diagonal ofẼ consists of Regge-Wheeler operators D s (up to constant, ωand B l -dependent factors), while the equivalence maps are rational differential operators of order at most 1 and poles only at r = 0, 2M .
(iii) The equivalence k, k ,k,k and homotopy h,h operators in (2.1) depend rationally on ω and B l , with poles only at ω = 0 and B l = 0, 2, respectively, with one exception. The exception is the case E = 2 e , where some of these operators have additional poles at the algebraically special frequencies ω = ±i A l 12M . In that case, using the Chandrasekhar transformation (Section 4.1) to replace the Regge-Wheeler operator D 2 by the Zerilli operator D + 2 inẼ removes the poles at the algebraically special frequencies at the price of introducing factors of (B l −2+3f 1 ) in some of the denominators in the equivalence diagram.
(iv) The equivalence diagram (2.1) has finite limit as M → 0.Ẽ becomes diagonal in that limit.
(vi) The rows and columns ofẼ can be permuted such that it remains upper triangular and block diagonal with respect to blocks of equal s of the Regge-Wheeler operators D s . The most general rational automorphism ofẼ takes the formẼA = AẼ, where A is a constant matrix that is block diagonal with respect to equal spin blocks and is upper triangular within each such block, with the single exception of the s = 0 block when E = 2 e .
In each of these cases, the upper triangular Regge-Wheeler system has the following form: (a) when E = 0 , theñ Proof . Section 5 gives a detailed presentation of the explicit formulas for the operators in the equivalence diagram (2.1), thus proving (ii), for each case being considered: (5.2) =⇒ (a), (5.5) =⇒ (b), (5.2.2) =⇒ (c), (5.10) =⇒ (d), (5.12) =⇒ (e). The rest of the theorem is simply a summary, collected here in a way convenient for future reference, of properties of these operators that can be gleaned from their explicit formulas.

Spectral problem
Proof . The corollary follows from a direct generalization to larger operator matrices of Remark 4.7, about the relative boundedness of off-diagonal matrix elements with respect to the Regge-Wheeler operators on the diagonal, which here follows from Theorem 2.1(v).
As point out in Remark 4.7, the novel results of Section 4 concerning the cokernel of the differential operator in (4.4) are crucial for guaranteeing the possibility of choosing the offdiagonal elements to be relatively bounded while simplifying an upper triangular Regge-Wheeler system. In practice, we have found that at the initial stages of this simplification, which is part of the general strategy of Section 3, the off-diagonal elements were generally not relatively bounded. Using only the results from [34] about only the solutions of (4.4), or other off-the-shelf tools from the literature, it would have been impossible to decide a priori whether those off-diagonal elements that cannot be chosen to be zero could be chosen to be relatively bounded.
From Corollary 2.2 one can essentially conclude that all the properties of the resolventẼ −1 (like its dependence on ω), and by equivalence hence also of E −1 , can be deduced from the properties of the resolvent D −1 s , which has been very well studied. In particular, we can conclude that the ω-spectrum of each E in Theorem 2.1 is purely real, as it is for D s [20], once the relevant function spaces are chosen to respect the equivalences that we have constructed. This information is an important starting point for the integral representation of the solutions of the corresponding Cauchy problem and of the large time asymptotics of such solutions [20]. It is also crucial for the explicit construction of propagators in quantum field theory [8,12,22].
We now make the above reasoning slightly more precise. When E is one of 0 , 1 o , 1 e , 2 o or 2 e , both it and its decoupled triangular formẼ are rational ODEs with a regular singular point at r = 2M and an irregular singular point at r = ∞. Hence, for each system, and at each singular point, using the methods of asymptotic analysis of ODEs [48], each solution can be uniquely identified by its asymptotics. Being bijective on solutions, the equivalence morphisms k,k of Theorem 2.1 must hence be bijective on these asymptotics. Moreover, since these morphisms are differential operators, they can be applied to the asymptotics directly, without knowledge of the full solution. This means that if we know the connection coefficients between the singular points at r = 2M and r = ∞ (essentially the transmission and reflection coefficients for incoming and outgoing waves) at frequency ω forẼ, the equivalence morphisms transfer them to the corresponding connection coefficients for E. This argument of course works only when the equivalence morphisms are well-defined, which excludes l = 0, 1 in some cases and ω = 0, as well as the algebraically special ω = ±i A l 12M in the case E = 2 e . But these excluded cases can be studied individually. For instance, for the algebraically special frequencies, it is sufficient to replace the Regge-Wheeler operator D 2 with the Zerilli operator D + 2 in˜ 2 e , as pointed out in Theorem 2.1(iii). For ω = 0 and l = 0, 1 one can get more information by studying the coefficients of the Laurent expansion of the equivalence morphisms. Thus, we arrive at Corollary 2.3. For ω = 0 and l ≥ 2, the equivalence morphisms k,k from Theorem 2.1 are bijective on solution asymptotics of each E,Ẽ pair near r = 2M and r = ∞. Hence, after bijectively identifying their solutions, the connection coefficients for E andẼ are equal.

Symmetries and potentials
Each of Section 5.1, 5.2 and 5.3 starts with separation of variables on the static spherically symmetric Schwarzschild black hole background, which turns spacetime partial differential operators into ordinary differential operators acting on radial mode functions. The time derivative i∂ t is replaced by the frequency parameter ω, while the spherical Laplacian D A D A is replaced by the angular momentum parameter −B l (shifted by a constant depending on the tensor rank of its argument), among other substitutions. Not every operator on radial mode functions comes from such a separation of variables, but it is straightforward to work out the rules for recognizing those operators that do. Thus, directly from the statement of Theorem 2.1, we can also draw conclusions about spacetime partial differential operators.
The first conclusion is about symmetries. A symmetry of a differential equation E[u] = 0 is a pair of differential operators S, S such that ES = S E (Section 3.1). Such symmetries are extremely important in understanding the structure of a differential equation [39].
The second conclusion is about a variational formulation [39]. A formally self-adjoint differential operator E = E * (Section 3.3) is the Euler-Lagrange equation of a variational principle u, E[u] dx. Conversely, the Euler-Lagrange equation of a variational principle is always formally self-adjoint. For linear differential equations the two notions are essentially equivalent. Having a variational formulation leads to the possibility of using Noether's theorem to convert symmetries to conservation laws. Even more interesting is when the same equation has more than one variational principle. This is the case for a self-adjoint equation   (b) By similar arguments, the morphisms k, k andk,k from (2.1) can be promoted to spacetime differential operators, possibly after being multiplied by a constant coefficient polynomial in B l and ω of sufficiently high degree. These spacetime differential operators will effect a morphism between the original harmonic gauge tensor wave equation and the triangular system of Regge-Wheeler wave equations from (a).
(c) The systems of Regge-Wheeler wave equations from (a) are self-adjoint on spacetime after composition with the Σ operator from Theorem 2.1(vii). Hence they possess a variational formulation.
Remark 2.5. Note that Corollary 2.4 cannot hold for the Zerilli form of˜ 2 e mentioned in Theorem 2.1(iii), because no spacetime differential operator has a radial mode decomposition with B l -dependent factors like (B l − 2 + 3f 1 ) in the denominator.
A Debye potential is a differential operator that maps solutions of a simple auxiliary PDE equation (typically a scalar wave equation with a potential) to solutions of another PDE of interest (typically a more complicated tensor wave equation), with the usual additional requirement that the image solutions are non-trivial (typically not pure gauge, in an appropriate sense). In the terminology of Section 3, a Debye potential a special kind of morphism between PDEs. Debye potentials for Maxwell and linearized Einstein equations, on both flat and some curved backgrounds (including black hole backgrounds), have been known for some time, but have also gained more attention recently [3,4,45,47]. In the case of black hole backgrounds, the Debye potentials for Maxwell and linearized Einstein equations typically give solutions in some version of radiation gauge [45,47], which is adapted to the algebraically special nature of these backgrounds, but may not be optimal for other purposes. The recent work [30] produced a Debye potential in a (non-local) generalized harmonic gauge. Of course, a solution in one gauge can be (at least locally) transformed into any other gauge, but it is by no means obvious when such a transformation can be done by a differential operator. With that in mind, we draw attention to the next corollary (first implicitly obtained in [9]): Proof . From Theorem 2.1, it is clear that the spin s = 1 components decouple from the other components in˜ 1 o and˜ 1 e , so they satisfy independent Regge-Wheeler equations. Hence, the corresponding columns of thek 1o andk 1e operators, once promoted to spacetime differential operators as in Corollary 2.4, are Debye potentials for the vector wave equation 1 . However, the composition properties of k 1o , k 1e with T 1o , T 2e (as computed in Sections 5.2.1 and 5.2.2) imply that the image of the Debye potential is annihilated by the spacetime operator T 1 . In other words, the Debye potentials actually generate solutions of Maxwell equations in harmonic gauge. For linearized Einstein equations, the discussion is completely analogous, starting with the Debye potentials for the Lichnerowicz wave equation generated by solutions of Regge-Wheeler equations of spin s = 2. In that case, the composition properties with T 2 and tr operators (as computed in Sections 5.3.1 and 5.3.2) imply that the image of these Debye potentials is both harmonic and trace free.
3 Formal properties of differential equations and operators

Equations and morphisms
For the purposes of this work, a (partial ) differential operator, say E, is always linear, with smooth coefficients, which we will write as E[u], where u is a possibly vector valued function. Later on we will even specialize to ordinary differential operators with rational coefficients, but all abstract statements expressed in terms of operator identities or commutative diagrams will be valid at the greater level of generality. Compositions of differential operators will be written E 1 • E 2 or just E 1 E 2 when no confusion could arise.
An operator E can always be thought of as a matrix of scalar differential operators acting on the components of u. By a scalar differential operator, we mean an operator that takes possibly vector valued functions into scalar valued functions (corresponding to a matrix with a single row). Scalars are real or complex numbers (for the sake of generality we will generally allow complex scalars). Differential operators could be of any order, including order zero, which just corresponds to multiplication by some matrix valued function. While the operators can be thought of as acting on smooth (vector valued) functions, we will be mostly concerned with composition identities among operators, and so we will not bother specifying precisely the domain or codomain of each operator. This information can always be deduced from the context (e.g., the size of the matrix). Also, for the abstract discussion below, it is also not necessary to fix the number of independent variables, but in the rest of this work we will be mostly concerned with applications to ordinary differential operators (i.e., acting on functions of a single independent variable).
Quick example: the operator E = r −2 ∂ r r 2 should be interpreted as E[u] = r −2 ∂ r r 2 u = ∂ r u + 2u/r, so alternatively it can be rewritten as E = ∂ r + 2/r.
Below, we state some basic definitions and results concerning differential equations and morphisms between them, which essentially correspond to differential operators that map solutions to solutions. The presentation is logically self-contained and does not extend far beyond what is needed in the rest of the paper. We generalize and streamline slightly the presentation previously given in [32,Section 2], which can be consulted for a more pedagogical exposition. For proper context, these ideas can be seen as simple special cases of concepts coming from the more general frameworks of D-modules [43, Section 10.5], the category of differential equations [35, Section VII.5] or homological algebra [49]. We will not delve into the precise connection with these larger frameworks, but it is useful to mention that all arrows/morphisms need to be reversed when our statements are interpreted in the language of D-modules.
Central objects of our attention are differential equations, like E[u] = 0. Sometimes we will refer just to the operator E as the equation itself and vice versa. This should not lead to any confusion. Technically, it becomes convenient to consider not just an equation E (0) , but also a set of Noether (or compatibility) identities that comes with it, namely an operator E such that E • E (0) = 0. One can equally consider higher stage Noether identities, for example an operator E such that E • E = 0, etc. Such a sequence of operators E (0) , E , E , . . . , E (n) , is also called a complex (a term from homological algebra [49]), but we might as well still refer to it as a differential equation, albeit equipped with the extra data of Noether and higher stage Noether identities. We get back to a notion of solutions via the idea of cohomology H (n) (E) := ker E (n) / im E (n−1) of the complex, which may be non-vanishing in any location and can be interpreted as the space of solutions of coincides with the usual notion of solutions. Note however that the precise spaces of solutions and cohomologies depend on the choice of the function space on which our differential operators act. The various transformations and constructions involving complexes that we introduce below are designed to preserve cohomologies, hence also solutions, once the function spaces have been chosen, which is the main reason that we abstract from solutions themselves and work only at the level of complexes of differential operators.
The next most important concept is that of a morphism between two differential equations. Referring again to terminology from homological algebra, a morphism for us will be a cochain map between the complexes corresponding to the two equations. In more detail, given two complexes E (0) , E , . . . and F (0) , F , . . ., a cochain map consists of differential operators k (0) , k , . . . with domains and codomains prescribed by the diagram in Figure 1(a), which also needs to be commutative (any sequence of operators starting and ending at the same point must compose to the same thing). We have denoted by bullets the appropriate function spaces, which are less important to us than the structure and properties of the differential operators acting between them. When needed, e.g., for the purposes of defining the cohomology or for defining a cochain map between complexes of different lengths, any complex can be freely extended by zero maps in both directions. This explains the extra zero maps appended and prepended to the complexes in Figure 1(a). But since this can be done implicitly, we will not show such extensions and rather use more compact diagrams like the rest of Figure 1. When a complex consists of a single operator E (E = E = · · · = 0, and we drop the now unnecessary superscript (0) ), the cohomology has the obvious identification H(E) = ker E/ im 0 = ker E and H (E) = ker 0/ im E = coker E. Once again, this highlights the interpretation of a complex as a generalization of a linear differential equation and of its cohomology as generalizing the space of solutions. It is not hard to check that the defining properties of a morphism (cochain map) k, k , . . . imply that it descends to well-defined maps in cohomology k (n) : In the simplest example of a single operator, this means that k : ker E → ker F maps solutions of E[u] = 0 to solutions of F [v] = 0 (also k : coker E → coker F in this case). An automorphism is a morphism from a given equation to itself, in which case the operator k mapping solutions of the equation to solutions is also known as a symmetry operator. A morphism is induced by a homotopy if there exists a homotopy, that is, a sequence of operators h (0) , h , . . . (again, extended by zero operators in either direction, as necessary) fitting into the morphism diagram in Figure 1(b), such that each horizontal arrow satisfies the formula which is indeed a morphism. It is a straightforward exercise to check that a morphism induced by a homotopy always descends to the zero map in cohomology. Two morphisms k induced by a homotopy. Two equations are equivalent up to homotopy if there exists morphisms between them that are mutually inverse up to homotopy (their compositions are equivalent to the identity morphisms up to homotopy). All the relevant operators can be exhibited in an equivalence diagram, like in Figure 1(c), where the solid arrows commute and the homotopy corrections enter the identities Equivalence up to homotopy, both for operators and equations, of course defines an equivalence relation.
Hence, as promised, in this case, a morphism induces a map from solutions to solutions, which is a familiar notion in PDE theory. Similarly, a symmetry (or automorphism) in our sense induces a map from solutions to solutions of a given equation. Extending our earlier example by an inverse morphismk up to homotopy implies a bijection between the solutions of But even more than that, an equivalence allows us to transfer the solutions of inhomogeneous problems between equations, , which can be checked by direct calculation. Essentially, our notion of morphism is the same as the usual PDE notion of a map between solution spaces induced by a differential operator, but equipped with extra structure. For some equations that are not in some sense pathological, the extra operators k (n) may be reconstructed (perhaps up to homotopy equivalence) from the knowledge of a single operator k (0) . Then, the above two notions actually coincide. See [32] for a bit more discussion on this point. Considering subtle cases when differences between the two notions actually arise is beyond the scope of this work. We will restrict ourselves to using the more structured notion of morphism, which is happens to be most convenient for our work.

Triangular decoupling strategy
In this section, we outline our general decoupling strategy. We basically summarize the more pedagogical discussion from [32]. There are two small differences. We use slightly different notation, more adapted to the later needs of Section 5. Also, we take advantage of having defined, in the preceding section, morphisms for complexes of possibly more than one differential operators. Part of the discussion in [32] was unnecessarily awkward because we insisted on using morphisms only with single differential operators (namely, the homotopy corrections on the right-hand side of (3.1) did not have the same uniform structure because of the truncation).
The input is a differential equation E[u] = 0, together with the following morphisms The end result is an equation in block upper triangular form with D D , D Φ and D T on the diagonal, where both the diagonal blocks and the off-diagonal ones have also been simplified as much as possible, keeping upper triangular structure, as shown in (3.4). For our purposes, we may refer to the degrees of freedom captured by the D operator as pure gauge modes, those by the Φ operator as gauge invariant modes, and those by the T operator as constraint violating modes. This terminology will be justified in how the strategy will be applied in Section 5 to Maxwell and linearized Einstein equations in harmonic gauge (see also the detailed discussion in [32]). At the abstract level, we have simply introduced convenient names to the degrees of freedom corresponding to the 3 sectors of the eventual block upper triangular decouplingĒ in (3.4). The first non-trivial step is to find an equivalence diagram, illustrated in (3.3a), between D D (the residual gauge equation) and the joint system E[u] = 0, Φ[u] = 0 (vanishing of gauge invariant modes), T [u] = 0 (gauge fixing condition). That is, we must complement the input differential operators already introduced in (3.2) by a bunch of new operators that are defined by their relation to the following diagram being an equivalence up to homotopy: supplemented by two more operatorsH Φ ,H T satisfying the relation If such a diagram does not exist, than our strategy cannot proceed directly. In this work, we will only consider cases where this step of the strategy succeeds.
There are a few potential points of failure, which should all be checked. First, the gauge-invariant and constraint violating modes should be compatible and independent from the pure gauge modes, namely Φ T • D ∼ 0 up to homotopy. This condition is captured by the following subset of identities from (3.3a): Strictly speaking, to justify the adjective gauge-invariant, we should have Φ • D = 0 (H T = 0), but here we are allowing Φ to mix with E and T , that is, terms that vanish on-shell or by the gauge-fixing condition. In practive, within this work we will always be able to make an initial choice of Φ such that H T = 0, and then update it to fit our needs while maintaining the above relations. Note that we have chosen to represent the E-Φ-T equation in diagram (3.3a) by a complex of length 2, with the second operator representing Noether identities for the first one. If we do not choose this complex correctly, then the construction of the equivalence diagram may be obstructed, for instance when some Noether identity for the first operator is not a consequence of the rows of the second operator. When such obstructions are absent, the existence of operators H Φ ,H T in (3.3b) actually follows from the relations already exhibited by (3.3a). In this work, we will only encounter cases where both sides of the equivalence diagram (3.3a) are correctly chosen, but to save space we simply postulate the factorizations (3.3b) separately.
Finally, another point of failure could be that the T [u] = 0 and Φ[u] = 0 conditions are simply not strong enough to reduce the number of remaining degrees of freedom to equal that of the D D equation. Again, in this work we will only encounter cases where this step succeeds. Next, using only operators that have already been introduced above, together with their defining properties, we can construct two more equivalence diagrams that are concatenated and explicitly shown in Figure 2. Composing the morphisms in Figure 2, in both directions, we obtain the much more compact equivalence with an upper triangular form, which we denote byẼ, Conversely, knowing the final equivalence diagram (3.4) of E with the decoupled formẼ, simply by placing individual operators in the appropriate places, we can reproduce the equivalence (3.3a) of the E-Φ-T system with D D and the relations (3.3b) The triangular systemẼ may still be rather complicated. So the strategy proceeds with two kinds of simplifications: simplifying the diagonal blocks, and simplifying the off-diagonal blocks (cf. [32,). The first kind is summarized in the following Lemma 3.3. If there exists an equivalence diagram between equations E 2 andĒ 2 , then it implies the following equivalence diagram for an upper triangular system, where E 2 is the middle block: , provided the extra identity kh =hk also holds (though it is only needed to verify the form of the barred homotopy operator on the right).
The proof is by direct calculation. Though, we should note that when E 2 has no nontrivial Noether identities (for example, in the ODE case, when it can be solved for all highest derivatives), then the extra identity on the homotopy operators is automatic. Indeed, then and since E 2 has no non-trivial Noether identities, we must have kh −hk = 0.
Note that Lemma 3.3 remains valid when either of the untransformed diagonal blocks E 1 or E 3 is zero dimensional (a 0 × 0 matrix, or just absent, in other words).
The second kind of simplification can be achieved by applying the next Lemma 3.4. For any operators δ and δ the following equivalence diagram holds: The proof is again by direct calculation. This Lemma implies that by a clever choice of δ and δ it might be possible to set∆ = 0, thus decoupling the two diagonal blocks, or at least set it some preferred canonical form. For instance, in the case when E 1 and E 2 are ODEs that can be solved for their higher derivatives, say respectively of orders n 1 and n 2 , then by an obvious choice of δ and δ the off-diagonal term∆ can be reduced to order min(n 1 , n 2 ).
Our decoupling strategy proceeds by applying the last two lemmas recursively, repartitioning our triangular equation into blocks as necessary, until no further simplification is possible or practical. In the successful applications of this strategy in Section 5 we stop when the systemẼ has been reduced to fully upper triangular form, where each operator on the diagonal is a spin-s Regge-Wheeler operator (Section 4), while off-diagonal terms are mostly zero, with those that are non-vanishing expressed in terms of a small number of canonical representatives. We dedicate the following section (Section 4) to giving an effective method of deciding when Lemma 3.4 can be used to cancel (or at least "maximally simplify") specific off-diagonal terms in a triangular system with Regge-Wheeler operators on the diagonal.

Formal adjoints
Given a local sesquilinear scalar paring −, − = dx −, − x , with −, − x the corresponding pointwise non-degenerate pairing, the formal adjoint E * of a differential operator E is defined by for smooth compactly supported test functions u and v. If an operator satisfies E * = E, then it is (formally) self-adjoint. As defined, formal adjoints always exist and are unique. In local coordinates they may be computed by the standard rules On a manifold with a metric, for tensor fields v and u, we define the paring v, u to be the integral with respect to the metric volume form of the metric contraction of v * (the complex conjugate of v) and u. We use this convention in Section 5, both for the spacetime metric 4 g µν on M = R 2 × S 2 and for the metric g ab on the radio-temporal R 2 -factor.
In Sections 4 and 5, we will also be dealing with ordinary differential operators with respect to the Schwarzschild radial variable r. In that context, we take the pairing v, u = ∞ 2M v * u dr, where u is a column vector valued function and v * is the conjugate transpose of a column vector valued function.
When the operators depend on spectral parameters, like ω and l, they are considered to be real, so that ω * = ω and l * = l. However, that does not prevent us from considering complex values of these parameters once all the adjoints have been computed.

Morphisms and extensions of Regge-Wheeler equations
In this section, we present the theory needed for simplifying the off-diagonal elements of rational ODE systems in upper triangular form. We first concisely summarize some necessary notions and results from [34] (where a more detailed and pedagogical treatment can be found) and then extend them. In this section, we will work exclusively with rational differential operators.
The following functions will appear in the Schwarzschild metric and its derivatives: Let also l, s and ω be real (or more generally complex) constants, with B l = l(l + 1) a convenient notation. Throughout this work, we use only ∂ r to denote differentiation, while primes are used only as decorations (to distinguish two objects from each other, like to operators k and k ). Define the (generalized) spin-s Regge-Wheeler operator with mass parameter M , angular momentum quantum number l and frequency ω by (following the notation conventions of Section 3.1) We will usually assume that l = 0, 1, 2, . . ., that s is a non-negative integer (in fact, only s = 0, 1, 2 will concern us), and that ω = 0. Of course, D s has rational coefficients in r and all of these parameters. Consider the following equivalence of upper triangular Regge-Wheeler systems: As discussed in Section 3, this system is reducible to diagonal form by the above equivalence (Lemma 3.4) if and only if the following equation is satisfied by some δ and ε with rational coefficients: For simplicity of notation, where it does not introduce any ambiguities, we set∆ = 0 below. By [34, Theorem 3.1], without loss of generality, we can consider this problem restricted to first order δ, and ∆, with δ uniquely determining ε and vice versa. For instance, if ∆ is of order higher than one, we can always reduce the order by writing ∆ = ∆ 1 + ∆ • D s 1 with ∆ 1 of order at most one, and absorb ∆ into ε (or analogously with δ).
Let us parametrize these operators as follows: where δ ± , ∆ ± are arbitrary rational functions and we have used the anti-commutator notation {g, h} = g•h+h•g for differential operators, namely {g 0 , ∂ r } = 2g 0 ∂ r +(∂ r g 0 ) for any function g 0 . Plugging this parametrization into (4.2) and comparing coefficients of powers of ∂ r , we find the equivalent ODE system The parametrization we have used for ∆ and δ is different from the one we used previously in [34], but it is more convenient for some simple technical reasons. If ∆ is parametrized by (∆ + , ∆ − ), then ∆ * is parametrized by (∆ + , −∆ − ). Moreover, given a solution δ of (4.4) parametrized by (δ + , δ − ), if we simultaneously replace ∆ by ∆ * and exchange s 0 with s 1 , then the solution to the new problem is parametrized in the same way but by (δ * + , −δ * − ). The last property is due to the formal self-adjointness of D * s = D s . With our new parametrization, it is also easy to check that the operator D(s 0 , s 1 ) in (4.4) is itself formally self-adjoint, and setting s 0 = s 1 puts it into diagonal form. The version of the decoupling equation used in [34] did not manifestly exhibit these properties.
Furthermore, the ODE coefficients in (4.4) are not just rational functions, but even more conveniently they are Laurent polynomials in r. The rationality of the coefficients means that we can use the results from [34] to systematically solve (4.4), or check that it has no solution. Actually, we will need to extend the results of [34] in an important way. When ∆ is such that no solution δ exists, we would like to find some canonical choice of∆ such that replacing the right-hand side by ∆ −∆ does make the equation solvable. That is, we would like to be able to pick some convenient representatives of the equivalence classes in the cokernel of D(s 0 , s 1 ) in (4.4). In fact, exploiting the fact that its coefficients are actually Laurent polynomials, we will show in Theorem 4.5 that fixing the locations of the poles of ∆, the cokernel of the operator in (4.4) is finite dimensional and a complete set of representatives can be explicitly selected.
First, let us recall some useful terminology and notions from [34], where further details and more references to the literature on this topic can be found. We will be working with 5 complex valued rational functions in r, C(r), polynomials in r and r −1 , C[r] and C r −1 , formal power series in r and r −1 , C[[r]] and C r −1 , as well as formal Laurent series in r. Let us distinguish between unbounded Laurent series C r, r −1 , bounded (from below) Laurent series C r −1 [[r]], bounded from above Laurent series C[r] r −1 and Laurent polynomials C[r, r −1 ]. Of course, we could also consider Laurent series centered at some other r = ρ = 0, but for convenience of notation whenever possible we will stick with ρ = 0.
By convention, matrix ordinary differential operators E[φ] act on column vectors φ with components either in C(r), or belonging to one of these classes of Laurent series. When E[φ] has bounded Laurent series coefficients, its action on bounded Laurent series produces bounded Laurent series (from above or below, respectively). Rational coefficients become bounded Laurent series upon appropriate expansion. For the action of E[φ] to be well-defined on unbounded Laurent series, its coefficients must be Laurent polynomials. In each case, let us call such an E a compatible differential operator.
We call a Laurent polynomial matrix S = S(r) Laurent unimodular if its inverse S(r) −1 is also Laurent polynomial (equivalently, the determinant of S is proportional to a single power of r and non-vanishing). We say that Laurent unimodular matrices S and T are respectively source and target leading multipliers of E when, after expanding all rational coefficients as bounded Laurent series, we have E S(r)φ n r n = T (r) E n φ n r n + r n O(r) , with the components of O(r) all in rC [[r]] and E n an r-independent matrix that is invertible for almost all n (i.e., all but finitely many). We call E n the leading characteristic matrix of E with respect to the given multipliers. Similarly, we say that S and T are respectively source and target trailing multipliers of E when, after expanding all rational coefficients as bounded from above Laurent series, we have E S(r)φ n r n = T (r) E n φ n r n + r n O r −1 , with the components of O r −1 all in r −1 C r −1 and E n an r-independent matrix that is invertible for almost all n. We call E n the trailing characteristic matrix of E with respect to the given multipliers.
Those integers n ∈ Z such that det E n = 0, which is a polynomial equation in n, are called (respectively leading or trailing) (integer) characteristic roots or exponents of E with respect to given multipliers S, T . We denote the set of such leading characteristic exponents byσ(E) and the set of such trailing characteristic exponents byσ(E), with implicit dependence on the S, T multipliers, of course.
If the leading (trailing) multipliers S, T of E are known, the formal adjoint operator E * obviously has S E * = T − * and T E * = S − * as corresponding multipliers, where A − * = (A * ) −1 is the inverse of the conjugate transpose matrix. Dependence on n in E n is due to the action of differential operators r∂ r r n = nr n . Due to the identity (r∂ r ) * = −(r∂ r + 1), using the above multipliers for E * allows us identify the characteristic exponents of E * asσ(E * ) = −σ(E) − 1 andσ(E * ) = −σ(E) − 1.
We will not dwell on when leading or trailing multipliers exist or on their uniqueness, but will just assume that they are given for any particular problem. Our focus, instead, is to use these data as a certificate that a concrete, finite dimensional linear algebra computation can prove some statement about a rational ODE. In practice, often S and T may be taken to be diagonal, with appropriately chosen powers of r on the diagonal. Otherwise, they could be determined by a recursive procedure similar to that used in the analysis of regular and irregular singularities for ODEs with meromorphic coefficients [48]. Related algorithms can be found in the recent papers [1,2].
For the inhomogeneous problem E[φ] = β, the solution φ may only have poles at a subset of the union of the poles of β and of the (regular or irregular) singular points of E[φ]. When the number of such singular points is finite, the knowledge of the leading and trailing multipliers at the singular points (including r = ∞) and at the poles of β is sufficient to reduce the problem of finding all rational solutions φ for rational β to a finite dimensional linear algebra problem [34,Theorem 2.4]. For our purposes, we now need to extend this result to characterize both the kernel and cokernel of E[φ] when acting on the different classes of Laurent series. There exist results in the literature [15,38] that are, on an abstract level, equivalent or even more general than what we present below. But, being written in the abstract language of D-modules, they would require substantial interpretation to be applicable to our very concrete situation.
Let us fix the following notation. We denote the kernel of E by ker E in the C r, r −1 class, by ker + E in the C r −1 [[r]] class, by ker − E in the C[r] r −1 class, and by ker −+ E in the C r, r −1 class. For the cokernel of E, we write coker E, coker + E, coker − E and coker −+ E, using the same conventions. In each case, the operator E[φ] must have a well-defined action on the Laurent series of the appropriate class. If we have in mind Laurent series centered at r = ρ instead of r = 0, then we write ker ρ and coker ρ with appropriate subscripts. Finally, we write ker ρ 1 ,ρ 2 ,... rat and coker ρ 1 ,ρ 2 ,... rat when we restrict our attention to rational functions with possible poles only at r = ρ 1 , ρ 2 , . . ..
Proof . For rational functions with prescribed poles, this is the statement of [34,Corollary 2.5]. For the other cases, the proof can be attacked with similar methods, so we merely sketch it. For ker + E, we can split φ = φ + Sφ + , using the leading multiplier S, where the powers of r in S −1 φ are bounded from below byň and from above byn satisfyingň σ(E) n, while the powers of r in φ + are only bounded from below byn. Then there are no obstructions to solving E[Sφ + ] = −E[φ ] for φ + order by order, since no characteristic exponents will appear by our conditions onn. Then finding all φ such that E[φ + Sφ + ] = 0 reduces to a finite dimensional linear algebra problem. The ker − E case is completely analogous. For ker −+ E, one needs to start with the splitting φ =Šφ − + φ + Sφ + , where nowŠ is the trailing multiplier and again proceed analogously.
To deal with cokernels, it is helpful to introduce the sesquilinear residue pairing where we take the residue at r = ρ (defined as the coefficient of (r − ρ) −1 ) of the product of a column vector φ and the conjugate transposed column vector α (r * = r and its powers are treated as real under conjugation). We may drop the subscript when ρ = 0, −, − = −, − 0 . The residue pairing is of course well-defined for rational functions, but it is also well-defined for the following pairs of Laurent series classes: (bounded from above and below, unbounded), (bounded from below, bounded from below), and (bounded from above, bounded from above). This pairing has several nice properties. It is easy to see that it is non-degenerate in either argument, for each possible combination of pairs Laurent series classes. If S is a Laurent unimodular matrix, then the pairing is obviously invariant under S − * α, Sφ = α, φ . The following lemma is slightly less obvious and will be particularly helpful. Proof . The last statement holds simply because α, φ * = φ, α . To prove the adjoint relation, we need to make use of the defining formula (4.5). First, note that it is obviously true that α, Cφ = C * α, φ , when C is a constant matrix, and that α, r n φ = r n α, φ , for any power of r. Then, due to the Leibniz identity α * ∂ r φ − (−∂ r α) * φ = ∂ r (α * φ) and the observation that Res 0 (∂ r w) = 0 for any Laurent series w, it is also true that α, ∂ r φ = −∂ r α, φ . Hence, the adjoitness property follows for all differential operators, since any operator can be written as a sum of compositions of these three kinds of operators.
We can now state and prove the following lemmas characterizing cokernels in the bounded Laurent series and Laurent polynomial cases. In each case, since E * * = E, we are free to exchange the roles of E and E * . Proof . Obviously, E[φ] is a compatible differential operator. We will consider Laurent series bounded from below (the bounded from above case is completely analogous).
Let L and R denote the spaces of left and right arguments of the residue pairing. Denote by L n ⊂ L and R n ⊂ R the subspaces consisting of Laurent series containing powers of r only equal to n or greater. Note that L ⊥ n = R −n and vice versa. By the invariance of the residue pairing, without loss of generality, we can assume that E[φ] has trivial source and target multipliers (they are just identity matrices). Then E[R n ] ⊂ R n for any n.  Proof . Clearly, E is a compatible differential operator, the induced pairing is well-defined and the non-degeneracy in the first argument follows directly from the non-degeneracy of the residue pairing. It remains to show non-degeneracy in the second argument.
Finally, we can package the above results in the following, for our purposes, conveniently phrased Theorem 4.5. Consider an ordinary differential operator E[φ] with Laurent polynomial coefficients with a finite number of singular points (with given multiplier matrices). Suppose that the leading multipliers and integer characteristic exponents are known at each point ρ ∈ C. The formal adjoint E * [α] will share the same properties.
(a) Given ρ ∈ C, there is a finite basis α ρ j spanning a subspace of ker ρ + E * and as many polynomials ψ ρ k in (r − ρ) −1 such that the residue pairing R ρ jk = Res ρ α ρ j * ψ ρ k is full rank and such that, for any rational ψ, there exist constants c k ρ and a Laurent polynomial φ ρ such that has no poles at r = ρ.
(b) There is a finite basis (α j ) for ker −+ E * and as many Laurent polynomials (ψ k ) such that R jk = Res ρ (α * j ψ k ) is full rank, for any Laurent polynomial ψ, there exist constants c k and a Laurent polynomial φ such that Proof . Part (a) follows from Proposition 4.1 and Lemma 4.3. The implication is that both ker ρ + E * and coker ρ + E are finite dimensional. Hence, the dimension of the subspace of coker ρ + E represented by polynomials in (r − ρ) −1 must also be finite. Let ψ ρ k be a set of polynomials in (r − ρ) −1 that represent a basis of this subspace. Then, by non-degeneracy of the residue pairing, we can pick a finite subset α k ρ satisfying the desired non-degeneracy condition. So by subtracting the appropriate linear combination k c k ρ ψ ρ k from ψ we can find a representativeψ ρ from the same equivalence class in coker ρ + E that has no poles at r = ρ. This argument gives us φ ρ that is merely bounded from below as a Laurent series, but once we have it, we can simply truncate it after some sufficiently high finite power or r.
Part (b) follows from Proposition 4.1 and Lemma 4.4 similarly. Once we know that ker −+ E * and coker E are finite dimensional, it remains only to pick a basis (α j ) for ker −+ E * and a finite set (ψ k ) whose representatives form a basis of coker E. for ψ and φ rational with potential pole locations fixed at ρ 1 , ρ 2 , . . ., is finite dimensional. The functions ψ k , ψ ρ 1 k , ψ ρ 2 k , . . . then constitute a complete set of independent representatives of the cokernel equivalence classes exactly when they satisfy the conditions in (a) and (b) of the theorem.
• ρ = ∞:σ ∞ (D(s 0 , s 1 )) = {0}, with • ρ ∈ {0, 2M, ∞}:σ ρ (D(s 0 , s 1 )) = {0, 1, 2}, with Notice that, besides r = 0 and r = ∞, the only singular point of D(s 0 , s 1 ) is r = 2M . In practice (in Section 5), we will only deal with rational functions that have poles at these three points, so we will restrict our attention only to them. Given the above data, it is possible to explicitly calculate (using computer algebra) the dimensions of the solution spaces ker 2M + D(s 0 , s 1 ) * and ker −+ D(s 0 , s 1 ) * , as well as their bases. Since these solutions are presented by infinite series, it is impractical to write them here explicitly. Instead, we will simply report their dimensions, while explicitly writing a choice corresponding cokernel representatives ∆ k , ∆ 2M k . In practice, the self-adjointness D(s 0 , s 1 ) * = D(s 0 , s 1 ) turns out to be extremely useful and saves a significant amount of effort since the same computer algebra code can be used for both operators.
For ∆ involving poles at points other than {0, ∞}, we can report the following.
For ∆ a Laurent polynomial, we can report the following. This data depends on the values of s 0 and s 1 , so we concentrate only on the cases that will be needed later (in Section 5).
Finally, we would like to record here another result that was actually already implicit in the examples presented in [34, Section 4.1], but not made explicit there. If we consider the Regge-Wheeler decoupling equation (4.4) with zero right-hand side, then solutions constitute morphisms from D s 1 φ 1 = 0 to D s 0 φ 0 = 0. Such a morphism, if it existed between different s 0 and s 1 it could be called a spin raising/lowering operator for radial modes. Such spin raising and lowering operators are known to exist for angular modes [44] and their existence for radial modes would have important applications. Unfortunately, explicit computation shows that the result is negative, at least when we restrict ourselves to operators with rational coefficients.

Equivalence of Regge-Wheeler and Zerilli equations
We have defined the spin-s Regge-Wheeler operator D s in (4.1). When s = 2, D 2 is closely related to the so-called Zerilli operator In fact, the equations D 2 φ = 0 and D + 2 φ = 0 are equivalent (in the sense of Section 3) for all values of ω except the algebraically special frequencies [16,36] satisfying This can be seen from the following equivalence diagram: where are the Chandrasekhar transformation operators [14,Section 4.26], [21,24]. The validity of the above equivalence diagrams may be checked by direct calculation, or can be seen more abstractly from the identities The necessity of the failure of the equivalence identities in (4.7) when α = 0 (due to division by α) can be explained as follows. When α = 0, the two operators simplify to D 2 = −C − f C + and D + 2 = −C + f C − . Reducing each of them to a first order system gives the equivalences Then, using the same machinery that we applied earlier to upper triangular Regge-Wheeler systems (Theorem 4.8 and other results), one can conclude that (a) there are no non-vanishing morphisms between the equations C + φ + = 0 and C − φ − = 0 and (b) neither of the upper triangular systems above can be decoupled to diagonal form. On the other hand, given (a), it is not hard to see that the existence of an equivalence D 2 ∼ D + 2 would imply that both of the above upper triangular systems decouple to the diagonal form C + 0 0 C − , which contradicts (b). Thus, no manipulation will be able to get rid of the poles at α = 0 from the equivalence in (4.7). Remark 4.9. Since both D 2 and D + 2 are formally self-adjoint, while α and M are real, taking formal adjoints of all the operators in the equivalence diagram (4.7) gives another equivalence diagram. The way that we have distributed various factors of α and the adjoint identity C + = C * − makes the new equivalence diagram identical to the original one, and hence the whole diagram formally self-adjoint. While we could have chosen to distribute various constant factors among the operators in (4.7) in different ways, the choice of no prefactor for D + 2 and the self-adjointness of the whole diagram dictates the coefficient α (12M ) 2 in front of D 2 .

Tensor wave equations on Schwarzschild
Consider the exterior Schwarzschild spacetime (M, 4 g) of mass M > 0, where M ∼ = R 2 × S 2 . When S 2 , Ω is the unit round sphere and (t, r), with −∞ < t < ∞ and 2M < r < ∞, are the Schwarzschild coordinates on the (R 2 , g) factor, the exterior Schwarzschild metric has the following (2 + 2)-warped product form: For convenience, we also define (consistently with Section 4) This is a vacuum spacetime the Einstein tensor 4 G µν = 4 R µν − 1 2 r R 4 g µν = 0 vanishes describing the asymptotically flat region outside the horizon of a static spherically symmetric black hole of mass M .
The goal of this section is to write down the vector and Lichnerowicz wave equations 1 v µ := 4 v µ = 0 and 2 p µν := 4 p µν − 2 4 R (µ λκ ν) p λκ = 0, on this spacetime and to reduced them to their radial mode equations, with respect to spherical and time harmonic modes. Eventually we will also show their equivalence to a simplified triangular system of Regge-Wheeler equations (4.1). For background, let us recall that the vector wave operator 1 can be obtained from the Maxwell equation by adding a term factoring through the harmonic (Lorenz) gauge fixing condition 4 ∇ ν v ν = 0, namely Similarly, the Lichnerowicz wave operator 2 can be obtained from the linearized Einstein equation by adding a term factoring through the harmonic (de Donder) gauge fixing condition 4 ∇ λ p νλ = 0, namely where p µν = p µν − 1 2 p λ λ 4 g µν is the trace reversal operation and 4Ġ µν [p] is the linearized Einstein tensor, defined by 4 G µν 4 g + p = 4Ġ µν [p] + O p 2 , keeping in mind that 4 G µν 4 g = 0 for the background Schwarzschild metric.
We will proceed by first using the covariant formalism of [37] for (2 + 2)-warped products and spherical harmonics, then by using Schwarzschild coordinates, then by presenting the ingredients needed to implement the general decoupling strategy from Section 3.2, and finally by giving the explicit equivalence. The intermediate calculations in the decoupling strategy are rather long, unenlightening and best performed by computer algebra. Thus, we present only the final result, but along the way giving some guiding comments on how it could be reproduced. The reader is referred to [32], where the intermediate steps were worked out in more detail for the vector wave operator 1 .
To start, let us briefly introduce the covariant 2+2 formalism from [37]. We use Greek indices, (µ, ν, . . .), for spacetime or M-tensors, lower or upper case Latin indices, (a, b, . . .) or (A, B, . . .), respectively, for tensors on the R 2 or S 2 factors. The tangent and cotangent subspaces of R 2 and S 2 are naturally identified with the appropriate subspaces of T M and T * M, allowing us to map tensors from its factors to M. Thus, any spacetime tensor can be decomposed into sectors with possibly mixed lower-upper case Latin indices, e.g., Note the extra factors of r, which make it convenient to raise and lower Greek indices with 4 g µν , while their sector components have their indices raised and lowered respectively by g ab and Ω AB . This convention extends to higher rank tensors by respecting tensor products. We define the following spacetime tensors by their decomposition into sectors: We denote respectively by ∇ a and D A the Levi-Civita connections on R 2 , g and S 2 , Ω . In terms of them and the decomposition into sectors, the spacetime Levi-Civita connection is determined by where r a = ∇ a r = (dr) a . Now, working in Schwarzschild (t, r) coordinates, we define t a = (∂ t ) a , which is a timelike Killing vector. We also introduce the notation t a = g ab t b = −f (dt) a and ε ab = 2(dt) [a (dr Then r a r a = f , t a t a = −f , and t a = −ε ab r b . The Levi-Civita connection ∇ a is determined by The respective Riemann tensors of the R 2 , g and S 2 , Ω factors are Finally, the spacetime Riemann tensor (computed in [33, equation (23)], for instance) is conveniently expressed as Scalar spherical harmonics Y lm , where l = 0, 1, 2, . . . and −l ≤ m ≤ l are respectively the angular momentum and magnetic quantum numbers, are normalized eigenfunctions of the Laplacian on the unit sphere S 2 , Ω , defined by the conditions where we have used the short-hand − = S 2 (−), with AB the volume form on S 2 , Ω . Where no confusion would arise, we shall regularly omit the l and m indices from Y lm and their coefficients. We use the convenient notations (already mentioned at the top of Section 2) B l := l(l + 1) and A l := (l − 1)l(l + 1)(l + 2) = B l (B l − 2).
By taking covariant derivatives of the Y scalar harmonics and combining them with Ω AB and AB , we can define also vector and tensor spherical harmonics: The above definitions are chosen so that the following convenient orthogonality and normalization properties hold: Note the consequent vanishing of Y 0m The tensorial spherical harmonics also satisfy the following eigenvalue and divergence identities: Under the parity inversion of the sphere, all the Y -harmonics transform as Y → (−) l Y and so are called even, while all the X-harmonics transform as X → (−) l+1 X and so are called odd [40]. Our conventions for spherical harmonics follow [37], whose Appendix A may be consulted for comparison with other conventions in the literature. When needed, we will also decompose sectorial tensors on the R 2 , g factor into coordinate components with respect to the Schwarzschild coordinates (t, r) with the conventions If we raise and lower indices with g ab , this convention implies the Remark 5.1. In the sequel, while working at the mode level, we allow ourselves the freedom of dividing by the expressions ω, B l or A l in the mode parameters. The expressions for various radial mode operators will be valid for any value of the mode parameters where they are welldefined, which might exclude ω = 0, l = 0 or l = 1, when dividing by ω or B l . These singular cases would need to be investigated separately. However, since all expressions will be rational, it would be sufficient to consider Laurent expansions in the mode parameters and just compare corresponding coefficients.
Harmonic dependence on time with frequency ω will refer to being proportional to e −iωt with time independent coefficients. Remark 5.2. We will use the convention that r and ω −1 will carry the dimensions of length, r −1 , ω and ∂ r will carry the dimensions of inverse length, while f , f 1 , B l and A l will be dimensionless. Clearly, by this convention, length dimension is additive under multiplication and composition of differential operators. For practical convenience, upon radial mode decomposition, we will try to make use of only dimensionless differential operators, by for instance inserting extra factors of ω, where needed. Remark 5.3. Whenever we need to take a formal adjoint of a differential operator, we will follow the conventions from Section 3.3.
Before proceeding to mode decomposing concrete equations, let us make some general remarks about the notation. If s µν... is a tensor field, then schematically we will denote its mode coefficients (omitting indices) as The s coefficients represent the even and odd multiplets obtained after the decomposition of a spacetime tensor into sectors, absorbing all angular dependence and spherical tensor indices into the spherical harmonics Y and X. The s coefficients represent the even and odd multiplets after decomposing the remaining sectorial tensors in v into coordinate components with respect to the Schwarzschild coordinates (t, r) on (R 2 , g), absorbing all time dependence into the harmonic factor e −iωt . When convenient, the s and s coefficients may be defined with extra r-, land ω-dependent normalization factors. Given a spacetime differential operator s = O[s] (in general s and s need not be of the same tensor rank), which is even under inversion of the S 2 , Ω factor, we will schematically denote its action on the mode coefficients as The O e and O o operators acting on the s coefficients are precisely the radial mode operator versions of the original spacetime operator O that we will eventually be working with. In principle, the intermediate stage of the decomposition involving s and s can be omitted from the presentation. But we believe that it is worth keeping them, because these expressions would be useful for obtaining the explicit radial mode equations for other coordinate systems on the R 2 , g factor, for instance the Eddington-Finkelstein coordinates instead of the Schwarzschild one. In addition, these intermediate expressions would be useful checkpoints for anyone trying to reproduce our calculations.

Scalar wave equation
Before proceeding to tensor wave equations, for illustration purposes, let us first handle the simplest case of the scalar wave equation z = 0 z := 4 z = 0. (5.1a) Following earlier conventions, we use the notation z = zY = zY e −iωt and z = z r 2 Y = z r 2 Y e −iωt . First, we separate out S 2 dependence via spherical harmonics: Then, we separate out the t coordinate via harmonic time dependence: The last formula defines the radial mode ODE for the scalar wave equation (5.1a). Note that, to economize notation, we have reused the symbol 0 , though it changes meaning when acting on different function spaces, which should be clear from context. Note that the final expression for 0 z(r) is formally self-adjoint, 0 * = 0 . That is no accident, since we have chosen our definitions such that the sesquilinear forms (with compactly supported y(r) and z(r)) all agree and are symmetric up to complex conjugation, because the original 0 operator is formally self-adjoint on spacetime. We will follow the same convention in the rest of the section. This means that some of the components of the radial mode equations in the tensor case might have spurious looking factors of B l or A l . Such factors arise from absorbing the norms of the corresponding tensor harmonics and are needed to maintain manifest formal self-adjointness of the radial mode ODEs.
The transformation In the sense of Section 3 we have an equivalence of 0 with D 0 , exhibited by the diagram Notice that the original operator 0 is dimensionless. We have chosen to include various factors of ω to keep all the other operators in the diagram dimensionless as well. We will follow the same convention throughout all further calculations.
Remark 5.4. Note also that the above equivalence diagram is self-adjoint in the following sense: 0 * = 0 ,˜ 0 * =˜ 0 , h * 0 = h 0 andh * 0 =h 0 , whilek 0 = (k 0 ) * andk 0 = (k 0 ) * . That is, taking formal adjoints of all operators in the diagram, we obtain the same diagram. The various factors of i and ω were chosen also to exhibit this property. In particular, the powers of ω in k 0 , k 0 and˜ 0 are not all independent. There is an obvious economy to self-adjoint equivalence diagrams: the morphism need only be specified in one direction, with the inverse morphism recovered by taking adjoints. In subsequent cases, we have succeeded in identifying a self-adjoint equivalences, but with slightly modified notion of self-adjointness for the decoupled triangular form. It is not obvious to us at the moment why we were successful in that sense and whether that is due to a fortuitous coincidence or a general pattern.

Vector wave equation
To present the radial mode equation for 1 v µ = 0, we follow the same pattern as in the case of the scalar wave equation (Section 5.1).
The tensor spherical harmonic decomposition of a covariant vector field is For conciseness, as before, we will omit the lm indices. Following the notation conventions given at the top of Section 5, we first separate out the spherical harmonics (bold coefficients) and then harmonic time dependence (fraktur coefficients). We will need to use two conventions to normalize the coefficients, which we distinguish by primes (recall that primes are decorations and not derivatives): With these conventions, again reusing the symbol 1 at each stage of mode separation of after separating out the S 2 dependence, and after separating out the time dependence. Our normalization conventions for the mode coefficients were chosen precisely to maintain manifest formal self-adjointness, 1 * e = 1 e and 1 * o = 1 o , at each stage of the mode separation. Next, we will introduce the operators needed to decouple radial modes into pure gauge and constraint violating modes, to set up the decoupling strategy (cf. Section 3.2). The key identities, analogous to E 0 T = T E 1 and E 1 D = D E 0 , hold already at the spacetime level: In addition, analogous to T D = H T E 0 , we have the composition identity Now we define the operators T 1 , T 1 , D 1 , D 1 and their successive mode decompositions: The extra iω factors in the normalizations have been chosen so that all operators acting on radial modes are dimensionless and the first two of the following commutative diagrams remain valid at each stage of mode separation: The last diagram holds at the radial mode level and shows the decoupling of gauge invariant modes (cf. Section 3.2). It decomposes into independent even and odd sectors as follows: We are now ready to give the full triangular decoupling of the even and odd sectors of the vector wave equation 1 v µ = 0.

Odd sector
This sector is structurally similar to the scalar wave case (Section 5.1). It is particularly simple; since it does not contain any gauge modes, it is not constrained by the harmonic gauge condition and the operators Φ 1o , Φ 1o are directly invertible. Thus, the transformation In diagrammatic form, we have This equivalence diagram is self-adjoint in the sense explained in Section 5.1 for the scalar wave equation.

Even sector
This sector is more complicated and is the first prototype for the abstract approach to triangular decoupling that we outlined earlier in Section 3. All that we need to feed into it are the even sectors of the diagrams (5.4) and the additional composition identities Since the implementation of our approach in this particular case was explained in detail in [32] (in particular, see the step-by-step discussion there in Section 3.2), we only state the final result, which has only been slightly adjusted for the purposes of this work (the end of the section shows how). The final decoupled triangular form is While˜ 1 e is obviously not formally self-adjoint, its equivalence with the formally self-adjoint 1 e survives in the existence of an operator Σ 1e effecting the equivalence between˜ 1 e and˜ 1 * The precise equivalence identities take the form while the remaining operators can be recovered from the identities where also h * 1e = h 1e ,h * 1e = Σ 1eh1e Σ 1e . Remark 5.6. These last identities embody the modified sense in which the equivalence diagram (5.2.2) is self-adjoint, in the sense discussed earlier in Section 5.1, but with a twist provided by the matrix Σ 1e . In fact, if we replace˜ 1 e by Σ 1e˜ 1 e , it is no longer upper triangular, but it is formally self-adjoint Σ 1e˜ 1 e * = Σ 1e˜ 1 e , and the twist by Σ 1e is no longer necessary. Since we place importance on the upper triangular form, we will not take advantage of this replacement.
The above upper triangular form of˜ 1 e allows us to classify all of its symmetries (or automorphisms in the sense of Section 3). The key result is the absence of non-vanishing morphisms between the Regge-Wheeler equations D s 0 and D s 1 , except the identity morphism when s 0 = s 1 (Theorem 4.8). This prevents the coupling of the D 0 and D 1 blocks by an automorphism. Also, the single non-vanishing (and non-removable) off-diagonal element in˜ 1 e prevents the exchange of the order of the D 0 blocks. Hence, any automorphism of˜ 1 e must also be upper triangular, with every non-vanishing matrix element proportional to the identity.
With the above logic in mind, the most general automorphism takes the form˜ 1 e A = A˜ 1 e , with  a 1 , b 1 , a 2 . It is invertible when a 1 , a 2 = 0.
With the above choice of k 1e and k 1e , up to homotopy, we have the following equivalences of operators: Finally, for the record, denoting by ˜ 1 e old , k 1e old and k 1e old by corresponding operators from [32], let us note the explicit relations

Lichnerowicz wave equation
To present the radial mode equation for 2 p µν = 0, we follow the same pattern as in the case of the scalar and vector wave equations (Sections 5.1 and 5.2). The tensor spherical harmonic decomposition of a symmetric covariant 2-tensor is For conciseness, as before, we will omit the lm indices. Following the notation conventions given at the top of Section 5, we first separate out the spherical harmonics (bold coefficients) and then harmonic time dependence (fraktur coefficients). We will need to use two conventions to normalize the coefficients, which we distinguish by primes (recall that primes are decorations and not derivatives): With these conventions, again reusing the symbol 2 at each stage of mode separation of p µν = 2 p µν , we have after separating out the S 2 dependence, and after separating out the time dependence the final result is The operators 2 e and 2 o defined in (5.6c) and (5.6d) are our radial mode equations.
Remark 5.7. Note that the final radial mode equations are manifestly self-adjoint, 2 * e = 2 e and 2 * o = 2 o . These expressions are also valid for all values of ω and l. When l = 0, Y l=0 A = X l=0 A = Y l=0 AB = X l=0 AB = 0, and when l = 1, Y l=1 AB = X l=1 AB = 0, meaning that their coefficients are not well-defined. By consistency, the corresponding components of 2 e and 2 o simply vanish.
Next, we will introduce the operators needed to decouple radial modes into pure gauge and constraint violating modes, to set up the decoupling strategy (cf. Section 3.2). The key identities, analogous to E 1 T = T E 2 , E 2 D = D E 1 and so on, hold already at the spacetime level: Some of the above radial mode operators are related by taking adjoints, namely For the trace operators again we return to the convention for scalars from Section 5.1: using z Y e −iωt = z Y = 4 g λκ p λκ =: tr[p] we get tr o := 0; and using 1 r 2 z Y e −iωt = 1 r 2 z Y = 4 g λκ p λκ =: tr [p ], we get The extra iω factors in the normalizations have been chosen so that all operators acting on radial modes are dimensionless and the first three of the following commutative diagrams remain valid at each stage of mode separation: The last diagram holds at the radial mode level and shows the decoupling of gauge invariant modes (cf. Section 3.2). It decomposes into independent even and odd sectors as follows. The odd sector is described by the following radial mode operators: Remark 5.8. The operator Φ 2o is indeed gauge invariant, Φ 2o • D 2o = 0. It is essentially the well-known odd sector Regge-Wheeler scalar [9, equation (2.15)], whose decoupling was first identified in [40].
In the even radial mode sector, we could choose D Φ 2e to be proportional to D 2 as well, but the corresponding Φ 2e and Φ 2e would be rather large and, at the same time, their complexity would obscure some important properties that we would like to highlight. Instead, we choose D Φ 2e to be a 2 × 2 first order system, with corresponding Φ 2e and Φ 2e that are simple 6 and satisfy the gauge invariance condition Φ 2e • D 2e = 0. The corresponding formulas are the even sector. More precisely, the equivalence is established by concatenating the following equivalence diagrams, where the middle vertical arrow corresponds to a self-adjoint first order reduction of D 2 : where the operators in the first square are given by and α (equation (4.6)) is the frequency dependent constant which vanishes at so-called algebraically special frequencies [16]. Notice that α appears in the denominator of q andq , thus creating poles at the algebraically special frequencies, the only poles in frequency other than ω = 0 that enter the equivalence morphisms with eventual triangular decoupling of 2 e (Section 5.3.2). The same constant α and the poles created by it appears in the Chandrasekhar transformations that relate the Regge-Wheeler D 2 and Zerilli D + 2 operators. We give a complete discussion of this relation in Section 4.1 and the necessity of these poles. Remark 5.9. Using the Zerilli operator D + 2 (and its first order reduction) instead of the Regge-Wheeler operator D 2 in (5.8) results in equivalence maps that are free of poles at α = 0, the algebraically special frequencies (in fact, free of all poles in frequency except ω = 0). We have checked that this change also eliminates the α = 0 from the equivalence morphisms with the triangular decoupled form in Section 5.3.2. However, a different undesirable property appears. As is well-known and can be explicitly seen from our discussion in Section 4.1, both the Zerilli operator D + 2 and the corresponding decoupling morphisms from 2 e contain r-dependent poles in the angular momentum quantum number l. This means that the corresponding radial mode operators cannot be easily interpreted as the mode separated form of differential operators on spacetime. On the contrary poles at α = 0 (since it is a constant independent of r and depending polynomially on ω and B l ) do not prevent us from finding such an interpretation in terms of differential operators on spacetime (Corollary 2.4). Also, besides being simpler, we already have a good understanding of D 2 with respect to Section 4, while the analogous treatment of D + 2 would have to be done separately. For these reasons, we refrain from using the Zerilli operator in this work.
We are now ready to give the full triangular decoupling of the even and odd sectors of the Lichnerowicz wave equation 2 p µν = 0.

Odd sector
This sector is structurally very similar to the even sector of the vector wave equation (Section 5.2.2). What we need to feed into the decoupling strategy of Section 3.2 are the odd sectors of the diagrams (5.7) and the additional composition identities The key equivalence diagram (3.3a) for the E 2o -Φ 2o -T 2o system is then obtained by following exactly the same steps as for the vector even sector (see [ While˜ 2 o is obviously not formally self-adjoint, its equivalence with the formally self-adjoint 2 o survives in the existence of an operator Σ 2o effecting the equivalence between˜ 2 o and˜ 2 * o , The precise equivalence identities take the form where the operatorsk 2o ,k 2o , h 2o ,h 2o arē while the remaining operators can be recovered from the identities where also When it comes to the properties of diagram (5.10) with respect to taking formal adjoints, Remark 5.6 applies here equally well, with appropriate transposition of notation.
The above upper triangular form of˜ 2 o allows us to classify all of its symmetries (or automorphisms in the sense of Section 3). The key result is the absence of non-vanishing morphisms between the Regge-Wheeler equations D s 0 and D s 1 , except the identity morphism when s 0 = s 1 (Theorem 4.8). This prevents the coupling of the D 1 and D 2 blocks by an automorphism. Also, the single non-vanishing (and non-removable) off-diagonal element in˜ 2 o prevents the exchange of the order of the D 1 blocks. Hence, any automorphism of˜ 2 o must also be upper triangular, with every non-vanishing matrix element proportional to the identity.
With the above logic in mind, the most general automorphism takes the form˜ 2 o A = A˜ 2 o , with parametrized by the 3 constants a 1 , b 1 , a 2 . It is invertible when a 1 , a 2 = 0.
With the above choice of k 2o and k 2o , up to homotopy, we have the following equivalences of operators:

Even sector
The even sector of the Lichnerowicz wave equation is by far the most complicated case of those we presented in this work. It requires the full arsenal of tools that we have outlined in , and the composition identities is clearly equivalent to the diagonal operator We use this equivalence to simplify the first diagram above before proceeding. 3. Next, we need to simplify the joint system 2 e p = 0, Φ 2e [p] = 0, tr e [p] = 0 and T 2e [p] = 0.
It is of mixed second, first and zeroth orders. We reduce it to mixed first and zeroth orders only, as follows. We use 2 e to eliminate second order derivatives from ∂ r T 2e , giving us six independent first order equations. Adding to the list the first component of Φ 2e , we get seven independent first order equations, from which we can solve for the first derivatives of the seven metric components h tr , j t , h tt , h rr , K, j r and G. Using these first order equations to eliminate all derivatives from ∂ r tr, we get a zeroth order equation. Including also tr itself and the second component of Φ 2e , we get three independent zeroth order equations.
4. Next, we use these three zeroth order equations to eliminate the h tr , h tt and h rr components completely. The choice of these three components is dictated by the simplicity of the denominators needed for the inversion. What remains is the following first order system for the remaining four metric components: 5. The above 4 × 4 first order system is then equivalent to the 2 × 2 second order D 0 0 0 B l D 1 . In one direction, the equivalence is given by the appropriately modified D 2e . Obviously, the 2 × 2 second order system can also be reduced to a 4 × 4 first order system, while homotopy equivalence reduces the morphism induced by D 2e to a zeroth order operator (a 4 × 4 matrix, in fact) between these two first order systems. This 4 × 4 matrix is invertible and its inverse, after undoing the order reductions, gives the equivalence in the other direction.
6. In the above steps, the various eliminations of components introduce divisions only by ω, r and f , except the final 4 × 4 matrix inversion, which introduces division by B l . Whenever encountering choices, e.g., which components to eliminate, we are guided by the desire not to introduce other factors into the denominators, in order to keep the intermediate expressions in the calculations more manageable.
7. From this point on, it is enough to continue blindly following the strategy from Section 3.2 to end up with a block upper triangular form of 2 e . At this stage, the diagonal blocks are not all proportional to Regge-Wheeler operators D s , but include also 0 , 1 e and D Φ 2e . Now we can use the results of Sections 5.1, 5.2.2 and the equivalence (5.8) to transform each of these blocks themselves to upper triangular form (using the transformation rules from Lemma 3.3), now with only Regge-Wheeler operators on the diagonals. 8. It remains now only to apply the methods of Section 4 to systematically simplify the off-diagonal components of the above upper triangular system. Practically speaking, it turns out to be convenient to first decouple the D 2 equations from the rest, then the D 1 equations from the remaining D 0 equations. Within each still coupled block, it is convenient to complete the simplification by going through successive super-diagonals.
The final decoupled triangular form (recalling the definition of α in (4.6)) of 2 e p even = 0 is Remark 5.10. At this point, it is worth noting that equation (5.11) is the first appearance of the full triangular decoupled form˜ 2 e in the literature. Previously, the original investigations in [9] had only produced the (φ 0 , φ 1 , φ 2 ) subsystem, without studying how it couples to the remaining degrees of freedom, or considering the maximal simplification of the off-diagonal couplings as we have done here. A similar remark can be made about˜ 2 o from equation (5.9).
While˜ 2 e is obviously not formally self-adjoint, its equivalence with the formally self-adjoint 2 e survives in the existence of an operator Σ 2e effecting the equivalence between˜ 2 e and˜ 2 * where the (rather lengthy) explicit formulas fork 2e ,k 2e and h 2e ,h 2e can be found in Appendix B, 7 while the remaining operators can be recovered from the identities where also h * 2e = h 2e ,h * 2e = Σ 2eh2e Σ 2e . (5.13) When it comes to the properties of diagram (5.12) with respect to taking formal adjoints, Remark 5.6 applies here equally well, with appropriate transposition of notation. Just in the case of˜ 2 o (Section 5.3.1), the above upper triangular form of˜ 2 e allows us to classify all of its symmetries (or automorphisms in the sense of Section 3). Following the same logic, the most general automorphism takes the form˜ 2 e A = A˜ 2 e , with parametrized by the 9 constants a 0 , b 0 , c 0 , d 0 , e 0 , g 0 , a 1 , b 1 , a 2 . It is invertible when a 0 , b 0 , a 1 , a 2 = 0. Note that it is allowed to deviate from upper triangular form by mixing the two inner D 0 blocks.
With this choice of k 2e and k 2e , up to homotopy, we have the following equivalences of operators: iωr 2 tr •k 2e ∼ 0 0 1 0 1 0 0 , Note that the adjointness properties of the equivalence diagram (5.12) do not uniquely fix the operators k 2e and k 2e . It turns out that a subfamily of (5.14), the symmetries of˜ 2 e , respects the adjointness properties of the equivalence when A −1 = Σ 2e A * Σ 2e . This constraint is satisfied for instance by the family a 0 = a 1 = b 0 = 1, d 0 = e 0 = 0, and g 0 = c 2 0 . The remaining c 0 parameter may be uniquely fixed by requiring that the operator equivalent to iωr 2 tr involves only the two inner spin-0 components φ 0 and χ 0 , and the operator equivalent to D 2e has the right sign. No other ambiguity remains after that.

Discussion
Inspired by previous more ad-hoc work [9,41] on decoupling the radial mode equations of the vector wave and Lichnerowicz wave equations on the Schwarzschild spacetime, in [32,34] we have initiated a research program of finding a systematic approach to the problem, as well as simplifying and extending the results. This work completes part of our program, which consists of successfully applying the developed methods to both the odd and even sectors of the Lichnerowicz wave equation on Schwarzschild. The results include explicit formulas for reducing the complicated mode separated radial equations to a triangular system of sparsely coupled spin-s Regge-Wheeler equations, by separating the degrees of freedom into pure gauge, gauge invariant, and constraint violating modes. Our systematically derived results extend those of [9,41] because these earlier works did not attempt to uncover the full triangular structure of the decoupled equations or their maximal simplification (in the context of rational ODEs).
Our main results are summarized in Theorem 2.1 in Section 2, where we also list a number of quick corollaries that indicate important applications of our results to the spectral analysis of electromagnetic and gravitational perturbations of the Schwarzschild black hole in harmonic gauges. One application that we intend to pursue in future work, by leveraging the well-understood spectral theory of Regge-Wheeler operators [20], is an explicit mode-level construction of the retarded/advanced Green functions for these perturbations. Such a construction is a first step toward classical stability analysis and the construction of quantum field theoretic propagators. The absence decoupling results, such as ours, have so far constituted a great obstacle in these applications.
An open question is whether our decoupling strategy could be applied to Kerr black holes, where so far the harmonic gauges have received very little attention. The main obstacle to a direct generalization seems to be the lack of a complete separation of the vector wave and Lichnerowicz wave equations on the Kerr background. Other than that, the basic ingredients p (10) = h 2 , Y (10) µν → 0 0 0 r 2 X AB .
Our conventions on the definition of the spherical harmonics Y , Y A , X A , Y AB and X AB are taken directly from [37] (MP), where they are compared to spin-weighted and pure spin harmonics in [37, Appendix A]. However, our normalizations for the coefficients differ by some factors of r [37, equations (4.1)-(4.5), (5.1)-(5.4)]: The conventions for the spherical harmonics used in [5, Appendix A] (BL) are different. In their notation [5, equation (8)], the trace-reversed metric perturbationp µν = p µν − 1 2 p λ λ 4 g µν is expanded as where they factored out some normalizing constants a (i)l [5, equations (9,10)] from their basis. The relation of their basis to ours is the following, which can be read off by comparing the corresponding coordinate formulas in [37, Appendix A] and [5, equations (9,10),Appendix A]:

B Operators for the Lichnerowicz even sector
Below, we give explicit formulas for the morphism operatorsk 2e ,k 2e and homotopy operators h 2e , h 2e that appear in the equivalence diagram (5.11) in Section 5.3.2. The morphism operators are first order. Hence, they can be written as 2e .
Let us denote by C andC the coefficient matrices of ∂ 2 r in 2 e and˜ 2 e , respectively. They are clearly diagonal and invertible. By comparing the highest order coefficients in the identity 2 ek2e =k 2e˜ 2 e , we find the relation Note that the 4th column of K is proportional to 1/α, which is reflected in the notation in a special way. It is easy to see (comparing the coefficients of highest order terms of the relevant identities) that the homotopy corrections h 2e andh 2e depend only on k (1) 2e andk (1) 2e coefficients. More explicitly, plugging in their dependence on K, we obtain h 2e = C −1 KCΣ 2e K * C −1 , andh 2e = Σ 2e K * C −1 K.