# Efficiently Separating the Source of the Teukolsky Equation

Andrew Spiers  School of Mathematical Sciences and School of Physics and
Astronomy, University of Nottingham, University Park, Nottingham NG7 2RD,
United Kingdom;
School of Mathematical Sciences and STAG Research Centre, University of
Southampton, Southampton, SO17 1BJ, United Kingdom

###### Abstract

Recent gravitational wave detections from black hole mergers have underscored
the critical role black hole perturbation theory and the Teukolsky equation
play in understanding the behaviour of black holes. The separable nature of
the Teukolsky equation has long been leveraged to study the vacuum linear
Teukolsky equation; however, as theory and measurements advance, solving the
sourced Teukolsky equation is becoming a frontier of research. In particular,
second-order calculations, such as in quasi-normal mode and self-force
problems, have extended sources. This paper presents a novel method for
efficiently separating the Teukolsky equation’s source analytically, a non-
trivial problem due to the angular and radial mixing of generic quantities in
Kerr spacetime. We provide a proof of concept demonstration of our method’s
accuracy, separating the Teukolsky source produced by the stress-energy tensor
of an ideal gas cloud surrounding a Kerr black hole. The detailed application
of our method is provided in an accompanying Mathematica notebook. Our
approach opens up a new avenue for efficient, accurate black hole perturbation
theory calculations with sources in both the time and frequency domain.

##  I  Introduction

Since the first detection of gravitational waves from a black hole binary [1]
in 2015, gravitational wave astronomy has been a rapidly growing field. The
current catalogue of binary detections [2] is centred on binaries where the
two compact objects have a similar mass. However, recently, the
LIGO/VIRGO/KAGRA collaboration detected a binary with mass-ratio of
$\varepsilon=\frac{\mu}{M}\sim 0.04$ [3]. Gravitational astronomy is now not
confined to the comparable mass regime; the so-called intermediate-mass-ratio
inspirals (IMRIs, with mass ratio $10^{-2}\gtrsim\varepsilon\gtrsim 10^{-4}$)
are increasingly more detectable. Additionally, millihertz frequency
gravitational waves will be detectable with space-based gravitational wave
interferometers, sensitive to lower mass ratios than ground-based detectors.
In the 2030s, three space-based gravitational wave interferometer missions are
planned: LISA [4], TianQin [5], and Taiji [6]. Extreme-mass-ratio inspirals
(EMRIs, $10^{-5}\gtrsim\varepsilon\gtrsim 10^{-8}$) are a key target for these
missions [7].

Numerical Relativity [8] has had great success in accurately modelling
comparable mass binaries. Progress is being made to extend the catalogue of NR
waveforms into the IMRI regime [9, 10]. Synergies with black hole perturbation
theory can also improve efficiency [11, 12]. However, these are challenging
tasks due to the disparate length scales of such systems. In the EMRI regime,
it is unrealistic to expect such methods to be capable of filling the
parameter space of initial conditions in the LISA mission time frame.

Alternatively, black hole perturbation theory provides an ideal method for
approximating the binary spacetime in the small mass-ratio limit. The
evolution of the binary can be modelled using the self-force approach [13]. In
Refs. [14, 15, 16], a two-timescale approximation [17, 18] was implemented for
self-force evolution of a binary, including up to (the dissipative) second-
perturbative-order self-force. This model was limited to a binary of two
Schwarzschild black holes in a quasi-circular inspiral. The waveforms were
produced for a range of mass ratios and showed outstanding agreement with
Numerical Relativity (before the inner-most stable circular orbit), even in
the $\varepsilon\sim\frac{1}{10}$ regime. These results show that self-force
waveforms have incredible potential for gravitational wave science. Crucially,
the dissipative piece of the second-order self-force is a necessary
contribution for high accuracy waveforms. There is also growing interest in
second-order quasi-normal mode calculations [19, 20, 21, 22, 23, 24, 25, 26,
27].

The recent second-order self-force results are limited to a Schwarzschild
background black hole, whereas the primary black hole of astrophysical EMRIs
is expected to have significant spin ($a\sim 0.9M$). Adding the linear in-spin
contributions perturbatively [28, 29] is straightforward. However, including
the full nonlinear spin of the primary object is one of the outstanding
problems in second-order self-force. Working non-linearly in spin involves the
background spacetime being a Kerr black hole [30]. The Kerr metric [30], in
Boyer–Lindquist coordinates [31], is

$\displaystyle ds^{2}$
$\displaystyle=-\Big{(}1-\frac{2Mr}{\Sigma}\Big{)}dt^{2}+\frac{\Sigma}{\Delta}dr^{2}+\Sigma
d\theta^{2}$ $\displaystyle\ \ \
+\Big{(}r^{2}+a^{2}+\frac{2Ma^{2}r\sin[\theta]}{\Sigma}\Big{)}\sin^{2}[\theta]d\phi^{2}$
$\displaystyle\ \ \ -\frac{4Mar\sin^{2}[\theta]}{\Sigma}dtd\phi,$ (1)

where $M$ is the black hole mass, $a$ is the angular momentum per unit mass,
$\Sigma=r^{2}+a^{2}\cos^{2}[\theta]$, and $\Delta=r^{2}-2Mr+a^{2}$. Kerr
spacetime is stationary and axially symmetric, but not spherically symmetric;
hence, there is radial (r) and polar angle ($\theta$) mixing in the metric
(1). This lack of symmetry makes the linearised Einstein field equations
generally non-separable in Kerr; attempting to reduce the linearised Einstein
field equations, a coupled set of PDEs, to a coupled set of ODEs by separation
of variables seems fruitless. However, a separable field equation does exist,
the Teukolsky equation [32, 33]. Solving the sourced Teukolsky equation
separably requires separating the angular and radial dependency in the source.
In this paper, we produce a formalism for efficiently separating general
sources of the Teukolsky equation analytically. Our method promises to
simplify and improve efficiency in solving the second-order Teukolsky equation
for second-order self-force and quasi-normal mode calculations.

Before presenting our formalism, we briefly introduce black hole perturbation
theory and the separability of the Teukolsky equation in the remainder of this
section. In section II, we present our method for decomposing the source. In
section III, we apply our method to a toy example to show it accurately
decomposes sources. Finally, in section IV, we make our concluding remarks.

###  I.1 Black hole perturbation theory and separability

Black hole perturbation theory involves solving the Einstein field equation
perturbatively around a background black hole spacetime. We write Einstein
field equation in natural units as

$\displaystyle G_{ab}[g_{ab}]=8\pi T_{ab},$ (2)

where $G_{ab}$ is the Einstein tensor, $g_{ab}$ is the metric, and $T_{ab}$ is
the stress-energy tensor. The metric and stress-energy tensor are expanded in
orders of the small parameter $\varepsilon$,

$\displaystyle g_{ab}$ $\displaystyle=g^{(0)}_{ab}+\varepsilon
h^{(1)}_{ab}+\varepsilon^{2}h^{(2)}_{ab}+\ldots+\varepsilon^{n}h^{(n)}_{ab}+\ldots,$
(3) $\displaystyle T_{ab}$ $\displaystyle=T^{(0)}_{ab}+\varepsilon
T^{(1)}_{ab}+\varepsilon^{2}T^{(2)}_{ab}+\ldots+\varepsilon^{n}T^{(n)}_{ab}+\ldots$
(4)

We take $g^{(0)}_{ab}$, the background metric, to be the Kerr metric.
Therefore, $T^{(0)}_{ab}=0$, as Kerr is a vacuum solution. Under the
expansions in Eqs. 3 and 4, the Einstein field equation, Eq. 2, can be
expressed as a hierarchical set of linear field equations in ascending order
of $\varepsilon$ [15, 34],

$\displaystyle\delta G_{ab}[h^{(1)}_{ab}]$ $\displaystyle=8\pi T^{(1)}_{ab},$
(5) $\displaystyle\delta G_{ab}[h^{(2)}_{ab}]$ $\displaystyle=8\pi
T^{(2)}_{ab}-\delta^{2}G_{ab}[h^{(1)}_{ab},h^{(1)}_{ab}],$ (6)
$\displaystyle\delta G_{ab}[h^{(3)}_{ab}]$ $\displaystyle=\ldots$ (7)

where $\delta G_{ab}$ is the linearised Einstein operator and
$\delta^{2}G_{ab}$ is the quadratic Einstein operator [34].

The abundant symmetry of Schwarzschild spacetime makes the linearized Einstein
operator, $\delta G_{ab}$, separable [35, 36]. However, the linearized
Einstein operator in Kerr spacetime is generally non-separable. It is,
therefore, remarkable that separable field equations exist, the Teukolsky
equations [32, 33],

$\displaystyle\mathcal{O}\psi_{0}^{(1)}={\mathcal{S}}_{0}T^{(1)}_{ab},\ \ \ \
\mathcal{O}^{\prime}\psi^{(1)}_{4}={\mathcal{S}}_{4}T^{(1)}_{ab},$ (8)

where $\psi_{0}^{(1)}={\mathcal{T}}_{0}h^{(1)}_{ab}$,
$\psi_{4}^{(1)}={\mathcal{T}}_{4}h^{(1)}_{ab}$, and ${\mathcal{T}}_{0}$,
${\mathcal{T}}_{4}$, ${\mathcal{S}}_{0}$, ${\mathcal{S}}_{4}$, $\mathcal{O}$,
and $\mathcal{O}^{\prime}$ are second order differential operators, see Ref.
[37].

It is not immediately obvious that Eq. 8 is separable; one must choose a
tetrad and convert to master Teukolsky equation form [32, 33, 37]. In the
Kinnersley tetrad [38], the master Teukolsky equations are

$\displaystyle\hat{\mathcal{O}}_{+2}[\psi_{0}^{(1)}]$
$\displaystyle=-\Sigma\mathcal{S}_{0}[8\pi T^{(1)}_{ab}]$ (9)
$\displaystyle\hat{\mathcal{O}}_{-2}[\zeta^{4}\psi_{4}^{(1)}]$
$\displaystyle=-\Sigma\zeta^{4}\mathcal{S}_{4}[8\pi T^{(1)}_{ab}],$ (10)

where $\hat{\mathcal{O}}_{s}$ is the spin $s$ master Teukolsky operator and

$\displaystyle\zeta:=r-ia\cos[\theta].$ (11)

The separability ansatz relies on the time dependence of the Teukolsky master
variable, $\psi_{0}^{(1)}$ or $\zeta^{4}\psi_{4}^{(1)}$, being expressed as
proportional to $e^{-i\omega t}$; for example, working in the frequency
domain. The spin $s$ master Teukolsky operator can then be expressed as

$\displaystyle\hat{\mathcal{O}}_{s}=\mathcal{R}_{s}+\Theta_{s},$ (12)

where [37], in Boyer–Lindquist coordinates,

$\displaystyle\mathcal{R}_{s}$
$\displaystyle=\bigg{(}\Delta^{-s}\frac{d}{dr}\Big{(}\Delta^{s+1}\frac{d}{dr}\Big{)}+\frac{K^{2}-2is(r-M)K}{\Delta},$
$\displaystyle\ \ \ +4is\omega r-{}_{s}\lambda_{\ell\mathtt{m}}\bigg{)},$ (13)

where $K:=(r^{2}+a^{2})\omega-a\mathtt{m}$, and

$\displaystyle\Theta_{s}$
$\displaystyle=\bigg{(}\frac{d}{d\chi}\Big{(}(1-\chi^{2})\frac{d}{d\chi}\Big{)}+a^{2}\omega^{2}\chi^{2}-\frac{(\mathtt{m}+s\chi)}{1-\chi^{2}}$
$\displaystyle\ \ \ -2as\omega\chi+s+A\bigg{)},$ (14)

where $\chi:=\cos[\theta]$ and
$A:={}_{s}\lambda_{\ell\mathtt{m}}+2am\omega-a^{2}\omega^{2}$ (and the
eigenvalues ${}_{s}\lambda_{\ell\mathtt{m}}$ depend on $a\omega$). The
eigenfunctions of $\Theta_{s}$ are spin-weighted spheroidal harmonics111In the
limit $a\omega\rightarrow 0$ the spin-weighted spheroidal harmonics reduce to
spin-weighted spherical harmonics. [32, 33],
${}_{s}S_{\ell\mathtt{m}}[\theta,\phi,a\omega]$; that is,

$\displaystyle\Theta_{s}\big{[}{}_{s}S_{\ell\mathtt{m}}\big{]}=0.$ (15)

More generally, the Teukolsky equation can be shown to be separable in any
principle-null direction aligned tetrad in Kerr without using coordinates [39,
37]. The Teukolsky equation can be written as

$\displaystyle\zeta\bar{\zeta}\mathcal{O}=\mathfrak{R}-\mathfrak{S},$ (16)

where $\mathfrak{R}$ and $\mathfrak{S}$ are operators that commute; hence, the
equation is separable.

Research has largely focused on solving the first-order vacuum Teukolsky
equation [40, 41, 42] or point-particle source problem [13, 37]. Hence,
previous research has generally not required efficiently separating the right-
hand side of the Teukolsky equation into modes of the separable Teukolsky
equation. As calculations move to second order, sources are unavoidable. The
second-order Teukolsky equation [43, 44, 34] is sourced by a quadratic
operator acting on the first-order metric perturbation. The reduced second-
order Teukolsky equation takes the form

$\displaystyle\mathcal{O}[\psi_{0L}^{(2)}]$
$\displaystyle=\mathcal{S}_{0}\big{[}8\pi
T^{(2)}_{ab}-\delta^{2}G_{ab}[h^{(1)}_{ab},h^{(1)}_{ab}]\big{]},$ (17)
$\displaystyle\mathcal{O}^{\prime}[\psi_{4L}^{(2)}]$
$\displaystyle=\mathcal{S}_{4}\big{[}8\pi
T^{(2)}_{ab}-\delta^{2}G_{ab}[h^{(1)}_{ab},h^{(1)}_{ab}]\big{]},$ (18)

where $\psi_{0L}^{(2)}={\mathcal{T}}_{0}h^{(2)}_{ab}$,
$\psi_{4L}^{(2)}={\mathcal{T}}_{4}h^{(2)}_{ab}$. Hence, the source is non-
compact, extending from the horizon of the primary to future null infinity.

