# Gauge-Independent Metric Reconstruction of Perturbations of Vacuum
Spherically-Symmetric Spacetimes

Michele Lenzi<EMAIL_ADDRESS>Institut de Ciències de l’Espai (ICE, CSIC),
Campus UAB, Carrer de Can Magrans s/n, 08193 Cerdanyola del Vallès, Spain
Institut d’Estudis Espacials de Catalunya (IEEC), Edifici Nexus, Carrer del
Gran Capità 2-4, despatx 201, 08034 Barcelona, Spain  Carlos F. Sopuerta
<EMAIL_ADDRESS>Institut de Ciències de l’Espai (ICE, CSIC), Campus
UAB, Carrer de Can Magrans s/n, 08193 Cerdanyola del Vallès, Spain  Institut
d’Estudis Espacials de Catalunya (IEEC), Edifici Nexus, Carrer del Gran Capità
2-4, despatx 201, 08034 Barcelona, Spain

###### Abstract

Perturbation theory of vacuum spherically-symmetric spacetimes (including the
cosmological constant) has greatly contributed to the understanding of black
holes, relativistic compact stars and even inhomogeneous cosmological models.
The perturbative equations can be decoupled in terms of (gauge-invariant)
master functions satisfying $1+1$ wave equations. In this work, building on
previous work on the structure of the space of master functions and equations,
we study the reconstruction of the metric perturbations in terms of the master
functions. To that end, we consider the general situation in which the
perturbations are driven by an arbitrary energy-momentum tensor. Then, we
perform the metric reconstruction in a completely general perturbative gauge.
In doing so, we investigate the role of Darboux transformations and Darboux
covariance, responsible for the isospectrality between odd and even parity in
the absence of matter sources and also of the physical equivalence between the
descriptions based on all the possible master equations. We also show that the
metric reconstruction can be carried out in terms of any of the possible
master functions and that the expressions admit an explicitly covariant form.

##  I Introduction

General Relativity (GR) is a classical theory of gravity that can accommodate
all the physical observations/experiments involving gravitational fields,
despite a better theoretical understanding of some phenomena, like dark matter
and dark energy, may be desirable. It also predicts the generation of
gravitational waves (GWs), whose existence has been confirmed by recent
observations made with the ground-based laser-interferometer detectors LIGO
and Virgo [1, 2, 3], and more recently by the LIGO-Virgo-KAGRA scientific
collaboration [4]. This new type of Astronomy, usually known as Gravitational
Wave Astronomy, constitutes one of the best tools we have to test gravity in
general and GR, specially in dynamical situations involving the strongest
possible gravitational field configurations corresponding to objects like
neutron stars and black holes (BHs) (see [5] for tests of GR with the first GW
event ever detected [6]). The expectation is that the precision of these tests
will improve as the sensitivity of the current detectors increases with the
different instrument upgrades. However, we may need to wait for third-
generation ground detectors, as Einstein Telescope [7] or Cosmic Explorer [8],
to be able to make groundbreaking discoveries in fundamental physics with GW
observations. Alternatively, we can also expect revolutionary science from
space-based GW detectors operating in the low-frequency band like LISA [9, 10,
11, 12], which is expected to be launched around 2035.

To a large extent, GW science requires the use of precise theoretical models
of the expected GW signals, the gravitational waveforms (see [13] for the case
of LISA) which, in turn, require the development of a number of exact,
approximate, and numerical techniques. Of particular relevance is relativistic
perturbation theory, which provides a good theoretical basis for the
description of the generation and propagation of GWs as well as for the
computation of associated physical quantities like GW energy and angular
momentum fluxes (see [14] for textbook accounts). When applied to BHs, the so-
called black hole perturbation theory (BHPT) can successfully describe
important physical phenomena such as scattering processes (we can compute the
black hole greybody factors) and their quasinormal oscillations (the BH
_ringdown_ in the context of the coalescence of binary BH systems). See [15,
16, 17, 18] for BHPT formalisms (and [19, 20] for D-dimensional maximally-
symmetric background spacetimes) and [21, 22, 23, 24] for reviews describing
perturbations of BHs and neutron stars. Applications of BHPT to gravitational
wave astronomy can be found in [25, 26, 27], and to fundamental physics in
[28, 29, 30].

BHPT is also the main tool to describe one of the most important sources of
GWs for space-based observatories like LISA, the so-called extreme-mass-ratio
inspirals (EMRIs), where a stellar-mass compact object is captured by a
supermassive BH and performs a long inspiral before it plunges into the big BH
(see [31, 32, 33, 34, 35] for details on the physics of and science from
EMRIs). The inspiral is driven by gravitational backreaction which, in the
context of BHPT, is described by the so-called self-force [36] (see [37, 38,
39, 40] for accounts on the self-force programme). The computation of the
self-force requires regularization (it becomes singular at the stellar-mass
compact object location) of the metric perturbations as it is made out of
their gradients and evaluated at the particle location. It also requires, in
principle, the knowledge of all the components of the metric perturbations.
This means that we need to have a good understanding of BHPT.

The first steps in BHPT for non-spinning BHs were given in the 1950s by Regge
and Wheeler [41]. Using the spherical symmetry of the BH background, they
decomposed the perturbations in scalar, vector, and tensor harmonics, and
showed that the harmonics decouple, and also that modes with different parity
decouple. Then, they managed to study the odd-parity perturbations of a
Schwarzschild BH [42] and were able to decouple the perturbative Einstein
equations for each harmonic mode in terms of a single complex function, the
_master_ function, which satisfies a wave type equation with a potential that
represents the response of the BH to the perturbations. It turns out that the
master functions carries the true degrees of freedom of the gravitational
field and hence, we can extract from it all the physical information encoded
in the perturbations, in particular the gravitational radiation emitted and
its energy-momentum content. As a consequence, it has been shown that master
functions are gauge invariant. The even-parity sector of the metric
perturbations of Schwarzschild was not put in the same status until the 1970s,
when Zerilli [43, 44] managed to decouple the perturbative equations and find
master functions and master equations with a different potential (see also
[45] and also [46, 47]).

In the case of spinning BHs, given that the Kerr metric was found in the 1960s
[48], a similar framework as in the non-spinning case was developed in the
1970s from the work of Teukolsky [49, 50]. The main difference is that while
in the Schwarzschild case the master functions are functions of the metric
perturbations and its first-order derivatives, in the Kerr case, they are
given in terms of gauge-invariant components of the perturbative Weyl
curvature tensor, which are expressed in a compact way if one uses the Newman-
Penrose formalism [51] with a basis of vectors adapted to the main null
principal directions of the Kerr Weyl tensor (it has a special algebraic
structure corresponding to the so-called type D in the Petrov classification
[52]).

Although master functions carry the physical information of the perturbations,
there are situations where we may need to know all the metric perturbations,
independently on whether they are gauge invariant or not. This is the case of
the self-force, which tells us how an object moves around a BH under its own
self-gravity. The self-force is a gauge-dependent quantity, but it is
necessary to derive the perturbed trajectory, from where we can compute gauge-
invariant gravitational wave polarizations (at the next perturbative order).
Therefore, we need in BHPT a procedure to recover the metric perturbations
from the master functions.

The goal of this paper is to generalize existing metric reconstruction
procedures for perturbations of vacuum spherically-symmetric backgrounds [53,
54, 55, 56, 57], including a cosmological constant (throughout the paper we
refer to these spacetimes as $\Lambda$-vacuum spacetimes). In our analysis, we
allow for energy-momentum distributions appearing at the perturbative level
(i.e., not affecting the background). We will provide metric reconstruction
formulae that are independent of the gauge (they are derived in a completely
general gauge). We will also consider the case where there is a general
energy-momentum tensor at the perturbative level (at first-order in the
perturbations), which is the situation we have when computing the
perturbations induced by a point-like object moving in the background geometry
and whose rest-mass can be considered small so that the modified gravitational
field can be treated in a perturbative way. Finally, we also consider the
metric reconstruction procedure in the framework of the symmetries of the
perturbative equations, studied by the authors in previous works. In
particular, we consider here the metric reconstruction procedure in the
context of what we called _Darboux covariance_ , a symmetry in the space of
possible master equations that tells us that the physics they describe is
equivalent.

The starting point of this work is a previous study [58] where we analyzed the
phase space of master functions and equations in $\Lambda$-vacuum background
spacetimes. To that end we assumed: (i) The master functions are linear
combinations of the metric perturbations and its first-order derivatives. (ii)
The coefficients of these linear combinations are time-independent, i.e. they
only depend on the radial area coordinate $r$. (iii) The master functions
satisfy a wave equations of the form:

$\left(\square_{2}-\Omega_{\ell}(r)\right)\Psi_{\ell m}=0\,,$ (1)

where $\square_{2}$ is the d’Alambertian operator constructed from the metric
on the two-dimensional Lorentzian manifold $M^{2}$ orthogonal to the 2-spheres
of symmetry (i.e. the orbits of the rotation group of Killing symmetries of
the background spacetime), and $\Omega(r)$ is essentially the potential which,
in the study of [58] is initially arbitrary and its form is only determined by
the perturbative Einstein equations. It is important to remark that the
operator in (1), acting on the master function $\Psi_{\ell m}$, is constructed
only from the background geometry.

Assumptions (i) and (ii) are shared by the known master functions, in
particular those introduced by Regge-Wheeler and by Zerilli. The first one is
a simplification that limits the space of master functions. The second one is
a consequence of the structure of the master functions: Since they are linear
in the perturbations, the coefficients of the linear combination have to be
functional of the background metric. The outcome of the study of [58] can be
summarized as follows: There are two branches of master equations: (i) The
standard branch, which contains master equations with the known potentials,
i.e. the Regge-Wheeler potential [41] for odd-parity perturbations and the
Zerilli potential [44] for even-parity perturbations. (ii) The Darboux branch,
where there is an infinite set of master equations with new potentials.

Some interesting consequences to remark from the study are: (i) All the
resulting master functions turn out to be gauge-invariant despite it was not
imposed. This is consistent with the fact that master functions encode the
true degrees of freedom of the gravitational perturbations. (ii) All the
master functions admit a fully covariant form with respect to the two-
dimensional Lorentzian metric. (iii) In the same way that the Regge-Wheeler
and Zerilli potentials coincide for the case of a maximally-symmetric
background, the equations for the potentials in the second branch also
coincide.

On the other hand, in [59] it was shown that, for a given harmonic, all the
master functions are physically equivalently, independently of the parity. The
main reason behind this it that all the master equations/functions can be
related by means of Darboux transformations (see also [60]). A direct
consequence is that Darboux transformations preserves the spectrum, that is,
they are isospectral transformations. More specifically, Darboux _covariance_
implies that all these master equations lead to the same reflection and
transmission coefficients, which is generalization of the result of Refs. [61,
62]. In [59] it was also shown that there is another type of isospectral
symmetries of the master equations, actually an infinite number of them,
associated with deformations along the flow of the Korteweg-de Vries (KdV)
equation [63]. The associated conservation laws produce an infinite set of
conserved quantities [64, 65, 66], the KdV integrals. These integrals were
also shown to be invariant under Darboux transformations. In [67], it was
shown that one can determine the reflection and transmission coefficients for
BH scattering only in terms of the KdV integrals via a moment problem. In
[68], practical computations are described for the case of the Schwarzschild
potentials and also for a Pöschl-Teller potential.

In this paper, we consider the problem of reconstructing all the metric
perturbations in terms of the solutions of the master equation in the case of
spherically-symmetric $\Lambda$-vacuum backgrounds spacetimes, without fixing
the perturbative gauge, i.e. in a gauge-independent way111Note the difference
between gauge independence and gauge invariance. The first one referes to
perturbative computations in a generic gauge, while the second is a property
of a perturbative object which does not change its functional form in terms of
metric perturbations under changes in the perturbative gauge.. This
reconstruction process is of great importance for the computation of the self-
force [37, 38, 39], and it can be very helpful for any other application
requiring the knowledge of the metric perturbations beyond the gauge-invariant
sector. This includes Hamiltonian formulations of the perturbations (see [69]
for a treatment of the odd-parity sector), in particular towards quantization
in a semiclassical context (see, e.g. [70, 71, 72]). In our metric
reconstruction we have included a generic energy-momentum tensor sourcing the
perturbations, generalizing in this way several results in our previous papers
[58, 59]. We also show a systematic way of carrying the reconstruction, which
may be useful for generalizing this to other situations where we may have
different (spherically-symmetric) background (for instance, for the treatment
of compact relativistic stars), or even a different theory of gravity.
Finally, we put the metric reconstruction process in the context of the hidden
symmetry of the space of master functions and equations studied in [59].
Actually, we explain how to realize the metric reconstruction in terms of
master functions different from the ones known in the broad literature on the
subject.

Plan of the paper. In Sec. II we introduce the main elements of BHPT that use
in this work: The perturbative Einstein equations; the background spacetimes;
the harmonic decomposition; and gauge invariant quantities. In Sec. III we
describe the space of master functions and equations introduced in [58], and
whose structure was studied in [59], extending the results to the case in
which we have an arbitrary energy-momentum tensor at the first perturbative
order. We also extend the study of Darboux covariance clarifying some open
questions in [59]. In Sec. IV, we describe a systematic procedure to carry out
the metric reconstruction and provide general results taking in account the
full space of master functions introduced previously. In Sec. V we discuss
these results, compare them with other formalisms in General Relativity and
consider possible applications and extensions to other situations. We have the
following appendices: In Appendix A we briefly discuss the family of
background spacetime geometries. In Appendix B we list some properties of
scalar, vector, and tensor spherical harmonics that are relevant for this
paper. In Appendix C, we list Einstein perturbative equations, including a
generic energy-momentum tensor contribution at first-order, both in terms of
the original metric perturbations and in terms of the gauge invariant metric
perturbations. Finally, in Appendix D we list the harmonic components of the
energy-momentum conservation equations.

Conventions. Through this paper we use units in which $c=8\pi G=1$, where $c$
is the speed of light and $G$ the Newtonian gravitational constant. We use a
definition of the Einstein tensor that includes the cosmological constant term
[see Eq. (2)]. Regarding indices: Greek letters are used for spacetime
indices; capital Latin letters are used for indices in the 2-sphere ($S^{2}$);
small-case Latin letters are used for indices in the Lorentzian two-
dimensional manifold $M^{2}$. The metric on $S^{2}$ is written as
$\Omega_{AB}$ while in $M^{2}$ it is $g_{ab}$. Covariant differentiation in
the four-dimensional background spacetime is either denoted by a semicolon or
by the nabla operator $\widehat{\nabla}$. Instead, in $S^{2}$ is denoted using
a vertical bar ($\Omega_{AB|C}=0$) and in $M^{2}$ by a colon ($g_{ab:c}=0$).
The antisymmetric covariant unit tensor associated with the volume form (Levi-
Civita tensor) in $S^{2}$ is denoted by $\epsilon_{AB}$, and the corresponding
one on the Lorentzian manifold $M^{2}$ is denoted by $\varepsilon_{ab}\,$.

##  II Black Hole Perturbation Theory in the non-rotating case

BHPT is a particular case of relativistic perturbation theory (see [73, 74]),
in which a physically realistic system can be described as a deviation from an
idealized situation. We assume the existence of a one-parameter family of
spacetimes, $(\mathcal{M}_{\lambda},g_{\lambda})$, with $\lambda=0$
corresponding to the background spacetime describing the idealized situation,
and construct the physical spacetime by Taylor expanding around $\lambda=0$.
In this way, the parameter $\lambda$ formally controls the magnitude of the
perturbations (see e.g. [73, 75, 76, 77]). Since in our case $\lambda$ is a
dummy parameter we ignore it through the paper.

In our case, the background spacetime is taken to be a $\Lambda$-vacuum
spherically symmetric spacetime. Then, the background spacetime metric,
$\widehat{g}_{\mu\nu}$, satisfies the following version of Einstein’s field
equations 222We use a hat to denote quantities associated with the background
spacetime, like $\widehat{Q}$.:

$\displaystyle\widehat{G}_{\mu\nu}=\widehat{R}_{\mu\nu}-\frac{1}{2}\widehat{g}_{\mu\nu}\widehat{R}+\Lambda\,\widehat{g}_{\mu\nu}=0\,,$
(2)

where $\widehat{R}_{\mu\nu}$ and $\widehat{G}_{\mu\nu}$ denote the Ricci and
Einstein tensors of the background metric respectively,
$\widehat{R}=\widehat{g}^{\mu\nu}\widehat{R}_{\mu\nu}$ is the background
scalar curvature, and $\Lambda$ is the cosmological constant. Since the
background is $\Lambda$-vacuum we have from Eq. (2):

$\widehat{R}_{\mu\nu}=\Lambda\widehat{g}_{\mu\nu}\,,\quad\widehat{R}=4\Lambda\,.$
(3)

After a correspondence between the physical (perturbed) and background
spacetimes is established (see Sec. II.3), we can write the physical spacetime
metric $g_{\mu\nu}$ as:

$g_{\mu\nu}=\widehat{g}_{\mu\nu}+h_{\mu\nu}\,,$ (4)

where $h_{\mu\nu}$ ($|h_{\mu\nu}|\ll|\widehat{g}_{\mu\nu}|$) are the metric
perturbations.

Let us now fix the notation we adopt. Given any quantity $Q$ in the perturbed
spacetime, a $\delta$ in front of it denotes its perturbative part, i.e.
$\delta Q=Q-\widehat{Q}$. For the example, in the case of the metric we have:
$h_{\mu\nu}=\delta g_{\mu\nu}=g_{\mu\nu}-\widehat{g}_{\mu\nu}$. One can then
find that the perturbation of the Christoffel symbols at first order can be
expressed in terms of the background covariant derivatives (denoted here by a
semicolon) of the metric perturbations as follows:

$\displaystyle\delta\Gamma_{\mu\nu}^{\rho}=\frac{1}{2}\widehat{g}^{\rho\sigma}\left(h_{\mu\sigma;\nu}+h_{\nu\sigma;\mu}-h_{\mu\nu;\sigma}\right)\,.$
(5)

