{'fraction_non_alphanumeric': 0.047193009993198344, 'fraction_numerical': 0.013707947470308167, 'mean_word_length': 4.045802534318902, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Gauge invariant perturbations of Schwarzschild black holes in horizon-penetrating coordinates

We derive a geometrical version of the Regge-Wheeler and Zerilli equations, which allows us to study gravitational perturbations on an arbitrary spherically symmetric slicing of a Schwarzschild black hole. We explain how to obtain the gauge-invariant part of the metric perturbations from the amplitudes obeying our generalized Regge-Wheeler and Zerilli equations and vice-versa. We also give a general expression for the radiated energy at infinity, and establish the relation between our geometrical equations and the Teukolsky formalism. The results presented in this paper are expected to be useful for the close-limit approximation to black hole collisions, for the Cauchy perturbative matching problem, and for the study of isolated horizons.


I. INTRODUCTION
In the last decade, perturbation theory for black holes has played, in several different ways, a key role in numerical and computational relativity. Already in the seventies it proved to be a very valuable tool to predict gravitational waveforms from processes such as a particle falling towards a black hole. Since the early nineties, due to the network of interferometric gravitational wave detectors in construction, there has been renewed interest in predicting waveforms for strong sources of gravitational waves such as black hole collisions. In particular, the first predictions using perturbation theory in this new era have been quite striking [1]. Some of the applications of perturbation theory in recent years involved computing the evolution for different conformally flat initial data describing black holes in the close limit in order to predict radiated energy and angular momentum [2], to provide both analytical understanding and benchmarking of full numerical results [3], or to quantify the amount of spurious radiation in conformally flat initial data [4] (see [5] for a general review). The usual Regge-Wheeler (RW) -Zerilli [6], [7] and Teukolsky [8] formalisms have also been extended to second order [9], a necessary step in providing first-order perturbations with their own "error bars" [10]. Other recent approaches use black hole perturbations to extend the computational domain in numerical simulations to the radiative zone via Cauchy-perturbative matching [11], or concentrate full numerical resources in the nonlinear regime and let perturbation theory take over in the late stage of black hole collisions [12].
All of the applications just mentioned, though diverse, have a common feature: they are limited to perturbations of Schwarzschild black holes in Schwarzschild coordinates, and Kerr black holes in Boyer-Lindquist coordinates. The reason for this is that, until very recently, most of the initial data typically used in numerical relativity were for maximal slicing, and thus reduced, in the various regimes where perturbation is used (far region, late times, initially close black holes, etc), precisely to the Schwarzschild and Kerr spacetimes in Schwarzschild and Boyer-Lindquist coordinates, respectively. In recent years, however, work was started on Kerr-Schild-type initial data [13], which are not maximal. Part of the motivation for introducing this new kind of initial data is to avoid the typical grid stretching that maximal slicings produce near the event horizon 1 , a stretching that eventually causes numerical simulations to crash 2 . One is then faced with the fact that in order to accommodate these new initial data, either for the close-limit approximation or for Cauchy-perturbative matching, a formalism is needed that allows perturbations in more general slicings than Schwarzschild and Boyer-Lindquist.
Another important motivation for having such a formalism in place is to study the recently developed isolatedhorizon formalism [17] in the perturbative regime. For such studies, one needs to be able to analyze a neighborhood of the background horizon, which necessitates the use of horizon-penetrating coordinates.
The two most used approaches to black hole perturbations have been the RW-Zerilli and the Teukolsky ones. Each of these methods has its own advantages and limitations: The Teukolsky formalism can be used for rotating black holes, but one cannot obtain the whole perturbed geometry but, rather, Ψ 4 or Ψ 0 (this is enough to compute radiation, though) 3 . The RW -Zerilli technique, on the other hand, provides the whole perturbed metric, but is limited to non rotating black holes.
The Teukolsky equation, in its original formulation, can, in fact, be used to describe perturbations around any Petrov type-D background, without relying on a particular choice of coordinates. Work has started very recently on the application of this to Kerr-Schild black hole perturbations [18]. This paper, in turn, develops an appropriate extension of the RW-Zerilli formalism to perturbations of a Schwarzschild black holes in arbitrary spherically symmetric coordinates. One can imagine a huge variety of applications of such an extension; here we have concentrated on the aspects of the formalism that we need in order to proceed with our main motivations. In order to generalize the RW-Zerilli formalism, we start from a perturbation formalism introduced by Gerlach and Sengupta [19] and derive two master equations which hold in any spherically symmetric coordinates of the background, but reduce to the equations obtained by Regge-Wheeler and Zerilli if one uses the standard Schwarzschild coordinates.
Our approach is organized as follows. In section II we present the basic formalism that decouples the field equations into the generalized RW and Zerilli ones. The special cases with total angular momentum l = 1 and l = 0 are treated carefully. In section III we work out a relation needed for Cauchy-perturbative matching, namely, the one between the RW and Zerilli functions and the ADM quantities. In section IV we establish the relation between the present formalism and the Teukolsky one, a relation that is desirable not only to compute the radiated energy and make contact with [18], but also from a conceptual point of view. Finally, in section V we comment on the properties of the RW and Zerilli equations, and on a numerical code that we have written to solve them. In order to establish the contact between the abstract formalism in the body of this paper and more direct applications, we give some explicit expressions in appendix A. Finally, in appendix B, we summarize some properties of spin-weighted spherical harmonics which are needed in section IV.

II. THE GENERALIZED RW AND ZERILLI EQUATION
In what follows, we assume that the background spacetime (M, g) can be represented as a product ofM = M/SO(3) and S 2 with metric Hereĝ = dΩ 2 is the standard metric on S 2 , andg and r denote the metric tensor and a positive function, respectively, defined on the two-dimensional pseudo-Riemannian orbit spaceM . In what follows, lower-case Latin indices refer to coordinates on (M ,g), while capital Latin indices refer to the coordinates ϑ and ϕ on (S 2 ,ĝ). Below, we will derive perturbation equations which do not depend explicitly on the background metric coefficients. In fact, we will only use the background equations which are given by the components of the Einstein tensor Here, N =g(dr, dr), andκ denotes the Gauss curvature of the metricg. A coordinate-invariant definition of the ADM mass is given by We can see by inspection that this M is the mass if Schwarzschild coordinates are used; on the other hand, M is defined in terms of scalars onM , and therefore M represents the ADM mass in any coordinates onM . Note that in a vacuum spacetime, equation (2) Since the background is spherically symmetric, it is convenient to expand the perturbed metric in spherical harmonics: where H ab denotes a tensor field, Q b and h b vector fields, and K, F and k scalar fields onM . Here, Y ≡ Y lm are the standard spherical harmonics, and S A = ( * dY ) A =ǫ B A∇ B Y and 2∇ (A S B) ≡∇ A S B +∇ B S A form a basis of odd-parity vector fields and symmetric tensor fields, respectively, on S 2 (See Appendix D of Ref. [21] for more details on spherical tensor harmonics.). We suppress the indices lm and the sum over these indices since the modes belonging to different pairs of lm decouple in the perturbation equation. The Y lm are normalized with respect to the standard metricĝ on S 2 , an exception being the cases l = 0 and l = 1: There, we choose the normalization such that Y 00 = 1, and S 2 Y 1mȲ 1m dΩ = 4π/3.
In what follows, it will also be convenient to use a coordinate-free notation for differential forms on (M ,g): * and d † ≡ * d * denote the Hodge dual and the co-differential operator, respectively, with respect tog. That is, * u a dx a =ǫ ab u a dx b , whereǫ andǫ denote the standard volume elements in (M ,g) and (S 2 ,ĝ), respectively 4 . For example, we havẽ A further simplification comes from the fact that a spherically symmetric metric is invariant under parity transformation x → −x. As a consequence, the above defined amplitudes decouple into two sets, one set transforming like Y (called scalar perturbations or even-parity perturbations) and the other set transforming like S = * dY (called vector perturbations or odd-parity perturbations) under parity transformations. In this sense, the amplitudes H ab , Q b , K and F have even parity while the amplitudes h b and k have odd parity.
A. The odd-parity sector We start with the simpler case of the odd-parity sector. The perturbations of g µν are parameterized in terms of a scalar field k and a one-form h = h a dx a , where k and h a depend on the coordinates x b only. Note that for l = 1,∇ (A S B) vanishes and k is not present. For l = 0, S A = 0 and there are no gravitational perturbations.

Coordinate-invariant amplitudes
A vector field X = X µ ∂ µ generating an infinitesimal coordinate transformation with odd parity is determined by Using the fact that to linear order, δg µν transforms like the Lie derivative of the background metric with respect to X, we find the following transformations Note that one can choose a gauge in which k = 0. This gauge, which is usually called the RW gauge, is unique.
For l ≥ 2, one can construct the coordinate-invariant one-form For l = 1, we will see that only the invariant two-form enters the perturbation equations.
In terms of these gauge-invariant quantities, the components of the linearized Einstein tensor are δG ab = 0, where the background equations have been used, and where here and in the following, λ ≡ (l − 1)(l + 2).

The master equation
The vacuum perturbations with odd parity are obtained from equation (4), which yields The usual way to derive the RW equation for l ≥ 2 from equation (6) is to decompose the one-form h (inv) with respect to Schwarzschild coordinates, h (inv) = h (inv) t dt+h (inv) r dr, and to use the integrability condition (5) This yields an equation for h (inv) r alone, which is then cast into a wave equation for the functionΦ = (1−2M/r)h (inv) r /r (Φ defined below). This can also be achieved in a coordinate-invariant way as follows: One uses the integrability conditiond † h (inv) = 0 to introduce the scalar potential Φ according to h (inv) = * d(rΦ) =ǫ ab∇ a (rΦ)dx b . Equation (6) may then be integrated to yield the following wave equation where the two-dimensional Laplacian of a function is∆Φ ≡ −d † dΦ =∇ a∇ a Φ. Here, the free constant in the potential Φ has been used to set the integration constant to zero. Equation (7) is the coordinate-invariant version of the RW equation. Indeed, we have not specified any coordinates on the orbit manifoldM . Using the coordinate-independent vacuum background equation 0 = r 2 G a a = 2(r∆r + N − 1), equation (7) finally assumes the form with ¿From equation (7), we also get the following relation which enables us to compute Φ from the gauge-invariant one-form h (inv) . For l = 1 equation (7) is immediately seen to admit the solution 1/r. Since λ = 0, we may also directly integrate equation (6). This yields * where 6J is a constant of integration. At this point, it is important to recall that the one-form h is not coordinateinvariant, but transforms according to h → h + r 2 d(f /r 2 ). This implies that the solution of the homogeneous part of the above equation is pure gauge. A special solution is As explicitly shown in [21], this describes the Kerr metric in Boyer-Lindquist coordinates in first order of the rotation parameter a = J/M . By equation (10), J is defined in a coordinate-invariant way. In summary, a general l = 1 perturbation is given by with f an arbitrary function on the orbit space.

B. The even-parity sector
The even-parity perturbations of g µν are parameterized by a symmetric tensor field H ab , a one-form Q b and two scalar fields K and G on the orbit spaceM , Here, the basis of symmetric tensors in δg AB is chosen to be orthogonal with respect to the inner product induced by g. Furthermore, one has∇ A∇B Y + 1 2 l(l + 1)ĝ AB Y = 0 for l = 0, 1; hence the amplitude G is not present in those cases. For l = 0, the amplitude Q b is also absent.

Coordinate-invariant amplitudes
An infinitesimal coordinate transformation with even parity is generated by a vector field X with where ξ a and f are a vector field and a function, respectively, onM . With respect to this, the metric perturbations transform according to Here and in the following, ξ b|a ≡∇ a ξ b denotes the covariant derivative with respect to the orbit metricg and v a ≡ r |a /r .
For l ≥ 2, one can construct the following set of coordinate-invariant amplitudes: where p a = Q a − 1 2 r 2 G |a . The (generalized) RW gauge is defined by choosing ξ a and f such that Q b and G vanish. We see that in this gauge, which is also unique, H (inv) ab and K (inv) coincide with H ab and K. For l = 1, there is no such simple choice of coordinate-invariant amplitudes, since G is not present in this case. Nevertheless, we can always chose the gauge such that Q b vanishes. One then remains with H ab and K, which are subject to the residual coordinate transformations as in (13) For l = 0, Q b and G are absent anyway and one can arrange the gauge such that K = 0. In summary, it is sufficient to derive the linearized Einstein equations for the perturbed metric where for l ≥ 2, H ab and K can be replaced by their coordinate-invariant counterparts defined in (14,15).

The master equation
The long but straightforward computation of the linearized Einstein tensor is given in [20]. The equations' structure becomes much more transparent if one first splits the two-tensor H ab into its trace and traceless part and then introduces the one-form whereĤ ab denotes the traceless part of H ab . A similar split is performed for the components of the Einstein tensor. As a result, the relevant components of the Einstein tensor define two scalars S and T and two one-forms U and V according to The vacuum field equations are then expressed in terms of the one-form C and the two scalars H =g ab H ab and K. The simplest equation, which is present only for l ≥ 2, gives hence H ab is traceless (For l = 0, 1, we can make use of the residual gauge freedom in order to impose H = 0. Residual coordinate transformations are then of the form (13) with ξ |a a = 0 and ξ a = −r 2 f |a for l = 1.). Using H = 0, the remaining equations reduce to where equation (18) is void for l = 0. We recall that for l ≥ 2, we should replace H ab by H (inv) ab and K by K (inv) in the above equations in order to give them a coordinate-invariant meaning. For l ≥ 1, we compute the component of U parallel to * dr: This motivates us to replace the one-form C with the one-form In terms of Z and K Einstein's equations become where we have defined Using the background equation 0 = r 2 G a a = 2(r∆r + N − 1) and N = 1 − 2M/r, one finds a 0 = 6M/r. In view of equation (21), we may introduce the scalar field ζ according to Z = dζ. Equation (24) may then be integrated to yield It is now clear how to obtain a single, second order differential equation for ζ: First, we eliminate∆K from equations (22) and (23). This gives Next, this equation is used to eliminate K in (25). Hence, (Note that for l = 1, this equation is equivalent to equation (26) and thus is also valid in that case.) Finally, we define the new scalar function Ψ by in order to remove the first order derivatives. This yields the Zerilli equation [7], where Before we discuss the special cases l = 0 and l = 1, we make two remarks: First, for l ≥ 2, the scalar field Ψ can be obtained from the Zerilli one-form Z using equation (27). The second point is that it is also possible to obtain the RW equation for the scalar Φ e = r 2d † Z. In fact, Chandrasekhar (see, e.g. [27]) has shown that the equations of RW and Zerilli for a Schwarzschild background are equivalent in the frequency domain. However, in the time domain, we were not able to express Z in terms of Φ e and its derivatives alone. For this reason, we will use the Zerilli equation in the even-parity sector and not the RW equation.
For l = 1, equation (27) reduces to However we recall that for l = 1, H ab , and hence also ζ are not defined in a coordinate-invariant manner. Under the residual gauge-freedom (13) with ξ a = −r 2 f |a and we find that Z transforms according to Since f is an arbitrary solution of equation (31), and since the equations (30) and (31) are equivalent, it is clear that every solution of (30) corresponds to pure gauge. In particular, we can choose the residual gauge in order for ζ to vanish. In this gauge K vanishes as well, as a consequence of equation (25). The even-parity sector is therefore empty for l = 1.
For l = 0, one can choose the gauge such that both H and K vanish. Then equations (17) and (19) yield which has the general solution Here, δM is a constant describing the variation of the ADM mass, and h is a function that only depends on r.
Comparing this with a residual gauge, which is generated by ξ a =ǫ ab k |b for a function k of r, we get showing that the function h(r) above corresponds to pure gauge. This can also be seen in a gauge-invariant way: Recall that for any spherically symmetric metric of the form (1) we defined the mass parameter M through 1−2M/r = N =g ab r |a r |b . Using the fact that (for Y=1) δg ab = H ab and δ(r 2 ) = r 2 K, we obtain 2δM = rr |a r |b H ab − rr |b (rK) |b + M K.
It can be checked that the RHS is indeed a gauge-invariant combination. On the other hand, for K = 0, equation (33) yields 2δM =g(rC, dr), as above.

C. Summary
In both the odd-and the even-parity sector, perturbations on any spherically symmetric vacuum background are described by a wave equation of the form where∆ is the Laplacian with respect to the orbit metricg and where the potential V depends on the ADM mass M , r and the angular momentum number l only. For l = 0 and l = 1, there are no dynamical modes. The only physical solutions in those cases are stationary, describing variation of the mass and angular momentum. The gauge-invariant part of the metric can be obtained from u and vice-versa. These relations are going to be made more precise in the next section.
Finally, we would like to mention that our gauge-invariant perturbation formalism has also been generalized to the case where matter fields are coupled to the metric [20,21]. In the case of Einstein-Maxwell, we were able to generalize the equations obtained by Moncrief [22]. However, as we have argued in a recent Letter [23], the perturbation formalism presented here fails to yield a wave equation of the form (34) with a symmetric potential, V = V T , when non-Abelian fields are coupled to the metric.

III. RELATION TO THE ADM QUANTITIES
As mentioned in the introduction, one of the motivations for the present work is Cauchy-perturbative matching. This amounts to matching numerically, at each time step, the variables used in a nonlinear code with the ones used in the perturbative regime (in our case the RW and Zerilli functions). The matching takes place at a timelike boundary. For this purpose we explicitly show the relation between the RW and Zerilli gauge-invariant potentials and the ADM quantities, namely, the three-metric and the extrinsic curvature. This does not restrict the formulation of Einstein's equations to be used in the nonlinear regime, since for a formulation other than the standard ADM (e.g. conformal ADM, or a hyperbolic formulation) the relevant quantities can be obtained from the three-metric and the extrinsic curvature, and vice-versa.
So our aim is to make explicit the relationship between the scalar fields Φ and Ψ satisfying the RW and Zerilli equations (8) and (28) and the components of the linearized 3-metric and extrinsic curvature. We will show in this section that -modulo gauge transformations -there is a one-to-one correspondence between δḡ ij , δK ij and the scalar amplitudes Φ,Φ ≡ ∂ t Φ, Ψ,Ψ. Furthermore, this correspondence involves no time-derivatives. For example, it is possible to expressΨ in terms of purely spatial quantities, i.e. δḡ ij , δK ij and their spatial derivatives only.
We assume that the full metric, satisfying the nonlinear field equations, has the ADM form where µ is a variational parameter, such that for µ = 0, the metric is spherically symmetric. With respect to the 2+2 split (1), the orbit metricg takes the formg where x is any radial coordinate, α and β ≡ β x are the background lapse and shift, respectively, and γ 2 ≡ḡ xx . The components of the extrinsic curvature are where a prime denotes differentiation with respect to x, and where we have also introduced the normal derivative The components of the linearized metric have the form where i = x, A. Note that we use perturbations of the coshift rather than the shift vector. This fact will turn out to be important when we try to express the ADM quantities in terms of the RW and Zerilli scalars. Similarly, the components of the linearized extrinsic curvature are given by

A. The odd-parity sector
In the odd-parity sector with l ≥ 2, the only non-vanishing perturbations can be parameterized according to δK xA = π 1 S A , δK AB = 2π 2∇(A S B) .

The potentials in term of the ADM quantities (l ≥ 2)
We want to express Φ andΦ in terms of the quantities b, h 1 , h 2 , π 1 , π 2 and their spatial derivatives. First, we observe that h t = b, h x = h 1 , and k = h 2 where h t , h x and k are the amplitudes introduced in (3). Therefore, we obtain Next, the equations (36) yield the relations Using the last of these relations, eq. (40), to re-express, in eq. (37), time derivatives of h 2 in terms of spatial quantities, the components of the gauge-invariant one-form h (inv) take the form where h . Next, one uses equation (9), trading time derivatives for spatial ones with the aid of (39), to find which is one of the formulae we were looking for. In order to obtain the time derivative of Φ, one uses the definition of Φ, i.e. h (inv) = * d(rΦ). This yields which one can solve forΦ. Using eqs. (43,41,42), the result iṡ 2. The ADM quantities from the potentials (l ≥ 2) On the other hand, given Φ andΦ, we obtain δg ij and δK ij in the following way: First, we compute h 3. The special case l = 1 For l = 1, the amplitudes h 2 and π 2 are absent. According to the analysis in the last section, the only physical solution is the Kerr mode. Using equation (10), one finds that the rotation parameter (the only gauge invariant for l = 1 ) can be extracted from the ADM quantities according to On the other hand, using (11), one finds where f parameterizes the gauge freedom, and where N = 1 − 2M/r = −(∂ 0 r) 2 /α 2 + r ′2 /γ 2 . The amplitude π 1 then follows from (52).

B. The even-parity sector
Here the perturbations are 1. The potential in terms of the ADM quantities (l ≥ 2) ¿From equations (35) one finds H tt = −2α a − β 2 h + 2βb 1 , H tx = b 1 , H xx = h, Q t = b 2 , Q x = q, while K and G agree with their definitions in the previous section. The expressions for the linearized curvature tensor, equation (36), yield At first sight it is not clear how the Zerilli one-form Z, defined in (20) can be expressed in terms of spatial amplitudes only, since from the definition of H (inv) ab one sees that second time derivatives of metric components can appear. However, it turns out that only the two-form ω ab = p b|a − p a|b , which contains no second time derivatives of h, q, K and G, appears in the Zerilli one-form. Using the fact that H (inv) ab is traceless as a consequence of the field equations, one obtains Now, using (54), it is easy to find Using the background equation v b|a = M r −3g ab − r −2 r |a r |b , one eventually obtains The Zerilli scalar ζ (and Ψ) can now be obtained from its definition, eq.(25), with and the one-form Z given by equations (57,58). On the other hand, the latter equations also give usζ fromζ = Z 0 + βZ x . Note that -as in the odd-parity sector -the scalars ζ andζ do not depend on the perturbed lapse nor on the perturbed shift. Finally, we see that for a Schwarzschild slicing where ∂ 0 r = 0, Φ andζ are linear combinations of the extrinsic curvature components only. These combinations precisely agree with the ones obtained in a perturbative approach on a static background in terms of curvature-based quantities [28].

The ADM quantities from the potentials (l ≥ 2)
If ζ andζ are known, equation (25) tells us how to obtain K (inv) andK (inv) . Next, the traceless part of H Here, p a and G parameterize the gauge freedom. The amplitudes π h , ... π G are obtained from this and (54).
3. The special case l = 1 For l = 1, according to the analysis of the last section, we have only gauge modes, and G and π G are absent. Therefore, the ADM-based amplitudes are obtained from the same equations as above, but where H (inv) ab and K (inv) are set to zero, and p a and G are replaced by ξ a and 2f , respectively, where ξ a and f parameterize the gauge transformation that brings us from the RW gauge to the actual gauge that one wants to use.

The special case l = 0
For l = 0, b 2 , q and π q are also absent. Using equation (33) and the relations (54), the perturbed mass parameter is found to be In order to obtain the perturbed three-metric and extrinsic curvature in terms of δM , one uses equation (32) which gives H ab = 4δM rN 2 r |a r |b − N 2g ab + gauge, and the ADM quantities are obtained in a similar way to above.

C. Gauge fixing vs choices of lapse and shift
We have shown above how to construct the three metric and extrinsic curvature from the potentials (and vice-versa), up to gauge freedom. In numerical simulations, however, usually one does not fix the gauge but rather chooses lapse and shift, perhaps as prescribed functions of spacetime ("exact lapse" or "exact shift") or as dynamical quantities coupled to the the three metric and/or extrinsic curvature ("live gauges"). In general this does not fix the gauge completely, which means that we have to relate the gauge freedom to the choice of lapse and shift. The properties of such relations depend on the details of how the lapse and shift are chosen, and it is therefore not possible to give a general discussion. These equations, for example, might be elliptic if some kind of minimal distortion is imposed, hyperbolic as in the case we discuss below, or of some other (perhaps unknown) type.
Here we will concentrate on a specific simple prescription, but it should be clear that other cases can be treated similarly. The case we are going to discuss is exact-coshift, exact-lapse; that is, the lapse and shift covector are arbitrary but a priori given functions on the orbit space.
We start with the odd-parity sector. The perturbed lapse is zero, and, for l ≥ 2, the perturbed coshift is given by the right hand side of equation (47). The function k parameterizes the gauge freedom, and it is thus given by the equation Since b is a given function, this equation can be solved for k, provided we supply initial data. Given any three-metric and extrinsic curvature at t 0 , the initial data for k is given by k(t 0 , x) = h 2 (t 0 , x). The treatment for l = 1 is similar. Now the gauge freedom is controlled by f , which can be related to the coshift by equation (53), rewritten as In the even-parity case with l ≥ 2 the gauge functions p t , p x and G are related to the lapse and shift by evolution equations which are straightforwardly obtained from equations (59,60,61). These evolution equations form a 3 × 3 coupled system, first order in space and in time,Ġ = 0 + l.o. (64) where l.o. stands for lower order terms. Initial data for the system (64,65,66) is given by the three-metric and the extrinsic curvature at some time t 0 and the formulas (55,56) for p t and p x . It is easy to see that equations (64,65,66) constitute a weakly hyperbolic (see, eg, [24]) system if β = 0 and a stricly hyperbolic system otherwise. That is, if these equations are written as u t = Au x + l.o., with u = (G, p x , p t ) T , the matrix A has three different real eigenvalues if β = 0; and a single degenerate real eigenvalue (zero) with only two independent eigenvectors if β = 0. The structure of the equations for l = 1 is the same, replacing p a , G by ξ a , 2f , respectively. Finally, for l = 0 the situation is similar but simpler: G does not appear, and the principal part of the evolution equations for two gauge quantities p x and p t is also given by (65, 66). As before, these equations are weakly hyperbolic if the background shift is zero, and strictly hyperbolic otherwise. If we use densitized lapse, as is usually done in hyperbolic formulations (see, e.g., [25]), the above system of equations is strongly hyperbolic even in the case where β = 0. In contrast to this, the system is ill posed if we use exact shift instead of coshift. For l = 0, this fact has already been noted in [26].

IV. RELATION TO THE TEUKOLSKY FORMALISM
In order to compare our perturbation equations with the Teukolsky equation for a non-rotating background, we introduce a NP null tetrad that is adapted to the spherically symmetric metric (1), i.e. l = l a dx a , n = n a dx a , m = m A dx A , where l and n form a null dyad ofg,g ab = −l a n b − l b n a , and m is a complex one-form such that Here and in the following, a bar denotes complex conjugation. Note that The only non-vanishing NP coefficients are ρ = − 1 r Dr, µ = 1 r ∆r, ǫ = 1 2 dl(l, n) = 1 2 l a n b∇ a l b , γ = 1 2 dn(l, n) = − 1 2 l b n a∇ a n b , Am B∇ AmB is a NP coefficient with respect to the dyad defined bym ≡ 1 r m. Here, D = l a∇ a and ∆ = n a∇ a . We also introduce, for later use, the angular derivative operatorδ =m A∇ A . From the NP vacuum equations (see, e.g. [29]), it then follows that all Weyl scalars but Ψ 2 vanish, Ψ 0 = Ψ 1 = Ψ 3 = Ψ 4 = 0. In terms of the invariant definition of M given in section 2, Ψ 2 can be expressed as In particular, the metric is of type D with repeated principal null vectors aligned with l a and n a .
In what follows, we study the decoupled equation derived by Teukolsky [8] governing linear fluctuations of Ψ 4 on any spherically symmetric vacuum background. To linear order, Ψ 0 and Ψ 4 are invariant with respect to both infinitesimal coordinate transformations and null tetrad rotations. The reason why we focus on Ψ 4 and not Ψ 0 is that we want to study outgoing radiation at null infinity, which is described by Ψ 4 (see [8]). However, by interchanging l with n and m withm in what follows, one easily obtains the corresponding results for Ψ 0 , describing ingoing radiation at the event horizon.
With respect to the chosen null tetrad, the pulsation operator acting on the linearized field Ψ (1) 4 splits into the sum of an orbital and an angular operator, Next, we compute the perturbed Weyl scalar Ψ (1) 4 : With respect to the background metric (1), one obtains Performing the multipole decomposition as described in section II, we obtain Ψ (1) 4 = n a n b∇ a h in the odd-parity sector and Ψ (1) in the even-parity sector. Using the definition of the derivative operatorδ and the NP coefficientα, one can check that in both parity sectors, the angular part is proportional to the spin-weighted spherical harmonics Y lm −2 defined in appendix B. Explicitly, we have where C 2 l = (l − 1)l(l + 1)(l + 2)/4. It remains to express h and H (inv) ab in terms of the scalar fields Φ and ζ: Using the definitions of the RW and Zerilli potentials Φ and ζ, as well as equations (20) and (25), we obtain =ǫ ab∇ a (rΦ), where we recall that a 0 = 6M/r. Eventually, we get which gives Ψ (1) 4 in terms of the RW and Zerilli potentials Φ lm and Ψ lm = ζ lm /(a 0 + λ) introduced in section 2. Here, we have re-introduced the indices l and m. In order for the metric perturbation to be real, we must haveΨ lm = Ψ l−m andΦ lm = Φ l−m (a bar denoting complex conjugation). Equation (69) (or its Fourier transform in time) is some kind of generalization of the Chandrasekar transformation (see [27], also [30]).
Provided that Φ lm and Ψ lm satisfy the RW and Zerilli equations (7) and (28), respectively, and using the vacuum NP equations and the commutation relations (D + 2(s+1)ǫ + q ρ) (∆ + 2sγ + p µ) − (∆ + 2(s−1)γ + p µ) (D + 2sǫ + q ρ) = p + q 2r 2 + 2(s + p + q)Ψ 2 , where s, p and q are arbitrary real numbers, one can show that indeed, Ψ 4 satisfies the Teukolsky equation (67). Finally, the total radiated energy per unit time can be obtained from where asymptotically flat coordinates and an asymptotically flat NP tetrad are chosen, and where u = t − r. In our case Ψ 4 = 0 on the background, so the radiated energy depends only quadratically on Ψ (1) 4 . Since the fields Φ lm and Ψ lm are scalars with respect to the background metricg, we can evaluate (70) using any asymptotically flat coordinates on the background. Using the fact that at infinity, ∆ = 1 2 (∂ t − ∂ r ) + O(r −1 ), a, γ, µ = O(r −1 ), and imposing the outgoing wave conditionΦ lm + Φ ′ lm = 0,Ψ lm + Ψ ′ lm = 0 at infinity, one arrives at (In the derivation, we have also used the orthogonality of the Y lm −2 andΨ lm = Ψ l−m ,Φ lm = Φ l−m .) As a consistency check, it is useful to note that this coincides with the usual well known result for Schwarzschild black holes in Schwarzschild coordinates 5 . Eq.(71), however, holds for any coordinates.

V. THE RW AND ZERILLI EQUATIONS AND NUMERICAL ISSUES
The RW and Zerilli equations are wave equations with exactly the same differential operator (the Laplacian on the orbit space); they differ only in the corresponding potentials, cf. eq. (34). This simplifies their analysis, since properties such as well-posedness do not depend on lower order terms (as the potentials are). We now discuss certain properties of the RW and Zerilli equations. In particular, we wish to note the fact that they are perfectly well defined as long as the background is regular (both in the coordinate and curvature sense). Provided the latter holds there is no pathology in the equations (or the solutions) at, for example, the event horizon.
We start writing these equations explicitly by introducing coordinates. We then express the whole metric as g total µν = g µν + δg µν , with the background metric given by g = (−α 2 + γ 2 β 2 )dt 2 + 2γ 2 βdtdx + γ 2 dx 2 + r 2 (dϑ 2 + sin 2 ϑdφ 2 ) The RW and Zerilli equations areZ where Z denotes either the RW or Zerilli functions. The coefficients c i are c 1 = 2β and the corresponding potentials are To make the hyperbolic character of these equations manifest, we introduce new variables y := Z ′ , w :=Ż, and write it as a first-order system for the "vector" u = (w, y, Z) T , i.e.,u = Au ′ + Bu, where in this case the principal part is In our case, c 2 1 + 4c 2 = 4α 2 γ −2 and, so, the eigenvectors of A are independent provided the background metric is locally well defined. Thus, the system is strongly hyperbolic, which is enough to prove well-posedness for the initial-boundary value problem if one gives boundary data for the characteristic modes that enter the domain [24]. For the close-limit evolution of black holes in horizon-penetrating coordinates, for example, one would put the inner boundary inside the black hole, check that the characteristic modes are indeed leaving the computational domain (i.e. that the eigenvalues of A are positive), and thus not put boundary conditions there ("excision"). At the outer boundary one would typically put zero boundary conditions for the ingoing modes.
One of the additional advantages of having a hyperbolic equation is that one can write codes that can a priori be shown to be convergent [32]. We have indeed written two such codes for the RW and Zerilli equations with an arbitrary background. One of them uses fourth-order centered differences in space and fourth order Runge-Kutta in time. It uses extrapolation at the inner boundary (assumed to be inside the black hole), and gives zero boundary data to the characteristic mode that enters the computational domain at the outer boundary. The other code is second-order; it also uses Runge-Kutta for time integration and centered differencing in space, but now needs some dissipation (one can prove that this scheme is unstable without dissipation, see [32]). In future work we will present numerical details of these codes applied to the close-limit collision of superposed Kerr-Schild and Painlevé-Gullstrand black holes. The perturbation for the four metric with l ≥ 2 is given by (Y φ = ∂ φ Y , etc.) δg ϑφ = k cos ϑY ϑ + sin −1 ϑY φφ − sin ϑY ϑϑ δg ϑt = 1 γα −γ 2 βrΦ + r(γ 2 β 2 − α 2 )Φ ′ + (−α 2 r ′ −ṙβγ 2 + r ′ γ 2 β 2 )Φ + rk − 2kṙ r Y φ sin ϑ δg φφ = 2k [cos ϑY φ − sin ϑY ϑφ ] δg φt = − 1 γα −γ 2 βrΦ + r(γ 2 β 2 − α 2 )Φ ′ + (−α 2 r ′ −ṙβγ 2 + r ′ γ 2 β 2 )Φ + rk − 2kṙ r sin ϑY ϑ It is straightforward to compute the linearized Ricci or Einstein tensor for the above perturbed metric and see that they are indeed annihilated if the master equation (73) holds.
The l = 1 components of the metric, on the other hand, are given by r + 2J γ 2 β(βr ′ −ṙ) − α 2 r ′ γα(2M − r) sin ϑY ϑ and it is also straightforward to check that this linearized metric satisfies the linearized vacuum equations.