We need a precise definition for separating the source of the Teukolsky
equation. Take the master Teukolsky equation with a generic source
$f[t,r,\theta,\phi]$,

$\displaystyle\hat{\mathcal{O}}_{s}[\psi]=f,$ (19)

for the spin-$s$ master Teukolsky variable $\psi$. In order to solve the
equation separably, $\psi$ and $f$ must be expressed in a mode decomposition
form:

$\displaystyle\psi[t,r,\theta,\phi]=\sum_{\ell=|s|}^{\infty}\sum_{m=-l}^{l}\psi_{\ell\mathtt{m}}[r]e^{-i\omega
t}{}_{s}S_{\ell\mathtt{m}}[\theta,\phi],$ (20) $\displaystyle
f[t,r,\theta,\phi]=\sum_{\ell=|s|}^{\infty}\sum_{m=-l}^{l}f_{\ell\mathtt{m}}[r]e^{-i\omega
t}{}_{s}S_{\ell\mathtt{m}}[\theta,\phi].$ (21)

Then, Eq. 19 reduces to

$\displaystyle\mathcal{R}_{s}[\psi_{\ell\mathtt{m}}]=f_{\ell\mathtt{m}}$ (22)

for each $\ell$ and $\mathtt{m}$.

As spin-weighted spheroidal-harmonics are an orthonormal set [45] of
eigenvectors on $S^{2}$, one can separate the source of the Teukolsky equation
of spin weight $s$, using

$\displaystyle f_{\ell\mathtt{m}}[r]e^{-i\omega
t}=\oint{}_{s}S^{*}_{\ell\mathtt{m}}f[t,r,\theta,\phi]d\Omega,$ (23)

where ${}_{s}S^{*}_{\ell\mathtt{m}}$ is the bi-orthogonal dual [46] of the
spin-weighted spheroidal harmonic which simplifies to the complex conjugate
when $\omega$ is real. $f_{\ell\mathtt{m}}[t,r]$ in Eq. 23 will then satisfy
Eq. 21. As spin-weighted spheroidal harmonics generally have no known closed
form, integrating Eq. (23) is challenging. One could integrate numerically,
but $f[t,r,\theta,\phi]$ is not generally separable in the radial and polar
angle coordinates. Therefore, the integral will need to be computed
numerically at each radial point on a grid for each $\ell$ and $\mathtt{m}$
mode. This method has been used to separate the source of the Teukolsky
equation in Ref. [27] for a quadratic quasi-normal mode calculation. However,
the inefficiency is problematic for self-force calculations where multiple
modes must be calculated, and second-order calculations must fill a four-
dimensional parameter space of initial data to create the waveform templates
for LISA data analysis.

Alternatively, one could avoid separating the source of the Teukolsky equation
by solving the Teukolsky equation in PDE form [47, 48]. However, numerical
methods for solving PDEs convert the PDE into a system of ODEs. Hence, it is
likely that leveraging the separability of the Teukolsky equation will develop
the most efficient algorithm for separating the Teukolsky equation into ODEs.

##  II Decomposing General functions in Kerr spacetime

Our formalism is designed to separate general sources in Kerr into the sum of
the product of radial functions and spin-weighted spheroidal harmonics.
Essentially, we separate the radial and polar angle dependency222Separating
the time and azimuthal angle dependencies is trivial due to the symmetries of
Kerr spacetime. and express the angular dependence in terms of spin-weighted
spheroidal harmonics.

An immediate simplification we can make is re-expanding the spheroidal
harmonics into an easier basis of functions to work with. Leavers method [49]
is the most accurate method of expressing spin-weighted spheroidal harmonics
but does not expand in terms of spin-weighted spherical harmonics. Press and
Teukolsky [50] were the first to expand spin-weighted spheroidal harmonics.
Ref. [51] then found a more efficient method, which has been implemented in
the Black Hole Perturbation Toolkit [52]. The expansion in terms of spherical
harmonics converges to Leaver’s method solutions and is invertible. Hence, for
the remainder of this paper, we assume any inputs that are generally expressed
in terms of spin-weighted spheroidal harmonics can be re-expanded in terms of
spin-weighted spherical harmonics. Similarly, it is sufficient to express the
source in terms of spin-weighted spherical harmonics as we can re-expand it in
terms of spin-weighted spheroidal harmonics. That is, we can simplify our goal
to expressing the source as

$\displaystyle
f[t,r,\theta,\phi]=\sum_{\ell=|s|}^{\infty}\sum_{\mathtt{m}=-\ell}^{\ell}f_{\ell\mathtt{m}}[r]e^{-i\omega
t}{}_{s}Y_{\ell\mathtt{m}}[\theta,\phi].$ (24)

where ${}_{s}Y_{\ell\mathtt{m}}[\theta,\phi]$ are the spin-weighted spherical
harmonics.

The rest of this section is dedicated to expressing the source of the
Teukolsky equation in the form of Eq. (24). To simplify this task, we use the
Newman–Penrose [53], GHP [54], and Held [55] formalisms. We also find there is
a preferred tetrad and coordinate scheme to express the source in Eq. (24)
form. In such a tetrad and coordinates, almost all the Kerr background
quantities are naturally expressed as single spin-weighted spherical
harmonics, and all angular derivatives become spin-weight raising/lowering
operators. One problematic background quantity is unavoidable, $\rho$ (and
$\bar{\rho}$), where

$\displaystyle\rho=-\frac{1}{\zeta}=\frac{-1}{(r-ia\cos(\theta)}.$ (25)

As the radial and angular dependence appears on the denominator, $\rho$ does
not naturally separate into a radial and angular function. In Sec. II.6, we
separate the radial and angular dependence in Eq. (25) using a Fourier
expansion. Combining these methods, we produce a complete formalism for
expressing the source as a convergent sum of spin-weighted spherical harmonics
(Eq. (24)).

###  II.1 Spin-weighted spherical harmonics

The spin-weighted spherical harmonics were defined in Ref. [56], and we follow
their conventions up to a minus sign in the definition of the spin raising and
lowering operators (similarly to Ref. [36]),

$\displaystyle\hat{\hbox{\ec\char 240\relax}}$
$\displaystyle=(\sin(\theta))^{-s}(\partial_{\theta}+i\csc(\theta)\partial_{\phi})(\sin(\theta))^{s}$
$\displaystyle=\partial_{\theta}+i\csc(\theta)\partial_{\phi}-s\cot(\theta),$
(26) $\displaystyle\hat{\hbox{\ec\char 240\relax}^{\prime}}$
$\displaystyle=(\sin(\theta))^{-s}(\partial_{\theta}-i\csc(\theta)\partial_{\phi})(\sin(\theta))^{s}$
$\displaystyle=\partial_{\theta}-i\csc(\theta)\partial_{\phi}-s\cot(\theta),$
(27)

when acting on a spin-$s$ weighted quantity. One can define the spin-$s$
spherical harmonics from the scalar harmonics,

$\displaystyle{}_{s}Y_{\ell\mathtt{m}}:$
$\displaystyle=\sqrt{\frac{(l-|s|)!}{(l+|s|)!}}(-1)^{s}\hat{\hbox{\ec\char
240\relax}}^{s}Y_{\ell\mathtt{m}},\ \ 0\geq s\geq\ell,$ (28) $\displaystyle:$
$\displaystyle=\sqrt{\frac{(l-|s|)!}{(l+|s|)!}}\hat{\hbox{\ec\char
240\relax}}^{\prime|s|}Y_{\ell\mathtt{m}},\ \ \ \ \ -\ell\geq s\geq 0.$ (29)

The following relations then hold [56],

$\displaystyle\hat{\hbox{\ec\char 240\relax}}\ {}_{s}Y_{\ell,\mathtt{m}}=$
$\displaystyle-\sqrt{(l-s)(l+s+1)}\ {}_{s+1}Y_{\ell,\mathtt{m}},$ (30)
$\displaystyle\hat{\hbox{\ec\char 240\relax}}^{\prime}\
{}_{s}Y_{\ell,\mathtt{m}}=$ $\displaystyle\sqrt{(l+s)(l-s+1)}\
{}_{s-1}Y_{\ell,\mathtt{m}},$ (31) $\displaystyle\hat{\hbox{\ec\char
240\relax}}^{\prime}\hat{\hbox{\ec\char 240\relax}}\
{}_{s}Y_{\ell,\mathtt{m}}=$ $\displaystyle-(l-s)(l+s+1)\
{}_{s}Y_{\ell,\mathtt{m}},$ (32)
$\displaystyle{}_{s}\overline{Y}_{\ell,\mathtt{m}}=$ $\displaystyle(-1)^{s+m}\
{}_{-s}Y_{l,-m}.$ (33)

The product of two spin-weighted spherical harmonics can be expressed as the
sum of spin-weighted spherical harmonics,

$\displaystyle{}_{s_{1}}Y_{\ell_{1},\mathtt{m}_{1}}{}_{s_{2}}Y_{\ell_{2},\mathtt{m}_{2}}=\sum_{\ell_{1}=s_{1},\mathtt{m}_{1}}^{\infty}\sum_{\ell_{2}=s_{2},\mathtt{m}_{2}}^{\infty}C^{\ell,\mathtt{m},s}_{\ell_{1},\mathtt{m}_{1},s_{1},\ell_{2},\mathtt{m}_{2},s_{2}}{}_{s}Y_{\ell,\mathtt{m}}.$
(34)

where
$C^{\ell,\mathtt{m},s}_{\ell_{1},\mathtt{m}_{1},s_{1},\ell_{2},\mathtt{m}_{2},s_{2}}$
is equivalent to a surface integral,

$\displaystyle
C^{\ell,\mathtt{m},s}_{\ell_{1},\mathtt{m}_{1},s_{1},\ell_{2},\mathtt{m}_{2},s_{2}}=\oint{}_{s}\bar{Y}_{\ell,\mathtt{m}}{}_{s_{1}}Y_{\ell_{1},\mathtt{m}_{1}}{}_{s_{2}}Y_{\ell_{2},\mathtt{m}_{2}}.$
(35)

As spin-weighted spherical harmonics are related to Wigner-D matrices, Eq. 35
can be evaluated algebraically using 3j symbols [57, 58, 36],

$\displaystyle
C^{\ell\mathtt{m}s}_{\ell^{\prime}\mathtt{m}^{\prime}s^{\prime}\ell^{\prime\prime}\mathtt{m}^{\prime\prime}s^{\prime\prime}}$
$\displaystyle=(-1)^{\mathtt{m}+s}\sqrt{\frac{(2\ell+1)(2\ell^{\prime}+1)(2\ell^{\prime\prime}+1)}{4\pi}}\times$
$\displaystyle\ \ \ \begin{pmatrix}\ell&\ell^{\prime}&\ell^{\prime\prime}\\\
s&-s^{\prime}&-s^{\prime\prime}\end{pmatrix}\begin{pmatrix}\ell&\ell^{\prime}&\ell^{\prime\prime}\\\
-\mathtt{m}&\mathtt{m}^{\prime}&\mathtt{m}^{\prime\prime}\end{pmatrix}\\!.$
(36)

###  II.2 The Newman–Penrose formalism

The Newman–Penrose (NP) formalism [53] utilises a basis of four null vectors,
called a tetrad, to express curvature quantities. The NP basis vectors are
labelled

$\displaystyle
e_{[a]}^{a}:=\\{e_{[1]}^{a},e_{[2]}^{a},e_{[3]}^{a},e_{[4]}^{a}\\}:=\\{l^{a},n^{a},m^{a},\bar{m}^{a}\\},$
(37)

where indices in square brackets are tetrad indices. The vectors $l^{a}$ and
$n^{a}$ are real, and $m^{a}$ is complex. Overbars denote a complex
conjugation. Conventionally, for a positive metric signature, the orthonormal
relationship of the tetrad takes the form

$l^{a}n_{a}=-1,m^{a}\bar{m}_{a}=1,$ (38)

As the tetrad is orthonormal, following Eq. (38), the metric can be expressed
as [53],

$g_{ab}=-2l_{(a}n_{b)}+2m_{(a}\bar{m}_{b)}.$ (39)

The NP formalism uses Ricci rotation coefficients to express the connection
[59],

$\gamma_{[c][a][b]}=e_{[c]}^{k}e_{[a]k;i}e_{[b]}^{i}.$ (40)

There are 24 independent components of Ricci rotation coefficients. In the NP
formalism, the components are expressed as 12 complex scalars, known as spin
coefficients:

$\displaystyle\kappa=-\gamma_{[3][1][1]},\ \tau=-$
$\displaystyle\gamma_{[3][1][2]},\ \sigma=-\gamma_{[3][1][3]},$
$\displaystyle\ \rho=-\gamma_{[3][1][4]},\pi=-$
$\displaystyle\gamma_{[2][4][1]},\ \nu=-\gamma_{[2][4][2]},$ $\displaystyle\
\mu=-\gamma_{[2][4][3]},$ $\displaystyle\ \lambda=-\gamma_{[2][4][4]},$
$\displaystyle\epsilon=-\frac{\gamma_{[2][1][1]}+\gamma_{[3][4][1]}}{2},$
$\displaystyle\ \gamma=-\frac{\gamma_{[2][1][2]}+\gamma_{[3][4][2]}}{2},$
$\displaystyle\beta=-\frac{\gamma_{[2][1][3]}+\gamma_{[3][4][3]}}{2},$
$\displaystyle\ \alpha=-\frac{\gamma_{[2][1][4]}+\gamma_{[3][4][4]}}{2}].$
(41)

In vacuum spacetimes, the Ricci curvature is zero. The Weyl tensor contains
the vacuum curvature. In the NP formalism, the ten degrees of freedom of the
Weyl tensor are expressed using five complex scalars, known as Weyl scalars:

$\displaystyle\psi_{0}=$ $\displaystyle C_{[1][3][1][3]},$ (42)
$\displaystyle\psi_{1}=$ $\displaystyle C_{[1][3][1][2]},$ (43)
$\displaystyle\psi_{2}=$ $\displaystyle C_{[1][3][4][2]},$ (44)
$\displaystyle\psi_{3}=$ $\displaystyle C_{[1][2][4][2]},$ (45)
$\displaystyle\psi_{4}=$ $\displaystyle C_{[2][4][2][4]}.$ (46)