While the Christoffel symbols of the background are not tensors with respect
to changes of the coordinates of the background spacetime, their perturbations
are since they are the subtraction of two Christoffel symbols. Then, we can
write the perturbations of the Riemann tensor in terms of the perturbed
Christoffel symbols:

$\delta
R^{\mu}_{\nu\rho\sigma}=\delta\Gamma^{\mu}_{\nu\sigma;\rho}-\delta\Gamma^{\mu}_{\nu\rho;\sigma}=2\delta\Gamma^{\mu}_{\nu[\sigma;\rho]}\,,$
(6)

and the perturbations of the Ricci tensor follows directly:

$\displaystyle\delta
R_{\mu\nu}=\delta\Gamma^{\rho}_{\mu\nu;\rho}-\delta\Gamma^{\rho}_{\rho\mu;\nu}\,.$
(7)

The Einstein tensor can be decomposed as:
$G_{\mu\nu}=\widehat{G}_{\mu\nu}+\delta G_{\mu\nu}$, and we can write $\delta
G_{\mu\nu}$ in terms of the previous perturbative objects as follows [see Eq.
(2)]:

$\delta G_{\mu\nu}=\delta R_{\mu\nu}-\frac{1}{2}\widehat{g}_{\mu\nu}\delta
R-\Lambda h_{\mu\nu}\,.$ (8)

Then, the perturbation of the Einstein tensor can finally be written as

$\displaystyle\delta G_{\mu\nu}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\widehat{\square}\,\bar{h}_{\mu\nu}-\widehat{R}^{\rho}_{\mu}{}^{\sigma}_{\nu}\bar{h}_{\rho\sigma}+\widehat{\nabla}_{(\mu}\mathcal{L}_{\nu)}$
(9) $\displaystyle-$
$\displaystyle\frac{1}{2}\widehat{g}_{\mu\nu}\left(\widehat{\nabla}_{\rho}\mathcal{L}^{\rho}\right)-\Lambda\bar{h}_{\mu\nu}\,,$

where $\widehat{R}^{\rho}_{\mu}{}^{\sigma}_{\nu}$ is the background Riemann
tensor, and we have introduced the trace-reversed metric perturbations:

$\bar{h}_{\mu\nu}=h_{\mu\nu}-\frac{1}{2}\widehat{g}_{\mu\nu}h\,,$ (10)

with $h$ being the trace of $h_{\mu\nu}$ with respect to the background
metric:
$h=\widehat{g}^{\mu\nu}h_{\mu\nu}=-\widehat{g}^{\mu\nu}\bar{h}_{\mu\nu}$. The
operator $\widehat{\square}$ is the d’Alambertian associated with the
background metric, thus:

$\widehat{\square}\,\bar{h}_{\mu\nu}=\bar{h}_{\mu\nu;\rho}{}^{;\rho}\,.$ (11)

Finally, in Eq. (9) we have introduced the quantity:

$\mathcal{L}_{\mu}=\widehat{g}^{\rho\sigma}\widehat{\nabla}_{\rho}\bar{h}_{\sigma\mu}\,,$
(12)

to isolate terms which vanish in the Lorentz gauge, i.e.
$\mathcal{L}_{\mu}=0$. Nevertheless, in this paper we work in a completely
general gauge.

In the following, are going to consider the situation in which the
perturbations are sourced by matter fields whose nature we do not specify and
hence, it is arbitrary. In other words, we consider a completely general
energy-momentum tensor $T_{\mu\nu}$, which only affects the metric
perturbations but not the background. In this way, the perturbative field
equations become:

$\delta G_{\mu\nu}=T_{\mu\nu}\,,$ (13)

with $\delta G_{\mu\nu}$ given by Eq. (9). Until here, we have only imposed
the background to be a $\Lambda$-vacuum solution [Eq. (2)]. In what follows,
we summarize the ingredients involved in the treatment of perturbations of
spherically-symmetric background spacetimes.

###  II.1 Structure of the Background Spacetime

Our background spacetimes are spherically-symmetric solutions of the Einstein
vacuum equations with a cosmological constant [see Eq. (2) and Appendix A].
They have a special geometric structure as they are the _warped_ product of
two manifolds: $M^{2}\times_{r}S^{2}$, where $S^{2}$ is the $2$-sphere, $r$ is
the radial area coordinate defining the warp factor, and $M^{2}$ is a two-
dimensional Lorentzian manifold. Therefore, the background metric is the
semidirect product of a Lorentzian metric on $M^{2}$, $g_{ab}$, and the unit
curvature metric on $S^{2}$, $\Omega_{AB}$:

$\displaystyle
ds^{2}=\widehat{g}_{\mu\nu}dx^{\mu}dx^{\nu}=g_{ab}dx^{a}dx^{b}+r^{2}\Omega_{AB}d\Theta^{A}d\Theta^{B}\,,$
(14)

where

$\displaystyle g_{ab}dx^{a}dx^{b}$ $\displaystyle=$
$\displaystyle-f(r)\,dt^{2}+\frac{dr^{2}}{f(r)}\,,$ (15)
$\displaystyle\Omega_{AB}d\Theta^{A}d\Theta^{B}$ $\displaystyle=$
$\displaystyle d\theta^{2}+\sin^{2}\theta d\varphi^{2}\,,$ (16)

where $(x^{a})=(t,r)$ and $(\Theta^{A})=(\theta,\varphi)$ are coordinates on
$M^{2}$ (Schwarzschild coordinates) and $S^{2}$ respectively, and $f(r)$ is a
function parameterising time translations (related to the redshift). Using the
background Einstein’s equations (2) one derive ordinary differential equations
(ODEs) for $f(r)$:

$\displaystyle rf^{\prime}+f+\Lambda r^{2}-1=0\,,$ (17) $\displaystyle
r\left(f^{\prime\prime}+2\Lambda\right)+2f^{\prime}=0\,.$ (18)

From these equations we can write the cosmological constant in terms of $f(r)$
and its derivatives in two different forms:

$\displaystyle\Lambda$ $\displaystyle=$
$\displaystyle\frac{1}{r^{2}}\left[1-\left(rf\right)^{\prime}\right]=-\frac{1}{2r^{2}}\left(r^{2}f^{\prime}\right)^{\prime}\,.$
(19)

Now, combining this expression for $\Lambda$ with Eq. (17) we can write

$\displaystyle\left[r\left(1-f-\frac{\Lambda}{3}r^{2}\right)\right]^{\prime}=0\quad\Rightarrow$
$\displaystyle\Rightarrow\quad
r\left(1-f-\frac{\Lambda}{3}r^{2}\right)=2M=\text{const.}\,,$ (20)

where $M$ is the spacetime mass. Actually, from this equation we obtain the
following well-known expression for $f(r)$:

$f(r)=1-\frac{2M}{r}-\frac{\Lambda}{3}r^{2}\,.$ (21)

See Appendix A for the different well-known metrics included.

On the Lorentzian two-dimensional spacetime $M^{2}$ we can introduce a basis
of vectors that is going to be useful in order to express some of the results
we are presenting in this work:

$\displaystyle r_{a}=r_{:a}\,,\qquad t^{a}=-\varepsilon^{ab}r_{b}\,.$ (22)

The first one, $r_{a}$, is the gradient of the area radial coordinate, an
invariant of the rotation group of symmetries. The second one is a timelike
Killing vector, which makes the spacetime static (it is irrotational and
hence, it generates orthogonal hypersurfaces). Moreover, they are orthogonal:

$t^{a}r_{a}=0\,,$ (23)

and their norms are:

$\displaystyle r^{a}r_{a}=f\,,\qquad t^{a}t_{a}=-f\,.$ (24)

On the other hand, it is convenient for future purposes to introduce the
d’Alambertian operator associated with the metric of $M^{2}$, i.e. $g_{ab}$:

$\square_{2}\Psi\equiv
g^{ab}\Psi_{:ab}=\frac{1}{\sqrt{-g}}\partial_{a}\left(\sqrt{-g}g^{ab}\partial_{b}\Psi\right)\,.$
(25)

In Schwarzschild coordinates we have the following expression for
$\square_{2}$:

$\square_{2}\phi=-\frac{1}{f}\frac{\partial^{2}\phi}{\partial
t^{2}}+\frac{\partial}{\partial r}\left(f\frac{\partial\phi}{\partial
r}\right)\,.$ (26)

###  II.2 Expansion of the Perturbations in Spherical Harmonics

The warped geometry of the background leads to the separation of solutions of
the wave equations whose operator is built with the background metric. In
particular, we can separate the dependence on the coordinates of $M^{2}$ from
the angular dependence on $S^{2}$ for the perturbative Einstein equations
(13). In practice, this is done by expanding the metric perturbations using
the eigenfunctions of the Laplacian, the scalar spherical harmonics, and the
vector and tensor objects we derive from them. In turn, we can distinguish
between even- and odd-parity harmonics depending on how they transform under
parity transformation given by: $(\theta,\phi)$ $\rightarrow$
$(\pi-\theta,\phi+\pi)$. If a given harmonic object $\mathcal{O}^{\ell m}$
transforms as: $\mathcal{O}^{\ell m}\rightarrow(-1)^{\ell}\mathcal{O}^{\ell
m}$ it is said to be of the even-parity type; while if it transforms as
$\mathcal{O}^{\ell m}\rightarrow(-1)^{\ell+1}\mathcal{O}^{\ell m}$ it is said
to be of the odd-parity type.

With this in mind, the scalar, vector and tensor spherical harmonics are:
Scalar spherical harmonics: $Y^{\ell m}$. As we have just mentioned, they are
eigenfunctions of the Laplace operator on $S^{2}$ (see Appendix B):

$\displaystyle\Omega^{AB}Y^{\ell m}_{|AB}=-\ell(\ell+1)Y^{\ell m}\,.$ (27)

Vector Spherical Harmonics: They are defined for $\ell\geqslant 1$, are given
by:

$\displaystyle Y^{\ell m}_{A}\equiv Y^{\ell m}_{|A}\qquad~{}$
$\displaystyle\quad\text{Even (polar) parity}\,,$ (28) $\displaystyle X^{\ell
m}_{A}\equiv-\epsilon_{A}{}^{B}Y^{\ell m}_{B}$ $\displaystyle\quad\text{Odd
(axial) parity}\,.$ (29)

Tensor Spherical Harmonics: The basis of symmetric $2$nd-rank tensor spherical
harmonics, which are defined for $\ell\geqslant 2$, are given by

$\displaystyle T_{AB}^{\ell m}\equiv Y^{\ell
m}\,\Omega_{AB}\qquad\qquad\qquad\quad$ $\displaystyle\quad\text{Even
parity}\,,$ (30) $\displaystyle Y^{\ell m}_{AB}\equiv Y_{|AB}^{\ell
m}+\frac{\ell(\ell+1)}{2}Y^{\ell m}\,\Omega_{AB}$
$\displaystyle\quad\text{Even parity}\,,$ (31) $\displaystyle X_{AB}^{\ell
m}\equiv X^{\ell m}_{(A|B)}\qquad\qquad\qquad\qquad$
$\displaystyle\quad\text{Odd parity}\,,$ (32)

with the following traces

$\Omega^{AB}T^{\ell m}_{AB}=2Y^{\ell m}\,,\quad\Omega^{AB}Y^{\ell
m}_{AB}=\Omega^{AB}X^{\ell m}_{AB}=0\,.$ (33)

Some orthogonality and differential properties of these spherical harmonics,
necessary to manipulate and separate the perturbed Einstein equations, can be
found in Appendix B.

The expansion of the metric perturbations in scalar, vector, and tensor
spherical harmonics leads to the decoupling of the equations, not only for
each harmonic $(\ell,m)$, but also for modes with different parity [15, 16].
Then, we write

$\displaystyle h_{\mu\nu}=\sum_{\ell,m}h^{\ell
m,\mathrm{odd}}_{\mu\nu}+h^{\ell m,\mathrm{even}}_{\mu\nu}\,,$ (34)

with

$h^{\ell m,\mathrm{odd}}_{\mu\nu}=\begin{pmatrix}0&h_{a}^{\ell m}X^{\ell
m}_{A}\\\\[5.69054pt] \ast&h_{2}^{\ell m}X^{\ell m}_{AB}\end{pmatrix}\,,$ (35)

and

$h^{\ell m,\mathrm{even}}_{\mu\nu}=\begin{pmatrix}h_{ab}^{\ell m}Y^{\ell
m}&\mbox{\rm\j}^{\ell m}_{a}Y^{\ell m}_{A}\\\\[5.69054pt]
\ast&r^{2}\left(K^{\ell m}T^{\ell m}_{AB}+G^{\ell m}Y^{\ell
m}_{AB}\right)\end{pmatrix}\,,$ (36)

where the asterisk tells us that the corresponding component is obtain form
the fact that the metric perturbations is a symmetric tensorial object.
Moreover, $K^{\ell m}$, $G^{\ell m}$ and $h_{2}^{\ell m}$ denote the scalar
perturbations; $h_{a}^{\ell m}$ and $j^{\ell m}_{a}$ the vector perturbations;
and $h_{ab}^{\ell m}$ the tensorial ones. All these quantities depend only on
the coordinates of $M^{2}$, i.e. on $\\{x^{a}\\}$.

Since we are considering the case in which the perturbations are sourced by a
distribution of energy and momentum, we also need to decompose the energy-
momentum tensor in spherical harmonics. We can follow closely the
decomposition of the metric perturbations as both are symmetric tensor. First,
we separate even- and odd-parity harmonics:

$\displaystyle T_{\mu\nu}=\sum_{\ell,m}T^{\ell
m,\mathrm{odd}}_{\mu\nu}+T^{\ell m,\mathrm{even}}_{\mu\nu}\,,$ (37)

The odd-parity sector has the following structure

$T^{\ell m,\mathrm{odd}}_{\mu\nu}=\begin{pmatrix}0&\mathcal{S}_{a}^{\ell
m}X^{\ell m}_{A}\\\\[5.69054pt] \ast&\mathcal{S}^{\ell m}X^{\ell
m}_{AB}\end{pmatrix}\,,$ (38)

where the odd-parity harmonic modes are defined as follows:

$\displaystyle\mathcal{S}^{a}_{\ell m}$ $\displaystyle=$
$\displaystyle\frac{r^{2}}{\ell(\ell+1)}\int_{S^{2}}d\Omega\,T^{aA}\bar{X}_{A}^{\ell
m}\,,$ (39) $\displaystyle\mathcal{S}_{\ell m}$ $\displaystyle=$
$\displaystyle
2r^{4}\frac{(\ell-2)!}{(\ell+2)!}\int_{S^{2}}d\Omega\,T^{AB}\bar{X}_{AB}^{\ell
m}\,,$ (40)

where $d\Omega$ is the two-sphere _volume_ element. The even-parity sector has
a different structure given by:

$T^{\ell m,\mathrm{even}}_{\mu\nu}=\begin{pmatrix}\mathcal{Q}_{ab}^{\ell
m}Y^{\ell m}&\mathcal{P}^{\ell m}_{a}Y^{\ell m}_{A}\\\\[5.69054pt]
\ast&r^{2}\left(\mathcal{T}^{\ell m}T^{\ell m}_{AB}+\mathcal{P}^{\ell
m}Y^{\ell m}_{AB}\right)\end{pmatrix}\,,$ (41)

and the harmonic components read:

$\displaystyle\mathcal{Q}^{ab}_{\ell m}$ $\displaystyle=$
$\displaystyle\int_{S^{2}}d\Omega\,T^{ab}\bar{Y}^{\ell m}\,,$ (42)
$\displaystyle\mathcal{P}_{\ell m}^{a}$ $\displaystyle=$
$\displaystyle\frac{r^{2}}{\ell(\ell+1)}\int_{S^{2}}d\Omega\,T^{aA}\bar{Y}_{A}^{\ell
m}\,,$ (43) $\displaystyle\mathcal{P}_{\ell m}$ $\displaystyle=$
$\displaystyle
2r^{2}\frac{(\ell-2)!}{(\ell+2)!}\int_{S^{2}}d\Omega\,T^{AB}\bar{Y}_{AB}^{\ell
m}\,,$ (44) $\displaystyle\mathcal{T}_{\ell m}$ $\displaystyle=$
$\displaystyle\frac{r^{2}}{2}\int_{S^{2}}d\Omega\,T^{AB}\bar{T}_{AB}^{\ell
m}\,.$ (45)

In summary, $(\mathcal{P}^{\ell m}(x^{a}),\mathcal{T}^{\ell
m}(x^{a}),\mathcal{S}^{\ell m}(x^{a}))$ are the scalar (polar) harmonics;
$(\mathcal{P}^{\ell m}_{a}(x^{b}),\mathcal{S}^{\ell m}_{a}(x^{b}))$ are the
vector harmonics; and $\mathcal{Q}_{ab}^{\ell m}(x^{c})$ are the tensor ones.

We can now insert the harmonic decomposition of the metric perturbations and
of the energy-momentum tensor into the perturbative Einstein equations (13) to
obtain the decoupled equations for each harmonic and parity (see Appendix C).
In addition, it is important to take into account the energy-momentum tensor
conservation equations, which are a consequence of the perturbative second
Bianchi identities:

$\widehat{\nabla}^{\mu}\delta
G_{\mu\nu}=0\quad\Rightarrow\quad\widehat{\nabla}^{\mu}T_{\mu\nu}=0\,.$ (46)

That is, the local energy-momentum conservation law is with respect to the
background Levi-Civita connection. Introducing here the harmonic decomposition
of $h_{\mu\nu}$, $T_{\mu\nu}$, and derived quantities, we obtain the decoupled
equations for each harmonic and parity (see Appendix D).

###  II.3 Gauge Invariance