Kerr spacetime has two principal null vectors (two pairs of principal null
vectors which coincide) [60, 61]. That is, Kerr is Petrov type D. A tetrad can
be chosen such that $l^{a}$ and $n^{a}$ are tangent to the principal null
directions and four of the Weyl scalars and four spin coefficients then vanish
[62, 61],

$\displaystyle\psi_{0}=0,\psi_{1}=0,\psi_{3}=0,\psi_{4}=0,$ (47)
$\displaystyle\kappa=0,\lambda=0,\nu=0,\sigma=0.$ (48)

These simplifications were used to help derive the Teukolsky equation [33].

The covariant derivative is also expressed using tetrad components in the NP
formalism:

$\displaystyle D\eta:=\eta_{|[1]}:=l^{a}\nabla_{a}\eta,\ \ \ \
\Delta\eta:=\eta_{|[2]}:=n^{a}\nabla_{a}\eta,$
$\displaystyle\delta\eta:=\eta_{|[3]}:=m^{a}\nabla_{a}\eta,\ \ \ \
\overline{\delta}\eta:=\eta_{|[4]}:=\overline{m}^{a}\nabla_{a}\eta.$ (49)

###  II.3 The GHP formalism

The GHP formalism builds on the NP formalism, helping to represent the
symmetry in principal null direction aligned tetrads in Petrov type D
spacetimes. While constraining $l^{a}$ and $n^{a}$ to point in the same
direction, two tetrad rotation degrees of freedom remain unconstrained. These
freedoms can be associated with spin and boost transformations and are
isomorphic to the group of multiplication by a complex number, $\vartheta$
[63]. To express the weight of a GHP quantity $f$ we use the notation
$f\circeq\\{p,q\\}$, meaning that under a spin and boost transformation

$f\xrightarrow{}\vartheta^{p}\bar{\vartheta}^{q}f\ \ \Longleftrightarrow\ \
f\circeq\\{p,q\\}.$ (50)

$p$ and $q$ are weights that can be equated to the spin ($s=\frac{1}{2}(p-q)$)
and boost ($b=\frac{1}{2}(p+q)$) weight. The $[b,s]$ weights of the tetrad
vectors are, $[1,0]$, $[-1,0]$, $[0,1]$, and $[0,-1]$ for $l^{a}$, $n^{a}$,
$m^{a}$, and $\bar{m}^{a}$ respectively. Similarly, the $\\{p,q\\}$ weights of
the tetrad vectors are $\\{1,1\\}$, $\\{-1,-1\\}$, $\\{1,-1\\}$, and
$\\{-1,1\\}$ respectively.

Additionally, there is a freedom to interchange $l^{a}$ and $n^{a}$. GHP
introduced a prime operation to represent this interchange. The prime
operation affects the spin and boost weight accordingly,
$f^{\prime}\circeq\\{-p,-q\\}$. Complex conjugation affects the weight as
follows: $\bar{f}\circeq\\{q,p\\}$.

Half of the NP spin coefficients are relabelled using the prime operation
notation,

$\displaystyle\kappa^{\prime}:=-\nu,\ \sigma^{\prime}:=-\lambda,\
\rho^{\prime}:=-\mu,$ $\displaystyle\tau^{\prime}:=-\pi,\
\beta^{\prime}:=-\alpha,\ \epsilon^{\prime}:=-\gamma.$ (51)

The GHP weights of the spin coefficients (and their primes) follow directly
from the weights of the tetrad vectors (using Eqs. (41)); that is,

$\displaystyle\kappa\circeq\\{3,1\\},\ \sigma\circeq\\{3,-1\\},\
\rho\circeq\\{1,1\\},\tau\circeq\\{1,-1\\}.$ (52)

The spin coefficients $\epsilon$, $\epsilon^{\prime}$, $\beta$, and
$\beta^{\prime}$ do not have well-defined weights. Similarly, the NP
derivative operators do not have well-defined weights. GHP found by combining
these poorly defined weight quantities, one can produce derivative operators
with well-defined weights,

$\displaystyle\hbox{\ec\char
222\relax}\eta=(D-p\epsilon-q\bar{\epsilon})\eta,\ \hbox{\ec\char
222\relax}^{\prime}\eta=(\Delta+p\epsilon^{\prime}+q\bar{\epsilon}^{\prime})\eta,$
$\displaystyle\hbox{\ec\char
240\relax}\eta=(\delta-p\beta+q\bar{\beta}^{\prime})\eta,\ \hbox{\ec\char
240\relax}^{\prime}\eta=(\bar{\delta}+p\beta^{\prime}-q\bar{\beta})\eta.$ (53)

Two clear advantages of the GHP formalism are that the equations are more
condensed than in NP form, and they offer a straightforward consistency check
by checking that the weights of an equation are consistent.

We will also make use of the GHP commutation relations (with Eqs. 47 and 48
imposed),

$\displaystyle\hbox{\ec\char 222\relax}^{\prime}\hbox{\ec\char 222\relax}\eta$
$\displaystyle=\hbox{\ec\char 222\relax}\hbox{\ec\char
222\relax}^{\prime}\eta-\eta\big{(}p(\psi_{2}+\tau\tau^{\prime})+q({\bar{\tau}}^{\prime}{\bar{\tau}}-\bar{\psi}_{2})\big{)}$
$\displaystyle\ \ \ -({\bar{\tau}}-\tau^{\prime})\hbox{\ec\char
240\relax}\eta-(\tau-\bar{\tau}^{\prime})\hbox{\ec\char
240\relax}^{\prime}\eta,$ (54) $\displaystyle\hbox{\ec\char
240\relax}\hbox{\ec\char 222\relax}\eta$ $\displaystyle=\hbox{\ec\char
222\relax}\hbox{\ec\char
240\relax}\eta-q\eta{\bar{\rho}}{\bar{\tau}}^{\prime}+{\bar{\tau}}^{\prime}\hbox{\ec\char
222\relax}\eta-{\bar{\rho}}\hbox{\ec\char 240\relax}\eta,$ (55)
$\displaystyle\hbox{\ec\char 240\relax}\hbox{\ec\char 222\relax}^{\prime}\eta$
$\displaystyle=\hbox{\ec\char 222\relax}^{\prime}\hbox{\ec\char
240\relax}\eta+p\rho^{\prime}\tau+\tau\hbox{\ec\char
222\relax}^{\prime}\eta-\rho^{\prime}\hbox{\ec\char 240\relax}\eta.$ (56)

We apply the commutation relations, Eqs. 54, 55 and 56 (and their complex
conjugates) to commute all ð and $\hbox{\ec\char 240\relax}^{\prime}$
derivatives to the right, and all $\hbox{\ec\char 222\relax}^{\prime}$
derivatives to the right of Þ. This simplifies equations and this ordering
will place all spin-raising and lowering operators to the right of radial
derivatives in section II.5.1.

###  II.4 The Held(–Kinnersley) Formalism

The Held formalism builds on the GHP formalism, producing an algorithmic
integration technique for a section of the GHP equations (the Ricci and
Bianchi identities in GHP form). The Held integration method has recently been
applied to solve black hole perturbation theory problems [44, 63, 64]. Here,
we review the main aspects of the Held formalism and clarify its restriction
to the Kinnersly tetrad [38]. In Boyer–Lindquist coordinates, the Kinnersley
tetrad is

$\displaystyle l^{a}$
$\displaystyle=\frac{1}{\Delta}\Big{\\{}r^{2}+a^{2},\Delta,0,a\Big{\\}},$ (57)
$\displaystyle n^{a}$
$\displaystyle=\frac{1}{2\Sigma}\Big{\\{}r^{2}+a^{2},-\Delta,0,a\Big{\\}},$
(58) $\displaystyle m^{a}$
$\displaystyle=\frac{1}{\sqrt{2}\bar{\zeta}}\Big{\\{}ia\sin[\theta],0,1,\frac{i}{\sin[\theta]}\Big{\\}},$
(59) $\displaystyle\bar{m}^{a}$
$\displaystyle=\frac{1}{\sqrt{2}\zeta}\Big{\\{}-ia\sin[\theta],0,1,\frac{-i}{\sin[\theta]}\Big{\\}},$
(60)

I also comment on the advantages of the Held formalism and when the
integration technique can be used.

The Held formalism leverages the Ricci identities (called “Field equations" in
[55]), Bianchi identities, and commutation relations to build a set of
equations for derivatives of background Kerr quantities. For example, in a
vacuum Petrov type-D spacetime, the $R_{[1][3][1][4]}$ Ricci identity [61]
gives

$\displaystyle\hbox{\ec\char 222\relax}\rho$ $\displaystyle=\rho^{2}.$ (61)

This equation is key to the Held integration technique as it allows many
complex differential equations to be solved analytically. For example, one may
solve the complex ordinary differential equation.

$\displaystyle\hbox{\ec\char 222\relax}A[\rho]=B^{\circ}\rho^{3},$ (62)

where $A[\rho]$ is an unknown function of $\rho$ and $B^{\circ}$ is a function
independent of $\rho$ ($\hbox{\ec\char 222\relax}B^{\circ}=0$, all variables
labelled with $\circ$ superscripts share this property and will be referred to
as Held scalars333While quantities labelled with a superscript $\circ$ are
commonly known as Held scalars, the notation was originally introduced by
Kinnersley [38].). Using Eq. 61 one can find the solution,

$\displaystyle A[\rho]=\frac{B^{\circ}\rho^{2}}{2}+C^{\circ},$ (63)

where $C^{\circ}$ is any function independent of $\rho$.

However, Held integration does not apply to general GHP equations. When an
alternative derivative operator (such as $\hbox{\ec\char 222\relax}^{\prime}$,
ð, or $\hbox{\ec\char 240\relax}^{\prime}$) acts on an unknown variable in a
differential equation, the Held method of integration is not helpful. For
example,

$\displaystyle\hbox{\ec\char 240\relax}D[\rho]=E^{\circ}\rho^{3},$ (64)

cannot be solved using Held integration. Within these limitations, the Held
integration method is useful for solving certain problems in black hole
perturbation theory [44, 63, 64].

To apply Held integration to more generic equations than Eq. (63), such as
when GHP derivatives act on spin coefficients, the Held formalism extracts the
$\rho$ dependency out of such quantities. Held began by modifying the
derivative operators of GHP. The Held derivative operators are Þ,
$\tilde{\hbox{\ec\char 222\relax}}^{\prime}$, $\tilde{\hbox{\ec\char
240\relax}}$, and $\tilde{\hbox{\ec\char 240\relax}}^{\prime}$. The new
operators are defined as [44, 63],

$\displaystyle\tilde{\hbox{\ec\char 222\relax}}^{\prime}$
$\displaystyle=\hbox{\ec\char 222\relax}^{\prime}-\bar{\tau}\hbox{\ec\char
240\relax}-\tau\hbox{\ec\char
240\relax}^{\prime}+\tau\bar{\tau}\left(\frac{p}{\bar{\rho}}+\frac{q}{\rho}\right)+\frac{1}{2}\left(\frac{q\bar{\Psi}_{2}}{\bar{\rho}}+\frac{p\Psi_{2}}{\rho}\right),$
(65) $\displaystyle\tilde{\hbox{\ec\char 240\relax}}$
$\displaystyle=\frac{1}{\bar{\rho}}\hbox{\ec\char
240\relax}+\frac{q\tau}{\rho},$ (66) $\displaystyle\tilde{\hbox{\ec\char
240\relax}}^{\prime}$ $\displaystyle=\frac{1}{\rho}\hbox{\ec\char
240\relax}^{\prime}+\frac{p\bar{\tau}}{\bar{\rho}},$ (67)

with GHP weights $\\{-1,-1\\}$, $\\{0,-2\\}$, and $\\{-2,0\\}$ respectively.
These derivative operators were chosen to satisfy the following commutation
relations [55],

$\displaystyle\Big{[}\hbox{\ec\char 222\relax},\tilde{\hbox{\ec\char
222\relax}}^{\prime}\Big{]}\eta^{\circ}=\Big{[}\hbox{\ec\char
222\relax},\tilde{\hbox{\ec\char
240\relax}}\Big{]}\eta^{\circ}=\Big{[}\hbox{\ec\char
222\relax},\tilde{\hbox{\ec\char 240\relax}}^{\prime}\Big{]}\eta^{\circ}=0,$
(68)

where the ∘ superscript notation denotes a quantity that is annihilated by Þ
[38]. Eq. (68) can be interpreted as $\tilde{\hbox{\ec\char
222\relax}}^{\prime}$, $\tilde{\hbox{\ec\char 240\relax}}$, and
$\tilde{\hbox{\ec\char 240\relax}}^{\prime}$ containing no explicit $\rho$
dependency444Held integration is still limited to equations where only  acts
on the unknown quantity; like Eq. 64, $\tilde{\hbox{\ec\char
240\relax}}D[\rho]=E^{\circ}\rho^{3}$ cannot be solved for $D[\rho]$ using
Held integration..

To begin extracting the $\rho$ dependency in the spin coefficients in Kerr,
Held used a relation for $\tau$ taken from Ref. [38],

$\displaystyle\tau=\rho\bar{\rho}\tau^{\circ},$ (69)

with $\tau^{\circ}\circeq\\{-1,-3\\}$. Eq. (69) only holds in the Kinnersley
tetrad as the analysis in Ref. [38] uses

$\displaystyle\epsilon=0.$ (70)

Hence, the Held formalism is in the Kinnersley tetrad. This claim can be
tested by inputting Carter tetrad quantities for
$\rho=-\frac{1}{\zeta}\sqrt{\frac{\Delta}{2\Sigma}}$ and
$\tau=-\frac{ia\sin[\theta]}{\zeta\sqrt{2\Sigma}}$ into Eq. (69). This would
give $\tau^{\circ}=-ia\sin[\theta]\frac{\bar{\zeta}}{\sqrt{\Delta}}$ which is
not annihilated by Þ. Therefore, the Held formalism is not applicable in the
Carter tetrad555The Held formalism being specific to the Kinnersley tetrad is
not made clear in Ref. [55]. In fact, Eq. (5.24) of Ref. [55] tries to use the
tetrad degrees of freedom despite the tetrad being fixed. Commonly, the Held
formalism in modern literature [44, 63, 64] only involves the analysis before
Eq. (5.24) in Ref. [55], so it is consistent. We propose the name
“Held–Kinnersley” formalism to make clear the applicability of the formalism
and give credit to Kinnersley for the work which Held built upon.. As the Held
formalism is specific to the Kinnersley tetrad, it is less flexible than the
GHP formalism and does not respect the symmetries of Kerr spacetime to the
same extent.

To extract the $\rho$ dependency from the other spin coefficients, Held
integrated various Ricci identities, Bianchi identities, and commutation
relations [55], finding

$\displaystyle\tau^{\prime}$ $\displaystyle=-\rho^{2}\bar{\tau}^{\circ},$ (71)
$\displaystyle\psi_{2}$ $\displaystyle=\Psi^{\circ}\rho^{3},$ (72)
$\displaystyle\rho^{\prime}$
$\displaystyle=\rho^{\prime\circ}\bar{\rho}-\frac{1}{2}\Psi^{\circ}\rho^{2}-\left(\hbox{\ec\char
240\relax}^{\prime}\bar{\tau}^{\circ}+\frac{1}{2}\Psi^{\circ}\right)\rho\bar{\rho}-\tau^{\circ}\bar{\tau}^{\circ}\rho^{2}\bar{\rho}.$
(73)

which defines the Held scalars $\Psi^{\circ}\circeq\\{-3,-3\\}$ and
$\rho^{\prime\circ}\circeq\\{-2,-2\\}$.

There is an additional Held scalar that can be derived from the spin
coefficient $\rho$,

$\displaystyle\Omega^{\circ}:=\frac{1}{{\bar{\rho}}}-\frac{1}{\rho},$ (74)

with $\Omega^{\circ}\circeq\\{-1,-1\\}$

At this point, we will deviate from the conventional Held formalism notation.
Due to the Held formalism being in the Kinnersley tetrad, we can express the
Held quantities in a coordinated form. In Boyer-Lindquist coordinates
($\\{t,r,\theta,\phi\\}$) [44, 63],

$\displaystyle\tau^{\circ}$
$\displaystyle=-\frac{ia\sin\theta}{\sqrt{2}},\qquad\rho^{\prime\circ}=-1/2,$
(75) $\displaystyle\qquad\Omega^{\circ}$ $\displaystyle=-2ia\cos\theta,\ \ \ \
\ \Psi^{\circ}=M.$ (76)

Clearly, $\rho^{\prime\circ}$ and $\Psi^{\circ}$ are coordinate invariant. In
the following, we replace $\rho^{\prime\circ}$ and $\Psi^{\circ}$ with $-1/2$
and $M$ respectively. Note, this relabelling does not have a well-defined spin
and boost weight, but this does not affect the form of Held equations because
$\hbox{\ec\char 222\relax}\Psi^{\circ}=\tilde{\hbox{\ec\char
222\relax}}^{\prime}\Psi^{\circ}=\tilde{\hbox{\ec\char
240\relax}}\Psi^{\circ}=\tilde{\hbox{\ec\char
240\relax}}^{\prime}\Psi^{\circ}=\hbox{\ec\char
222\relax}\rho^{\prime\circ}=\tilde{\hbox{\ec\char
222\relax}}^{\prime}\rho^{\prime\circ}=\tilde{\hbox{\ec\char
240\relax}}\rho^{\prime\circ}=\tilde{\hbox{\ec\char
240\relax}}^{\prime}\rho^{\prime\circ}=0$ (from the Ricci and Bianchi
identities [55]).

Held used the remaining Ricci and Bianchi identities (and commutation
relations) to derive simplifications for the Held derivatives acting on the
Kerr–Held scalars.

$\displaystyle\tilde{\hbox{\ec\char 222\relax}}^{\prime}\rho$
$\displaystyle=-\frac{\rho^{2}}{2}-\frac{\rho^{2}M}{2}(\rho+{\bar{\rho}})-\rho^{3}\bar{\rho}\tau^{\circ},$
$\displaystyle\tilde{\hbox{\ec\char 240\relax}}\rho$
$\displaystyle=\rho^{2}\tau^{\circ},\ \tilde{\hbox{\ec\char
240\relax}}^{\prime}\rho=-\rho^{2}{\bar{\tau}}^{\circ},$
$\displaystyle\tilde{\hbox{\ec\char 222\relax}}^{\prime}\tau^{\circ}$
$\displaystyle=0,\ \tilde{\hbox{\ec\char 240\relax}}\tau^{\circ}=0,\
\tilde{\hbox{\ec\char
240\relax}}^{\prime}\tau^{\circ}=\frac{1}{2}\bar{\Omega}^{\circ},$
$\displaystyle\tilde{\hbox{\ec\char 222\relax}}^{\prime}\Omega^{\circ}$
$\displaystyle=0,\ \tilde{\hbox{\ec\char
240\relax}}\Omega^{\circ}=2\tau^{\circ},\ \tilde{\hbox{\ec\char
240\relax}}^{\prime}\Omega^{\circ}=-2{\bar{\tau}}^{\circ}.$ (77)

####  II.4.1 Held derivatives and spin-raising/lowering operators

A further convenience of the Held formalism is that $\tilde{\hbox{\ec\char
240\relax}}$ and $\tilde{\hbox{\ec\char 240\relax}}^{\prime}$ relate to the
spin-raising and lowering operators of the spin-weighted spherical harmonics.
Next, we comment on how this relationship is not unique to the Held formalism
and how the Kinnersley tetrad is one of a select class of tetrads where the
GHP derivatives ð and $\hbox{\ec\char 240\relax}^{\prime}$ relate to the spin-
raising and lowering operators.

In the Kinnersley tetrad, it is not challenging to extract the spin-raising
and lowering operators from the GHP derivatives ð and $\hbox{\ec\char
240\relax}^{\prime}$ (Eq. (53)) when expressed in coordinate form. First,
express the GHP derivatives ð and $\hbox{\ec\char 240\relax}^{\prime}$ in NP
form using Eq. 53. From Eq. 53, we see the relevant NP quantities are
$\delta$, $\beta$, $\alpha$, and their complex conjugates. In the Kinnersley
tetrad and Boyer–Lindquist coordinates, Eqs. 57, 58, 59 and 60 give

$\displaystyle\delta$
$\displaystyle=\frac{1}{\sqrt{2}\bar{\zeta}}\Big{(}ia\sin[\theta]\partial_{t}+\partial_{\theta}+\frac{i}{\sin[\theta]}\partial_{\phi}\Big{)},$
(78) $\displaystyle\beta$
$\displaystyle=\frac{\cot[\theta]}{2\sqrt{2}\bar{\zeta}},\ \
-\alpha=\beta^{\prime}=\frac{\cot[\theta]}{2\sqrt{2}\zeta}-\frac{ia\sin[\theta]}{\sqrt{2}\zeta^{2}}.$
(79)

Noting the common factor of $\frac{1}{\sqrt{2}\bar{\zeta}}$ in Eqs. 78 and 79,
and using Eq. 53, it is straightforward to show the GHP ð relates to the spin
raising operator $\hat{\hbox{\ec\char 240\relax}}$,

ð $\displaystyle=\frac{1}{\sqrt{2}\bar{\zeta}}\hat{\hbox{\ec\char
240\relax}}+\frac{ia\sin[\theta]\partial_{u}}{\sqrt{2}\bar{\zeta}}+\frac{iq\sin[\theta]}{\sqrt{2}\bar{\zeta}^{2}}.$
(80)

The Kinnersley tetrad is not the only tetrad where GHP ð and $\hbox{\ec\char
240\relax}^{\prime}$ relate to spin-raising and lowering operators. This can
be seen by examining the Type-III tetrad rotations [61] (under which the
direction of the vectors $l^{a}$ and $n^{a}$ are unchanged and, therefore,
remain aligned with the principle null directions),

$\displaystyle l^{a}\rightarrow A^{-1}l^{a},\ n^{a}\rightarrow An^{a},\
m^{a}\rightarrow e^{i\theta}m^{a},\ \bar{m}^{a}\rightarrow
e^{-i\theta}\bar{m}^{a},$ (81)

where $A$ and $\theta$ are real functions. The relevant spin coefficients (see
Eq. 53) transform as [61],

$\displaystyle\beta$ $\displaystyle\rightarrow
e^{i\theta}\beta+\frac{1}{2}ie^{i\theta}\delta\theta-\frac{1}{2}A^{-1}e^{i\theta}\delta
A,$ (82) $\displaystyle\beta^{\prime}$ $\displaystyle\rightarrow
e^{-i\theta}\beta^{\prime}-\frac{1}{2}ie^{-i\theta}\bar{\delta}\theta+\frac{1}{2}A^{-1}e^{-i\theta}\bar{\delta}A,$
(83)

Take a tetrad related to the Kinnersley tetrad where $\theta$ and $A$ satisfy
$\delta\theta=\delta A=0$; then $\beta$, $\bar{\beta}^{\prime}$, and $m^{a}$
are only re-scaled by a factor of $e^{i\theta}$, so the relation of the GHP
derivatives ð and $\hbox{\ec\char 240\relax}^{\prime}$ to the spin raising and
lowering operators remains trivial. The Hartle–Hawking tetrad [65, 66] is one
such example, where $A=\frac{2(r^{2}+a^{2})}{\Delta}$ and $\theta=0$. However,
this condition does not hold for all tetrads, such as the Carter tetrad. To
transform from the Kinnersley to the Carter tetrad
$A=\sqrt{\frac{2\Sigma}{\Delta}}$ and
$\theta=-i\ln\Big{[}\frac{\bar{\zeta}}{\sqrt{\Sigma}}\Big{]}$. Hence,

$\displaystyle
m^{a}_{C}=\frac{1}{\sqrt{2\Sigma}}\Big{\\{}ia\sin[\theta],0,1,\frac{i}{\sin[\theta]}\Big{\\}}$
(84)

and the relevant spin coefficients in the Carter tetrad are

$\displaystyle\beta_{C}=\beta^{\prime}_{C}=-\frac{i}{\zeta}\frac{a+ir\cos[\theta]}{2\sin[\theta]\sqrt{2\Sigma}}.$
(85)

Therefore, in the Carter tetrad, the $\cot[\theta]$ terms in $\beta$ and
$\beta^{\prime}$ do not have the same coefficient as the $\theta$ derivative
in $\delta$, so they are not eliminated by the conversion to a spin-raising
operator. Also, the $\frac{a}{\sin[\theta]}$ terms will remain after
expressing the $\theta$ derivatives in terms of raising and lowering
operators. $\cot[\theta]$ and $\frac{1}{\sin[\theta]}$ are singular at the
poles, and when expressed as a sum of spin-weighted spherical harmonics, the
sum does not converge. Therefore, there are preferred tetrads for expressing
Kerr quantities such that they are smooth on the 2-sphere: the Kinnersley
tetrad and tetrads related to the Kinnersley tetrad by Eq. 81 with
$\delta\theta=\delta A=0$.

While the Kinnersley tetrad in Boyer–Lindquist coordinates avoids
singularities in ð and $\hbox{\ec\char 240\relax}^{\prime}$, we can see from
Eqs. 57 and 53 that Þ contains $\frac{1}{\Delta}$ terms. These terms are
singular at the horizons as $\Delta=0$ for

$\displaystyle r_{\pm}=M\pm\sqrt{M^{2}-a^{2}},$ (86)

the outer and inner horizons, respectively. To avoid singularities at the
horizon and make Þ and $\tilde{\hbox{\ec\char 222\relax}}^{\prime}$ as simple
as possible, we work in Kerr–Newman coordinates [44].The outgoing Kerr–Newman
coordinates, $\big{\\{}u,r,\theta,\phi_{*}\big{\\}}$, are related to
Boyer–Lindquist coordinates as follows, [44, 67]

$\displaystyle u$
$\displaystyle=t-r_{*}=t-r-\frac{r_{+}^{2}+a^{2}}{r_{+}-r_{-}}\ln\Big{[}\frac{r-r_{+}}{r_{+}}\Big{]}$
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
+\frac{r_{-}^{2}+a^{2}}{r_{+}-r_{-}}\ln\Big{[}\frac{r-r_{-}}{r_{-}}\Big{]},$
(87) $\displaystyle\phi_{*}$
$\displaystyle=\phi-\frac{a}{r_{+}-r_{-}}\ln\Big{[}\frac{r-r_{+}}{r-r_{-}}\Big{]},$
(88)

In these coordinates [44, 63],

$\displaystyle l^{a}$ $\displaystyle=\Big{\\{}0,1,0,0\Big{\\}},$ (89)
$\displaystyle n^{a}$
$\displaystyle=\frac{1}{\Sigma}\Big{\\{}r^{2}+a^{2},-\frac{\Delta}{2},0,a\Big{\\}},$
(90) $\displaystyle m^{a}$
$\displaystyle=\frac{1}{\sqrt{2}\bar{\zeta}}\Big{\\{}ia\sin[\theta],0,1,\frac{i}{\sin[\theta]}\Big{\\}},$
(91) $\displaystyle\bar{m}^{a}$
$\displaystyle=\frac{1}{\sqrt{2}\zeta}\Big{\\{}-ia\sin[\theta],0,1,\frac{-i}{\sin[\theta]}\Big{\\}};$
(92)

Using Eqs. 89 and 70,

Þ $\displaystyle=\partial_{r},$ (93)

When acting on Held scalars [63]

$\displaystyle\tilde{\hbox{\ec\char 222\relax}}^{\prime}$
$\displaystyle=\partial_{u}.$ (94)

When $\tilde{\hbox{\ec\char 222\relax}}^{\prime}$ acts on a quantity which is
not a Held scalar (i.e., not annihilated by Þ), its coordinate form is more
complicated and dependent on the quantities spin- and boost-weight (which we
address in section II.5).