A key relation in relativistic perturbation theory [73, 75] is represented by
Eq. (4), where the perturbed metric is written as the sum of a background
metric and the perturbations. As we already mentioned in Sec. II, this
relation requires the prescription of a correspondence (a diffeomorphism)
between the background and the physical spacetimes. Given one of the possible
(infinite) correspondences, we can pull back the physical metric and related
tensorial structure into the background tensorial structure. Therefore, the
relation established by Eq. (4) depends on the choice of correspondence.
Fixing the correspondence between background and physical spacetimes is what
we refer to as fixing the perturbative gauge. To understand this better,
consider a point $p$ in the physical spacetime and two different
diffeomorphisms to the background (two different gauges), say $\phi_{1}$ and
$\phi_{2}$, such that $q_{1}=\phi_{1}(p)$ and $q_{2}=\phi_{2}(p)$ are two
points in the background spacetime (in principle they are different unless
$\phi_{1}=\phi_{2}$). Since they are diffeomorphisms, we can write
$q_{2}=\phi_{2}\circ\phi^{-1}_{1}(q_{1})$. The composed mapping,
$\Phi=\phi_{2}\circ\phi^{-1}_{1}$ is what we call a gauge transformation, from
gauge $\phi_{1}$ to gauge $\phi_{2}$, and acts by changing the reference point
in the background to which we associate the perturbation. In perturbation
theory we assume that the gauge tranformation induces a motion on the
background point such that, when written in a given coordinate system, the
coordinate version of the gauge transformation can be Taylor-expanded in
$\lambda$ (see the beginning of Sec. II), so that it looks, at first-order in
$\lambda$, as follows (see, e.g. [75]):

$x^{\mu}\quad\longrightarrow\quad x^{\prime\mu}={x}^{\mu}+\lambda\xi^{\mu}\,,$
(47)

where ${x}^{\mu}$ and $x^{\prime\mu}$ are the coordinates of two points
$q_{2}$ and $q_{1}$ respectively, and the vector field $\xi^{\mu}$ is the
first-order generator of the gauge transformation (see [76, 77] for the
analogous equation in multi-parameter perturbation theory). In what follows,
we absorb the perturbative parameter $\lambda$ into the generator $\xi^{\mu}$,
and assume it is small in the same way as we do with the metric perturbations,
i.e. $|\xi^{\mu}|\ll|\widehat{g}_{\mu\nu}|$.

The gauge transformation in Eq. (47) generates the following transformation of
the metric perturbations:

$h_{\mu\nu}\quad\longrightarrow\quad{h^{\prime}}_{\mu\nu}=h_{\mu\nu}-2\,\xi_{(\mu;\nu)}\,.$
(48)

Given that we have split the metric perturbations and energy-momentum tensor
in harmonics, we should also split the gauge transformation generator
$\xi^{\mu}$. The even-parity sector harmonic $(\ell,m)$ of the gauge generator
has the form:

$\xi^{\ell m,\mathrm{even}}_{\mu}=\left(\alpha^{\ell m}_{a}(x^{b})Y^{\ell
m}\,,r^{2}\beta^{\ell m}(x^{a})Y^{\ell m}_{A}\right)\,,$ (49)

and the corresponding odd-parity harmonic mode

$\xi^{\ell m,\mathrm{odd}}_{\mu}=\left(0\,,r^{2}\gamma^{\ell m}(x^{a})X^{\ell
m}_{A}\right)\,.$ (50)

Note that there are three gauge functions for even-parity perturbations and
one for the odd-parity ones.

Introducing the multipolar decomposition of the metric perturbations and the
gauge vector into Eq. (48), we obtain the transformation rules for even-parity
metric perturbations under general perturbative gauge transformations:

$\displaystyle h^{\prime\ell m}_{ab}$ $\displaystyle=$ $\displaystyle h^{\ell
m}_{ab}-2\,\alpha^{\ell m}_{(a:b)}\,,$ (51)
$\displaystyle\mbox{\rm\j}^{\prime\ell m}_{a}$ $\displaystyle=$
$\displaystyle\mbox{\rm\j}^{\ell m}_{a}-\left(\alpha^{\ell
m}_{a}+r^{2}\beta^{\ell m}_{:a}\right)\,,$ (52) $\displaystyle K^{\prime\ell
m}$ $\displaystyle=$ $\displaystyle K^{\ell m}+\ell(\ell+1)\beta^{\ell
m}-2\frac{r^{:a}}{r}\alpha^{\ell m}_{a}\,,$ (53) $\displaystyle G^{\prime\ell
m}$ $\displaystyle=$ $\displaystyle G^{\ell m}-2\beta^{\ell m}\,.$ (54)

and the ones for odd-parity metric perturbations:

$\displaystyle h^{\prime\ell m}_{a}$ $\displaystyle=$ $\displaystyle h^{\ell
m}_{a}-r^{2}\gamma^{\ell m}_{:a}\,,$ (55) $\displaystyle h^{\prime\ell m}_{2}$
$\displaystyle=$ $\displaystyle h^{\ell m}_{2}-2\,r^{2}\gamma^{\ell m}\,.$
(56)

It turns out that there are combinations of the metric perturbations and its
derivatives that are invariant under gauge transformations. In the case of
even-parity metric perturbations there are four independent gauge-invariant
quantities, which can be written as:

$\displaystyle\tilde{h}_{ab}$ $\displaystyle=$ $\displaystyle
h_{ab}-\kappa_{a:b}-\kappa_{b:a}\,,$ (57) $\displaystyle\tilde{K}$
$\displaystyle=$ $\displaystyle
K+\frac{\ell(\ell+1)}{2}G-2\frac{r^{a}}{r}\kappa_{a}\,,$ (58)

where

$\displaystyle\kappa_{a}$ $\displaystyle=$
$\displaystyle\mbox{\rm\j}_{a}-\frac{r^{2}}{2}G_{:a}\,,\quad
r_{a}=r_{:a}~{}\Rightarrow~{}r^{a}=g^{ab}r_{b}\,.$ (59)

In the case of odd-parity metric perturbations there are two independent
gauge-invariant quantities:

$\tilde{h}_{a}=h_{a}-\frac{1}{2}h_{2:a}+\frac{r_{a}}{r}h_{2}\,,$ (60)

We can also analyze the changes in the components of the energy-momentum
tensor $T_{\mu\nu}$ under the gauge transformation of Eq. (47). In analogy
with what happens with the components of the metric perturbations [see Eq.
(48)], we can show that the components of $T_{\mu\nu}$ transform as follows:

$T_{\mu\nu}\quad\longrightarrow\quad{T^{\prime}}_{\mu\nu}=T_{\mu\nu}-\pounds_{\xi}\widehat{T}_{\mu\nu}\,,$
(61)

where $\pounds_{\xi}$ denotes Lie differentiation (see, e.g. [14, 78]) with
respect to the vector field generating the gauge transformation, $\xi^{\mu}$,
and $\widehat{T}_{\mu\nu}$ is the background energy-momentum tensor, which in
our case is identically zero. And hence it follows that the components of the
energy-momentum tensor are gauge-invariant. This is exactly what the well-
known Stewart-Walker lemma [73] tell us, and we could have just invoked it.
But this discussion gives more relevance to the fact that we introduce the
energy-momentum distribution already at first perturbative order.

##  III The Space of Master Functions and Equations

The basic idea leading to master equations and master functions is to find
linear functionals of the metric perturbations (for each harmonic and parity
mode with $\ell\geq 2$), the master functions, containing only the gauge-
invariant metric perturbations themselves and their derivatives

$\Psi_{\ell m}=\Psi_{\ell m}[r,\vec{H}^{\ell m},\partial_{a}\vec{H}^{\ell
m},\partial^{2}_{ab}\vec{H}^{\ell m}\,,\cdots]\,,$ (62)

where

$\vec{H}^{\ell m}=(\tilde{h}^{\ell m}_{ab},\tilde{K}^{\ell m};\tilde{h}^{\ell
m}_{a})\,,$ (63)

so that the (vacuum) perturbative Einstein equations reduce to wave equations
for the master functions of the form:

$\left(\square_{2}-\frac{1}{f}V_{\ell}\right)\Psi_{\ell m}=0\,,$ (64)

where $V_{\ell}(r)$ is the $\ell$-dependent potential, which will be in
general different for even- and odd-parity modes. From a more physical point
of view, the master functions capture the true degrees of freedom of the
gravitational field at the perturbative level, and hence they must be gauge
invariant as reflected in the previous equations.

In [58], the full space of master functions and equations, for
$\Lambda$-vacuum spherically-symmetric spacetimes, was studied in a systematic
way. All the possible master equations were found by assuming that the master
functions are linear combinations, with coefficients depending only on $r$, of
the metric perturbations and their first-order derivatives. Two branches of
possible pairs of potentials/master functions, $\\{(V,\Psi)\\}$ were
identified:

(i) The standard branch. We call it standard because it contains a single
potential for each parity, and these potentials are the known ones: The Regge-
Wheeler potential, $V_{\mathrm{RW}}$, for odd-parity perturbations [41] and
the Zerilli potential, $V_{\mathrm{Z}}$, for even-parity perturbations [44].

(ii) The Darboux branch. This is a completely new branch for both parities,
where we find an infinite set of potentials and master functions,
$\\{(V,\Psi)\\}$. In particular, the set of allowed potentials is determined
by a non-linear second-order ODE.

There are some relevant comments about this construction: First, in Eq. (62)
there is no explicit dependence in time. As we mentioned before, this is a
consequence of the existence of a timelike Killing vector in the background,
$K^{\mu}=(t^{a},0)$. Second, in the study of [58] we did not impose gauge
invariance [in the ansatz for the master function we included all the metric
perturbations, in contrast to Eq. (62), where only gauge invariant ones are
considered], it came as a logical result of the analysis. We can again invoke
the Stewart-Walker lemma and apply it to the perturbed Einstein equations,
which are gauge invariant as the background Einstein tensor vanishes
identically. Since the master equations are linear combinations of the
perturbed Eintein equations (and their derivatives), they must also be gauge
invariant. And since they only contain the master function, it has to be gauge
invariant too.

The structure of the space of master functions and equations was analyzed in
Ref. [59] (see also [67]), where we found that all the master equations are
connected by _Darboux transformations_ , i.e. a particular set of derivative
transformations preserving the spectrum of the system (see Sec. III.3).
Therefore, there are infinite possible descriptions of the perturbations of
$\Lambda$-vacuum spherically-symmetric spacetimes in terms of master
equations, all of them with the same spectrum. This hidden symmetry of the
master equations has been called Darboux covariance [59], since it actually
represent a covariance in the description of gravitational perturbations.

In this section, we generalize the work of Ref. [58] to include generic matter
sources that enter at first-order, i.e. without affecting the background
spacetime, which is still a $\Lambda$-vacuum solution. The generalization goes
along the lines of the derivations of [58], in the sense that the same
combinations of the perturbed Einstein equations lead to master equations with
the same operators acting on the master functions but with a source term that
is build as linear combinations of the energy-momentum tensor $T_{\mu\nu}$ and
its derivatives. That is, we obtain master equations of the form:

$\left(\square_{2}-\frac{1}{f}V_{\ell}\right)\Psi_{\ell m}=\mathcal{F}_{\ell
m}\,,$ (65)

where the source term $\mathcal{F}_{\ell m}$ has the following structure:

$\mathcal{F}_{\ell m}=\mathcal{F}_{\ell m}[r,g_{ab},\vec{\mathcal{J}}_{\ell
m},\partial_{a}\vec{\mathcal{J}}_{\ell m}]\,,$ (66)

where

$\vec{\mathcal{J}}_{\ell m}=(\mathcal{Q}^{\ell m}_{ab},\mathcal{P}^{\ell
m}_{a},\mathcal{T}_{\ell m},\mathcal{P}_{\ell m};\mathcal{S}^{\ell
m}_{a},\mathcal{S}_{\ell m})\,,$ (67)

are the energy-momentum tensor harmonic components, which are gauge invariant
as we have argued in the previous section.

These master equations, using Schwarzschild coordinates, $(x^{a})=(t,r)$, and
transforming to the radial tortoise coordinate

$\frac{dx}{dr}=\frac{1}{f(r)}\,,$ (68)

can be written as a wave equation in one spatial dimension with a potential
and a source term:

$\left(-\partial^{2}_{t}+L_{V}\right)\Psi_{\ell m}=f\,\mathcal{F}_{\ell m}\,,$
(69)

where the operator $L_{V}$ is a Schrödinger-type operator

$L_{V}=\partial^{2}_{x}-V_{\ell}\,.$ (70)

In the rest of the paper, for the sake of simplicity, we omit the harmonic
numbers in the notation.

###  III.1 Master Functions and Equations in the Standard Branch

The standard branch is characterized by master equations determined by the
known potentials (Regge-Wheeler and Zerilli). The most general master function
is a linear combination of two master functions:

${}_{S}\Psi^{\mathrm{even}}_{\mathrm{odd}}=\mathcal{C}_{1}\Psi^{\mathrm{ZM}}_{\mathrm{CPM}}+\mathcal{C}_{2}\Psi^{\mathrm{NE}}_{\mathrm{RW}}\,,$
(71)

where $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ are two arbitrary constants. We
present them in a explicitly gauge-invariant and covariant (with respect to
coordinates changes in $M^{2}$) form. In the odd-parity case, the two master
functions can be taken to be the well-known Regge-Wheeler [41] master function

$\Psi_{\mathrm{RW}}=\frac{r^{a}}{r}\tilde{h}_{a}\,.$ (72)

and the Cunningham-Price-Moncrief master function [79, 80, 81] (see also [82,
83]):

$\Psi_{\mathrm{CPM}}=\frac{2r}{(\ell-1)(\ell+2)}\varepsilon^{ab}\left(\tilde{h}_{b:a}-\frac{2}{r}r_{a}\tilde{h}_{b}\right)\,.$
(73)

In the absence of matter sources ($T_{\mu\nu}=0\Rightarrow\mathcal{F}=0$), the
Regge-Wheeler master function is exactly the time derivative of the
Cunningham-Price-Moncrief master function (see e.g. [83]). With matter sources
it is modified in the following way

$t^{a}\Psi_{{\mathrm{CPM}}\,:a}=2\,\Psi_{\mathrm{RW}}-\frac{4r}{(\ell+2)(\ell-1)}r^{a}\mathcal{S}_{a}\,.$
(74)

In the standard branch of the odd-parity case the potential is the Regge-
Wheeler potential [41]:

$\displaystyle
V_{\mathrm{RW}}(r)=f\left\\{\Lambda+\frac{\ell(\ell+1)+3(f-1)}{r^{2}}\right\\}\,.$
(75)

The source term of the master equation (65)] depends on the choice of master
function within the family given in Eq. (71). For the Regge-Wheeler master
function it is:

$\displaystyle\mathcal{F}^{\mathrm{RW}}=\frac{1}{r}\left[r^{a}\left(\mathcal{S}_{:a}-2\mathcal{S}_{a}\right)+\frac{1-3f-\Lambda
r^{2}}{r}\mathcal{S}\right]\,,$ (76)

while for the Cunningham-Price-Moncrief master function it is:

$\displaystyle\mathcal{F}^{\mathrm{CPM}}=\frac{4r}{(\ell+2)(\ell-1)}\varepsilon^{ab}\mathcal{S}_{a:b}\,.$
(77)

In the even-parity case, the two master functions can be taken to be: (i) The
Zerilli-Moncrief master function $\Psi_{\mathrm{ZM}}$ [43, 44, 45], introduced
first by Zerilli [44] and later by Moncrief [45] (see also [84, 83]):

$\Psi_{\mathrm{ZM}}=\frac{2r}{\ell(\ell+1)}\left\\{\tilde{K}+\frac{2}{\lambda}\left(r^{a}r^{b}\tilde{h}_{ab}-rr^{a}\tilde{K}_{:a}\right)\right\\}\,,$
(78)

where

$\displaystyle\lambda(r)$ $\displaystyle=$ $\displaystyle
rf^{\prime}-2\left(f-1\right)+(\ell+2)(\ell-1)$ (79) $\displaystyle=$
$\displaystyle(\ell+2)(\ell-1)-\Lambda r^{2}-3\left(f-1\right)\,.$

(ii) The following less known master function (using the name given in [58]):

$\Psi_{\mathrm{NE}}(t,r)=\frac{1}{\lambda(r)}t^{a}\left(r\tilde{K}_{:a}-\tilde{h}_{ab}r^{b}\right)\,.$
(80)

This master function was previously unknown to us. But it turns out that a
non-covariant form in the Regge-Wheeler gauge was used in Ref. [85], where the
authors found it more convenient for second-order perturbative calculations in
the context of the close-limit approximation [86, 87, 88, 89]. Moreover, more
recently it has also appeared as the time derive of the Zerilli-Moncrief
master function [90]. In [58] we provided the explicitly gauge-invariant and
covariant form for perturbations of $\Lambda$-vacuum spherically-symmetric
spacetimes.

In the absence of matter sources, the master function $\Psi_{\mathrm{NE}}$ is
exactly the time derivative of the Zerilli-Moncrief master function (see [90]
and [58] for the covariant version). When considering matter sources, we find
that it is modified with the addition of harmonic components of the energy-
momentum tensor in the following way:

$t^{a}\Psi_{\mathrm{ZM}:a}=2\,\Psi_{\mathrm{NE}}+\frac{4r^{2}}{\ell(\ell+1)\lambda(r)}\mathcal{Q}_{ab}t^{a}r^{b}\,.$
(81)

where $\lambda(r)$ was introduced in Eq. (79). In the standard branch of the
even-parity case the potential is the Zerilli potential [43]

$V_{\mathrm{Z}}(r)=f(r)\frac{\lambda^{3}(r)-2\Lambda
r^{2}\left[\lambda(r)-(\ell+2)(\ell-1)\right]^{2}+2\left(\ell+2\right)^{2}\left(\ell-1\right)^{2}\left(\ell^{2}+\ell+1\right)}{3r^{2}\lambda^{2}(r)}\,.$
(82)

As in the odd-parity sector, the source term of the master equation (65)
depends on the choice of master function. For the Zerilli-Moncrief master
function we find the following source term:

$\displaystyle\mathcal{F}^{\mathrm{ZM}}$ $\displaystyle=$
$\displaystyle-2r\mathcal{P}+\frac{4}{\lambda}r^{a}\mathcal{P}_{a}+\frac{16r}{\ell(\ell+1)\lambda}\left[1-\frac{(\ell+2)(\ell-1)}{2\lambda}\right]Q_{ab}t^{a}t^{b}+\frac{2r}{\ell(\ell+1)}\left[1+\frac{2(r^{2}\Lambda-1+2f)}{\lambda}\right]\mathcal{Q}^{a}_{a}$
(83) $\displaystyle+$
$\displaystyle\frac{4r^{2}}{\ell(\ell+1)\lambda}t^{c}\varepsilon^{ab}\mathcal{Q}_{cb:a}\,,$

which is equivalent to the covariant expression in [83] after considering the
energy-momentum conservation equations (see Appendix D). For the other even-
parity master function, $\Psi_{\mathrm{NE}}$ we find the following expression
for the source term

$\displaystyle\mathcal{F}^{\mathrm{NE}}$ $\displaystyle=$
$\displaystyle\left[\frac{3\ell(\ell+1)}{r\lambda}-\frac{2}{r}+\frac{8f}{r\lambda}\left(1\frac{(\ell+2)(\ell-1)}{2\lambda}-1\right)\right]t^{a}\mathcal{P}_{a}+\frac{6}{\lambda^{2}}\left(f-1+\frac{\Lambda
r^{2}}{3}\right)\mathcal{Q}_{ab}t^{a}r^{b}$ (84) $\displaystyle-$
$\displaystyle
rt^{a}\mathcal{P}_{:a}+\frac{2f}{\lambda}\varepsilon^{ab}\mathcal{P}_{a:b}-\frac{r}{\lambda}r^{c}\varepsilon^{ab}\mathcal{Q}_{cb:a}\,.$

This completes the description of master equations in the standard branch.

###  III.2 Master Functions and Equations in the Darboux Branch

The Darboux branch contains a family of potentials that are different from the
ones of the standard branch, i.e. the Regge-Wheeler and Zerilli potentials.
This family of potentials is characterized by the fact that they have to
satisfy the following non-linear ODE:

$\left(\frac{\delta{V}_{,x}}{\delta{V}}\right)_{,x}+2\left(\frac{V^{\mathrm{RW/Z}}_{,x}}{\delta{V}}\right)_{,x}-\delta{V}=0\,,$
(85)

which is valid for both parities: For the even-parity case one has to select
the Zerilli potential, $V_{\mathrm{Z}}$, in the second term, while for the
odd-parity case we have to select the Regge-Wheeler potential,
$V_{\mathrm{RW}}$. Moreover, we have introduced the following definition:

$\displaystyle\delta V$ $\displaystyle=$ $\displaystyle
V^{\mathrm{odd/even}}-V^{\mathrm{RW/Z}}\,,$ (86)

The master functions in the Darboux branch are also different from the ones in
the standard branch and they depend on the particular potential, solution of
Eq. (85), that we consider. In order to set the ground for extending the
analysis of [59] to the presence of perturbative sources (see Sec. III.3 for
this), it is convenient to revisit the results of Ref. [58] for the Darboux
branch, clarifying an issue with the derivation made there. The Darboux master
functions presented in [58] were obtained from a systematic procedure and they
appear as a linear combination of master functions with two arbitrary
constants. However, one of these constants (the one that multiples a master
function from the standard branch) should have been set to zero by one of the
conditions of the systematic analysis that was not used [See, in Ref. [58],
equation (86) for the odd-parity case and equation (145) for the even-parity
one]. Then, in the Darboux branch, the most general master function can be
written as:

${}_{D}\Psi^{\mathrm{even}}_{\mathrm{odd}}=\mathcal{C}\left(\Sigma^{\mathrm{even}}_{\mathrm{odd}}\,\Psi^{\mathrm{ZM}}_{\mathrm{CPM}}+\Phi^{\mathrm{even}}_{\mathrm{odd}}\right)\,,$
(87)

with $\mathcal{C}$ being the only arbitrary constant. The functions
$\Phi^{\mathrm{even}}_{\mathrm{odd}}$ are gauge-invariant linear combinations
of the metric perturbations and their first-order derivatives, but only the
combination with $\Psi^{\mathrm{ZM}}_{\mathrm{CPM}}$ in Eq. (87) is a true
master function. In the odd-parity case, $\Phi^{\mathrm{odd}}$ is given by:

$\displaystyle\Phi_{\mathrm{odd}}(x^{a})=-\varepsilon^{ab}\tilde{h}_{a:b}\,.$
(88)

while in the even-parity case $\Phi_{\mathrm{even}}$ is:

$\displaystyle\Phi_{\mathrm{even}}(x^{a})=\tilde{K}\,.$ (89)

On the other hand, in Eq. (87)
$\Sigma^{\mathrm{even}}_{\mathrm{odd}}=\Sigma^{\mathrm{even}}_{\mathrm{odd}}(x)$
is a function that contains the integral of the potential. In the odd-parity
case it is given by:

$\displaystyle\Sigma_{\mathrm{odd}}(x)=\frac{1}{2}\left[\frac{1-f(r)+\Lambda\,r^{2}}{2\,r}+\int^{x}_{x_{0}}dx^{\prime}V^{\mathrm{odd}}(x^{\prime})\right]\,,$
(90)

and in the even-parity case it is

$\displaystyle\Sigma_{\mathrm{even}}(x)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[\frac{(\ell+2)(\ell-1)-\lambda(r)}{2\,r}-\mathcal{C}_{\rm
E}\right.$ (91) $\displaystyle+$
$\displaystyle\left.\int^{x}_{x_{0}}dx^{\prime}V^{\mathrm{even}}(x^{\prime})\right]\,,$

where $\mathcal{C}_{\rm E}$ is a constant proportional to the cosmological
constant and the mass:

$\displaystyle\mathcal{C}_{\rm E}=\frac{4M\Lambda}{(\ell+2)(\ell-1)}\,.$ (92)

The lower limit of integration appearing in these equations, $x_{0}$, depends
on the particular choice of background spacetime (see Appendix A, where the
notation is introduced). For the BH background cases (Sch, SchdS, SchAdS),
where $M\neq 0$, it is $x_{0}\rightarrow-\infty$. And for the maximally
symmetric cases (M, dS, AdS), where $M=0$, it is $x_{0}=0$. The different
values is a consequence of the different properties of the tortoise coordinate
for the different background cases.

Finally, the general master equation in the Darboux branch is:

$\displaystyle\left(-\partial^{2}_{t}+\partial^{2}_{x}-V^{\mathrm{even/odd}}_{\ell}\right)_{D}\Psi^{\mathrm{even}}_{\mathrm{odd}}=f\,\mathcal{F}^{\mathrm{even}}_{\mathrm{odd}}\,,$
(93)

and the source terms $\mathcal{F}^{\mathrm{even}}_{\mathrm{odd}}$, associated
to the master functions in the Darboux branch, are discussed in the next
section, where we write them in a way that avoids lengthy and cumbersome
expressions.

###  III.3 Darboux Covariance of the Metric Perturbations

We are now going to study the hidden symmetry that joins the elements of the
space of master functions and equations, and which is realized in the form of
Darboux transformations [91, 92, 93]. In particular, we study whether Darboux
covariance is still present when we have matter sources and what it requires
to accommodate it, extending in this way the study of Ref. [59]. To that end,
let us first review what exactly we mean by a Darboux transformation in the
context of master equations. Then, let us consider two master wave equations
of the type given in Eq. (69), that is

$\left(-\partial^{2}_{t}+L_{v}\right)\Phi=\sigma\,,\quad\left(-\partial^{2}_{t}+L_{V}\right)\Psi=S\,.$
(94)

The novelty here is that we are considering non-homogeneous equations with
sources $\sigma$ and $S$. The study of Darboux transformation for the
homogeneous case $S=\sigma=0$, in the context of the space of master functions
and equations, was carried out in [59] where it was introduced as a mapping
between the two pairs of master functions and potentials [see Eq. (94) above]:

$\displaystyle(v,\Phi)\rightarrow(V,\Psi):\quad\left\\{\begin{array}[]{lcl}\Psi&=&\Phi_{,x}+g\,\Phi\,,\\\\[5.69054pt]
V&=&v+2\,g_{,x}\,,\end{array}\right.$ (97)

where $g$ is the Darboux transformation generating function. In the case where
there are no sources, $S=\sigma=0$, the mapping between master equations
requires that the following consistency condition on $g$ is fulfilled:

$\displaystyle g_{x}-g^{2}+v=\mathcal{C}_{R}\,,$ (98)

where $\mathcal{C}_{R}$ is an arbitrary integration constant, as the original
consistency condition is the derivative of this equation. This is a Riccati
equation that restricts the set of possible Darboux transformation generating
functions $g$. When the source terms in Eq. (94) are present, we find that an
additional condition, on the sources $\sigma$ and $S$, is required

$\displaystyle S=\sigma_{,x}+g\sigma\,,$ (99)

This new relation tells us how the source terms have to transform in order for
the two master equations in Eq. (94) to be Darboux related. This condition can
taken as part of the definition of the Darboux transformation when source
terms are present. From these definitions and equations it is clear that the
inverse Darboux transformation is generated by $-g$ (see Ref. [93] for details
on the invertibility of Darboux transformations).

The link to our infinite collection of master functions and equations is based
on a mapping between the Darboux generating function and the potentials. This
comes from the Darboux transformation of the potential [Eq. (97)] and from the
derivative of the Riccati equation [Eq. (98)] for the Darboux transformation
generating function. Indeed, these two equations can be manipulated to obtain
the following expressions for $g_{,x}$ and $g$ in terms of the potentials:

$g_{,x}=\frac{V-v}{2}\,,\quad g=\frac{(V+v)_{,x}}{2(V-v)}\,.$ (100)

It is straighforward to see that the consistency between these two equations
for $g$ and $g_{,x}$, that is, asking that the derivative of the second
equation is equal to the first one, yields

$\left(\frac{\delta V_{,x}}{\delta V}\right)_{,x}+2\left(\frac{v_{,x}}{\delta
V}\right)_{,x}-\delta V=0\,,$ (101)

where $\delta V=V-v\,$. But, if we make the correspondence
$v\longleftrightarrow V^{\mathrm{RW/Z}}$, this is precisely Eq. (85). In other
words, Eq. (101) coincides with the condition that any potential belonging to
the Darboux branch should satisfy [58]. In the absence of sources, this is
enough to establish that all master equations in the Darboux branch are
connected with the standard branch via Darboux transformations. The
corresponding Darboux generating function can be obtained by integrating the
first expression in Eq. (100), i.e.

$\displaystyle
g^{\mathrm{Z}\rightarrow\mathrm{even}}_{\mathrm{RW}\rightarrow\mathrm{odd}}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\int_{x}dx^{\prime}\,\left(V^{\mathrm{even}}_{\mathrm{odd}}-V^{\mathrm{Z}}_{\mathrm{RW}}\right)$
(102) $\displaystyle=$
$\displaystyle\frac{1}{2}\int^{x}_{x_{0}}dx^{\prime}\,\left(V^{\mathrm{even}}_{\mathrm{odd}}-V^{\mathrm{Z}}_{\mathrm{RW}}\right)+\mathcal{C}_{\mathrm{S}\rightarrow\mathrm{D}}\,.$

The choice of finding $g$ in integral form, instead of the derivative one
given by the second expression in Eq. (100), is motivated by the fact that the
master functions in the Darboux branch (87) depend on integral functions of
the potential. At first glance, it seems like there is an arbitrary
integration constant, $\mathcal{C}_{\mathrm{S}\rightarrow\mathrm{D}}$, in this
definition. However, this constant is fixed by requiring that the expression
for $g$ in the second relation in Eq. (100) has to be equal to the integral of
the expression for $g_{,x}$ in the first one. That is, we can write

$\displaystyle\frac{(V+v)_{,x}}{(V-v)}=\int_{x}dx^{\prime}\left(V-v\right)\,.$
(103)

Something similar happens between the odd- and even-parity sectors in the
standard branch. Indeed, we can introduce a Darboux generating function that
transforms the Regge-Wheeler potential into the Zerilli in the same way, i.e.

$\displaystyle g^{\mathrm{RW}\rightarrow\mathrm{Z}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int^{x}_{x_{0}}dx^{\prime}\,\left(V^{\mathrm{Z}}-V^{\mathrm{RW}}\right)+\mathcal{C}_{\mathrm{RW}\rightarrow\mathrm{Z}}\,.$
(104)

The constant $\mathcal{C}_{\mathrm{RW}\rightarrow\mathrm{Z}}$ is fixed by Eq.
(103) and reads:

$\displaystyle\mathcal{C}_{\mathrm{RW}\rightarrow\mathrm{Z}}=\alpha-\frac{\mathcal{C}_{\rm
E}}{2}\,,$ (105)

where $\alpha$ is given by:

$\displaystyle\alpha=\frac{1}{12M}\frac{(\ell+2)!}{(\ell-2)!}\,,$ (106)

and $\mathcal{C}_{\rm E}$ is given by Eq. (92). We will go back to this
generating function between odd- and even-parity standard branches in Sec.
III.4, where we analyze the consequences of Darboux transformation when
considered in the frequency domain (the typical arena for this kind of
transformations [93]).

Before getting deeper into the general details of the Darboux transformation
between master functions, let us discuss a particular feature of it that only
applies in the case of the maximmally-symmetric background spacetimes
(Minkowski, de Sitter and anti-de Sitter). By looking back at Eqs. (104),
(105) and (106), it is clear that the Darboux generating function is not
defined in the limit $M\to 0$, corresponding to the maximmally-symmetric
backgrounds. This is reflected in the fact that, in this case, the even- and
odd-parity potentials (75) and (82) in the standard branch, are identical and
reduce to:

$V^{\rm odd}_{\ell}=V^{\rm even}_{\ell}=f\frac{\ell(\ell+1)}{r^{2}}\,.$ (107)

If we now look at the equation that potentials in the Darboux branch have to
satisfy, Eq. (85), it is a second-order non-linear ODE which is not trivial to
solve analytically. However, it turns out that in the maximally-symmetric case
it admits some simple solutions. Indeed, if we start from the potential of the
standard branch, Eq. (107), we obtain the following two solutions:

$\displaystyle V_{\ell+}=f\frac{(\ell+1)(\ell+2)}{r^{2}}\,,\quad
V_{\ell-}=f\frac{\ell(\ell-1)}{r^{2}}\,.$ (108)

The Darboux generating functions that produces the transformations to these
potentials are:

$g_{+}=-\frac{\ell+1}{r}\,,\qquad g_{-}=\frac{\ell}{r}\,.$ (109)

One can check that they satisfy the Riccati equation (98) and the consistency
condition (103) . As we can see, the first one, $g_{+}$, correspond to
shifting the harmonic number $\ell$ to $\ell+1$, while the second one,
corresponds to the shift: $\ell\rightarrow\ell-1$. This means that, in the
case of maximally-symmetric background spacetimes, the Darboux transformation
acts as a sort of ladder operator, increasing or decreasing by one the value
of the angular momentum number $\ell$. This is a very interesting result,
since it means that for these backgrounds, Darboux transformation not only
connects solutions of the master equations with the same harmonic numbers
$(\ell,m)$ but also solutions corresponding to different number $\ell$. In the
case where the master equations do not have sources, the solutions are in
principle independent of $m$, so we can connect all of them.

Another point of view consists in realizing that the Darboux transformations
generated by the functions in Eq. (109) do not change the shape of the
potential. To be more specific, a Darboux transformation is said to change the
potential in a _shape invariant_ way when it induces in the potential [see Eq.
(97)] a transformation that can be written in the form (see, e.g. [94, 95] for
a detailed exposition):

$V(x;f(\vec{\lambda}))=v(x;\vec{\lambda})+R(\vec{\lambda})\,,$ (110)

where the vector $\vec{\lambda}$ contains parameters on which the potential
depends. In our case, it is just $\ell$ so that $f(\vec{\lambda})\equiv
f(\ell)=\ell+1$ and $R(\vec{\lambda})\equiv 0$ in the case in which we raise
$\ell$ by one. A very interesting consequence of this property is that the
spectrum can be reconstructed in an analytical way. However, in our case, this
property only happens for the case of maximmally-symmetric backgrounds.

We now get back to the general case, but still without sources
($T_{\mu\nu}=0$), and we are going to show that the master functions of the
Darboux branch can be written in a form that shows their explicit connection
with Darboux transformations. In doing so, we are going to point to some
subtle aspects that were not addressed in [59] and which are very relevant for
extending all this formalism to the case with sources ($T_{\mu\nu}\neq 0$).
First, as it was already noticed by Chandrasekhar [62], the (Darboux)
transformation between odd- and even-parity master functions is a formal
transformation in the sense that it connects the master equations or,
equivalently, the potentials and the solutions of the master equations. But it
is not a transformation that can be formulated at level of the metric
perturbations, which is obvious because there is no possibility of
transforming odd-parity metric perturbations into even-parity ones. Following
this line of thought, we can see that things are different in the Darboux
branch of master functions and equations, as we are going to see that the
master functions come from explicit Darboux transformations of the standard
branch master functions. To see this, let us first consider the even-parity
sector. The first step is to realize that the gauge-invariant function
$\Phi_{\mathrm{even}}$ can be written as the following transformation of the
Zerilli-Moncrief master function:

$\Phi_{\mathrm{even}}=\Psi_{\mathrm{ZM},x}+{g}_{\ast}^{\mathrm{even}}\Psi_{\mathrm{ZM}}\,,$
(111)

where

${g}_{\ast}^{\mathrm{even}}=\frac{3M}{2r^{2}}+\frac{\mathcal{C}_{\rm
E}}{2}-\frac{1}{2}\int^{x}_{x_{0}}dx^{\prime}\,V^{\mathrm{Z}}\,.$ (112)

This transformation looks like a Darboux transformation but it is actually not
because the function ${g}_{\ast}^{\mathrm{even}}$ does not satisfy the
consistency condition (98) for a Darboux transformation generating function.
Nevertheless, if we introduce this transformation into the even-parity version
of Eq. (87) we find (for simplicity we set $\mathcal{C}=1$):

${}_{D}\Psi_{\mathrm{even}}=\Psi_{\mathrm{ZM},x}+g^{\mathrm{even}}\Psi_{\mathrm{ZM}}\,,$
(113)

where $g^{\mathrm{even}}$ is given by the expression in Eq. (102), but with
$\mathcal{C}_{\mathrm{S}\rightarrow\mathrm{D}}=0$, i.e.

$\displaystyle
g^{\mathrm{even}}=\frac{1}{2}\int^{x}_{x_{0}}dx^{\prime}\,\left(V^{\mathrm{even}}-V^{\mathrm{Z}}\right)\,,$
(114)

which satisfies the consistency condition by construction. This makes Eq.
(113) a true Darboux transformation, which holds at the level of the metric
perturbations by virtue of the perturbative Einstein equations.

The situation in the odd-parity sector is quite similar. The function
$\Phi_{\mathrm{odd}}$ can be written as the following transformation of the
Cunningham-Price-Moncrief master function:

$\Phi_{\mathrm{odd}}=\Psi_{\mathrm{CPM},x}+{g}_{\ast}^{\mathrm{odd}}\Psi_{\mathrm{CPM}}\,,$
(115)

where

${g}_{\ast}^{\mathrm{odd}}=-\frac{M}{2r^{2}}-\frac{r\Lambda}{3}-\frac{1}{2}\int^{x}_{x_{0}}dx^{\prime}\,V^{\mathrm{RW}}\,.$
(116)

Again, Eq. (115) is not a Darboux transformation for the same reason as
before, but we can be inserted it into Eq. (87) to find (again, for
simplicity, we set $\mathcal{C}=1$):

${}_{D}\Psi_{\mathrm{odd}}=\Psi_{\mathrm{CPM},x}+g^{\mathrm{odd}}\Psi_{\mathrm{CPM}}\,,$
(117)

where $g^{\mathrm{odd}}$ is given by Eq. (102), again with
$\mathcal{C}_{\mathrm{S}\rightarrow\mathrm{D}}=0$:

$\displaystyle
g^{\mathrm{odd}}=\frac{1}{2}\int^{x}_{x_{0}}dx^{\prime}\,\left(V^{\mathrm{odd}}-V^{\mathrm{RW}}\right)\,.$
(118)

It is important to mention that the systematic procedure of Ref. [58] contains
the necessary consistency conditions to fix
$\mathcal{C}_{\mathrm{S}\rightarrow\mathrm{D}}=0$ for both the odd- and even-
parity sectors. This seems to happen because the consistency condition in the
systematic procedure is imposed at the level of the allowed potentials, so
that the generating function defined from them automatically satisfies the
consistency condition. Therefore, the master function in the Darboux branches
are actual Darboux transformations of the standard branch master functions.

Darboux covariance provides another key piece of information about the
infinite landscape of master equations: It helps understanding why the Darboux
branch contains only master functions and equations that result from Darboux
transformations of the Cunningham-Price-Moncrief and Zerilli-Moncrief master
functions, $\Psi_{\mathrm{CPM}}$ and $\Psi_{\mathrm{ZM}}$, and their
associated potentials (which completely determine the master equations in the
absence of matter sources). Indeed, the Darboux transformations of the other
master functions in the standard branch, the Regge-Wheeler master function in
the odd-parity case, $\Psi_{\mathrm{RW}}$, and the other even-parity master
function [Eq. (80)] in the even-parity case, $\Psi_{\mathrm{NE}}$, are not
present at all in the Darboux branch. Actually, one could in principle
consider master functions that originate from $\Psi_{\mathrm{RW}}$ and
$\Psi_{\mathrm{NE}}$, that is,

$\displaystyle\Psi^{\mathrm{NE}}_{\mathrm{RW},x}+g^{\mathrm{even}}_{\mathrm{odd}}\Psi^{\mathrm{NE}}_{\mathrm{RW}}\,.$
(119)

In principle, these master functions, with their associated potentials, can
also constitute a valid description of the perturbations. However, one can
show that the master functions in Eq. (119) contain second-order derivatives
of the metric perturbations which cannot be eliminated by using Einstein
perturbed equations [see Appendix C]. But this is against the restrictions we
imposed on the space of master functions in [58], where we only considered
master functions that are linear combinations of the metric perturbations and
their first-order derivatives. This shows once more that the Cunningham-Price-
Moncrief and Zerilli-Moncrief master functions seem to play a special role in
BHPT.

Let us now focus on the extension to the case with sources generated by the
existence of a first-order energy-momentum tensor $T_{\mu\nu}$. Here the
situation is different and Eqs. (113) and (117) needs corrections that depend
on the harmonic components of the energy-momentum tensor. The relation between
the master functions in the standard and Darboux branches in the even-parity
sector is now

${}_{D}\Psi_{\mathrm{even}}=\Psi_{\mathrm{ZM},x}+g^{\mathrm{even}}\Psi_{\mathrm{ZM}}-\frac{4r^{2}}{\ell(\ell+1)\lambda}t^{a}t^{b}Q_{ab}\,,\qquad$
(120)

and in the case of the odd-parity sector:

${}_{D}\Psi_{\mathrm{odd}}=\Psi_{\mathrm{CPM},x}+g^{\mathrm{odd}}\Psi_{\mathrm{CPM}}+\frac{4r}{(\ell+2)(\ell-1)}t^{a}S_{a}\,.\qquad$
(121)

The source terms appearing in these equations come from the modification of
the transformations in Eqs. (111) and (115) between the functions
$\Phi^{\mathrm{even}}_{\mathrm{odd}}$ and $\Psi^{\mathrm{ZM}}_{\mathrm{CPM}}$
in the presence of matter. Clearly, these equations do not correspond to
Darboux transformations. Even more, one can check that the sources
$\mathcal{F}^{\mathrm{even}}_{\mathrm{odd}}$ are not Darboux transformations
of the standard branch sources $\mathcal{F}^{\mathrm{ZM}}_{\mathrm{CPM}}$
either [see Eq. (99) for the transformation of the sources of the master
equations under a Darboux transformation]. At first sight, one can then think
that Darboux covariance may be broken, but it turns out that we can rearrange
the variables in such a way that it can be preserved. Indeed, let us define,
for the even-parity sector the following function

${}_{D}\check{\Psi}_{\mathrm{even}}$ $\displaystyle=$
${}_{D}\Psi_{\mathrm{even}}+\frac{4r^{2}}{\ell(\ell+1)\lambda}t^{a}t^{b}Q_{ab}$
(122) $\displaystyle=$
$\displaystyle\Psi_{\mathrm{ZM},x}+g^{\mathrm{even}}\Psi_{\mathrm{ZM}}\,,$

where $g^{\mathrm{even}}$ is given in Eq. (114). And for the odd-parity case,
let us define

${}_{D}\check{\Psi}_{\mathrm{odd}}$ $\displaystyle=$
${}_{D}\Psi_{\mathrm{odd}}-\frac{4r}{(\ell+2)(\ell-1)}t^{a}S_{a}$ (123)
$\displaystyle=$
$\displaystyle\Psi_{\mathrm{CPM},x}+g^{\mathrm{odd}}\Psi_{\mathrm{CPM}}\,,$

where $g^{\mathrm{odd}}$ is given in Eq. (118). The master functions,
${}_{D}\check{\Psi}_{\mathrm{even}}$ and ${}_{D}\check{\Psi}_{\mathrm{odd}}$,
have the form of Darboux transformations of the standard branch master
functions and, more importantly, they satisfy master equations with same
potential as the original master functions, ${}_{D}\Psi_{\mathrm{even}}$ and
${}_{D}\Psi_{\mathrm{odd}}$, but with a different source term. That is,

$\displaystyle\left(-\partial^{2}_{t}+\partial^{2}_{x}-V^{\mathrm{even/odd}}_{\ell}\right)_{D}\check{\Psi}^{\mathrm{even}}_{\mathrm{odd}}=f\,\check{\mathcal{F}}^{\mathrm{even}}_{\mathrm{odd}}\,,\qquad$
(124)

where the expressions for the new source terms,
$\check{\mathcal{F}}^{\mathrm{even}}_{\mathrm{odd}}$, are given by

$\displaystyle\check{\mathcal{F}}^{\mathrm{even}}=\mathcal{F}^{\mathrm{even}}-\left(\square_{2}-\frac{1}{f}V^{\mathrm{even}}\right)\frac{4r^{2}}{\ell(\ell+1)\lambda}t^{a}t^{b}Q_{ab}\,,\qquad$
(125)

and

$\displaystyle\check{\mathcal{F}}^{\mathrm{odd}}=\mathcal{F}^{\mathrm{odd}}+\left(\square_{2}-\frac{1}{f}V^{\mathrm{odd}}\right)\frac{4r}{(\ell+2)(\ell-1)}t^{a}S_{a}\,.\qquad$
(126)

These new source terms are now related to the standard branch ones by the
Darboux condition (99), i.e.

$\displaystyle\check{\mathcal{F}}^{\mathrm{even}}$ $\displaystyle=$
$\displaystyle\mathcal{F}^{\mathrm{ZM}}_{,x}+g^{\mathrm{even}}\mathcal{F}^{\mathrm{ZM}}\,,$
(127) $\displaystyle\check{\mathcal{F}}^{\mathrm{odd}}$ $\displaystyle=$
$\displaystyle\mathcal{F}^{\mathrm{CPM}}_{,x}+g^{\mathrm{odd}}\mathcal{F}^{\mathrm{CPM}}\,.$
(128)

Then, in the presence of matter sources in the master equations, the pairs
$(\check{\Psi},\check{\mathcal{F}})$ are the true Darboux transformations of
the standard branch master functions/equations. Therefore, thanks to an
adequate rearrangement of the master variables, Darboux covariance is
preserved in the presence of matter sources.

###  III.4 Darboux Transformations and Isospectrality

One of the most important consequences of Darboux covariance in vacuum is
isospectrality. This can be seen very well in the frequency domain, for
instance by just considering a single frequency mode of the master function

$\Psi(t,r)=e^{ikt}\,\psi(k;x)\,,$ (129)

so that the master equation (69) becomes a time-independent Schrödinger
equation333In this subsection we assume there are no matter sources driving
the perturbations, since we are precisely interested in the (quasi)normal free
oscilations.:

$L_{V}\psi=-k^{2}\psi\,,$ (130)

where $L_{V}$ is given in Eq. (70) and $k$ is the frequency of the mode. It is
simple to check that a Darboux transformation [see Eq. (97)] not only
transforms this time-independent Schrödinger equation into another one, but it
does so in such a way that the eigenvalue $-k^{2}$ remains unchanged. In the
asymptotically-flat case this leads to the isospectrality of quasinormal modes
for all the master equations in our space. It also implies that the
transmission and reflection coefficients for scattering states are the same.
However, the isospectrality between odd- and even-parity quasinormal modes is
not guaranteed when $\Lambda\neq 0$ (see, e.g. [96, 97]). Moreover, it may
well be that isospectrality is a distinctive characteristic of General
Relativity since it is broken in other theories, opening the door to perform
fundamental physics tests with advanced gravitational-wave detectors [7, 10,
11]. The loss of isospectrality in other theories of gravity has been shown,
for instance, in [98, 99, 100, 101, 102] for non-rotating BHs (in particular
[103] in the context of loop quantum gravity) and [104] for the case of
spinning BHs, both in alternative theories of gravity. On the other hand, the
loss of isospectrality has been recently discussed in the context of the study
of the instability of quasinormal modes [105, 106] within General Relativity.

In this context, a remarkable fact, already pointed out in [107, 61] (see also
[59]) is that the Regge-Wheeler and Zerilli potentials can be derived from a
superpotential in the following way:

$V^{\mathrm{RW}}_{\mathrm{Z}}=\mp g_{,x}+g^{2}+k_{0}^{2}\,,$ (131)

where $g$ is the superpotential, which in the case of Schwarzschild looks as

$\displaystyle g(r)=ik_{0}+\frac{6\,Mf(r)}{\lambda(r)r^{2}}=ik_{0}+W(r)\,,$
(132)

and $k_{0}$ is the so-called algebraically special frequency, given by

$ik_{0}=\alpha=\frac{(\ell+2)(\ell+1)\ell(\ell-1)}{12\,M}\,,$ (133)

where $\alpha$ was introduced before in Eq. (106). In Eq. (132), $W(r)$ is a
potential used by Heading and Chandrasekhar to express the Regge-Wheeler and
Zerilli potentials in another compact form:

$V^{\mathrm{RW}}_{\mathrm{Z}}=\mp W_{,x}+\alpha W+W^{2}\,.$ (134)

The two potentials, $W$ and $g$, differ by just a constant, but $g$ plays a
special role as it is the generator of the Darboux transformation between of
the two parities of the standard branch. It also corresponds to what, in the
context of supersymmetric (SUSY) quantum mechanics, is called a
superpotential. The analogy with SUSY is actually exact if we absorb the
$k^{2}_{0}$ term in Eq. (131) into the eigenvalue in Eq. (130). For more
details on the role of Darboux transformations in SUSY quantum mechanics see
[108, 109, 110, 111]. It is important to notice that $g$ automatically
satisfies the Riccati consistency condition (98) with
$\mathcal{C}_{R}=k_{0}^{2}$, which is something that highlights the difference
between the Darboux tranformation as originally defined [92, 91] (see also
[93]) and Darboux covariance: Darboux transformations were originally
introduced in the frequency domain and they are generated from a particular
solution of the time-independent Schrödinger equation (130), typically an
eigenfunction with eigenvalue $k_{0}$, to connect two equations with same
eigenvalue $k$; instead, Darboux covariance is a symmetry of the full master
equation that can be used both in the time domain and in the frequency domain
and does not rely on having a special solution (like a bound state). This why
we need the consistency condition (98), which in the original Darboux
transformation is satisfied by construction.

The eigenfunction associated with the eigenvalue $k_{0}$ and $g$ in Eq. (132),
is usually called a algebraically special mode but, sometimes, in the context
of scattering theory, it is usually known as an antibound state. In
Schwarzschild spacetime it can be written as

$\psi_{0}=\frac{\lambda(r)}{2}e^{-ik_{0}x}\,,$ (135)

so that the function Darboux generating function, $g(r)$, in Eq. (132), is
obtained as the logarithmic derivative of $\psi_{0}$, i.e.

$g(r)=-(\ln{\psi_{0}})_{,x}\,.$ (136)

However, is quite straightforward to check that it is also a special solution
even with different asymptotics, i.e. in Schwarzschild-de
Sitter/Schwarzschild-Anti de Sitter spacetimes (see also [96, 97] for a nice
discussion of the extension of this formalism to $\Lambda$-vacuum spacetimes
and the special role of this solution in the _pole-skipping_ method to
reconstruct the quasi-normal mode spectrum).

##  IV Reconstruction of the Metric Perturbations: General Case

With all the elements of BHPT presented in the previous section, we are now
going to tackle the problem of reconstructing the metric perturbations in
terms of the master functions and perturbative gauge functions. Having the
full knowledge of the metric perturbations is of interest for some tasks, like
the computation of the gravitational self-force or in perturbative Hamiltonian
formulations.

Here, we show how to carry out the metric reconstruction for any background
within the family of the spherically-symmetric $\Lambda$-vacuum spacetimes,
not just for the case of the Schwarzschild BH. We also consider a general
situation in which the perturbations are driven by matter sources that are
generated by a general energy-momentum distribution $T_{\mu\nu}$ that enters
at the first perturbative order (it does not affect the background spacetime).
As we have done with the previous developments, we are going to carry out the
derivations in an arbitrary gauge and the expressions we obtain are valid in
an arbitrary coordinate system of the two-dimensional Lorentzian spacetime
$M^{2}$. That is, we are going to show that the metric reconstruction can be
in principle carried out in an arbitrary gauge. In order words, the metric
perturbations are completely determined once the master functions are found
and a perturbative gauge is fixed. This is of great relevance and gives a new
perspective into BHPT computations. It also reflects the fact that certain
combinations of the metric perturbations (and their derivatives) are gauge-
invariant. But this does not mean that the gauge-invariant quantities are the
true degrees of freedom of the gravitational field in the perturbative regime.
The true degrees of freedom are represented by the master functions which, as
is well known, determine all the physically relevant quantities about the
dynamics of the perturbations: Gravitational waveforms; Gravitational-waves
fluxes of energy; Gravitational-wave fluxes of angular momentum;
Gravitational-wave fluxes of linear momentum; etc. Another way of looking at
this fact, is to realize that while the construction of gauge-invariant
quantities from the metric perturbations is a process that can be carried out
independently of the Einstein equations, the reduction to the true dynamical
degrees of freedom (the master functions) requires the use of the
(perturbative) Einstein field equations. One illustrative example is the
reduction to the two gravitational wave polarizations in the case of the
linear theory (perturbations around Minkowski spacetime; see e.g. [14]).

Another important feature of our procedure for the metric reconstruction is
that it is systematic, thus we expect that it can be applied to other
scenarios: different backgrounds within spherical symmetry, different theories
of gravity, etc. After the derivation of the main results we discuss their
implications and scope in the light of Darboux covariance.

To start with, let us count metric perturbations. At the level of the
spacetime metric perturbations, we have ten independent metric perturbations
and the generator of the perturbative gauge transformations has four
components, therefore we can build six independent gauge-invariant metric
perturbations. At the level of the decomposition of the metric perturbations
into spherical harmonics, we have to distinguish the two different parity
sectors. In the odd-parity sector we have three independent metric
perturbations and just one gauge generator component and hence, there are 2
gauge invariant metric perturbations. In the even-parity sector we have seven
independent metric perturbations, while there are three components of the
gauge generator vector field, thus we can construct four gauge invariant
metric perturbations. It is important to clarify, that here gauge-invariant
metric perturbations means combinations of the metric perturbations and their
derivatives that are invariant under a general gauge transformation.

In the case of perturbations driven by matter sources, the components of the
energy-momentum tensor will also appear in the metric reconstruction. In the
vacuum case, the reconstruction can be carried out in two steps. In the first
step, we identify the pure gauge metric components, those that can be canceled
by a gauge transformation, and then we find how the rest of the components of
the metric perturbations can be written in terms of the gauge and the gauge-
invariant metric perturbations [Eqs. (57) and (58) for the even-parity case
and Eq. (60) for the odd-parity case]. The second step consists in writing the
gauge-invariant metric perturbations in terms of the master functions of the
standard branch.

###  IV.1 Metric Reconstruction in the Odd-Parity Sector

We start with the odd-parity sector as it involves less independent metric
perturbations, just three. Then, following the plan outlined above, let us
look at how the odd-parity metric perturbations transform [see Eqs. (55) and
(56)]. From the transformation in Eq. (56) we see that $h_{2}$ is a pure gauge
component of the metric perturbations, i.e. we can make it to have any value,
in particular we can be put it to zero by a gauge transformation, which is
what is done in the Regge-Wheeler gauge [41]. Therefore, we can write the
metric perturbations in the odd-parity case as follows:

$\displaystyle h_{a}$ $\displaystyle=$
$\displaystyle\tilde{h}_{a}+\frac{r^{2}}{2}\mathcal{G}^{\mathrm{odd}}_{:a}\,,$
(137) $\displaystyle h_{2}$ $\displaystyle=$ $\displaystyle
r^{2}\mathcal{G}^{\mathrm{odd}}\,,$ (138)

where $\mathcal{G}^{\mathrm{odd}}$ is the only odd-parity gauge function. This
completes the first step. The second step is to write the gauge-invariant
metric perturbations, $\tilde{h}_{a}$, in terms of the two independent master
functions of the standard branch: $\Psi_{\mathrm{RW}}$ and
$\Psi_{\mathrm{CPM}}$. Taking into account that the first is related to the
time derivative of the second [see Eq. (74)], we can always eliminate one of
them if we prefer. Later we will comment on how to recover the metric
perturbations in terms of the master functions of the Darboux branch.

There is a systematic way of carrying out the metric reconstruction of
$\tilde{h}_{a}$. It consists in writing the gauge-invariant metric
perturbations as a linear combination of the master functions and their first
derivatives. In this sense, the procedure is similar to the method used in
[58] to study the space of master functions and equations. To avoid
ambiguities we have to take into account that the Regge-Wheeler master
function is the time derivative of the Cunningham-Price-Moncrief master
function [see Eq. (74)], which is a consequence of the Einstein perturbative
field equations. Then, the most general combination is:

$\displaystyle\tilde{h}_{a}$ $\displaystyle=$
$\displaystyle\left(A_{1}t_{a}+A_{2}r_{a}\right)\Psi_{\mathrm{RW}}+\left(A_{3}t_{a}+A_{4}r_{a}\right)\Psi_{\mathrm{CPM}}$
(139) $\displaystyle+$
$\displaystyle\left(A_{5}t_{a}+A_{6}r_{a}\right)t^{b}\Psi_{{\mathrm{RW}}:b}+\left(A_{7}t_{a}+A_{8}r_{a}\right)r^{b}\Psi_{{\mathrm{RW}}:b}$
$\displaystyle+$
$\displaystyle\left(A_{9}t_{a}+A_{10}r_{a}\right)r^{b}\Psi_{{\mathrm{CPM}}:b}$
$\displaystyle+$
$\displaystyle\left(A_{11}t_{a}+A_{12}r_{a}\right)t^{b}\mathcal{S}_{b}+\left(A_{13}t_{a}+A_{14}r_{a}\right)r^{b}{S}_{b}$
$\displaystyle+$
$\displaystyle\left(A_{15}t^{a}+A_{16}r^{a}\right)\mathcal{S}\,.$

There are some comments in order. First, as we can see, there are no terms
proportional to $t^{a}\Psi_{{\mathrm{CPM}}:a}$ as we mentioned that this
quantity is essentially $\Psi_{\mathrm{RW}}$ by virtue of Eq. (74). Second,
since we are going to use the perturbative Einstein equations, it is important
to introduce the corresponding matter terms of the odd-parity sector, the last
two lines in Eq. (139). Then, all we have to do now is to introduce into this
combination the expressions for the master functions $\Psi_{\mathrm{ZM}}$ and
$\Psi_{\mathrm{CPM}}$ [Eqs. (72) and (73) respectively]. It is then a matter
of using the perturbative Einstein field equations [see Appendix C] and look
for the coefficients $A_{I}$ ($I=1,\ldots,16$). One can see that there is a
unique solution that corresponds to

$\displaystyle\tilde{h}_{a}$ $\displaystyle=$
$\displaystyle\frac{r}{f}\Psi_{\mathrm{RW}}r_{a}-\frac{r^{b}}{2f}\left(r\Psi_{\mathrm{CPM}}\right)_{:b}t_{a}$
(140) $\displaystyle-$
$\displaystyle\frac{2r^{2}}{(\ell+2)(\ell-1)f}(t^{b}\mathcal{S}_{b})t_{a}\,.$

As we have mentioned before, this expression can be rewritten entirely in
terms of $\Psi_{\mathrm{CPM}}$ and the odd-parity source terms:

$\displaystyle\tilde{h}_{a}=\frac{1}{2}\varepsilon_{a}{}^{b}\left(r\Psi_{\mathrm{CPM}}\right)_{:b}+\frac{2r^{2}}{(\ell+2)(\ell-1)}\mathcal{S}_{a}\,.$
(141)

Note that if we would like to express $\tilde{h}_{a}$ only in terms of
$\Psi_{\mathrm{RW}}$, we would end up with time integrals, which may be more
difficult to handle in practical situations. This shows again the preferred
role of $\Psi_{\mathrm{CPM}}$ over $\Psi_{\mathrm{RW}}$.

To sum up, combining equations (137) and (138) with Eq. (140) we can write the
reconstruction of the odd-parity metric perturbations $(h_{a},h_{2})$ as
follows:

$\displaystyle h_{a}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\varepsilon_{a}{}^{b}\left(r\Psi_{\mathrm{CPM}}\right)_{:b}+\frac{2r^{2}}{(\ell+2)(\ell-1)}\mathcal{S}_{a}$
(142) $\displaystyle+$
$\displaystyle\frac{r^{2}}{2}\mathcal{G}^{\mathrm{odd}}_{:a}\,,$
$\displaystyle h_{2}$ $\displaystyle=$ $\displaystyle
r^{2}\mathcal{G}^{\mathrm{odd}}\,.$ (143)

###  IV.2 Metric Reconstruction in the Even-Parity Sector

The first step is, again, to identify the components of the metric
perturbations that are pure gauge. To that end, we have to look at the
structure of the gauge transformation, which in the even-parity case is given
by Eqs. (51)-(54). From Eqs. (52) and (54), we can immediately see that the
even-parity metric perturbations $\mbox{\rm\j}_{a}$ and $G$ can be canceled by
a gauge transformation. Then, in the even-parity case, we can write the metric
perturbations as follows:

$\displaystyle h_{ab}$ $\displaystyle=$
$\displaystyle\tilde{h}_{ab}+2\mathcal{H}^{\mathrm{even}}_{(a:b)}\,,$ (144)
$\displaystyle\mbox{\rm\j}_{a}$ $\displaystyle=$
$\displaystyle\mathcal{G}^{\mathrm{even}}_{a}\,,$ (145) $\displaystyle K$
$\displaystyle=$
$\displaystyle\tilde{K}-\frac{\ell(\ell+1)}{2}\mathcal{G}^{\mathrm{even}}+2\frac{r^{a}}{r}\mathcal{H}^{\mathrm{even}}_{a}\,,$
(146) $\displaystyle G$ $\displaystyle=$
$\displaystyle\mathcal{G}^{\mathrm{even}}\,.$ (147)

where

$\mathcal{H}^{\mathrm{even}}_{a}=\mathcal{G}^{\mathrm{even}}_{a}-\frac{r^{2}}{2}\,\mathcal{G}^{\mathrm{even}}_{:a}\,,$
(148)

and here, $(\mathcal{G}^{\mathrm{even}}_{a},\mathcal{G}^{\mathrm{even}})$ are
the three even-parity gauge functions that we have in this case.

The second step consists in writing the gauge-invariant metric perturbations,
$\tilde{h}_{ab}$ and $\tilde{K}$, in terms of the master functions:
$\Psi_{\mathrm{ZM}}$ and $\Psi_{\mathrm{NE}}$, that is, we would like to have:
$\tilde{h}_{ab}=\tilde{h}_{ab}[r,\Psi_{\mathrm{ZM}},\Psi_{\mathrm{NE}}]$ and
$\tilde{K}=\tilde{K}[r,\Psi_{\mathrm{ZM}},\Psi_{\mathrm{NE}}]$. We are going
to use the same systematic procedure as in the odd-parity sector of the
perturbations. First, we need to construct a linear combination of the master
functions and their first derivatives that matches the gauge-invariant metric
perturbations. To prevent ambiguities in the construction we have to take into
account that the master function $\Psi_{\mathrm{NE}}$ is the time derivative
of the Zerilli-Moncrief master function [see Eq. (81)], which is a consequence
of the Einstein perturbative field equations. Let us start with the simplest
gauge-invariant even-parity metric perturbation, i.e. $\tilde{K}$. In this
case, the most general combination is:

$\displaystyle\tilde{K}$ $\displaystyle=$ $\displaystyle
B_{1}\Psi_{\mathrm{ZM}}+B_{2}\Psi_{\mathrm{NE}}+B_{3}r^{a}\Psi_{{\mathrm{ZM}}:a}+B_{4}t^{a}\Psi_{{\mathrm{NE}}:a}$
(149) $\displaystyle+$ $\displaystyle B_{5}r^{a}\Psi_{{\mathrm{NE}}:a}$
$\displaystyle+$ $\displaystyle
B_{6}\mathcal{T}+B_{7}\mathcal{P}+B_{8}t^{a}\mathcal{P}_{a}+B_{9}r^{a}\mathcal{P}_{a}$
$\displaystyle+$ $\displaystyle
B_{10}\mathcal{Q}_{ab}t^{a}t^{b}+B_{11}\mathcal{Q}_{ab}t^{a}r^{b}+B_{12}\mathcal{Q}_{ab}r^{a}r^{b}\,.$

Notice that we have not included a term proportional to
$t^{a}\Psi_{{\mathrm{ZM}}:a}$ since this, by virtue of Eq. (81), is equal to
$\Psi_{\mathrm{NE}}$. As in the treatment of the odd-parity sector, we have
introduced terms proportional to the even-parity components of the energy-
momentum tensor. They act as counter-terms as we have to use Einstein’s
perturbative equations [see Appendix C] and the components of the energy-
momentum tensor will appear and need to be compensated. Then, we can determine
the coefficients $B_{I}$ ($I=1,\ldots,12$) by introducing the expressions for
the master functions $\Psi_{\mathrm{ZM}}$ and $\Psi_{\mathrm{NE}}$ in terms of
the gauge-invariant metric perturbations [Eqs. (78) and (80) respectively]
into Eq. (149) and using the Einstein perturbative equations. As in the odd-
parity sector, by forcing the vanishing of the coefficients of the remaining
metric perturbations and their derivatives and of components of the energy-
momentum tensor, we find that there is a unique solution for the coefficients
$B_{I}$ so that we obtain the following expression for $\tilde{K}$:

$\displaystyle\tilde{K}$ $\displaystyle=$
$\displaystyle\left(\frac{\ell(\ell+1)}{2r}-W(r)\right)\Psi_{\mathrm{ZM}}+r^{a}\Psi_{{\mathrm{ZM}}:a}$
(150) $\displaystyle-$
$\displaystyle\frac{4r^{2}}{\ell(\ell+1)\lambda}\mathcal{Q}_{ab}t^{a}t^{b}\,,$

where $W(r)$ is the superpotential shown in Eq. (132). Note that this
expression for $\tilde{K}$ is consistent with the expressions we derived in
the study of the Darboux covariance of the metric perturbations [See Eqs.
(111) and (120) in Sec. III.3, and Eq. (89) for the definition of
$\Phi_{\mathrm{even}}(x^{a})$].

Let us now follow the same procedure for the reconstruction of the other gauge
invariant even-parity metric perturbations, i.e. $\tilde{h}_{ab}$. The
expression that we need to consider for $\tilde{h}_{ab}$ in terms of the
master functions and their first-order derivatives, as well as of the even-
parity components of the energy-momentum tensor, has the following structure:

$\displaystyle\tilde{h}_{ab}$ $\displaystyle=$
$\displaystyle\left(D_{1}t_{a}t_{b}+2D_{2}t_{(a}r_{b)}+D_{3}r_{a}r_{b}\right)\Psi_{\mathrm{ZM}}$
(154) $\displaystyle+$
$\displaystyle\left(D_{4}t_{a}t_{b}+2D_{5}t_{(a}r_{b)}+D_{6}r_{a}r_{b}\right)\Psi_{\mathrm{NE}}$
$\displaystyle+$
$\displaystyle\left(D_{7}t_{a}t_{b}+2D_{8}t_{(a}r_{b)}+D_{9}r_{a}r_{b}\right)r^{c}\Psi_{{\mathrm{ZM}}:c}$
$\displaystyle+$
$\displaystyle\left(D_{10}t_{a}t_{b}+2D_{11}t_{(a}r_{b)}+D_{12}r_{a}r_{b}\right)t^{c}\Psi_{{\mathrm{NE}}:c}$
$\displaystyle+$
$\displaystyle\left(D_{13}t_{a}t_{b}+2D_{14}t_{(a}r_{b)}+D_{15}r_{a}r_{b}\right)r^{c}\Psi_{{\mathrm{NE}}:c}$
$\displaystyle+$ $\displaystyle\left(\begin{array}[]{l}\mbox{General
combination of $T_{\mu\nu}$}\\\ \mbox{components. 21 more coefficients:}\\\
\mbox{$D_{16}\,$,\ldots, $D_{36}\,$.}\end{array}\right)_{ab}\,,\quad$

where again, to avoid ambiguities, we have not included terms proportional to
$t^{a}\Psi_{{\mathrm{ZM}}:a}$. The last line refers to a general combination
of the components of the energy-momentum tensor that we have not made explicit
for the sake of brevity and that will act as counter-terms when we use
Einstein’s perturbative equations. The only difference with the previous case,
the reconstruction of $\tilde{K}$, is that now we have a tensorial equation
and the number of coefficients to be determined is significantly higher:
$D_{I}$ ($I=1,\ldots 36$). The final result can be written in the following
form [where we have used the relationship between and $\Psi_{\mathrm{NE}}$ and
the time derivative of $\Psi_{\mathrm{ZM}}$ in Eq. (81)]:

$\displaystyle\tilde{h}_{ab}$ $\displaystyle=$
$\displaystyle\left[-\frac{rV_{\mathrm{Z}}}{2f^{2}}\Psi_{\mathrm{ZM}}+\frac{1}{\lambda
f^{3/2}}r^{c}\left(\frac{r\lambda}{\sqrt{f}}r^{d}\Psi_{{\mathrm{ZM}}:d}\right)_{:c}\right]\left(t_{a}t_{b}+r_{a}r_{b}\right)-\frac{4}{\lambda
f^{3/2}}r^{c}\left(\frac{r\lambda}{\sqrt{f}}\Psi_{\mathrm{NE}}\right)_{:c}\;t_{(a}r_{b)}$
(155) $\displaystyle+$
$\displaystyle\left\\{\frac{2r^{2}}{\ell(\ell+1)f^{2}\lambda^{2}}\left[\left(\lambda+3f\right)^{2}-2f\left(5\lambda+3-\Lambda
r^{2}+\frac{3f}{2}\right)\right]\mathcal{Q}_{cd}t^{c}t^{d}-\frac{4r^{3}}{\ell(\ell+1)f^{2}\lambda}r^{c}\left(\mathcal{Q}_{cd}t^{c}t^{d}\right)_{:c}\right\\}\left(t_{a}t_{b}+r_{a}r_{b}\right)$
$\displaystyle-$
$\displaystyle\frac{4r^{2}}{f^{2}\lambda}\left[\mathcal{Q}_{cd}t^{c}r^{d}+\frac{2f}{r}t^{c}\mathcal{P}_{c}\right]t_{(a}r_{b)}+\frac{2r^{2}}{f}\mathcal{P}t_{a}t_{b}\,.$

This expression deserves several comments. First of all, in the same way as
all the previous metric reconstructions that we have carried out, it has two
parts. A gravitational part that can be expressed only in terms of the master
functions and their derivatives along $r^{a}$, and a matter part that can be
expressed only in terms of the components of the energy-momentum tensor and
its first-order derivatives. This is quite relevant and useful for different
purposes. In particular, in the case of solving the master equations in the
time domain, that is, evolving Eq. (69) using a Cauchy evolution from initial
data on an initial time slice. Then, we can always reconstruct the metric
perturbations at any time sclice by taking the values of the master functions
and their spatial (radial) derivatives on that slice. A particularly
interesting feature of the reconstruction of $\tilde{h}_{ab}$ [Eq. (155)] is
that it contains second-order radial derivatives of the master function
$\Psi_{\mathrm{ZM}}$. This means we can use the master equation (69) to change
this expression, but this would bring second-order time derivatives of the
same master function and the first term would retain the same structure but
with a different coefficient.

Another remarkable fact of the metric reconstruction of $\tilde{h}_{a}$, Eq.
(140), is that it uses both master functions, $\Psi_{\mathrm{ZM}}$ and
$\Psi_{\mathrm{NE}}$. This change, as we already mentioned in the analysis of
the odd-parity sector, if we use Eq. (81) to express the metric reconstruction
in terms of $\Psi_{\mathrm{ZM}}$ only, but at the price of also having time
derivatives of this master function. An expression in terms of
$\Psi_{\mathrm{NE}}$ would unavoidably involve time integrals.

Finally, a quite remarkable fact about the gravitational part of the
reconstruction of $\tilde{h}_{ab}$ is that it has only two independent
components, namely one corresponding to $t_{a}t_{b}+r_{a}r_{b}$, and the other
one corresponding to $2t_{(a}r_{b)}$. Not only that, both components are
clearly traceless, and hence, it turns out that in the absence of matter
sources, so it is $\tilde{h}_{ab}$:

$g^{ab}\tilde{h}_{ab}=0\,.$ (156)

This is obviously a consequence of the perturbative Einstein equations for the
even-parity sector, but it is not completely obvious that it has to be this
way. In the general case, the trace of $\tilde{h}_{ab}$ is given by the last
term in its expression [Eq. (155)], that is:

$g^{ab}\tilde{h}_{ab}=-2r^{2}\,\mathcal{P}\,.$ (157)

This equation is one of the components of the perturbative Einstein field
equations [see Eq. (202) in Appendix C].

To sum up, using the metric reconstructions of Eqs. (150) and (155), the final
form of the metric reconstruction for even-parity metric perturbations is:

$\displaystyle h_{ab}$ $\displaystyle=$
$\displaystyle\left[-\frac{rV_{\mathrm{Z}}}{2f^{2}}\Psi_{\mathrm{ZM}}+\frac{1}{\lambda
f^{3/2}}r^{c}\left(\frac{r\lambda}{\sqrt{f}}r^{d}\Psi_{{\mathrm{ZM}}:d}\right)_{:c}\right]\left(t_{a}t_{b}+r_{a}r_{b}\right)-\frac{4}{\lambda
f^{3/2}}r^{c}\left(\frac{r\lambda}{\sqrt{f}}\Psi_{\mathrm{NE}}\right)_{:c}\;t_{(a}r_{b)}+2\mathcal{H}^{\mathrm{even}}_{(a:b)}$
(158) $\displaystyle+$
$\displaystyle\left\\{\frac{2r^{2}}{\ell(\ell+1)f^{2}\lambda^{2}}\left[\left(\lambda+3f\right)^{2}-2f\left(5\lambda+3-\Lambda
r^{2}+\frac{3f}{2}\right)\right]\mathcal{Q}_{cd}t^{c}t^{d}-\frac{4r^{3}}{\ell(\ell+1)f^{2}\lambda}r^{c}\left(\mathcal{Q}_{cd}t^{c}t^{d}\right)_{:c}\right\\}\left(t_{a}t_{b}+r_{a}r_{b}\right)$
$\displaystyle-$
$\displaystyle\frac{4r^{2}}{f^{2}\lambda}\left[\mathcal{Q}_{cd}t^{c}r^{d}+\frac{2f}{r}t^{c}\mathcal{P}_{c}\right]t_{(a}r_{b)}+\frac{2r^{2}}{f}\mathcal{P}t_{a}t_{b}\,,$
$\displaystyle\mbox{\rm\j}_{a}$ $\displaystyle=$
$\displaystyle\mathcal{G}^{\mathrm{even}}_{a}\,,$ (159) $\displaystyle K$
$\displaystyle=$
$\displaystyle\left(\frac{\ell(\ell+1)}{2r}-W(r)\right)\Psi_{\mathrm{ZM}}+r^{a}\Psi_{{\mathrm{ZM}}:a}-\frac{\ell(\ell+1)}{2}\mathcal{G}^{\mathrm{even}}+2\frac{r^{a}}{r}\mathcal{H}^{\mathrm{even}}_{a}-\frac{4r^{2}}{\ell(\ell+1)\lambda}\mathcal{Q}_{ab}t^{a}t^{b}\,,$
(160) $\displaystyle G$ $\displaystyle=$
$\displaystyle\mathcal{G}^{\mathrm{even}}\,.$ (161)

###  IV.3 Metric Reconstruction from Darboux Branch Master Functions

An obvious question in this construction is whether we can do the metric
reconstruction shown before but, instead of in terms of the master functions
of the standard branch, in terms of the master functions of the Darboux
branch. The answer is positive, but given that it involves longer computations
than in the case of the standard branch case that we have fully developed
here, we are only going to show that is possible to do it and the steps one
would need to follow. The main element is precisely the results we have for
the standard branch. We are going to see how we can use them to get to the
metric reconstruction in terms of the Darboux branch master functions. As we
have mentioned in Sec. III.3, the Darboux transformation is reversible in such
a way that the function generating the inverse Darboux transformation is just
$-g$, being $g$ the Darboux generating function (100). Therefore, we can just
exploit this property to invert Eqs. (122) and (123) to obtain:

$\displaystyle\Psi_{\mathrm{CPM}}$ $\displaystyle=$
${}_{D}\check{\Psi}_{\mathrm{odd},x}-g^{\mathrm{odd}}\,_{D}\check{\Psi}_{\mathrm{odd}}\,,$
(162) $\displaystyle\Psi_{\mathrm{ZM}}$ $\displaystyle=$
${}_{D}\check{\Psi}_{\mathrm{even},x}-g^{\mathrm{even}}\,_{D}\check{\Psi}_{\mathrm{even}}\,.$
(163)

Then, metric reconstruction in terms of ${}_{D}\check{\Psi}_{\mathrm{odd}}$
and ${}_{D}\check{\Psi}_{\mathrm{even}}$ can found in a straightforward way by
substituting these inverted Darboux transformations into the standard branch
metric reconstruction given in Eqs. (142)-(143) and (158)-(161). Notice
however that, since the Darboux transformations of $\Psi_{\mathrm{RW}}$ and
$\Psi_{\mathrm{NE}}$ [see Eq. (119) and the discussion around it in Sec.
III.3] do not belong to the space of master functions and equations that we
have considered, we can do this only with the metric reconstruction expressed
entirely in terms of $\Psi_{\mathrm{CPM}}$ and $\Psi_{\mathrm{ZM}}$.

##  V Discussion

In this paper, building on previous works presented in [58, 59, 67, 68], we
have dealt with the problem of reconstructing the metric perturbations of
$\Lambda$-vacuum spherically-spacetimes, generalizing previous results in the
literature in several ways: (i) Our analysis includes all the spherically-
symmetric background spacetimes described in Appendix A. (ii) We have included
in the analysis the presence of matter sources, at the first perturbative
order, by considering a completely general energy-momentum tensor. (iii) We
have performed the metric reconstruction in a way that is completely
independent of the choice of perturbative gauge. (iv) All the results we
present have been put in the context of the space of possible master functions
and equations, taking into account that all the elements of this space are
connected by a hidden symmetry, Darboux covariance, which is the key
ingredient to show that all the possible perturbative formulations that this
space allows are physically equivalent. Moreover, we have shown how to perform
the metric reconstruction in terms of any of the possible master functions in
the space we described. (v) All the expressions we present are covariant with
respect to changes of coordinates in $M^{2}$ (and of course, also in $S^{2}$).

The emerging point of view is that, in principle (up to considerations of
continuity and differentiability of the matter sources), we can always work
with the master equations and functions. The master functions represent the
true degrees of freedom of the gravitational field (remember the background is
not dynamical in the sense that the degrees of freedom of the gravitational
field are shut down). Then, it is in the reconstruction of the metric
perturbations where we can specify the particular perturbative gauge. Once the
gauge-invariant variables, which are reconstructed from the master functions,
are recovered, we have to solve equations/conditions that depend on the
specific choice of perturbative gauge.

Despite all these generalizations, this was not the only aim of this work. As
it is well-known, the master functions represent the true degrees of freedom
of the gravitational field described by the perturbations. Indeed, in terms of
these gauge-invariant quantites we can compute any physical observable
quantity that these perturbations can describe. In particular, the emission of
gravitational waves and their energy-momentum content (energy, angular
momentum, and linear momentum fluxes). Then, the possibility of isolating the
true degrees of freedom in BHPT, gives us some crucial advantages with respect
to the full theory, where several attempts have been described in the
literature as it can be the first step towards quantization. One remarkable
attempt is the introduction of the well-known Arnowitt-Deser-Misner formalism
(ADM) [112], which constitutes a Hamiltonian formalism for General Relativity
and the mathematical basis for the development of a large part of the area of
Numerical Relativity. Despite one cannot isolate the gravitational degrees of
freedom, there are some analogies with BHPT that can be more appreciated in
the framework we have presented here. In the ADM formulation and related ones,
one evolves the intrinsic metric and extrinsic curvature of a three-
dimensional initial Cauchy surfaces (encoding the physical gravitational
information of the system under consideration), and we have four gauge
functions (a scalar, the _lapse_ , and a three-dimensional vector, the
_shift_) that are freely specifiable. In this sense, a particularly
interesting formulation of Numerical Relativity for our purpose of comparison
is the Generalized Harmonic Gauge (GHG) approach introduced by Pretorius [113,
114], which is behind the Numerical Relativity breakthrough of solving the
binary black hole problem (BBH; see also [115, 116] for different Numerical
Relativity formalism that also managed to solve the BBH problem shortly
afterwards; these formulations are based on the formulations of Baumgarte &
Shapiro [117] and Shibata & Nakamura [118], the so-called BSSN formalism, an
extension of the ADM formalism). In the GHG formulation, we can obtain (non-
linear) wave-type equations for a set of metric components and have four
freely specifiable gauge functions. In analogy, in this paper it becomes clear
that in BHPT we can proceed to solve the master equations for the odd- and
even-parity perturbations, which are valid in any gauge, and then, afterwards,
we can reconstruct the full metric perturbations in any perturbative gauge of
interest for the particular problem under consideration. This provides a very
interesting perspective for BHPT computations, where we can always solve
master equations and, afterwards, find the gauge functions that corresponds to
the particular gauge we may be interested in for our physical situation. We
plan to illustrate, in a follow-up study, how this would work for the most
relevant gauge choices used in the literature

Finally, the work presented in this paper can be the basis for developing
applications of BHPT to different types of problems, from gravitational wave
astronomy to fundamental questions of gravitation. This is particularly useful
in cases where either a complete knowledge of the metric perturbations is
required, or when we need to perform computations in a non-standard gauge, in
particular in a gauge different from the Regge-Wheeler one, where all the
gauge functions are chosen to vanish identically. Another extension of this
work is to broaden its scope. By one hand, one can study possible extensions
within General Relativity: For instance, going to second- and higher-order
perturbations within spherical symmetry (see [119, 120, 56] for systematic
treatments of non-linear perturbations), or going beyond spherical symmetry,
in particular to generalize the present framework to the Kerr metric, whose
relevance and importance is well known. On the other hand, there is the
possibility of extending these results to other theories of gravity, where the
study of perturbations of spherically-symmetric systems may be of interest
[121, 102]. A particular example in this line, where a such a perturbative
scheme has been used, is in Ref. [122], where in the context of loop quantum
gravity/cosmology the authors study quantum effects in BH interiors.

###### Acknowledgements.

ML and CFS are supported by contracts PID2019-106515GB-I00 and
PID2022-137674NB-I00 (MCIN/AEI/10.13039/501100011033) and 2017-SGR-1469
(AGAUR, Generalitat de Catalunya). ML is also supported by Juan de la Cierva
contract FJC2021-047289-I funded by program MCIN/AEI/10.13039/501100011033
(Spanish Ministry of Science and Innovation) and by NextGenerationEU/PRTR
(European Union). This work was also partially supported by the program Unidad
de Excelencia María de Maeztu CEX2020-001058-M (Spanish Ministry of Science
and Innovation). We have used the Mathematica software [123] to develop
mathematical tools to manipulate Einstein’s equations for perturbations of
$\Lambda$-vacuum spherically-symmetric spacetimes (see Appendix C) and also to
check most of the computations that appear in this paper.

##  Appendix A The Background Spacetime Metric

Vacuum spherically symmetric solutions are locally determined by Birkhoff’s
uniqueness theorem [124] (published before by Jebsen [125, 126]; see also
[127, 128]). The theorem can be generalized to the case of a non-vanishing
cosmological constant [129] (see [130] for details), what we call a
$\Lambda$-vacuum spacetime.

The possible solutions have to be locally isometric either to one of the
Schwarzschild-de Sitter (SchdS [131]) and Schwarzschild-anti-de Sitter
(SchAdS) solutions or to the Nariai spacetime [132, 133], which can be seen as
the limit of SchdS when the cosmological and event horizons coincide. We do
not consider this case here as it may require a particular treatment.

This family of metrics includes, in the zero mass limit ($M\rightarrow 0$),
the well-knwon maximally-symmetric solutions of Einstein equations: Minkowski
flat spacetime (M; $\Lambda=0$), de Sitter (dS; $\Lambda>0$), and anti-de
Sitter (AdS; $\Lambda<0$).

All these metrics are determined by the function $f(r)$ in Eq. (21). In the
case of the Schwarzschild spacetime [42] (found independently also by Droste
[134]) we have

$f_{\mathrm{Sch}}=1-\frac{r_{s}}{r}\,,$ (164)

where we have introduced the Schwarzschild radius: $r_{s}=2GM/c^{2}=2M\,$. In
the case of de Sitter and anti-de Sitter spacetimes we have

$f_{\mathrm{dS}}=1-\frac{r^{2}}{L^{2}}\,,\qquad
f_{\mathrm{AdS}}=1+\frac{r^{2}}{L^{2}}$ (165)

where we have introduced the length scale $L$, which is related to the
absolute value of the cosmological constant by

$L=\sqrt{\frac{3}{|\Lambda|}}\,,$ (166)

##  Appendix B Spherical Harmonics

Adopting the conventions of [135, 136], the scalar spherical harmonics are
defined by

$Y^{\ell
m}(\theta,\varphi)=\sqrt{\frac{2\ell+1}{4\pi}\frac{\left(\ell-m\right)!}{\left(\ell+m\right)!}}\,P^{m}_{\ell}(\cos{\theta})e^{im\varphi}\,,$
(167)

where $P^{m}_{\ell}$ are the associated Legendre polynomials, which can be
introduced in the following form

$P^{m}_{\ell}(x)=\frac{(-1)^{m}}{2^{\ell}\ell!}(1-x^{2})^{m/2}\frac{d^{\ell+m}}{dx^{\ell+m}}(x^{2}-1)^{\ell}\,.$
(168)

An important property of the scalar spherical harmonics is:

$\bar{Y}^{\ell m}=(-1)^{m}Y^{\ell-m}\,.$ (169)

From the scalar spherical harmonics, we can introduce vector and tensor
spherical harmonics by the definition shown in Eqs. (29)-(32).

###  B.1 Orthogonality properties of Spherical Harmonics

Given the expansion of a scalar/vector/tensor in the corresponding spherical
harmonics, we can extract the coefficients of each mode by using the following
orthogonality properties:

$\displaystyle\int_{S^{2}}d\Omega\,Y^{\ell
m}\bar{Y}^{\ell^{\prime}m^{\prime}}$ $\displaystyle=$
$\displaystyle\delta^{\ell\ell^{\prime}}\delta^{mm^{\prime}}\,,$ (170)
$\displaystyle\int_{S^{2}}d\Omega\,\Omega^{AB}Y^{\ell
m}_{A}\bar{Y}^{\ell^{\prime}m^{\prime}}_{B}$ $\displaystyle=$
$\displaystyle\ell(\ell+1)\delta^{\ell\ell^{\prime}}\delta^{mm^{\prime}}\,,$
(171) $\displaystyle\int_{S^{2}}d\Omega\,\Omega^{AB}X^{\ell
m}_{A}\bar{X}^{\ell^{\prime}m^{\prime}}_{B}$ $\displaystyle=$
$\displaystyle\ell(\ell+1)\delta^{\ell\ell^{\prime}}\delta^{mm^{\prime}}\,,$
(172) $\displaystyle\int_{S^{2}}d\Omega\,\Omega^{AB}Y^{\ell
m}_{A}\bar{X}^{\ell^{\prime}m^{\prime}}_{B}$ $\displaystyle=$ $\displaystyle
0\,,$ (173) $\displaystyle\int_{S^{2}}d\Omega\,\Omega^{AC}\Omega^{BD}T^{\ell
m}_{AB}\bar{T}^{\ell^{\prime}m^{\prime}}_{CD}$ $\displaystyle=$ $\displaystyle
2\delta^{\ell\ell^{\prime}}\delta^{mm^{\prime}}\,,$ (174)
$\displaystyle\int_{S^{2}}d\Omega\,\Omega^{AC}\Omega^{BD}Y^{\ell
m}_{AB}\bar{Y}^{\ell^{\prime}m^{\prime}}_{CD}$ $\displaystyle=$
$\displaystyle\frac{(\ell+2)!}{2(\ell-2)!}\delta^{\ell\ell^{\prime}}\delta^{mm^{\prime}}\,,$
(175) $\displaystyle\int_{S^{2}}d\Omega\,\Omega^{AC}\Omega^{BD}X^{\ell
m}_{AB}\bar{X}^{\ell^{\prime}m^{\prime}}_{CD}$ $\displaystyle=$
$\displaystyle\frac{(\ell+2)!}{2(\ell-2)!}\delta^{\ell\ell^{\prime}}\delta^{mm^{\prime}}\,,$
(176) $\displaystyle\int_{S^{2}}d\Omega\,\Omega^{AC}\Omega^{BD}Y^{\ell
m}_{AB}\bar{X}^{\ell^{\prime}m^{\prime}}_{CD}$ $\displaystyle=$ $\displaystyle
0\,,$ (177) $\Omega^{AC}\Omega^{BD}Y^{\ell
m}_{AB}Y^{\ell^{\prime}m^{\prime}}_{CD}=\Omega^{AC}\Omega^{BD}X^{\ell
m}_{AB}T^{\ell^{\prime}m^{\prime}}_{CD}=0\,.$ (178)

###  B.2 Differential properties of Spherical Harmonics

Given that we have to introduce the expansion of the metric perturbations in
spherical harmonics into the Einstein perturbative equations, it is important
to consider the differential equations they satisfy. In the case of the scalar
harmonics we have Eq. (27) and the definitions of the other harmonics (vector
and tensors) in terms of them. On the other hand, the polar vector harmonics
$Y^{\ell m}_{A}$ satisfy the following differential identities:

$\displaystyle\Omega^{AB}Y^{\ell m}_{A|B}$ $\displaystyle=$
$\displaystyle-\ell(\ell+1)Y^{\ell m}\,,$ (179)
$\displaystyle\Omega^{BC}Y^{\ell m}_{B|CA}$ $\displaystyle=$
$\displaystyle-\ell(\ell+1)Y^{\ell m}_{A}\,,$ (180)
$\displaystyle\Omega^{BC}Y^{\ell m}_{A|BC}$ $\displaystyle=$
$\displaystyle\left[1-\ell(\ell+1)\right]Y^{\ell m}_{A}\,,$ (181)

and the axial vector harmonics $X^{\ell m}_{A}$ satisfy similar differential
identities:

$\displaystyle\Omega^{AB}X^{\ell m}_{A|B}$ $\displaystyle=$ $\displaystyle
0\,,$ (182) $\displaystyle\Omega^{BC}X^{\ell m}_{A|BC}$ $\displaystyle=$
$\displaystyle\left[1-\ell(\ell+1)\right]X^{\ell m}_{A}\,,$ (183)
$\displaystyle\Omega^{BC}X^{\ell m}_{B|AC}$ $\displaystyle=$ $\displaystyle
X^{\ell m}_{A}\,.$ (184)

On the other hand, the (symmetric) polar tensor harmonics $T^{\ell m}_{AB}$
and $Y^{\ell m}_{AB}$ satisfy the following differential identities:

$\displaystyle\Omega^{BC}T^{\ell m}_{BC|A}$ $\displaystyle=$ $\displaystyle
2\,Y^{\ell m}_{A}\,,$ (185) $\displaystyle\Omega^{BC}T^{\ell m}_{AB|C}$
$\displaystyle=$ $\displaystyle Y^{\ell m}_{A}\,,$ (186)
$\displaystyle\Omega^{CD}T^{\ell m}_{AB|CD}$ $\displaystyle=$
$\displaystyle-\ell(\ell+1)T^{\ell m}_{AB}\,,$ (187)
$\displaystyle\Omega^{CD}T^{\ell m}_{CD|AB}$ $\displaystyle=$ $\displaystyle
Y^{\ell m}_{AB}-\ell(\ell+1)T^{\ell m}_{AB}\,.$ (188)
$\displaystyle\Omega^{BC}Y^{\ell m}_{BC|A}$ $\displaystyle=$ $\displaystyle
0\,,$ (189) $\displaystyle\Omega^{BC}Y^{\ell m}_{AB|C}$ $\displaystyle=$
$\displaystyle-\frac{(\ell+2)(\ell-1)}{2}Y^{\ell m}_{A}\,,$ (190)
$\displaystyle\Omega^{CD}Y^{\ell m}_{AB|CD}$ $\displaystyle=$
$\displaystyle\left[4-\ell(\ell+1)\right]Y^{\ell m}_{AB}\,.$ (191)

Finally, the (symmetric) axial tensor harmonics $X^{\ell m}_{AB}$ satisfy the
following differential identities:

$\displaystyle\Omega^{BC}X^{\ell m}_{BC|A}$ $\displaystyle=$ $\displaystyle
0\,,$ (192) $\displaystyle\Omega^{BC}X^{\ell m}_{AB|C}$ $\displaystyle=$
$\displaystyle-\frac{(\ell+2)(\ell-1)}{2}X^{\ell m}_{A}\,,$ (193)
$\displaystyle\Omega^{CD}X^{\ell m}_{AB|CD}$ $\displaystyle=$
$\displaystyle\left[4-\ell(\ell+1)\right]X^{\ell m}_{AB}\,.$ (194)

##  Appendix C Equations for the Harmonic Components of the Metric
Perturbations

Following [83] and [58], the perturbative Einstein equations [see Eq. (13)]
can be decomposed in its spherical harmonic components. For a single
$(\ell,m)$-harmonic, the structure of $\delta G_{\mu\nu}$ is:

$\displaystyle\delta G^{\ell m}_{ab}(x^{c},\Theta^{A})$ $\displaystyle=$
$\displaystyle\mathcal{E}^{\ell m}_{ab}(x^{c})\;Y^{\ell m}(\Theta^{A})\,,$
(195) $\displaystyle\delta G^{\ell m}_{aA}(x^{b},\Theta^{B})$ $\displaystyle=$
$\displaystyle\mathcal{E}^{\ell m}_{a}(x^{b})\;Y^{\ell
m}_{A}(\Theta^{B})+\mathcal{O}^{\ell m}_{a}(x^{b})\;X^{\ell
m}_{A}(\Theta^{B})\,,$ (196) $\displaystyle\delta G^{\ell
m}_{AB}(x^{a},\Theta^{C})$ $\displaystyle=$ $\displaystyle\mathcal{E}^{\ell
m}_{T}(x^{a})\;T^{\ell m}_{AB}(\Theta^{C})+\mathcal{E}^{\ell
m}_{Y}(x^{a})\;Y^{\ell m}_{AB}(\Theta^{C})+\mathcal{O}^{\ell
m}_{X}(x^{a})\;X^{\ell m}_{AB}(\Theta^{C})\,.$ (197)

In this paper, we are considering the presence of matter source at first
perturbative order, described by an arbitrary energy-momentum tensor
$T_{\mu\nu}$ whose decomposition in spherical harmonics in given in Eqs.
(37)-(45). Then, by looking at these decompositions, the perturbative Einstein
equations are equivalent to impose the vanishing of the following harmonic
components, namely: (a) Polar-parity sector: $\mathcal{EE}^{\ell
m}_{ab}\equiv\mathcal{E}^{\ell m}_{ab}-\mathcal{Q}^{\ell m}_{ab}=0$,
$\mathcal{EE}^{\ell m}_{a}\equiv\mathcal{E}^{\ell m}_{a}-\mathcal{P}^{\ell
m}_{a}=0$, $\mathcal{EE}^{\ell m}_{T}\equiv\mathcal{E}^{\ell
m}_{T}-r^{2}\mathcal{T}^{\ell m}=0$, $\mathcal{EE}^{\ell
m}_{Y}\equiv\mathcal{E}^{\ell m}_{Y}-r^{2}\mathcal{P}^{\ell m}=0$. (b) Odd-
parity sector: $\mathcal{OO}^{\ell m}_{a}\equiv\mathcal{O}^{\ell
m}_{a}-\mathcal{S}^{\ell m}_{a}=0$, $\mathcal{OO}^{\ell
m}_{X}\equiv\mathcal{O}^{\ell m}_{X}-\mathcal{S}^{\ell m}=0$. The result is
[for simplicity we drop the harmonic numbers and these equations coincide with
the components of the perturbative Einstein equations up to global factors]:

$\displaystyle\mathcal{EE}_{ab}$ $\displaystyle:$ $\displaystyle
h_{c<a:b>}{}^{:c}-\frac{1}{2}h_{<ab>:c}{}^{:c}-\frac{1}{2}h_{:<ab>}+\frac{r^{c}}{r}\left(2h_{c<a:b>}-h_{<ab>:c}\right)-K_{:<ab>}-\frac{2}{r}r_{<a}K_{:b>}$
(198)
$\displaystyle+\frac{\ell(\ell+1)}{2r^{2}}\left(h_{<ab>}-2\,\mbox{\rm\j}_{<a:b>}\right)-\Lambda
h_{<ab>}-\frac{1}{2}g_{ab}\left[\frac{r^{c}}{r}\left(2h_{cd}{}^{:d}-h_{:c}\right)-K_{:c}{}^{:c}-4\frac{r^{c}}{r}K_{:c}+\frac{\ell(\ell+1)}{2r^{2}}h\right.$
$\displaystyle\left.+\frac{2}{r^{2}}\left(r^{c}r^{d}+rr^{:cd}\right)h_{cd}-\frac{\ell(\ell+1)}{r^{2}}\left(\mbox{\rm\j}_{c}{}^{:c}+2\frac{r^{c}}{r}\mbox{\rm\j}_{c}\right)+\left(\frac{\ell(\ell+1)}{r^{2}}-\frac{2}{r^{2}}(r^{c}r_{c}+rr_{:c}{}^{:c})-2\Lambda\right)K\right.$
$\displaystyle\left.+\frac{(\ell+2)(\ell+1)\ell(\ell-1)}{2r^{2}}G\right]=\mathcal{Q}_{ab}\,,$
$\displaystyle\mathcal{EE}_{a}$ $\displaystyle:$
$\displaystyle\frac{1}{2}h_{ab}{}^{:b}-\frac{1}{2}h_{:a}+\frac{r_{a}}{2r}h+\frac{1}{2}\mbox{\rm\j}^{c}_{:ca}-\frac{1}{2}\mbox{\rm\j}_{a:c}{}^{c}-\frac{r_{a}}{r}\mbox{\rm\j}^{b}_{:b}+\frac{r^{b}}{r}\mbox{\rm\j}_{b:a}+\left[\left(\frac{{}^{2}R}{4}-\Lambda\right)g_{ab}-\frac{r_{:ab}}{r}-\frac{r_{a}r_{b}}{r^{2}}\right]\mbox{\rm\j}^{b}-\frac{1}{2}K_{:a}$
(199) $\displaystyle-\frac{(\ell+2)(\ell-1)}{4}G_{:a}=\mathcal{P}_{a}\,,$
$\displaystyle\mathcal{OO}_{a}$ $\displaystyle:$
$\displaystyle\frac{1}{2}h^{c}_{:ca}-\frac{1}{2}h_{a:c}{}^{c}-\frac{r_{a}}{r}h^{b}_{:b}+\frac{r^{b}}{r}h_{b:a}+\left[\left(\frac{{}^{2}R}{4}-\Lambda\right)g_{ab}-\frac{r_{:ab}}{r}-\frac{r_{a}r_{b}}{r^{2}}\right]h^{b}+\frac{\ell(\ell+1)}{2r^{2}}h_{a}$
(200)
$\displaystyle-\frac{(\ell+2)(\ell-1)}{4r^{2}}\left(h_{2:a}-\frac{2r_{a}}{r}h_{2}\right)=\mathcal{S}_{a}\,,$
$\displaystyle\mathcal{EE}_{T}$ $\displaystyle:$ $\displaystyle
h_{ab}{}^{:ab}-h_{:a}{}^{:a}+\frac{r^{b}}{r}\left(2h_{ba}{}^{:a}-h_{:b}\right)+\left(\frac{\ell(\ell+1)}{2r^{2}}-\Lambda\right)h-\frac{\ell(\ell+1)}{r^{2}}\mbox{\rm\j}^{a}_{:a}-K_{:a}{}^{:a}-2\frac{r^{a}}{r}K_{:a}=-2\mathcal{T}\,,$
(201) $\displaystyle\mathcal{EE}_{Y}$ $\displaystyle:$
$\displaystyle\mbox{\rm\j}^{a}_{:a}-\frac{1}{2}r^{2}G_{:a}{}^{:a}-rr^{a}G_{:a}-\left(rr^{:a}_{:a}+r^{a}r_{a}-1+\Lambda
r^{2}\right)G-\frac{1}{2}h=r^{2}\mathcal{P}\,,$ (202)
$\displaystyle\mathcal{OO}_{X}$ $\displaystyle:$ $\displaystyle
h^{a}_{:a}-\frac{1}{2}h_{2:a}{}^{:a}+\frac{r^{a}}{r}h_{2:a}+\frac{1}{r^{2}}\left(1-2r^{a}r_{a}-\Lambda
r^{2}\right)h_{2}=\mathcal{S}\,,$ (203)

where in these equations we have used the trace-free notation for tensors,
which in two-dimensional case of $M^{2}$ reads: For any symmetric tensor
$A_{ab}$, we define: $A_{<ab>}=A_{ab}-(1/2)g_{ab}A\,$, with $A=g^{ab}A_{ab}$
being the trace.

Given that the Einstein tensor (including the cosmological constant) is gauge-
invariant at first-order in the perturbations (using the Stewart-Walker lemma
and taking into account that it vanishes in the background), we can use the
decomposition of the metric perturbations and rewrite them in terms of purely
gauge invariant quantities. In this way the equations simplify significantly.
For instance, Eq. (202) becomes Eq. (157).

##  Appendix D Equations for the Harmonic Components of the Energy-Momentum
Tensor

The energy-momentum conservation equations come from the second Bianchi
identities [see Eq. (46)]. The spherical-harmonic components of these
equations contain three polar-parity components and a single odd-parity
component, and they are given by

$\displaystyle\mathcal{Q}_{ab}{}^{:b}$ $\displaystyle=$
$\displaystyle-\frac{2r^{b}}{r}\mathcal{Q}_{ab}+\frac{\ell(\ell+1)}{r^{2}}\mathcal{P}_{a}+\frac{2r_{a}}{r}\mathcal{T}\,,$
(204) $\displaystyle\mathcal{P}_{a}{}^{:a}$ $\displaystyle=$
$\displaystyle-2\frac{r^{a}}{r}\mathcal{P}_{a}+\frac{(\ell+2)(\ell-1)}{2}\mathcal{P}-\mathcal{T}\,,$
(205) $\displaystyle\mathcal{S}_{a}{}^{:a}$ $\displaystyle=$
$\displaystyle-2\frac{r^{a}}{r}\mathcal{S}_{a}+\frac{(\ell+2)(\ell-1)}{2r^{2}}\mathcal{S}\,.$
(206)

## References

  * Abbott  _et al._ [2016a] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. 116, 061102 (2016a), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ [2019] B. P. Abbott _et al._ (LIGO Scientific, Virgo), GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs, Phys. Rev. X 9, 031040 (2019), arXiv:1811.12907 [astro-ph.HE] .
  * Abbott  _et al._ [2024] R. Abbott _et al._ (LIGO Scientific, VIRGO), GWTC-2.1: Deep extended catalog of compact binary coalescences observed by LIGO and Virgo during the first half of the third observing run, Phys. Rev. D 109, 022001 (2024), arXiv:2108.01045 [gr-qc] .
  * Abbott  _et al._ [2023] R. Abbott _et al._ (KAGRA, VIRGO, LIGO Scientific), GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo during the Second Part of the Third Observing Run, Phys. Rev. X 13, 041039 (2023), arXiv:2111.03606 [gr-qc] .
  * Abbott  _et al._ [2016b] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Tests of general relativity with GW150914, Phys. Rev. Lett. 116, 221101 (2016b), [Erratum: Phys.Rev.Lett. 121, 129902 (2018)], arXiv:1602.03841 [gr-qc] .
  * Abbott  _et al._ [2016c] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Properties of the Binary Black Hole Merger GW150914, Phys. Rev. Lett. 116, 241102 (2016c), arXiv:1602.03840 [gr-qc] .
  * Sathyaprakash  _et al._ [2012] B. Sathyaprakash _et al._ , Scientific Objectives of Einstein Telescope, Class. Quant. Grav. 29, 124013 (2012), [Erratum: Class.Quant.Grav. 30, 079501 (2013)], arXiv:1206.0331 [gr-qc] .
  * Evans  _et al._ [2021] M. Evans _et al._ , A Horizon Study for Cosmic Explorer: Science, Observatories, and Community, (2021), arXiv:2109.09882 [astro-ph.IM] .
  * Amaro-Seoane  _et al._ [2017] P. Amaro-Seoane _et al._ (LISA), Laser Interferometer Space Antenna, (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Barausse  _et al._ [2020] E. Barausse _et al._ , Prospects for Fundamental Physics with LISA, Gen. Rel. Grav. 52, 81 (2020), arXiv:2001.09793 [gr-qc] .
  * Arun  _et al._ [2022] K. G. Arun _et al._ (LISA), New horizons for fundamental physics with LISA, Living Rev. Rel. 25, 4 (2022), arXiv:2205.01597 [gr-qc] .
  * Colpi  _et al._ [2024] M. Colpi _et al._ , LISA Definition Study Report, (2024), arXiv:2402.07571 [astro-ph.CO] .
  * Afshordi  _et al._ [2023] N. Afshordi _et al._ (LISA Consortium Waveform Working Group), Waveform Modelling for the Laser Interferometer Space Antenna, (2023), arXiv:2311.01300 [gr-qc] .
  * Misner  _et al._ [1973] C. W. Misner, K. Thorne, and J. A. Wheeler,  _Gravitation_ (W. H. Freeman & Co., San Francisco, 1973).
  * Gerlach and Sengupta [1979] U. H. Gerlach and U. K. Sengupta, Gauge invariant perturbations on most general spherically symmetric space-times, Phys. Rev. D 19, 2268 (1979).
  * Gerlach and Sengupta [1980] U. H. Gerlach and U. K. Sengupta, Gauge invariant coupled gravitational, acoustical, and electromagnetic modes on most general spherical space-times, Phys. Rev. D 22, 1300 (1980).
  * Sarbach and Tiglio [2001] O. Sarbach and M. Tiglio, Gauge invariant perturbations of Schwarzschild black holes in horizon-penetrating coordinates, Phys. Rev. D 64, 084016 (2001), gr-qc/0104061 .
  * Clarkson and Barrett [2003] C. A. Clarkson and R. K. Barrett, Covariant perturbations of Schwarzschild black holes, Class. Quant. Grav. 20, 3855 (2003), arXiv:gr-qc/0209051 .
  * Mukohyama [2000] S. Mukohyama, Gauge invariant gravitational perturbations of maximally symmetric space-times, Phys. Rev. D 62, 084015 (2000), arXiv:hep-th/0004067 .
  * Kodama and Ishibashi [2003] H. Kodama and A. Ishibashi, A Master equation for gravitational perturbations of maximally symmetric black holes in higher dimensions, Prog. Theor. Phys. 110, 701 (2003), arXiv:hep-th/0305147 [hep-th] .
  * Nollert [1999] H.-P. Nollert, Quasinormal modes: the characteristic ‘sound’ of black holes and neutron stars, Class. Quantum Grav. 16, 159 (1999).
  * Kokkotas and Schmidt [1999] K. D. Kokkotas and B. G. Schmidt, Quasinormal modes of stars and black holes, Living Rev. Rel. 2, 2 (1999), arXiv:gr-qc/9909058 .
  * Andersson and Kokkotas [2001] N. Andersson and K. D. Kokkotas, The r-mode instability in rotating neutron stars, Int. J. Mod. Phys. D 10, 381 (2001), arXiv:gr-qc/0010102 .
  * Kokkotas and Schutz [1992] K. D. Kokkotas and B. F. Schutz, W-modes: A New family of normal modes of pulsating relativistic stars, Mon. Not. Roy. Astron. Soc. 255, 119 (1992).