$m^{a}$ and $\bar{m}^{a}$ in Eqs. 91 and 92 take the same form as in
Boyer–Lindquist coordinates (as do the spin-coefficients as they are scalars
which depend on only $r$ and $\theta$). The Held derivatives
$\tilde{\hbox{\ec\char 240\relax}}$ and $\tilde{\hbox{\ec\char
240\relax}}^{\prime}$ then relate similarly to the spin raising and lowering
operators,

$\displaystyle\tilde{\hbox{\ec\char 240\relax}}$
$\displaystyle=\frac{-1}{\sqrt{2}}\left(\partial_{\theta}+i\csc\theta\partial_{\phi_{*}}+ia\sin\theta\partial_{u}-\frac{1}{2}(p-q)\cot\theta\right)$
$\displaystyle=\frac{-1}{\sqrt{2}}\hat{\hbox{\ec\char
240\relax}}-\frac{ia\sin[\theta]\partial_{u}}{\sqrt{2}},$ (95)
$\displaystyle\tilde{\hbox{\ec\char 240\relax}}^{\prime}$
$\displaystyle=\frac{-1}{\sqrt{2}}\left(\partial_{\theta}-i\csc\theta\partial_{\phi_{*}}-ia\sin\theta\partial_{u}+\frac{1}{2}(p-q)\cot\theta\right)$
$\displaystyle=\frac{-1}{\sqrt{2}}\hat{\hbox{\ec\char
240\relax}}^{\prime}+\frac{ia\sin[\theta]\partial_{u}}{\sqrt{2}},$ (96)

with the subtlety that the spin raising and lowering operators now act on
${}_{s}Y_{\ell\mathtt{m}}[\theta,\phi_{*}]$ rather than
${}_{s}Y_{\ell\mathtt{m}}[\theta,\phi]$. The $\partial_{u}$ derivatives are
trivial as the background is time-independent and we assume the perturbations
have a time dependency proportional to $e^{-i\omega t}$ (converting the time
coordinate dependence is straightforward, $e^{-i\omega
t}=e^{-i\omega(u+r_{*})}$). Sections II.4.1 and 96 are similar to
Chandrasekhar’s operators
$\,{}_{-s}\mathcal{L}_{\omega}{}^{\dagger},\,{}_{+s}\mathcal{L}_{\omega}{}$ up
to a trivial normalization [61].

###  II.5 Half-Held–Kinnersley formalism

Here, we are not interested in the Held formalism for its integration method.
Rather, we are interested in how it separates quantities into an explicitly
$\rho$ dependent piece and a $\rho$ independent piece. We shall maintain this
property while making a small deviation from the conventions of Held,
producing a half-Held–Kinnersley formalism. As expressing
$\tilde{\hbox{\ec\char 222\relax}}^{\prime}$ in coordinate form is not
straightforward, we use the derivative operators Þ, $\hbox{\ec\char
222\relax}^{\prime}$, $\tilde{\hbox{\ec\char 240\relax}}$, and
$\tilde{\hbox{\ec\char 240\relax}}^{\prime}$ instead. Our formalism maintains
the advantageous Held $\tilde{\hbox{\ec\char 240\relax}}$ and
$\tilde{\hbox{\ec\char 240\relax}}^{\prime}$ derivatives, which are spin-
raising and -lowering operators of the spin-weighted spherical harmonics up to
time derivatives. As Þ and $\tilde{\hbox{\ec\char 222\relax}}^{\prime}$ do not
contain $\theta$ derivatives, all $\theta$ derivatives also appear as spin-
raising and lowering operators in the half-Held–Kinnersley formalism.

The relations for $\hbox{\ec\char 222\relax}^{\prime}$ acting on Held
quantities can be derived from the Ricci and Bianchi identities:

$\displaystyle\hbox{\ec\char 222\relax}^{\prime}\rho$
$\displaystyle=-\frac{\rho}{2}(M(2\rho^{2}+\rho{\bar{\rho}}+{\bar{\rho}}^{2})+2(-\rho\rho^{\prime\circ}+\rho^{2}{\bar{\rho}}\tau^{\circ}{\bar{\tau}}^{\circ}$
$\displaystyle\ \ \ +{\bar{\rho}}^{3}\tau^{\circ}{\bar{\tau}}^{\circ},$ (97)
$\displaystyle\hbox{\ec\char 222\relax}^{\prime}\Omega^{\circ}$
$\displaystyle=\frac{1}{2\rho{\bar{\rho}}}((\rho-{\bar{\rho}})(M(\rho^{2}+{\bar{\rho}}^{2})$
$\displaystyle\ \ \
+2(\rho^{3}+2\rho^{2}{\bar{\rho}}+2\rho{\bar{\rho}}^{2}+{\bar{\rho}}^{3})\tau^{\circ}{\bar{\tau}}^{\circ}),$
(98) $\displaystyle\hbox{\ec\char 222\relax}^{\prime}\tau^{\circ}$
$\displaystyle=\frac{1}{2{\bar{\rho}}}(\tau^{\circ}({\bar{\rho}}^{2}(-4\rho^{\prime\circ}+3M{\bar{\rho}})+4\rho^{3}{\bar{\rho}}\tau^{\circ}{\bar{\tau}}^{\circ}$
$\displaystyle\ \ \
+2\rho{\bar{\rho}}(\rho^{\prime\circ}-\Omega^{\circ}\rho^{\prime\circ}{\bar{\rho}}+{\bar{\rho}}^{2}\tau^{\circ}{\bar{\tau}}^{\circ})$
$\displaystyle\ \ \
+\rho^{2}(\rho^{\prime\circ}(2-2\Omega^{\circ}{\bar{\rho}})+{\bar{\rho}}(M+2{\bar{\rho}}\tau^{\circ}{\bar{\tau}}^{\circ})))),$
(99) $\displaystyle\hbox{\ec\char 222\relax}^{\prime}M$
$\displaystyle=\frac{3M}{2}(M(\rho^{2}+{\bar{\rho}}^{2})+2(-\rho^{\prime\circ}{\bar{\rho}}+\rho(\rho^{\prime\circ}-\Omega^{\circ}\rho^{\prime\circ}{\bar{\rho}})$
$\displaystyle\ \ \
+\rho^{3}\tau^{\circ}{\bar{\tau}}^{\circ}+\rho^{2}{\bar{\rho}}\tau^{\circ}{\bar{\tau}}^{\circ})),$
(100) $\displaystyle\hbox{\ec\char 222\relax}^{\prime}\rho^{\prime\circ}$
$\displaystyle=\frac{\rho^{\prime\circ}}{2({\bar{\rho}}+\Omega^{\circ}\rho{\bar{\rho}})}\Big{(}M\big{(}3{\bar{\rho}}^{3}+\rho^{2}{\bar{\rho}}(1-3\Omega^{\circ}{\bar{\rho}})$
$\displaystyle\ \ \
+\rho^{3}(4-2\Omega^{\circ}{\bar{\rho}})+\rho{\bar{\rho}}^{2}(-4+3\Omega^{\circ}{\bar{\rho}})\big{)}+2\rho\big{(}$
$\displaystyle\ \ \
\Omega^{\circ}\rho^{\prime\circ}{\bar{\rho}}({\bar{\rho}}-\rho+\Omega^{\circ}\rho{\bar{\rho}})+\tau^{\circ}{\bar{\tau}}^{\circ}(\rho^{3}-\Omega\rho^{3}{\bar{\rho}}$
$\displaystyle\ \ \
+\rho(-1+\Omega^{\circ}\rho){\bar{\rho}}^{2}+4(1+\Omega^{\circ}\rho){\bar{\rho}}^{3})\big{)}\Big{)}.$
(101)

Equation (101) appears to be problematic as the denominator
$({\bar{\rho}}+\Omega^{\circ}\rho{\bar{\rho}})$ introduces mixing between the
angular and radial coordinates. However, by partially introducing coordinates
using Eqs. 25, 75 and 76, it simplifies to

$\displaystyle\hbox{\ec\char 222\relax}^{\prime}\rho^{\prime\circ}=$
$\displaystyle\frac{-\rho{\bar{\rho}}}{2}\Big{(}2r+a^{2}(\rho+{\bar{\rho}})+r^{2}(\rho+{\bar{\rho}})-2M\big{(}1+r(\rho+{\bar{\rho}})\big{)}\Big{)}.$
(102)

####  II.5.1 Converting to coordinates

To convert expressions in the half-Held–Kinnersley formalism into Kerr–Newmann
coordinates ($\\{u,r,\theta,\phi_{*}\\}$), we first express the derivatives Þ
and $\hbox{\ec\char 222\relax}^{\prime}$ in NP form. From Eq. (53) we see this
re-introduces the derivatives $D$ and $\Delta$, and the spin-coefficients
$\epsilon$ and $\epsilon^{\prime}$. In the Kinnersley tetrad, in Kerr–Newman
coordinates,

$\displaystyle\epsilon$ $\displaystyle=0,$ (103)
$\displaystyle\epsilon^{\prime}$
$\displaystyle=\frac{\rho{\bar{\rho}}}{2}(r-M+(a^{2}-2Mr+r^{2})\rho),$ (104)
$\displaystyle D$ $\displaystyle=\partial_{r},$ (105) $\displaystyle\Delta$
$\displaystyle=\rho{\bar{\rho}}\Big{(}(r^{2}+a^{2})\partial_{u}-\frac{(r^{2}-2Mr+a^{2})}{2}\partial_{r}+a\partial_{\phi_{*}}\Big{)}.$
(106)

After substituting in Eqs. 75 and 76, the resulting coordinate expression
presents all the $\rho$ and ${\bar{\rho}}$ dependence explicitly666One may
encounter $\partial_{i}\rho$ For $i=\\{u,r,\phi_{*}\\}$ terms, it is trivial
to show $\partial_{r}\rho=\rho^{2}$,
$\partial_{u}\rho=\partial_{\phi_{*}}\rho=0$, and similarly for
${\bar{\rho}}$..

####  II.5.2 Mode decomposition

We assume the inputs ($A$) are expressed as a convergent sum of spin-weighted
spherical harmonics with coefficients that can depend on $r$, $\rho$, and
${\bar{\rho}}$:

$\displaystyle
A=\sum_{\ell=|s|}^{\infty}\sum_{\mathtt{m}=-\ell}^{\ell}A_{\ell\mathtt{m}}[r,\rho,{\bar{\rho}}]\
{}_{s}Y_{\ell\mathtt{m}}[\theta,\phi_{*}]\ e^{i\omega u},$ (107)

where $A_{\ell\mathtt{m}}[r,\rho,{\bar{\rho}}]$ are polynomials with positive
and finite powers of $\rho$ and ${\bar{\rho}}$; that is,
$A_{\ell\mathtt{m}}[r,\rho,{\bar{\rho}}]=\rho
A^{1,0}_{\ell\mathtt{m}}[r]+{\bar{\rho}}A^{0,1}_{\ell\mathtt{m}}[r]+\rho{\bar{\rho}}A^{1,1}_{\ell\mathtt{m}}[r]+\ldots+\rho^{i}{\bar{\rho}}^{j}A^{i,j}_{\ell\mathtt{m}}[r]$
for some finite $i$ and $j$. Note that $\phi_{*}$ derivatives, like time
derivatives, are trivial as ${}_{s}Y_{\ell\mathtt{m}}\sim
e^{i\mathtt{m}\phi_{*}}$; that is,

$\displaystyle\partial_{\phi_{*}}A=i\mathtt{m}A.$ (108)

While we can use Eq. Eqs. 75 and 76 to express the Held scalars in coordinate
form, for a mode decomposition, it is more useful to express them in terms of
spin-weighted spherical harmonics,

$\displaystyle\tau^{\circ}=-\sqrt{\frac{4\pi}{3}}a\ {}_{1}Y_{10},$
$\displaystyle\ \ {\bar{\tau}}^{\circ}=-\sqrt{\frac{4\pi}{3}}a\
{}_{-1}Y_{10},$ $\displaystyle\Omega^{\circ}=-4i$
$\displaystyle\sqrt{\frac{\pi}{3}}a\ {}_{0}Y_{10};$ (109)

where we have used

$\displaystyle
Z[\theta]=\sum_{\ell=s,0}^{\infty}c^{\ell,0,s}{}_{s}Y_{\ell,0},$ (110)

with

$\displaystyle c^{\ell,0,s}=\oint{}_{s}\bar{Y}_{\ell,0}Z[\theta].$ (111)

The resulting source expression will have products of many spin-weighted
spherical harmonics. To combine the spin-weighted spherical harmonics, we use
Eq. 34. The resulting expression for the source will be of the form

$\displaystyle
f[t,r,\theta,\phi]=\sum_{\ell=|s|}^{\infty}\sum_{\mathtt{m}=-\ell}^{\ell}\tilde{f}_{\ell\mathtt{m}}[r,\rho,{\bar{\rho}}]e^{-i\omega
t}{}_{s}Y_{\ell\mathtt{m}}[\theta,\phi_{*}].$ (112)

Where $\tilde{f}_{\ell\mathtt{m}}[r,\rho,{\bar{\rho}}]$ is a polynomial with
positive and finite powers of $\rho$ and ${\bar{\rho}}$; that is,
$\tilde{f}_{\ell\mathtt{m}}[r,\rho,{\bar{\rho}}]=\rho\tilde{f}^{1,0}_{\ell\mathtt{m}}[r]+{\bar{\rho}}\tilde{f}^{0,1}_{\ell\mathtt{m}}[r]+\rho{\bar{\rho}}\tilde{f}^{1,1}_{\ell\mathtt{m}}[r]+\ldots+\rho^{i}{\bar{\rho}}^{j}\tilde{f}^{i,j}_{\ell\mathtt{m}}[r]$
for some finite $i$ and $j$. Next, we must express $\rho$ and $\rho$ in terms
of spin-weighted spherical harmonics.

###  II.6 $\rho$ expansion

So far, our method does not completely express quantities separably as factors
where $\rho$ and ${\bar{\rho}}$ are present; see Eq. (25). The mixing between
the radial and polar angle coordinate in the denominator of Eq. (25) makes
separating the $r$ and $\theta$ dependency non-trivial. In this subsection, we
show the $r$ and $\theta$ dependency in $\rho^{i}{\bar{\rho}}^{k}$ (for any
positive integer $i$ and $k$) can be separated using a Fourier
expansion777Alternatively, one can separate the dependency using a two-point
Taylor series [68, 69, 70]; however, we found the Fourier series approximation
more accurate..

As Eq. (25) is a trigonometric function, a Fourier expansion seems like a
natural approach. Eq. (25) is periodic in $\theta$
($\rho[r,\theta]=\rho[r,\theta+2n\pi]$ for any $r$ and $\theta$, where $n$ is
an integer). Also, Eq. (25) is complex and non-singular for $r\geq r_{-}$
(where $r_{-}$ is the radius of the inner horizon). To assess the real and
imaginary behaviour of Eq. (25), one can write

$\displaystyle\rho=\frac{-(r+ia\cos[\theta])}{\Sigma},$ (113)

as

$\displaystyle\frac{1}{\Sigma}=\rho{\bar{\rho}}=\frac{1}{(r^{2}+a^{2}\cos^{2}[\theta])}.$
(114)

All factors of $\rho^{i}{\bar{\rho}}^{k}$ for $i\neq k$ can also be written in
a manner where the denominator is real: for $i\geq k$,

$\displaystyle\rho^{i}{\bar{\rho}}^{k}=\frac{(-(r+ia\cos[\theta]))^{i-k}}{\Sigma^{i}},$
(115)

and for $k\leq i$,

$\displaystyle\rho^{i}{\bar{\rho}}^{k}=\frac{(-(r-ia\cos[\theta]))^{k-i}}{\Sigma^{k}}.$
(116)

This form is useful as the Fourier expansion of the numerators in Eqs. (113),
(115) and (116) are trivial. Hence, only the Fourier expansion of the
denominator, which is real, is required.

A common denominator factor of $\Sigma$ is desirable to minimise the number of
Fourier expansions required. Using Eq. (113), we can express Eq. (112) as

$\displaystyle
f[t,r,\theta,\phi]=\frac{1}{\Sigma^{i}}\sum_{\ell=|s|}^{\infty}\sum_{\mathtt{m}=-\ell}^{\ell}f_{\ell\mathtt{m}}[r]e^{-i\omega
t}{}_{s}Y_{\ell\mathtt{m}}[\theta,\phi_{*}].$ (117)

for some $i$, which has collected all of the factors of $\rho$ and
${\bar{\rho}}$ outside of the sum in the form of $\frac{1}{\Sigma^{i}}$.
Therefore, we require only one Fourier expansion, that of
$\frac{1}{\Sigma^{i}}$. For example, take the function

$\displaystyle
f=A^{\circ}\rho+B^{\circ}\rho^{3}{\bar{\rho}}^{2}+C^{\circ}{\bar{\rho}}^{4}.$
(118)

To extract the same $\rho{\bar{\rho}}$ dependency, we note that the highest
power of either $\rho$ or ${\bar{\rho}}$ is ${\bar{\rho}}^{4}$. Hence, we
write this function as,

$\displaystyle f$
$\displaystyle=\rho^{4}{\bar{\rho}}^{4}(\frac{A^{\circ}}{\rho^{3}{\bar{\rho}}^{4}}+\frac{B^{\circ}}{\rho{\bar{\rho}}^{2}}+\frac{C^{\circ}}{\rho^{4}}).$
(119)

Using Eq. (25) to replace the inverse powers of $\rho$ and ${\bar{\rho}}$
gives

$\displaystyle
f=\frac{1}{\Sigma^{4}}(-(r-ia\cos[\theta])^{3}(r+ia\cos[\theta])^{4}A^{\circ}$
$\displaystyle-(r-ia\cos[\theta])(r+ia\cos[\theta])^{2}B^{\circ}+(r-ia\cos[\theta])^{4}C^{\circ}).$
(120)

We can use $\cos[\theta]=\sqrt{\frac{4\pi}{3}}\ {}_{0}Y_{10}$ and Eqs. (34)
and (II.1) to combine the angular dependence with that in $A^{\circ}$,
$B^{\circ}$, and $C^{\circ}$, allowing us to express $f$ in the form of
Eq.(117).

####  II.6.1 Fourier series

Figure 1:  Top: the function $\Sigma^{-6}$ is plotted in blue (for $a=0.9$ and
$r=r_{+}\approx 1.446M$), and the Fourier expansion (for $k=8$) is over-
plotted as a dashed yellow line. Bottom: the fractional disagreement between
the two plots is shown. There is less than $0.5\%$ error, and the error
decreases for increasing radius (see Fig. 2). Figure 2:  The fractional error,
$|\frac{\Sigma^{-6}}{F[\Sigma^{-6},k=12]}-1|$, for different digits of
precision used (green line is 8, red line is 32, and blue line is 128 digits
of precision), for $a=0.9$ and $\theta=\frac{\pi}{4}$. The plots show a clear
dependence of the error on the number of digits of precision used for large
$r$. This problem does not occur for $k<2i$ when approximating $\Sigma^{-i}$.

A Fourier series expresses periodic functions using the orthogonality of the
trigonometric functions $\sin[n\theta]$ and $\cos[n\theta]$. Truncated to
order $k$, over an interval of $2\pi$, the Fourier series of a function
$A[\theta]$ can be expressed as

$\displaystyle
F[A[\theta],k]=\frac{a_{0}}{2}+\sum_{n=1}^{k}a_{n}\cos[n\theta]+\sum_{n=1}^{k}b_{n}\sin[n\theta],$
(121)

where

$\displaystyle a_{0}$
$\displaystyle=\frac{1}{\pi}\int_{-\pi}^{\pi}A[\theta]d\theta,$ (122)
$\displaystyle a_{n}$
$\displaystyle=\frac{1}{\pi}\int_{-\pi}^{\pi}A[\theta]\cos[n\theta]d\theta,$
(123) $\displaystyle b_{n}$
$\displaystyle=\frac{1}{\pi}\int_{-\pi}^{\pi}A[\theta]\sin[n\theta]d\theta.$
(124)

$\frac{1}{\Sigma^{i}}$ is an even periodic function over $[-\pi,\pi]$; hence,
$b_{n}=0$. Also, as $\Sigma^{i}$ is a function of only even powers of
$\cos[\theta]$, $a_{2n+1}=0$ for any positive integer $n$. Note that
$\frac{1}{\Sigma^{i}}$ is not singular in the domain $r\in[r_{+},\infty)$ so
we expect our Fourier series to converge.

We now assess the accuracy of a truncated Fourier series. For example, take
the expansion of $\Sigma^{-6}$ to order $k=8$,

$\displaystyle F[$
$\displaystyle\Sigma^{-6},8]=\frac{1}{256r^{11}(a^{2}+r^{2})^{11/2}}\Big{(}(a^{2}+2r^{2})(63a^{8}$
$\displaystyle+224a^{6}r^{2}+352a^{4}r^{4}+256a^{2}r^{6}+128r^{8})-6a^{2}(21a^{8}$
$\displaystyle+112a^{6}r^{2}+240a^{4}r^{4}+256a^{2}r^{6}+128r^{8})\cos[2\theta]$
$\displaystyle+42a^{4}(a^{2}+2r^{2})(3a^{4}+8a^{2}r^{2}+8r^{4})\cos[4\theta]$
$\displaystyle-14a^{6}(9a^{4}+32a^{2}r^{2}+32r^{4})\cos[6\theta]$
$\displaystyle+126a^{8}(a^{2}+2r^{2})\cos[8\theta]\Big{)},$ (125)

calculated using Mathematica [71]. Eq. (II.6.1) accurately approximates
$\Sigma^{-6}$, as can be seen in Fig. 1. For $a=0.9$, $r=1.43M$, and $k=8$,
the error is $<0.5\%$. For larger $r$, the approximation becomes more accurate
(similarly to the blue line in Fig. 2). The most inaccurate region is for
small $r$ and near the poles. For $a=0.9$, the outer black hole horizon is
$r_{+}\approx 1.43M$ (Eq. 86). For higher values of $a$, the Fourier series is
less accurate near the horizon. Considering the accuracy requirements of LISA
and current estimates on the spins of supermassive black holes, the expansion
used in Fig. 1 is probably sufficient.

To achieve arbitrarily high accuracy, one can increase the number of terms in
the Fourier expansion, $k$. However, when attempting to approximate
$\Sigma^{n}$ with a Fourier series with $k\geq 2n$, errors can occur for large
$r$ if sufficiently high precision is not used, as shown in Fig. 2. We examine
the simplest example, $\Sigma^{-1}$, to understand why large errors are
encountered when using low precision. The Fourier series of $\Sigma^{-1}$ with
$k=2$ is

$\displaystyle
F[\Sigma^{-1},2]=\frac{a^{2}-2\big{(}a^{2}+2r(r-\sqrt{a^{2}+r^{2}})\big{)}\cos[2\theta]}{ra^{2}\sqrt{a^{2}+r^{2}}}$
(126)

Naively, the Fourier series appears to scale as $\sim r^{0}$ for large $r$,
whereas $\frac{1}{\Sigma}\sim r^{-2}$. However,
$r-\sqrt{a^{2}+r^{2}}\rightarrow 0$ for large $r$, so this contribution is
suppressed, so the Fourier series scales as $r^{-2}$, as expected. However,
resolving the suppression of such terms numerically requires high precision.
With insufficient precision, an error $\propto r^{2}$ is incurred.

In practice, using such high precision will slow calculations down.
Nonetheless, for EMRI models for LISA data analysis, the accuracy requirements
are likely sufficiently modest that Fourier expansions with $k<2i$ can be
used; that is, high precision evaluation will likely not be required.

###  II.7 Combining the spin-weighted spherical harmonics

Finally, we can replace $\rho^{i}{\bar{\rho}}^{i}=\Sigma^{-i}$ factor in Eq.
(117) with the Fourier expansion truncated to some finite order. To convert
the trigonometric quantities in the Fourier series with spin-weighted
spherical harmonics, we can use Eqs. 110 and 111. By again using Eqs. (34) and
(II.1), we can combine the spin-weighted spherical harmonics to express the
source as

$\displaystyle
f[t,r,\theta,\phi]=\sum_{\ell=|s|}^{\infty}\sum_{\mathtt{m}=-\ell}^{\ell}\hat{f}_{\ell\mathtt{m}}[r]e^{-i\omega
t}{}_{s}Y_{\ell\mathtt{m}}[\theta,\phi_{*}].$ (127)

Using the inversion of the spin-weighted spherical harmonic expansion of the
spin-weighted spheroidal harmonics [51, 52], one can re-expand Eq. 127 in
terms of spin-weighted spheroidal harmonics,

$\displaystyle
f[t,r,\theta,\phi]=\sum_{\ell=|s|}^{\infty}\sum_{\mathtt{m}=-\ell}^{\ell}f_{\ell\mathtt{m}}[r]e^{-i\omega
t}{}_{s}S_{\ell\mathtt{m}}[\theta,\phi_{*}].$ (128)

This completes our formalism as we have separated the source into the spin-
weighted spheroidal harmonic modes of the Teukolsky equation. The radial
Teukolsky equation (22) can now be solved with $f_{\ell\mathtt{m}}[r]$ in
(128). Our formalism is summarised in Fig. 3.

Take the source ($f$) of the master Teukolsky equation (9); E.g.
$f=-\Sigma\mathcal{S}_{0}[X[A_{1},A_{2},...,A_{n}]]$,for a general operator
$X$ acting on input tensors $A_{1},A_{2},...,A_{n}$.Express the inputs
$A_{1},A_{2},...,A_{n}$ in the form of Eq. 107,a spin-weighted spherical
harmonic mode decomposition, up to factors of $\rho$ and
${\bar{\rho}}$.Express $f$ in GHP form using sections II.2 and II.3.Apply the
commutation relations, Eqs. 54, 55 and 56 (and their complex
conjugates).Convert to half-Held form:Replace the spin coefficients (and
$\psi_{2}$) with their Held form using Eqs. 69, 71, 72, 73 and 74;Convert the
GHP derivatives to half-Held form using Eqs. 66 and 67;Simplify the expression
using Eqs. (II.5) to (100) and Eq. (77); Convert to coordinate form:Replace
Held scalars with their coordinate form using Eq. 109.Convert to
$\tilde{\hbox{\ec\char 240\relax}}$ and $\tilde{\hbox{\ec\char
240\relax}}^{\prime}$ to spin raising and lowering operators using Sections
II.4.1 and 96. Apply Eqs. 30 and 31, no spin raising and lowering operators
should remain.Replace Þ and $\hbox{\ec\char 222\relax}^{\prime}$ with their
coordinate form using Eq. 53 followed by Eqs. 70, 104, 105 and 106. Collect
the factors of $\rho$ and ${\bar{\rho}}$. Replace $\rho^{i}{\bar{\rho}}^{j}$
with the appropriate Fourier series; e.g., Section II.6.1 for
$\rho^{6}{\bar{\rho}}^{6}$. Express all the remaining trigonometric quantities
as spin-weighted spherical harmonics using Eqs. 110 and 111. Combine the spin-
weighted spherical harmonics using Eqs. 34 and II.1. The resulting expression
will be of the form Eq. 127. Re-expand the source in terms of spin-weighted
spheroidal harmonics using Refs. [51, 52].The source is now separated into
modes of the Teukolsky equation, Eq. 128. One can now solve the radial
Teukolsky equation, Eq. 22.  Figure 3:  Summary of our method for separating
the source of the Teukolsky equation.

##  III Example: First-order Teukolsky equation source with an extended source

In this section, we present a toy-model application of our method to show that
it provides a high-accuracy expression of the source of the Teukolsky
equation. The accompanying Mathematica notebook [72] provides the
implementation of this example for reference. This notebook was built upon
other publicly available notebooks [73, 74].

While the initial motivation for our method is separating the source of the
second-order Teukolsky equation, it also applies to the simpler case of the
first-order Teukolsky equation. The source of the spin $-2$ first-order
Teukolsky equation is

$\displaystyle f[t$ $\displaystyle,r,\theta,\phi]={\mathcal{S}}_{4}[T_{ab}]=$
$\displaystyle\frac{1}{2}(\hbox{\ec\char
240\relax}^{\prime}-\bar{\tau}-4\tau^{\prime})\Big{(}(\hbox{\ec\char
222\relax}^{\prime}-2\bar{\rho}^{\prime})T_{n\bar{m}}-(\hbox{\ec\char
240\relax}^{\prime}-\bar{\tau})T_{nn}\Big{)}$
$\displaystyle+\frac{1}{2}(\hbox{\ec\char
222\relax}^{\prime}-\bar{\rho}^{\prime}-4\rho^{\prime})\Big{(}(\hbox{\ec\char
240\relax}^{\prime}-2\bar{\tau})T_{n\bar{m}}-(\hbox{\ec\char
222\relax}^{\prime}-\bar{\rho}^{\prime})T_{nn}\Big{)}.$ (129)

For our example, we take Eq. 129 as our source and will express it in the form
of Eq. 127. Following Fig. 3, we will expand the stress-energy tensor in terms
of spin-weighted spherical harmonics, similarly to Eq. 107. We also express
Eq. 129 in our Half-Held formalism as outlined in Sec. II.5; please see the
accompanying Mathematica notebook [72] for the step-by-step application and
the resulting formulas.

We require a stress-energy tensor with an extended radial profile. The stress-
energy tensor of a perfect fluid is

$\displaystyle T^{ab}=(\varrho+p)u^{a}u^{b}+pg^{ab}_{(0)},$ (130)

where $\varrho$ is the density, $p$ is the pressure, $u^{a}$ is the four-
velocity. To produce a simple example for applying our method, we assume $p=0$
and

$\displaystyle\varrho=Pe^{-\frac{\Lambda}{r}},$ (131)

where $P$ and $\Lambda$ are constants. Similarly, we assume the fluid is
stationary in the rest-frame of the black hole,

$\displaystyle u^{a}=\Big{\\{}1,0,0,0\Big{\\}}.$ (132)

Using Eqs. 130, 132, 58, 60 and 1,

$\displaystyle T_{nn}$
$\displaystyle=\frac{\varrho\Delta^{2}}{4\Sigma^{2}}=\frac{\varrho\Delta^{2}}{4}\rho{\bar{\rho}},$
(133) $\displaystyle T_{n\bar{m}}$ $\displaystyle=\frac{-\varrho\Delta
ia\sin[\theta]}{2\sqrt{2}\Sigma\zeta}=\frac{\varrho\Delta
ia\sin[\theta]}{2\sqrt{2}}\rho^{2}{\bar{\rho}},$ (134) $\displaystyle
T_{\bar{m}\bar{m}}$ $\displaystyle=\frac{-\varrho
a^{2}\sin^{2}[\theta]}{2\zeta^{2}}=\frac{-\varrho
a^{2}\sin^{2}[\theta]}{2}\rho^{2},$ (135)

where we have extracted the power of $\rho$ and ${\bar{\rho}}$. Using Eqs. 110
and 111, we can decompose Eqs. 133, 134 and 135 into spin-weighted spherical
harmonics,

$\displaystyle T_{nn}$
$\displaystyle=\rho{\bar{\rho}}\sqrt{\pi}\frac{\varrho\Delta^{2}}{2}{}_{0}Y_{00},$
(136) $\displaystyle T_{n\bar{m}}$
$\displaystyle=\rho^{2}{\bar{\rho}}\Big{(}-i\sqrt{\frac{\pi}{3}}\varrho
a\Delta\Big{)}{}_{-1}Y_{10},$ (137) $\displaystyle T_{\bar{m}\bar{m}}$
$\displaystyle=\rho^{2}\Big{(}-2\sqrt{\frac{2\pi}{15}}\frac{\varrho
a^{2}}{2}\Big{)}{}_{-2}Y_{20}.$ (138)

The simplicity of our stress-energy tensor results in the mode decomposition
of its components being that of a single mode for each component.

Inserting Eqs. 136, 137 and 138 into Eq. 129 allows us to extract all the
$\rho$ and ${\bar{\rho}}$ factors, expressing the source in the form of Eq.
112. The highest factor is $\rho^{6}{\bar{\rho}}^{4}$. We can pull out a
common factor of $\rho^{6}{\bar{\rho}}^{4}$ using Eq. 25 and combine the spin-
weighted spherical harmonics using Eqs. 34 and II.1 to express

$\displaystyle
f[t,r,\theta,\phi]=\rho^{6}{\bar{\rho}}^{4}\sum_{\ell=|s|}^{\infty}\sum_{\mathtt{m}=-\ell}^{\ell}f_{\ell\mathtt{m}}[r]{}_{-2}Y_{\ell\mathtt{m}}[\theta,\phi_{*}].$
(139)

Note there is no $e^{-i\omega t}$ dependence, and the spin-weighted spheroidal
harmonics are equivalent to spin-weighted spherical harmonics because the
source is stationary ($\omega=0$).

To complete our decomposition method requires expressing
$\rho^{6}{\bar{\rho}}^{4}$ in a separable form and combining it into the mode
summation. To do so, we can use the Fourier expansion (with $k=8$) for
$\Sigma^{-6}$ in Section II.6.1, as
$\rho^{6}{\bar{\rho}}^{4}=\Sigma^{-6}(-r-ia\cos[\theta])^{2}$. The resulting
expression, after again applying Eqs. 34 and II.1, is of the form

$\displaystyle
f[t,r,\theta,\phi]=\sum_{\ell=|s|}^{\infty}\sum_{\mathtt{m}=-\ell}^{\ell}\hat{f}_{\ell\mathtt{m}}[r]{}_{-2}Y_{\ell\mathtt{m}}[\theta,\phi_{*}].$
(140)

The explicit form of $\hat{f}_{\ell\mathtt{m}}[r]$ can be found in the
accompanying Mathematica notebook.

In order to validate our expression for Eq. 139, we compare it to the four-
dimensional expression for the source, Eq. 129. It is straightforward to
calculate

$\displaystyle f[t,r,\theta,\phi]=$
$\displaystyle\frac{Pa^{2}e^{-\frac{r}{\Lambda}}\pi(2\Lambda+ia\cos{\theta}+r)(a^{2}+r(-2M+r))^{2}\sin{\theta}^{2}}{2\Lambda^{2}(a\cos{\theta}+ir)^{4}(ia\cos{\theta}+r)^{3}},$
(141)

due to the simplicity of the stress-energy tensor using Eqs. 53, 57, 58, 59
and 60 (see Ref. [37] for expressions of the spin-coefficients in the
Kinnersley tetrad). We compare our four-dimensional source to the sum in Eq.
139 in Fig. 4 near the horizon ($r=r_{+}\approx 1.446M$). In the top plot, the
two expressions overlap each other. The bottom plot shows that the error is
less than $0.5\%$ and takes the same form as the error in Fig. 1, as expected.

Near the horizon is the most inaccurate zone for our decomposition. This can
be seen in Fig. 5, which plots the radial profile of the error for $a=0.9$ and
$\theta=0.01$. Additional precision errors can be encountered near the horizon
if insufficient precision is used (less than 32 digits of precision). This is
due to the Kinnersley tetrad being singular at the horizon. Hence, when
expressing regular quantities using the Kinnersley tetrad near the horizon,
cancellations of large numbers are required. This problem could be avoided by
using a regular tetrad near the horizon, such as the Hartle–Hawking tetrad
[65, 66].

Figure 4:  Top: the source of the Teukolsky equation for a stationary ideal
gas cloud is plotted in blue (calculated using our decomposition method) and
red (calculated explicitly) for $a=0.9$, $\Lambda=1$, $P=1$, and
$r=r_{+}\approx 1.446M$. The lines overlay each other, showing our
decomposition method is accurate. Bottom: the fractional disagreement between
the two plots is shown; note the error is the same as Fig. 1, as expected.
Figure 5:  The fractional error in the source of the Teukolsky equation for a
stationary ideal gas cloud when calculated using our decomposition method. For
$a=0.9$, $\Lambda=1$, $P=1$, and $\theta=0.01$. The error is less than $0.5\%$
and decreases as one moves away from the horizon as expected.

##  IV Conclusion

In this paper, we produced a formalism for decomposing the source of the
Teukolsky equation into spin-weighted spheroidal harmonics analytically. Our
formalism leverages the Held-formalism’s [55] ability to extract the $\rho$
and ${\bar{\rho}}$ dependency in background Kerr quantities. The remaining
Kerr quantities exhibit a straightforward separation of variables in the
Kinnersley tetrad [38]. We then used a truncated Fourier series to separate
the angular and radial dependency in $\rho$ and ${\bar{\rho}}$. An order $k=8$
Fourier series has a less than $0.5\%$ error for approximating $\Sigma^{-6}$.
This error should be sufficiently small for second-order self-force
calculations, as the error rapidly decreases as one moves away from the
horizon, and second-order effects are small compared to first-order effects.
Additionally, the error can be decreased by increasing the order of the
expansion.

In section III, we demonstrated our method with a toy example: separating the
source of the first-order Teukolsky equation for a stress-energy tensor of a
pressure-less ideal gas. The calculation is publicly available in a
Mathematica notebook [72], which explicitly goes through each step of applying
our formalism, as summarised in Fig. 3. The calculations confirm the accuracy
of our formalism, as we show the error in our source is the same as the error
in our expansion for $\rho$ and ${\bar{\rho}}$, as expected.

There is an alternative method for separating the source of the second-order
Teukolsky equation. Using the orthogonality of the spin-weighted spheroidal
harmonics [45], as used in Ref. [27] for a quadratic quasi-normal mode
calculation. However, our method provides the efficiency advantage of being
analytic, avoiding the need to integrate at each radial point on a grid. This
efficiency saving may be crucial for the more involved second-order self-force
calculations in Kerr.

Implementing our method for the first second-order self-force calculations in
Kerr is currently in progress. Our method could also be used for second-order
(quadratic) quasi-normal mode calculations in Kerr [19, 20, 21, 22, 23, 24,
25, 26, 27]. Additionally, the formalism could be applied to solve Teukolsky
equations in theories of gravity beyond general relativity, which generally
have extended sources [75, 76, 77, 78, 79].

This paper focuses on applying our method in the frequency domain, but the
formalism is also applicable in the time domain. In the time domain,
decomposing into spin-weighted spherical harmonics reduces the 2+1 PDE
Teukolsky equation into a coupled set of 1+1 PDEs with only next-to-next-to
coupling [80, 81, 82, 83, 84]. Hence, our formalism could help produce the
first second-order self-force calculations in the time domain.

###### Acknowledgements.

AS acknowledges the partial support of a Royal Society University Research
Fellowship Enhancement Award and partial support from the STFC Consolidated
Grant no. ST/V005596/1. AS would like to thank Adam Pound, Zach Nasipak, and
Sam Dolan for helpful discussions. AS would like to thank Laura Sberna and
Alexander Grant for their comments on earlier manuscripts of this work.
ChatGPT 3.5 was used to help format equations and figures in this paper. This
work makes use of the Black Hole Perturbation Toolkit.

## References

  * Abbott  _et al._ [2016] B. P. Abbott, R. Abbott, T. Abbott, M. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. Adhikari, _et al._ , Observation of gravitational waves from a binary black hole merger, Physical review letters 116, 061102 (2016).
  * Abbott  _et al._ [2021a] R. Abbott, T. Abbott, S. Abraham, F. Acernese, K. Ackley, A. Adams, C. Adams, R. Adhikari, V. Adya, C. Affeldt, _et al._ , Population properties of compact objects from the second LIGO–Virgo gravitational-wave transient catalog, The Astrophysical journal letters 913, L7 (2021a).
  * Abbott  _et al._ [2021b] R. Abbott, T. Abbott, F. Acernese, K. Ackley, C. Adams, N. Adhikari, R. Adhikari, V. Adya, C. Affeldt, D. Agarwal, _et al._ , GWTC-3: compact binary coalescences observed by LIGO and Virgo during the second part of the third observing run, arXiv preprint arXiv:2111.03606 (2021b).
  * et al. [2017] K. D. et al., Lisa: A proposal in response to the esa call for l3 mission concepts, (2017).
  * Luo  _et al._ [2016] J. Luo, L.-S. Chen, H.-Z. Duan, Y.-G. Gong, S. Hu, J. Ji, Q. Liu, J. Mei, V. Milyukov, M. Sazhin, _et al._ , Tianqin: a space-borne gravitational wave detector, Classical and Quantum Gravity 33, 035010 (2016).
  * Ruan  _et al._ [2020] W.-H. Ruan, Z.-K. Guo, R.-G. Cai, and Y.-Z. Zhang, Taiji program: Gravitational-wave sources, International Journal of Modern Physics A 35, 2050075 (2020).
  * 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] .
  * Baumgarte and Shapiro [2010] T. W. Baumgarte and S. L. Shapiro,  _Numerical relativity: solving Einstein’s equations on the computer_ (Cambridge University Press, 2010).
  * Lousto and Healy [2020] C. O. Lousto and J. Healy, Exploring the small mass ratio binary black hole merger via zeno’s dichotomy approach, Physical Review Letters 125, 191102 (2020).
  * Rosato  _et al._ [2021] N. Rosato, J. Healy, and C. O. Lousto, Adapted gauge to small mass ratio binary black hole evolutions, Physical Review D 103, 104068 (2021).
  * Wittek  _et al._ [2023] N. A. Wittek _et al._ , Worldtube excision method for intermediate-mass-ratio inspirals: Scalar-field model in 3+1 dimensions, Phys. Rev. D 108, 024041 (2023), arXiv:2304.05329 [gr-qc] .
  * Dhesi  _et al._ [2021] M. Dhesi, H. R. Rüter, A. Pound, L. Barack, and H. P. Pfeiffer, Worldtube excision method for intermediate-mass-ratio inspirals: Scalar-field toy model, Phys. Rev. D 104, 124002 (2021), arXiv:2109.03531 [gr-qc] .
  * Barack and Pound [2018] L. Barack and A. Pound, Self-force and radiation reaction in general relativity, Reports on Progress in Physics 82, 016904 (2018).
  * Pound  _et al._ [2020] A. Pound, B. Wardell, N. Warburton, and J. Miller, Second-order self-force calculation of gravitational binding energy in compact binaries, Physical Review Letters 124, 021101 (2020).
  * Wardell  _et al._ [2021] B. Wardell, A. Pound, N. Warburton, J. Miller, L. Durkan, and A. L. Tiec, Gravitational waveforms for compact binaries from second-order self-force theory, arXiv preprint arXiv:2112.12265 (2021).
  * Warburton  _et al._ [2021] N. Warburton, A. Pound, B. Wardell, J. Miller, and L. Durkan, Gravitational-wave energy flux for compact binaries through second order in the mass ratio, Physical review letters 127, 151102 (2021).
  * Kevorkian and Cole [1996] J. Kevorkian and J. Cole, The method of multiple scales for ordinary differential equations, in  _Multiple Scale and Singular Perturbation Methods_ (Springer, 1996) pp. 267–409.
  * Flanagan and Hinderer [2008] E. E. Flanagan and T. Hinderer, Two-timescale analysis of extreme mass ratio inspirals in Kerr spacetime Orbital motion, Physical Review D 78, 064028 (2008).
  * Loutrel  _et al._ [2021] N. Loutrel, J. L. Ripley, E. Giorgi, and F. Pretorius, Second-order perturbations of Kerr black holes: Formalism and reconstruction of the first-order metric, Physical Review D 103, 104017 (2021).
  * Ripley  _et al._ [2021] J. L. Ripley, N. Loutrel, E. Giorgi, and F. Pretorius, Numerical computation of second-order vacuum perturbations of Kerr black holes, Physical Review D 103, 104018 (2021).
  * Ioka and Nakano [2007] K. Ioka and H. Nakano, Second-and higher-order quasinormal modes in binary black-hole mergers, Physical Review D 76, 061503 (2007).
  * Nakano and Ioka [2007] H. Nakano and K. Ioka, Second-order quasinormal mode of the Schwarzschild black hole, Physical Review D 76, 084007 (2007).
  * Cheung  _et al._ [2023] M. H.-Y. Cheung, V. Baibhav, E. Berti, V. Cardoso, G. Carullo, R. Cotesta, W. Del Pozzo, F. Duque, T. Helfer, E. Shukla, _et al._ , Nonlinear effects in black hole ringdown, Physical Review Letters 130, 081401 (2023).
  * Mitman  _et al._ [2023] K. Mitman, M. Lagos, L. C. Stein, S. Ma, L. Hui, Y. Chen, N. Deppe, F. Hébert, L. E. Kidder, J. Moxon, _et al._ , Nonlinearities in black hole ringdowns, Physical Review Letters 130, 081402 (2023).
  * Lagos and Hui [2023] M. Lagos and L. Hui, Generation and propagation of nonlinear quasinormal modes of a Schwarzschild black hole, Phys. Rev. D 107, 044040 (2023), arXiv:2208.07379 [gr-qc] .
  * Redondo-Yuste  _et al._ [2023] J. Redondo-Yuste, G. Carullo, J. L. Ripley, E. Berti, and V. Cardoso, Spin dependence of black hole ringdown nonlinearities, (2023), arXiv:2308.14796 [gr-qc] .
  * Ma and Yang [2024] S. Ma and H. Yang, The excitation of quadratic quasinormal modes for Kerr black holes, (2024), arXiv:2401.15516 [gr-qc] .
  * [28] J. Mathews, B. Wardell, A. Pound, N. Warburton, J. Thompson, L. Durkan, and J. Miller, Gravitational waveforms for compact binaries from second-order self-force theory: spinning binaries, in preparation.
  * Mathews  _et al._ [2022] J. Mathews, A. Pound, and B. Wardell, Self-force calculations with a spinning secondary, Physical Review D 105, 084031 (2022).
  * Kerr [1963] R. P. Kerr, Gravitational field of a spinning mass as an example of algebraically special metrics, Physical review letters 11, 237 (1963).
  * Boyer and Lindquist [1967] R. H. Boyer and R. W. Lindquist, Maximal analytic extension of the kerr metric, Journal of mathematical physics 8, 265 (1967).
  * Teukolsky [1972] S. A. Teukolsky, Rotating black holes: Separable wave equations for gravitational and electromagnetic perturbations, Physical Review Letters 29, 1114 (1972).
  * Teukolsky [1973] S. A. Teukolsky, Perturbations of a rotating black hole. 1. fundamental equations for gravitational electromagnetic and neutrino field perturbations, Astrophys. J. 185, 635 (1973).
  * Spiers  _et al._ [2023a] A. Spiers, A. Pound, and J. Moxon, Second-order Teukolsky formalism in Kerr spacetime: Formulation and nonlinear source, Phys. Rev. D 108, 064002 (2023a), arXiv:2305.19332 [gr-qc] .
  * Martel and Poisson [2005] K. Martel and E. Poisson, Gravitational perturbations of the Schwarzschild spacetime: a practical covariant and gauge-invariant formalism, Physical Review D 71, 104003 (2005).
  * Spiers  _et al._ [2023b] A. Spiers, A. Pound, and B. Wardell, Second-order perturbations of the schwarzschild spacetime: practical, covariant and gauge-invariant formalisms, arXiv preprint arXiv:2306.17847 (2023b).
  * Pound and Wardell [2021] A. Pound and B. Wardell, Black hole perturbation theory and gravitational self-force, arXiv preprint arXiv:2101.04592 (2021).
  * Kinnersley [1969] W. Kinnersley, Type D vacuum metrics, Journal of Mathematical Physics 10, 1195 (1969).
  * Aksteiner [2014] S. Aksteiner, Geometry and analysis on black hole spacetimes, (2014).
  * 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 .
  * Berti  _et al._ [2009] E. Berti, V. Cardoso, and A. O. Starinets, Quasinormal modes of black holes and black branes, Class. Quant. Grav. 26, 163001 (2009), arXiv:0905.2975 [gr-qc] .
  * Konoplya and Zhidenko [2011] R. A. Konoplya and A. Zhidenko, Quasinormal modes of black holes: From astrophysics to string theory, Rev. Mod. Phys. 83, 793 (2011), arXiv:1102.4014 [gr-qc] .
  * Campanelli and Lousto [1999] M. Campanelli and C. O. Lousto, Second order gauge invariant gravitational perturbations of a Kerr black hole, Physical Review D 59, 124022 (1999).
  * Green  _et al._ [2020] S. R. Green, S. Hollands, and P. Zimmerman, Teukolsky formalism for nonlinear Kerr perturbations, Classical and Quantum Gravity 37, 075001 (2020).
  * Breuer  _et al._ [1977] R. A. Breuer, M. P. Ryan, and S. Waller, Some properties of spin-weighted spheroidal harmonics, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 358, 71 (1977).
  * London [2023] L. T. London, Biorthogonal harmonics for the decomposition of gravitational radiation. i. angular modes, completeness, and the introduction of adjoint-spheroidal harmonics, Physical Review D 107, 044056 (2023).
  * Harms  _et al._ [2013] E. Harms, S. Bernuzzi, and B. Brügmann, Numerical solution of the 2+ 1 teukolsky equation on a hyperboloidal and horizon penetrating foliation of kerr and application to late-time decays, Classical and Quantum Gravity 30, 115013 (2013).
  * Panosso Macedo [2020] R. Panosso Macedo, Hyperboloidal framework for the Kerr spacetime, Class. Quant. Grav. 37, 065019 (2020), arXiv:1910.13452 [gr-qc] .
  * Leaver [1986] E. W. Leaver, Solutions to a generalized spheroidal wave equation: Teukolsky’s equations in general relativity, and the two-center problem in molecular quantum mechanics, Journal of mathematical physics 27, 1238 (1986).
  * Press and Teukolsky [1973] W. H. Press and S. A. Teukolsky, Perturbations of a rotating black hole. ii. dynamical stability of the kerr metric, Astrophysical Journal, Vol. 185, pp. 649-674 (1973) 185, 649 (1973).
  * Hughes [2000] S. A. Hughes, Evolution of circular, nonequatorial orbits of kerr black holes due to gravitational-wave emission, Physical Review D 61, 084004 (2000).
  * [52] Black Hole Perturbation Toolkit, (bhptoolkit.org).
  * Newman and Penrose [1962] E. Newman and R. Penrose, An approach to gravitational radiation by a method of spin coefficients, Journal of Mathematical Physics 3, 566 (1962).
  * Geroch  _et al._ [1973] R. Geroch, A. Held, and R. Penrose, A space-time calculus based on pairs of null directions, Journal of Mathematical Physics 14, 874 (1973).
  * Held [1974] A. Held, A formalism for the investigation of algebraically special metrics. i, Communications in Mathematical Physics 37, 311 (1974).
  * Goldberg  _et al._ [1967] J. N. Goldberg, A. J. MacFarlane, E. T. Newman, F. Rohrlich, and E. G. Sudarshan, Spin-s spherical harmonics and ð, Journal of Mathematical Physics 8, 2155 (1967).
  * Hecht [2012] K. T. Hecht,  _Quantum mechanics_ (Springer Science & Business Media, 2012).
  * Shiraishi [2013] M. Shiraishi,  _Probing the early universe with the CMB scalar, vector and tensor bispectrum_ (Springer Science & Business Media, 2013).
  * Wald [1984] R. M. Wald,  _General Relativity_ (The University of Chicago Press, 1984).
  * Petrov [2000] A. Z. Petrov, The classification of spaces defining gravitational fields, General Relativity and Gravitation 32, 1665 (2000).
  * Chandrasekhar [1992] S. Chandrasekhar,  _The Mathematical Theory of Black Holes_ (Clarendon Press Oxford, 1992) reprint: 2009.
  * Goldberg and Sachs [2009] J. Goldberg and R. Sachs, Republication of: A theorem on Petrov types, General Relativity and Gravitation 41, 433 (2009).
  * Toomani  _et al._ [2021] V. Toomani, P. Zimmerman, A. Spiers, S. Hollands, A. Pound, and S. R. Green, New metric reconstruction scheme for gravitational self-force calculations, Classical and Quantum Gravity 39, 015019 (2021).
  * Price  _et al._ [2007] L. R. Price, K. Shankar, and B. F. Whiting, On the existence of radiation gauges in Petrov type II spacetimes, Classical and Quantum Gravity 24, 2367 (2007).
  * Hawking and Hartle [1972] S. W. Hawking and J. Hartle, Energy and angular momentum flow into a black hole, Communications in mathematical physics 27, 283 (1972).
  * Poisson [2004] E. Poisson, Absorption of mass and angular momentum by a black hole: time-domain formalisms for gravitational perturbations, and the small-hole or slow-motion approximation, Physical Review D 70, 084044 (2004).
  * Cooperstock and Richardson [1992] F. Cooperstock and S. Richardson, Energy localization and the kerr-newman metric, in  _General Relativity And Relativistic Astrophysics-Proceedings Of The 4th Canadian Conference_ (World Scientific, 1992) p. 110.
  * Estes and Lancaster [1966] R. Estes and E. Lancaster,  _Two-point Taylor series expansions_ , Tech. Rep. (NASA, 1966).
  * Hablützel [2020] R. Hablützel, The easiest way towards multi-point taylor expansion, (2020), https://underthemath.wordpress.com/2020/06/12/polynomial-division-revisited/.
  * O’Boyle  _et al._ [2022] M. F. O’Boyle, C. Markakis, L. J. G. Da Silva, R. P. Macedo, and J. A. V. Kroon, Conservative evolution of black hole perturbations with time-symmetric numerical methods, arXiv preprint arXiv:2210.02550 (2022).
  * [71] W. R. Inc., Mathematica, Version 13.1, champaign, IL, 2022\.
  * Spiers [2024] A. Spiers, Teukolsky equation source decomposer, https://github.com/DrAndrewSpiers/Teukolsky-equation-source-decomposer (2024).
  * Gómez-Lobo [2009] A. G.-P. Gómez-Lobo, The n.p. and g.h.p. formalisms., http://www.xact.es/Documentation/English/PublicNPGHP.nb (2009).
  * Spiers [2023] A. Spiers, NP and GHP Formalisms for second-order Teukolsky equations, https://github.com/DrAndrewSpiers/NP-and-GHP-Formalisms-for-2nd-order-Teukolsky (2023).
  * Hussain and Zimmerman [2022] A. Hussain and A. Zimmerman, Approach to computing spectral shifts for black holes beyond kerr, Physical Review D 106, 104018 (2022).
  * Li  _et al._ [2023] D. Li, P. Wagle, Y. Chen, and N. Yunes, Perturbations of spinning black holes beyond general relativity: Modified teukolsky equation, Physical Review X 13, 021029 (2023).
  * Spiers  _et al._ [2023c] A. Spiers, A. Maselli, and T. P. Sotiriou, Measuring scalar charge with compact binaries: High accuracy modelling with self-force, (2023c), arXiv:2310.02315 [gr-qc] .
  * Wagle  _et al._ [2023] P. Wagle, D. Li, Y. Chen, and N. Yunes, Perturbations of spinning black holes in dynamical chern-simons gravity i. slow rotation equations, arXiv preprint arXiv:2311.07706 (2023).
  * D’Addario  _et al._ [2023] G. D’Addario, A. Padilla, P. M. Saffin, T. P. Sotiriou, and A. Spiers, Ringdowns for black holes with scalar hair: the large mass case, arXiv preprint arXiv:2311.17666 (2023).
  * Barack [1999] L. Barack, Late time decay of scalar, electromagnetic, and gravitational perturbations outside rotating black holes, Physical Review D 61, 024026 (1999).
  * Barack and Giudice [2017] L. Barack and P. Giudice, Time-domain metric reconstruction for self-force applications, Physical Review D 95, 104033 (2017).
  * Dolan [2013] S. R. Dolan, Superradiant instabilities of rotating black holes in the time domain, Physical Review D 87, 124026 (2013).
  * Long and Barack [2021] O. Long and L. Barack, Time-domain metric reconstruction for hyperbolic scattering, Physical Review D 104, 024014 (2021).
  * Markakis  _et al._ [2023] C. Markakis, S. Bray, and A. Zenginoğlu, Symmetric integration of the 1+ 1 teukolsky equation on hyperboloidal foliations of kerr spacetimes, arXiv preprint arXiv:2303.08153 (2023).