# Natural polynomials for Kerr quasi-normal modes

Lionel London  King’s College London, Strand, London WC2R 2LS, United Kingdom
Michelle Gurevich  King’s College London, Strand, London WC2R 2LS, United
Kingdom

###### Abstract

We present a polynomial basis that exactly tridiagonalizes Teukolsky’s radial
equation for quasi-normal modes. These polynomials naturally emerge from the
radial problem, and they are “canonical” in that they possess key features of
classical polynomials. Our canonical polynomials may be constructed using
various methods, the simplest of which is the Gram-Schmidt process. In
contrast with other polynomial bases, our polynomials allow for Teukolsky’s
radial equation to be represented as a simple matrix eigenvalue equation that
has well-behaved asymptotics and is free of non-physical solutions. We expect
that our polynomials will be useful for better understanding the Kerr
quasinormal modes’ properties, particularly their prospective spatial
completeness and orthogonality. We show that our polynomials are closely
related to the confluent Heun and Pollaczek-Jacobi type polynomials.
Consequently, our construction of polynomials may be used to tridiagonalize
other instances of the confluent Heun equation. We apply our polynomials to a
series of simple examples, including: (1) the high accuracy numerical
computation of radial eigenvalues, (2) the evaluation and validation of
quasinormal mode solutions to Teukolsky’s radial equation, and (3) the use of
Schwarzschild radial functions to represent those of Kerr. Along the way, a
potentially new concept, “confluent Heun polynomial/non-polynomial duality”,
is encountered and applied to show that some quasinormal mode separation
constants are well approximated by confluent Heun polynomial eigenvalues. We
briefly discuss the implications of our results on various topics, including
the prospective spatial completeness of Kerr quasinormal modes.

##  I Introduction

The results of single black hole (BH) perturbation theory have proved widely
useful to gravitational wave astronomy Vishveshwara (1970); Abbott  _et al._
(2016a, b, 2021a, 2021b, 2019); Hughes  _et al._ (2019); Khan  _et al._
(2016); Husa  _et al._ (2016); Pompili  _et al._ (2023). Most famously, the
aftermath of binary black hole (BBH) coalescence is expected to be a perturbed
Kerr BH that rings down, radiating away its perturbation, predominantly in the
form of damped oscillatory modes Press (1971); Press and Teukolsky (1973);
Teukolsky (1973); Teukolsky and Press (1974); Leaver (1985); Berti  _et al._
(2006a); London  _et al._ (2014); Hamilton  _et al._ (2023). The related
ringdown radiation is of practical relevance to many aspects of gravitational
wave astronomy, from the construction of signal models, to the ongoing search
for new physics London  _et al._ (2014); Hamilton  _et al._ (2023); Hughes
_et al._ (2019); Khan  _et al._ (2016); Husa  _et al._ (2016); Pompili  _et
al._ (2023); Carullo  _et al._ (2018); Baibhav  _et al._ (2023); Carullo  _et
al._ (2019); Cotesta  _et al._ (2022); Mishra  _et al._ (2010); Bhagwat  _et
al._ (2020); Abbott  _et al._ (2016b). Concurrently, ringdown radiation has
inspired varied activity in gravitational wave theory, from efforts to
understand BH nonlinearity, to emerging developments in time-dependent
perturbation theory Lagos and Hui (2022); Green  _et al._ (2023); Cannizzaro
_et al._ (2023); Sberna  _et al._ (2022); Redondo-Yuste  _et al._ (2023).
These efforts, among many others, aim to prepare gravitational wave science
for a future in which high-precision ringdown measurements may point the way
to a more complete understanding of gravity Arun  _et al._ (2022); Maselli
_et al._ (2017); Berti  _et al._ (2016); Barausse  _et al._ (2015); Klein  _et
al._ (2016); Bailes  _et al._ (2021); Vallisneri  _et al._ (2015); Karnesis
_et al._ (2022); Arun  _et al._ (2022). There is, therefore, good reason to
think that a deeper practical understanding of linearly perturbed Kerr BHs is
of use.

The present work is concerned with answering a number of questions regarding
the mathematical structure and practical use of gravitational wave ringdown.
Ringdown is comprised of Quasi-Normal Modes (QNMs), which are the analog of
normal modes for dissipative systems Kokkotas and Schmidt (1999); Nollert
(1999); Berti  _et al._ (2009a). Here, they will be treated under the
Teukolsky formalism, where the field equations are separable, and where the
QNMs are the joint point spectra of Teukolsky’s radial and angular equations
Press and Teukolsky (1973); Teukolsky (1973); Teukolsky and Press (1974);
Leaver (1985); Berti  _et al._ (2006a); Cook and Zalutskiy (2014). In this
context, the present work addresses the following questions about Kerr QNMs.

### What is the nature of the overtone label?

– Not long after the pioneering work of Teukolsky, Press, and many others on
linear perturbations of the Kerr geometry, Leaver lead the development of
analytic methods for computing the QNM’s functional dependence (i.e. in space
and time) and mode numbers Leaver (1985, 1986, 1986). Since, there has been
tremendous progress in understanding the mathematical structure of Kerr QNMs
Berti  _et al._ (2006a, 2003, b); Berti and Cardoso (2006); Yang  _et al._
(2013, 2012); Zimmerman  _et al._ (2015); Zimmerman and Mark (2016); Maggiore
(2008); Andersson (1995, 1997). For example, analysis of the gravitational
field equations confers that QNMs are described by four mode numbers, one for
each spacetime dimension Berti  _et al._ (2006a); Leaver (1985). In any
asymptotically valid spacetime coordinates, these are $\tilde{\omega}_{{\ell
mn}}$ (time), $\ell$ (polar angle), $m$ (azimuthal angle), and $n$ (radius).
The polar mode number, $\ell$, is unavoidably linked to the eigenvalues of
Jacobi polynomials London (2023a); Fackerell and Crossman (2008); London
(2023b). In this sense, the Jacobi polynomials are natural to the QNM problem.
Regarding the azimuthal angle $\phi$, this is similarly true of the monomials,
$\zeta^{m}$, where $\zeta=e^{i\phi}$. In contrast, the radial mode number,
most commonly referred to as the “overtone” index111The use of “overtone” is
known to be somewhat of a misnomer; e.g. acoustic overtones have higher
frequencies compared to a fundamental resonant frequency. In contrast, QNM
overtones are known to have slightly lower frequencies (but higher damping
rates) compared to the $n=0$, or fundamental, QNMs., has evaded rudimentary
explanation.

Here, the overtone index is shown to originate from an order parameter
inherent to the confluent Heun polynomials Ronveaux and Arscott (1995); DLMF ;
Fiziev (2009a); Hortacsu (2012); Dariescu  _et al._ (2021); Magnus  _et al._
(2021). Concurrently, it is observed perhaps for the first time that the
confluent Heun polynomials exhibit a kind of “polynomial/non-polynomial
duality”, where finite sets of confluent Heun polynomials are necessarily
paired with an infinite set of non-polynomial functions.

### Is a simple spectral solution to Teukolsky’s radial equation possible?

– It was pointed out in Ref. London (2023a) (hereafter Paper I) that there
exists a scalar product for which Teukolsky’s radial equation is self-adjoint,
and that this scalar product enables the direct application of Sturm-Liouville
theory Pinchover and Rubinstein (2005); Abramowitz and Stegun (1964); Arfken
_et al._ (2013); Courant and Hilbert. (1954). The scalar product’s weight
function coincides with the Green’s function approach to QNMs pioneered by
Leaver, and the deformed Jacobi weight known for Pollaczek–Jacobi
polynomialsLondon (2023a); Leaver (1986); Chen and Dai (2010). One use of the
scalar product is to represent Teukolsky’s radial equation as a simple linear
matrix equation Axler (2015).

While many studies have developed matrix representations of Teukolsky’s radial
problem, the lack of an appropriate scalar product has caused complications,
such as the need to numerically integrate the equations, or appearance of non-
physical modes Chung  _et al._ (2023a); Blázquez-Salcedo  _et al._ (2023);
Chung  _et al._ (2023b); Ghojogh  _et al._ (2023); Zhu  _et al._ (2023). In
contrast, the scalar product presented in Paper I was hypothesized to result
in matrix representations that are explicitly symmetric. Further, it was
postulated that there may exist a basis of classical-like polynomials that
enable the representation of Teukolsky’s radial equation as a symmetric and
tridiagonal matrix. Such matrices are known to have number of favorable
properties; e.g., their eigensolutions contain no non-physical modes, and they
are well developed methods to study their asymptotic structure Teschl (2000);
Kuijlaars  _et al._ (2003).

Here, it is shown there do indeed exist unique classical-like polynomials that
exactly tridiagonalize Teukolsky’s radial equation. Therefore, the resulting
method of solution to the QNM’s radial problem is spectral in the simplest
possible way: widely available linear algebra packages may be used to
explicitly diagonalize the problem, simultaneously yielding eigenvalues and
eigenfunctions. The polynomials in question will be referred to as “canonical
confluent Heun polynomials”; numerical examples are provided.

### Is deeper analytic understanding possible?

– Leaver’s algorithm for computing the QNMs’ functions and mode numbers has
influenced many studies, particularly those of the QNM frequencies,
$\tilde{\omega}_{{\ell mn}}$ Cook and Zalutskiy (2014); Magaña Zertuche  _et
al._ (2021); Berti  _et al._ (2009b). Concurrently, methods such as the WKB
and eikonal approximations have supported advances in the analytic
understanding of QNM’s behavior in certain limits, e.g. very high BH spins,
and/or very large values of $\ell$ Zimmerman  _et al._ (2015, 2015); Yang  _et
al._ (2012). Recently, there has been progress in understanding the pseudo-
spectral stability of the QNM frequency spectrum Jaramillo  _et al._ (2021);
Boyanov  _et al._ (2023); Destounis and Duque (2023); Courty  _et al._ (2023);
Torres (2023). Amid this activity, a number of related questions have remained
open. These include: Do the QNMs’ radial dependence have an underlying
structure that is polynomial? Is it heuristically accurate to think of the
radial mode number, $n$, as being the order of an underlying polynomial
object, in the same way that $\ell$ are related to the order of a Jacobi
polynomial?

Here, the QNMs’ radial functions are shown to be generally non-polynomial in
nature, but potentially well approximated by a confluent Heun polynomial of
order given by $n+1$, or in at least one case ($n=8$), of an order much
greater. A simple procedure for quantifying the approximate polynomial nature
of a QNM’s radial function is presented.

### Are the QNMs spatially complete?

– There has been a significant amount of research into the question of QNM
completeness222 As most commonly defined, a function sequence is complete
(i.e. is a basis) if any piecewise continuous function on the same domain can
be exactly expressed as a linear combination of sequence elements. It is a
standard result in Sturm-Liouville theory and functional analysis that self-
adjoint differential operators have complete eigenfunctions.. Some references
appear at tension regarding whether the radial field equation is self-adjoint,
and therefore whether or not QNMs radial functions may be complete. This
apparent tension arises from the fact it is insufficient to a label a
differential equation as non-self-adjoint, without additional statements about
the coordinate domain, boundary conditions, and related bilinear form (e.g.
inner or scalar product). For example, Refs. Jaramillo  _et al._ (2021); Leung
_et al._ (1994); Nollert (1999) describe analogs of the radial field equation
as non-self-adjoint; however, akin to Refs. Wald (2008); Beyer (1995); Green
_et al._ (2023), a key result of Paper I was the construction of a radial
scalar product such that the equation is self-adjoint.

Here, this scalar product is used to show that each QNM’s radial function may
be a member of an orthonormal basis. This new result places our understanding
of the QNM’s radial functions on the same footing as their angular functions
Cook and Zalutskiy (2014); O’Sullivan and Hughes (2014); London (2023b); Berti
and Klein (2014). In particular, it was shown in Ref. London (2023b) that each
QNM’s angular functions is a member of a basis, and that the properties of all
such angular bases lead to the conclusion that the angular functions, namely
the spheroidal harmonics with different QNM frequencies, are isomorphically
equivalent to the spherical harmonics, and therefore complete. The present
work takes a first step towards an analogous conclusion for the QNM’s radial
functions: it is shown that change-of-basis matrices between Kerr and
Schwarzschild are approximately diagonal, implying that they may be
isomorphically equivalent. This communicates that, for a large range of BH
spins, the radial functions of Kerr may be efficiently represented using those
of Schwarzschild.

This result and its consequences may be particularly relevant for the
formulation of a time-dependent perturbation theory for BHs, which presently
requires the rather strong assumption that QNMs are spatially complete Green
_et al._ (2023); Cannizzaro  _et al._ (2023); Redondo-Yuste  _et al._ (2023);
Sberna  _et al._ (2022).

##  II Scope and organization

### Scope

– The focus of this article will be on linear order curvature perturbations of
isolated Kerr BHs, as framed by Teukolsky, Press and Leaver Teukolsky (1973);
Press and Teukolsky (1973); Teukolsky and Press (1974); Leaver (1985, 1986).
While this article discusses only gravitational perturbations (i.e. those with
spin-weight $s=\pm 2$), all calculations present make no explicit assumption
about spin-weight, and in most cases have been verified to explicitly work
with other physically relevant spin weights, namely $|s|\leq 2$. The authors
do not aspire to the language of formal mathematical rigor; in place of
theorems and lemmas, arguments are presented as needed, and the most technical
of these are left to references. Numerical results will use QNM frequencies
due to their physical relevance; however, this article’s results use a
framework that applies to any single frequency. The concurrent treatment of
all QNM frequencies is left to future work. The unit convention, $G=c=1$, will
be used throughout.

### Organization

– This article is organized as follows. In Sec. III, a schematic background of
single BH perturbations is provided, culminating in the definition of the
QNM’s radial eigenvalue problem, and a brief review of Paper I. In Sec. IV,
the QNM’s radial eigenvalue problem is treated with confluent Heun
polynomials. This section concludes with a discussion of the confluent Heun
polynomials pros and cons. In Sec. V, this article’s central result, namely
the canonical confluent Heun polynomials, are introduced. In Sec. VI, the
canonical confluent Heun polynomials are shown to tridiagonalize Teukolsky’s
radial equation. In Sec. VII, numerical results are presented to demonstrate
the canonical confluent Heun polynomials’ end-to-end use. In Sec. VIII, the
results of the preceding sections are discussed in the context of past and
future work.

##  III Preliminaries

### Observables and perturbative quantities

– A gravitational wave observatory will detect a linear combination of the two
strain polarizations present in Einstein’s General Relativity (GR), $h_{+}$
and $h_{\times}$ Blanchet (2014). We may parameterize each using Boyer-
Lindquist coordinates, $\\{t,r,\theta,\phi\\}$ where, asymptotically far from
a radiating source, $t$ may be defined to measure time elapsed for an
observer, $r$ the luminosity distance, and $(\theta,\phi)$ the usual spherical
polar angles Boyer and Lindquist (2004). Since all $\theta$ dependence is in
the form of $u=\cos(\theta)$, $u$ will henceforth be used Leaver (1985).

From the perspective of Teukolsky’s equation, it is convenient to work with
the curvature scalar $\psi_{4}$ which is related to the strain polarizations
via

$\displaystyle h\;$ $\displaystyle=\;h_{+}\,-\,i\,h_{\times}\;$ (1a)
$\displaystyle=\;\int_{-\infty}^{t}\int_{-\infty}^{t^{\prime}}\psi_{4}(t^{\prime\prime},r,u,\phi)\;dt^{\prime\prime}dt^{\prime}\;.$
(1b)

In Eq. (1b), $\psi_{4}$ is the fifth scalar of the Weyl tensor in the Newman-
Penrose formalism Newman and Penrose (1962). It is a field of spin-weight
$s=-2$, and encodes all relevant information about gravitational radiation
Newman and Penrose (1962, 1966); Ruiz  _et al._ (2008); Merlin  _et al._
(2016). For a single isolated BH, the rescaled Weyl scalar

$\displaystyle\psi\;$ $\displaystyle=\;-(r-iau)^{-4}\;\psi_{4}\;.$ (2)

satisfies Teukolsky’s equation Teukolsky (1973); Press and Teukolsky (1973);
Teukolsky and Press (1974). Schematically, if we denote the second order
linear partial differential operator within Teukolsky’s equation as
$\mathcal{L}_{tru\phi}$, then for isolated black holes,

$\displaystyle\mathcal{L}_{tru\phi}\;\psi\;=\;0\;.$ (3)

Since $\mathcal{L}_{tru\phi}$ has partial derivatives on each spacetime
coordinate, the solutions $\psi$ are waves Teukolsky (1973).

### Quasinormal Modes

– The QNMs of an isolated Kerr BH comprise a sequence of solutions to Eq. (3)
that are ingoing waves at the BH horizon and purely outgoing radiation towards
spatial infinity Teukolsky (1973); Leaver (1985). These conditions on a QNM’s
phase velocity are known to constrain the functional form of $\psi$ by
defining the allowed radial behavior asymptotically near the event horizon and
spatial infinity Cook and Zalutskiy (2014); Press and Teukolsky (1973); Leaver
(1985).

An individual QNM has spatiotemporal dependence given by the product of four
functions,

$\displaystyle\psi\;\propto\;R(r)\,S(u)\,e^{-i\tilde{\omega}t}\,e^{-im\phi}\;.$
(4)

In Eq. (4), $R(r)$ will be referred to as simply a “radial function” and
$S(u)$ a spheroidal harmonic Teukolsky (1973). Both are implicitly
parameterized by the field’s spin weight, $s$, as well as the azimuthal and
temporal mode numbers, $m$ and $\tilde{\omega}$ respectively. For given values
of $\tilde{\omega}$ and $m$, the study of possible solutions for each $R(r)$
and $S(u)$ result in two additional mode numbers, $n$ and $\ell$, where $\ell$
is known to be closely related to the order of Jacobi polynomials London
(2023a); Fackerell and Crossman (2008); London (2023b). One goal of the
present work is to provide an analogous mathematical understanding of the
colloquially named “overtone” number, $n$.

In Eq. (4), $\tilde{\omega}\;=\;\omega\,-\,i/\tau$ is complex valued QNM or
“ringdown” frequency, and $\tau$ must be positive for time domain stability of
the BH Press and Teukolsky (1973); Leaver (1985); Whiting (1989).
Consequently, in the time domain, each QNM is a damped sinusoid with central
frequency $\omega$, and damping time $\tau$. Teukolsky’s equation constrains
the functional form of each QNM, and the details of the perturber constrain
its overall excitation amplitude Leaver (1986); Berti and Cardoso (2006).

Furthermore, application of Eq. (4) to Teukolsky’s equation yields that Eq.
(3) separates (with separation constant $A$) into

$\displaystyle\mathcal{L}{{}_{u}}\,S(u)\;$ $\displaystyle=\;-A\,S(u)\;,$ (5a)
$\displaystyle\mathcal{L}{{}_{r}}\,R(r)\;$ $\displaystyle=\;+A\,R(r)\;.$ (5b)

Henceforth, Eq. (5a) and Eq. (5b) will be referred to as Teukolsky’s angular
and radial equation, respectively. They are eigenvalue equations, where $-A$
and $+A$ are the respective complex valued333The complex nature of the
eigenvalues is a consequence of the fact that Eqs. (5a) and (5b) have complex
valued terms. It will be seen in Eq. (14) that this coincides with a non-
conjugate symmetric scalar product. eigenvalues London (2023a, b); Fiziev
(2009b); Leaver (1986). Therefore, $S(u)$ and $R(r)$ are the respective
eigenfunctions. In Eq. (5), $\mathcal{L}_{u}$ and $\mathcal{L}_{r}$ are second
order linear differential operators that will be referred to as Teukolsky’s
angular and radial operators, respectively. Each implicitly depends on $s$,
$m$ and $\tilde{\omega}$, as well as the BH mass $M$, and spin
$a=|\vec{J}|/M$.

Considered on their own, i.e. outside of the explicit physical context
described thus far where $\tilde{\omega}$ is assumed to be a QNM frequency,
Eqs. (5a) and (5b) are known to have a discrete space of eigenfunctions, which
may be labeled in ${\ell^{\prime}}$ and ${n^{\prime}}$, respectively Leaver
(1985); Berti  _et al._ (2007). Their eigenvalues implicitly inherit these
labels,

---
Figure 1:  Distributions of radial (black dots) and angular eigenvalues (open
circles) for $s=-2$, BH dimensionless spin $a/M=0.7$, and the
$(\ell,m,n)=(4,1,3)$ QNM with $M\tilde{\omega}_{413}=0.8545-0.6305i$. Since
both eigenvalue distributions have been computed with a QNM frequency, they
have one coincident value corresponding to
$(\ell^{\prime},n^{\prime})=(\ell,n)=(4,3)$. Radial egenvalues shown have been
computed using the results of the current work, and verified using Leaver’s
method Leaver (1985).

$\displaystyle\mathcal{L}^{\prime}{{}_{u}}\,S_{{\ell^{\prime}}}(u)\;$
$\displaystyle=\;-A_{{\ell^{\prime}}}\,S_{{\ell^{\prime}}}(u)\;,$ (6a)
$\displaystyle\mathcal{L}^{\prime}{{}_{r}}\,R_{{n^{\prime}}}(r)\;$
$\displaystyle=\;+A_{{n^{\prime}}}\,R_{{n^{\prime}}}(r)\;.$ (6b)

Equations (6a) and (6b) hold for any frequency parameter
$\tilde{\omega}^{\prime}$ that is independent of $\ell^{\prime}$ and
$n^{\prime}$. Since $\mathcal{L}{{}_{u}}$ and $\mathcal{L}{{}_{r}}$ depend on
the frequency parameter, they are also primed in Eq. (6). It is only for the
QNM frequencies that $A_{\ell^{\prime}}=A_{n^{\prime}}=A$, and Eq. (5) is
recovered. Since $\ell^{\prime}$ and $n^{\prime}$ are independent, for any one
QNM frequency, there will be only one value of $\ell^{\prime}$ and one value
of $n^{\prime}$ such that $A_{\ell^{\prime}}=A_{n^{\prime}}$. These special
labels are those for which $\ell^{\prime}=\ell$ and $n^{\prime}=n$. For these
reasons, the QNM frequencies and their separation constants are often labeled
with spatial mode numbers, i.e. $\tilde{\omega}_{{\ell mn}}{}$ and $A_{{\ell
mn}}$ respectively Berti  _et al._ (2006a); London  _et al._ (2014); Berti
_et al._ (2016); Magaña Zertuche  _et al._ (2021). As a result, the QNMs are a
non-trivial infinite sequence, wherein each element satisfies Eq. (6) with a
nominally distinct frequency parameter,
$\tilde{\omega}^{\prime}=\tilde{\omega}_{{\ell mn}}$ London (2023b).

An example of the angular and radial eigenvalue distributions is shown in
Figure 1. There, the typical qualitative features of each distribution are
shown: the angular eigenvalues, $A_{\ell^{\prime}}$, will typically be spaced
horizontally along the real axis, and the radial eigenvalues,
$A_{n^{\prime}}$, will vary most strongly in their imaginary part. Since the
case shown has been constructed using a QNM frequency, there will be a single
coincident eigenvalue, which is the one of physical interest (i.e. the one
compatible with the construction of Eq. 5).

### Notation

– For convenience of notation, QNM frequencies will henceforth be referred to
as $\tilde{\omega}_{{\ell mn}}$, and $\tilde{\omega}$ will refer to an
arbitrary frequency parameter. The related radial eigenvalue relation will be
referred to simply as

$\displaystyle\mathcal{L}{{}_{r}}\,R(r)\;$ $\displaystyle=\;A\,R(r)\;,$ (7)

where it is implicitly understood that the eigenvalues, $A$, and
eigenfunctions, $R(r)$, are elements of that equation’s discrete spectra. A
transformed version of the radial equation will be encountered shortly. At
that time, the same notational strategy will be used.

### A radial scalar product

– The angular equation, Eq. (6a), is well-known to be of the Sturm-Liouville
type (i.e. $\mathcal{L}_{u}$ is self-adjoint) London (2023b, a); Courant and
Hilbert. (1954); Lax (2002). Consequently, its eigenfunctions,
$S_{\ell^{\prime}}(u)$, are complete and orthogonal with respect to a scalar
product, ${\langle{a}\,|\,{b}\rangle}$, defined by the integral of $a(u)b(u)$
on $u\in[-1,1]$. The application of this angular scalar product is fundamental
to the use of spherical and spheroidal harmonics in gravitational wave theory
and data analysis Blanchet (2014); Khan  _et al._ (2016); García-Quirós  _et
al._ (2020).

In contrast, the radial equation, Eq. (6b), is known to not explicitly be of
the Sturm-Liouville type. This simply means that, in practice, the application
of Sturm-Liouville theory must begin with the construction of a weighted
scalar product, under which $\mathcal{L}_{r}$ is self-adjoint Courant and
Hilbert. (1954); Pinchover and Rubinstein (2005); Kristensson (2010). Below,
we provide a brief overview of the construction and evaluation of such a
radial scalar product. We refer the reader to Paper I for a complete
pedagogical development.

Using the compactified radial coordinate,

$\displaystyle\xi\;$ $\displaystyle=\;(r-r_{+})/(r-r_{-})\;,$ (8)

and applying the QNMs’ asymptotic boundary conditions via an appropriate
rescaling,

$\displaystyle R(r(\xi))\;=\;\mu(\xi)f(\xi)\;,$ (9)

the radial problem is transformed so that the radial function of interest is
now $f(\xi)$, and differential operator $\mathcal{L}_{r}$ becomes

$\displaystyle\mathcal{L}{{}_{\xi}}\;$
$\displaystyle=\;(\text{C}_{0}+\text{C}_{1}(1-\xi))$ (10a)
$\displaystyle+\;\left(\text{C}_{2}+\text{C}_{3}(1-\xi)+\text{C}_{4}(1-\xi)^{2}\right)\partial_{\xi}$
(10b) $\displaystyle+\;\xi(\xi-1)^{2}\partial_{\xi}^{2}\;.$ (10c)

The operator $\mathcal{L}_{\xi}$ may be compactly written by defining the
operator $\mathcal{D}{}_{\xi}$ to be the terms within Eqs. (10b) and (10c),

$\displaystyle\mathcal{L}_{\xi}\;=\;(\text{C}_{0}+\text{C}_{1}(1-\xi))\,+\,\mathcal{D}_{\xi}\;.$
(11)

Note that in Eq. (8), $\xi$ maps $r\in[0,\infty)$ to $\xi\in[0,1)$, thereby
conpactifying the radial coordinate. In Eq. (9), $\mu(\xi)$ is comprised of
factors related to the singular exponents of $\mathcal{L}_{r}$. Its explicit
form is given in Eq. (12) of Paper I. In Eq. (10), $\text{C}_{0}$ through
$\text{C}_{4}$ are complex valued constants determined by the BH mass $M$, the
BH spin $a$, the mode numbers $\\{\tilde{\omega},m,l,n\\}$, and the field’s
spin weight $s$:

$\displaystyle\text{C}_{0}\;$
$\displaystyle=\;-2am\tilde{\omega}-2i\tilde{\omega}(-\delta+M(2s+1))$ (12a)
$\displaystyle\quad\;\;+\tilde{\omega}^{2}(\delta+M)(\delta+7M)\;,$
$\displaystyle\text{C}_{1}\;$
$\displaystyle=\;8M^{2}\tilde{\omega}^{2}+s(4iM\tilde{\omega}-1)+6iM\tilde{\omega}-1$
(12b)
$\displaystyle\quad\;\;-(4M\tilde{\omega}+i)\frac{\left(am-2M^{2}\tilde{\omega}\right)}{\delta}\;,$
$\displaystyle\text{C}_{2}\;$ $\displaystyle=\;4i\delta\tilde{\omega}\;,$
(12c) $\displaystyle\text{C}_{3}\;$
$\displaystyle=\;-2(s+1)+4i\tilde{\omega}(M-\delta)\;,$ (12d)
$\displaystyle\text{C}_{4}\;$
$\displaystyle=\;s+3-6iM\tilde{\omega}+\frac{i\left(am-2M^{2}\tilde{\omega}\right)}{\delta}\;,$
(12e) $\displaystyle\delta\;$ $\displaystyle=\;\sqrt{M^{2}-a^{2}}\;.$ (12f)

Lastly, note that none of $\text{C}_{0}$ through $\text{C}_{4}$ become zero
when $a=0$. As a result, Eq. (11) is of the confluent Heun type for both
Schwarzschild and Kerr London (2023a).

Together, Eqs. (8-12) allow the radial eigenvalue relation (Eq. 7) to be
rewritten as

$\displaystyle\mathcal{L}{{}_{\xi}}\,f(\xi)\;$ $\displaystyle=\;A\,f(\xi)\;.$
(13a) The equivalent statement in standard bra-ket notation is
$\displaystyle\mathcal{L}{{}_{\xi}}\,|{f}\rangle\;$
$\displaystyle=\;A\,|{f}\rangle\;.$ (13b)

In Eq. (13b), and henceforth, the use of bra-ket notation is consistent with
the fact that $f(\xi)$ is a vector in space whose bilinear product is the
radial scalar product of interest. The practical use of kets, e.g.
$|{\Phi}\rangle$, and bras, $\langle{\Phi}|$, will be demonstrated in the
context of the scalar product.

A scalar product between two smooth functions of $\xi$, $a(\xi)$ and $b(\xi)$,
is defined as

$\displaystyle{\langle{\mathrm{a}}\,|\,{\mathrm{b}}\rangle}\;=\;\int_{0}^{1}\,\mathrm{a}(\xi)\,\mathrm{b}(\xi)\;\text{W}(\xi)\;d\xi\;,$
(14)

where $\text{W}(\xi)$ is the product’s weight function. In turn,
$\text{W}(\xi)$ is defined such that the operator, $\mathcal{L}_{\xi}$, is
equal to its adjoint, ${\mathcal{L}}^{\dagger}_{\xi}$, with respect to the
scalar product:

$\displaystyle{\langle{a}\,|\,{\mathcal{L}_{\xi}b}\rangle}\;$
$\displaystyle=\;{\langle{{\mathcal{L}}^{\dagger}_{\xi}a}\,|\,{b}\rangle}\;$
(15) $\displaystyle=\;{\langle{{\mathcal{L}}_{\xi}a}\,|\,{b}\rangle}\;.$

Note that there is no conjugation in Eq. (14). As a result, Eq. (15)’s adjoint
operator, ${\mathcal{L}}^{\dagger}_{\xi}$, does not involve complex
conjugation.

By combining the definition of $\mathcal{L}_{\xi}$ (Eq. 10) with the
requirement of self-adjointness (Eq. 15) and then applying integration by
parts, it may be shown that

$\displaystyle\text{W}(\xi)\;$
$\displaystyle=\;\xi^{\text{B}_{0}}(1-\xi)^{\text{B}_{1}}e^{\frac{\text{B}_{2}}{1-\xi}}\;,$
(16)

where constants $\text{B}_{0}$ through $\text{B}_{2}$ are,

$\displaystyle\text{B}_{0}\;$
$\displaystyle=\;\text{C}_{2}+\text{C}_{3}+\text{C}_{4}-1$ (17a)
$\displaystyle\text{B}_{1}\;$ $\displaystyle=\;-\text{C}_{2}-\text{C}_{3}-2$
(17b) $\displaystyle\text{B}_{2}\;$ $\displaystyle=\;\text{C}_{2}\;.$ (17c)

In practice, if the series expansion of two functions, $\text{a}(\xi)$ and
$\text{b}(\xi)$, exist and are known,

$\displaystyle\text{a}(\xi)\;=\;\sum_{j=0}\,\text{a}_{j}\,\xi^{j}\quad\text{
and }\quad\text{b}(\xi)\;=\;\sum_{k=0}\,\text{b}_{k}\,\xi^{k}\;,$ (18)

then the scalar product may be evaluated as

$\displaystyle{\langle{\text{a}}\,|\,{\text{b}}\rangle}\;=\;\sum_{j,k}\;\text{a}_{j}\,\text{b}_{k}\,\left<\xi^{j+k}\right>\;.$
(19)

In Eq. (19), $\left<\xi^{j+k}\right>$ are called monomial moments Chihara
(2011). Using the definition of the weight function, it has been observed
(e.g. Cho  _et al._ (2009); Leaver (1986); Abramowitz and Stegun (1964)) that
each monomial moment $\left<\xi^{p}\right>$ can be evaluated using the Gamma
function, $\Gamma(x)$, and the Tricomi confluent hypergeometric function,
$U(x,y,z)$,

$\displaystyle\left<\xi^{p}\right>\;=\;e^{\text{B}_{2}}\,\Gamma($
$\displaystyle\text{B}_{0}+p+1)$
$\displaystyle\times\;U(1+\text{B}_{0}+p,-\text{B}_{1},-\text{B}_{2})\;.$ (20)

Equations (19) and (III) enable the scalar product to be evaluated for complex
values of $\text{B}_{0}$, $\text{B}_{1}$ and $\text{B}_{2}$ via the standard
analytic continuations of $\Gamma(x)$ and $U(x,y,z)$ London (2023a).

### Confluent Heun polynomials

– In many areas of physics, polynomial solutions to a system’s equations of
motion play an important role in describing that system’s natural modes
Sakurai (1993); Cook and Zalutskiy (2014); London (2023b); Fackerell and
Crossman (2008). For the QNM’s radial problem (Eq. 13), one approach to
understanding polynomial solutions begins by studying the differential part of
$\mathcal{L}_{\xi}$, namely

$\displaystyle\mathcal{D}{{}_{\xi}}\;$
$\displaystyle=\;\left(\text{C}_{2}+\text{C}_{3}(1-\xi)+\text{C}_{4}(1-\xi)^{2}\right)\partial_{\xi}$
(21) $\displaystyle\hskip 14.22636pt+\;\xi(\xi-1)^{2}\partial_{\xi}^{2}\;.$

This operator defines polynomials that satisfy a two-parameter eigenvalue
problem,

$\displaystyle\mathcal{D}{{}_{\xi}}\;|{y_{pk}}\rangle\;$
$\displaystyle=\;(\lambda_{pk}\;+\;\mu_{p}\,\xi)\;|{y_{pk}}\rangle\;,$ (22a)
where $\displaystyle\mu_{p}\;$ $\displaystyle=\;p\,(\,p+\text{C}_{4}-1\,)\;.$
(22b)

In Eq. (22a), $\lambda_{pk}$ is the eigenvalue, and the linear factor,
$\mu_{p}$, is required for the polynomials, $|{y_{pk}}\rangle$, to terminate
after $p+1$ terms London (2023a); Ronveaux and Arscott (1995). For each
polynomial order $p$, there are $p+1$ solutions for $\lambda_{pk}$ that allow
for termination DLMF ; London (2023a). The label $k$ heuristically orders
these solutions according to the imaginary part of, $\text{Im}\,\lambda_{pk}$
London (2023a). Since Eq. (22a) is a confluent Heun equation,
$|{y_{pk}}\rangle$ comprise a class, or multiplex, of order-$p$ confluent Heun
polynomials Ronveaux and Arscott (1995).

The ostensibly unusual (two parameter) nature of Eq. (22a) is a consequence of
$\mathcal{D}_{\xi}$’s polynomial coefficient functions (Eq. 21). For example,
when acting upon any polynomial of order $p$,
$\xi(1-\xi)^{2}\partial_{\xi}^{2}$ will generally result in a polynomial of
order $p+1$. The same is true of
$\left(\text{C}_{2}+\text{C}_{3}(1-\xi)+\text{C}_{4}(1-\xi)^{2}\right)\partial_{\xi}$.
Therefore, $\mathcal{D}_{\xi}$ generally increases the order of polynomials by
one, which precludes the usual eigenrelation, where order is preserved.

A consequence of the multi-parameter eigenrelation is that the confluent Heun
polynomials do not obey exactly the same kind of orthogonality found in the
classical polynomials. This may be seen by noting that, like
$\mathcal{L}_{\xi}$, $\mathcal{D}_{\xi}$ is self-adjoint with respect to the
scalar product,

$\displaystyle{\langle{y_{pk^{\prime}}}\,|\,{\mathcal{D}_{\xi}y_{pk}}\rangle}\;=\;{\langle{\mathcal{D}_{\xi}y_{pk^{\prime}}}\,|\,{y_{pk}}\rangle}\;.$
(23)

Applying the eigenvalue relationship, Eq. (22a), to both sides of Eq. (23)
yields

$\displaystyle(\lambda_{pk}-\lambda_{p^{\prime}k^{\prime}})\,$
$\displaystyle{\langle{y_{p^{\prime}k^{\prime}}}\,|\,{y_{pk}}\rangle}$
$\displaystyle=\;(\mu_{p^{\prime}}-\mu_{p})\,\langle{y_{p^{\prime}k^{\prime}}}|\xi\;|{y_{pk}}\rangle\;.$
(24)

For self-adjoint and single-parameter eigenvalue problems, the right hand side
of Eq. (III) is identically zero, resulting in the well-known orthogonality of
eigenfunctions with unique eigenvalues Axler (2015); Sakurai (1993). Since
$\mu_{p}$ is not generally zero, only the polynomials of like order (i.e.
$\mu_{p^{\prime}}-\mu_{p}=0$) are orthogonal,

$\displaystyle{\langle{y_{pk^{\prime}}}\,|\,{y_{pk}}\rangle}\;$
$\displaystyle=\;\delta_{k^{\prime}k}\;.$ (25)

In Eq. (25), $|{y_{pk}}\rangle$ are normalized, i.e. scaled such that
${\langle{y_{pk}}\,|\,{y_{pk}}\rangle}=1$ London (2023a).

It was observed in Paper I that $\lambda_{pk}$ are non-degenerate across
physically relevant values of BH mass and spin. Consequently, each $p+1$
dimensional vector space of polynomials $y_{pk}$ is expected to be a basis,
meaning that the identity operator for order-$p$ polynomials may be written
as,

$\displaystyle{\mathbb{I}}_{p}\;=\;\sum_{k=0}^{p}\,|{y_{pk}}\rangle\langle{y_{pk}}|\;.$
(26)

However, since each basis for order $p$ polynomials spans all previous bases,
the total set of confluent Heun polynomials satisfying Eq. (22a) is redundant,
or rather overcomplete Christensen (2003).

##  IV The radial problem with Confluent Heun Polynomials

The properties of confluent Heun polynomials have three implications for each
QNM’s radial problem (Eq. 13b).

The first is that the radial problem, when represented in a basis of order-$p$
confluent Heun polynomials, may written simply, with a deformation parameter,

$\displaystyle\alpha_{p}\;$ $\displaystyle=\;\mu_{p}-\text{C}_{1}\;,$ (27)

that determines the size of all off-diagonal terms. Concurrently, off diagonal
terms are proportional to
${\langle{y_{pk^{\prime}}}\,|\,{\xi\,|y_{pk}}\rangle}$, meaning that, for non-
zero $\alpha_{p}$, the matrix representation of the radial problem may not be
approximately diagonal in the sense that elements do not necessarily decrease
far away from $k^{\prime}=k$.

The second implication is that, due to properties of $\alpha_{p}$, a kind of
duality exists in the radial problem’s solution space. Most simply, if
$\alpha_{p}=0$ and $p$ is a non-negative integer, then $p+1$ solutions will be
polynomials, while the remaining infinity of solutions will be infinite
series.

The final implication is that the physical QNMs will not generally be
approximately polynomial in nature. Along with the aforementioned duality,
this implies that there will typically be two, sometimes overlapping, sectors
of solution space. It happens that this concurrently highlights practical
issues that arise when attempting to approximate a function with an arbitrary
fixed-order polynomial basis.

While these implications provide an instructive reference for the basic
properties of non-polynomial solutions to the radial equation, they also
motivate the development of a basis comprised of uniquely ordered polynomials.
That is the topic of the next section. For now, the radial problem with
confluent Heun polynomials is to be discussed in more detail.

### The radial problem

– Any order-$p$ confluent Heun polynomial basis may be used to approximate the
radial problem,

$\displaystyle{\mathbb{I}}_{p}\,\mathcal{L}{{}_{\xi}}\,{\mathbb{I}}_{p}\,|{f}\rangle\;$
$\displaystyle=\;A\,{\mathbb{I}}_{p}\,|{f}\rangle\;,$ (28a)
$\displaystyle\hat{L}{{}_{p}}^{(\text{cH})}\,\vec{\text{f}}\;$
$\displaystyle=\;A\,\vec{\text{f}}\;.$ (28b)

In Eq. (28a), insertion of the identity operator (Eq. 26) signals that e.g.
$|{f}\rangle$ will be represented in the basis of order-$p$ confluent Heun
polynomials. In Eq. (28b), Eq. (26) has been applied, and it has been
recognized that the equation is equivalent to a matrix eigenvalue equation,
where $\hat{L}{{}_{p}}^{(\text{cH})}$ is a $(p+1)\times(p+1)$ matrix with
elements,

$\displaystyle
L^{(\text{cH})}_{p,\,k^{\prime}k}\;=\;{\langle{y_{pk^{\prime}}}\,|\,{\mathcal{L}_{\xi}\,|\,y_{pk}}\rangle}\;,$
(29)

and $\vec{\text{f}}$ has elements ${\langle{y_{pk}}\,|\,{f}\rangle}$.

In many cases, $\hat{L}{{}_{p}}^{(\text{cH})}$ will have many off-diagonal
terms. This may be understood by recalling that
$\mathcal{L}_{\xi}=(\text{C}_{0}+\text{C}_{1}[1-\xi])+\mathcal{D}_{\xi}$ (i.e.
Eq. 11), and then applying the confluent Heun polynomial eigenvalue
relationship to
${\langle{y_{pk^{\prime}}}\,|\,{\mathcal{L}_{\xi}\,|\,y_{pk}}\rangle}$,
yielding

$\displaystyle{\langle{y_{pk^{\prime}}}\,|\,{\mathcal{L}_{\xi}\,|\,y_{pk}}\rangle}\;=\;$
$\displaystyle(\text{C}_{0}+\text{C}_{1}+\lambda_{pk})\,\delta_{k^{\prime}k}\,$
$\displaystyle+\,\alpha_{p}\,{\langle{y_{pk^{\prime}}}\,|\,{\xi\,|\,y_{pk}}\rangle}\;,$
(30a)

In Eq. (30a), $\alpha_{p}$ is defined by Eq. (27).

If $\alpha_{p}=0$, Eqs. (28b) and (30a) yield that the QNM separation constant
is simply given by

$\displaystyle A_{pk}\;=\;{\text{C}_{0}}+{\text{C}_{1}}+\lambda_{pk}\;.$ (31)

If a QNM’s radial functions is well approximated by a single confluent Heun
polynomial, then Eq. (31) will be similarly approximate.

When $\alpha_{p}\neq 0$, the off-diagonal elements of
$\hat{L}{{}_{p}}^{(\text{cH})}$ will be proportional to
${\langle{y_{pk^{\prime}}}\,|\,{\xi\,|\,y_{pk}}\rangle}$. It was found in
Paper I that $|{\langle{y_{pk^{\prime}}}\,|\,{\xi\,|\,y_{pk}}\rangle}|$ do not
typically decrease in value as $|k^{\prime}-k|$ increases. Consequently,
$\hat{L}{{}_{p}}^{(\text{cH})}$ will typically not be non-band-diagonal when
$\alpha_{p}\neq 0$.

### Polynomial / non-polynomial duality

– From inspection of Eq. (30a), it is clear that when $\alpha_{p}=0$,
$\hat{L}{{}_{p}}^{(\text{cH})}$ is diagonal. This implies that
$\hat{L}{{}_{p}}^{(\text{cH})}$ may be approximately diagonal when
$|\alpha_{p}|\ll 1$.

This concept may be refined by defining pseudo-polynomial orders,
${p_{\star}^{+}}$ and ${p_{\star}^{-}}$, such that
$\alpha_{{p_{\star}^{\pm}}}$ is zero:

$\displaystyle\alpha_{{p_{\star}^{\pm}}}\;$
$\displaystyle=\;\mu_{{p_{\star}^{\pm}}}-\text{C}_{1}$ (32a)
$\displaystyle=\,{p_{\star}^{\pm}}\,(\,{p_{\star}^{\pm}}+\text{C}_{4}\,-\,1\,)\,-\;\text{C}_{1}=\;0\;.$
(32b)

In Eq. (32b), $\mu_{{p_{\star}^{\pm}}}$ has been rewritten using Eq. (22b).
Solving for ${p_{\star}^{\pm}}$, one finds that

$\displaystyle{p_{\star}^{\pm}}=\frac{1}{2}\left[(1-\text{C}_{4})\pm\sqrt{4\text{C}_{1}+(1-\text{C}_{4})^{2}}\right]\;.$
(33)

In Eq. (33), since ${p_{\star}^{+}}$ and ${p_{\star}^{-}}$ always come as a
pair, they are duals of each other. Henceforth, ${p_{\star}^{+}}$ wll be
referred to as the upper pseudo-polynomial order, and ${p_{\star}^{-}}$ will
be referred to as the lower pseudo-polynomial order.

$n$ |  $\quad\quad\quad\quad\quad{p_{\star}^{+}}$ |  $\quad\quad\quad\quad\quad{p_{\star}^{-}}$
---|---|---
$0$ |  $\quad 1.3879+0.5964i$ |  $\quad-0.6768+2.1304i$
$1$ |  $\quad 2.1725+0.5415i$ |  $\quad-0.0230+2.0846i$
$2$ |  $\quad 2.9791+0.4394i$ |  $\quad 0.6490+1.9996i$
$4$ |  $\quad 4.6187+0.1537i$ |  $\quad 2.0152+1.7615i$
$6$ |  $\quad 6.2586+0.0743i$ |  $\quad 3.3816+1.6954i$
$8$ |  $\quad 15.2263+-1.9170i$ |  $\quad 10.8539+0.0362i$
Table 1:  Select pseudo-polynomial orders, ${p_{\star}^{+}}$ and
${p_{\star}^{-}}$, for BH spin of $a/M=0.7$, and QNM indices
$(s,\ell,m)=(-2,2,2)$. Parameters for select QNM overtone numbers, $n$, are
shown.

For the moment, we are concerned with the implications of this duality on
solutions to the radial problem (Eq. 28), the QNMs radial functions, and the
ultimate practicality of confluent Heun polynomial representations for non-
polynomial problems.

Equation (33) implies a number of simple relationships between $\alpha_{p}$,
$\mu_{p}$, ${p_{\star}^{+}}$ and ${p_{\star}^{-}}$. For example, Eq. (33) may
used to show that

$\displaystyle{p_{\star}^{+}}+{p_{\star}^{-}}\;$
$\displaystyle=\;1-{\text{C}_{4}}\;.$ (34)

It also follows that ${\text{C}_{1}}=-{p_{\star}^{+}}{p_{\star}^{-}}$. In
terms of physical parameters (Eq. 12), ${p_{\star}^{+}}$ and ${p_{\star}^{-}}$
are

$\displaystyle{p_{\star}^{+}}\;$ $\displaystyle=\;-1\,+\,4iM\tilde{\omega}\;,$
(35a) $\displaystyle{p_{\star}^{-}}\;$
$\displaystyle=\;{p_{\star}^{+}}\,-s\,-\,2iM\tilde{\omega}\,-\,\frac{i\left(am-2M^{2}\tilde{\omega}\right)}{\delta}\;.$
(35b)

Lastly, Eq. (33) allows $\alpha_{p}$ to be written as

$\displaystyle\alpha_{p}\;=\;(p-{p_{\star}^{+}})(p-{p_{\star}^{-}})\;.$ (36)

From Eq. (34), it may be seen that if ${p_{\star}^{+}}$ is a polynomial order,
then ${p_{\star}^{-}}$ typically will not be unless ${\text{C}_{4}}$ is a
sufficiently negative integer. In such cases, the solution space of Eq. (28)
will be comprised of ${p_{\star}^{+}}+1$ polynomial solutions, and an infinity
of non-polynomial solutions. In this way, there exists a confluent Heun
polynomial/non-polynomial duality.

For QNMs, it is expected that neither ${p_{\star}^{+}}$ nor ${p_{\star}^{-}}$
will be polynomial powers, and that the proximity of each to an integer will
vary with different values of the overtone index, $n$. In particular,
solutions to the QNM radial problem may only be closely associated with
confluent Heun polynomials when $\tilde{\omega}_{{\ell mn}}$ are such that
${p_{\star}^{+}}$ or ${p_{\star}^{-}}$ is close to a positive integer in the
sense that,

$\displaystyle|\text{Im}\,{p_{\star}^{\pm}}|\;\ll\;\text{Re}\,{p_{\star}^{\pm}}\;,$
(37a) and
$\displaystyle{p_{\star}^{\pm}}-\lfloor\text{Re}\,{p_{\star}^{\pm}}\rfloor\;\ll\;1\;.$
(37b)

In Eq. (37), $\lfloor x\rfloor$ is the integral part of a real number (i.e.
$\lfloor 3.1416\rfloor=3$).

Table (1) provides examples of ${p_{\star}^{\pm}}$ for select QNM cases with
mode numbers $(\ell,m)=(2,2)$. There it is observed that the QNM’s
${p_{\star}^{+}}$ only very roughly tracks the overtone index, being
approximately off by one for e.g. $n=2$, but off by significantly less than
one for $n=6$. Notably, $n=8$ departs from this trend, having a
$\text{Re}\,{p_{\star}^{+}}$ that is just under two times the overtone index.
In all cases, it is observed that $n<{p_{\star}^{+}}$. Further, the case of
the fundamental QNM, $n=0$, demonstrates a general observation, namely that
the $n=0$ QNMs are typically not close to confluent Heun polynomial systems in
that their ${p_{\star}^{\pm}}$ are not close to integers.

### Polynomial noise

– Since the QNMs ${p_{\star}^{\pm}}$ are typically complex valued, it is not
guaranteed that using an order-$p$ confluent Heun polynomial basis with e.g.
$p=\lfloor\mathrm{max}(|{p_{\star}^{+}}|,|{p_{\star}^{-}}|)\rfloor$ will yield
an accurate representation of the radial problem. Worse, simply increasing the
order of the confluent Heun polynomial basis (i.e. increasing $p$) does not
necessarily result in a more precise representation. This may be understood as
resulting from a basic limitation of approximation with a fixed-order
polynomial bases.

As a somewhat extreme but plausible example, consider the order-$50$ confluent
Heun polynomial approximation of the number $1$,

$\displaystyle|{1}\rangle\;=\;\sum_{k=0}^{50}\,|{y_{50,k}}\rangle{\langle{y_{50,k}}\,|\,{1}\rangle}\;.$
(38)

Each polynomial, $|{y_{50,k}}\rangle$, has $51$ coefficients, and on the right
hand side of Eq. (38), values of ${\langle{y_{50,k}}\,|\,{1}\rangle}$ must
conspire such that only constant terms survive the sum (as is required since
$1$ is a constant). In practice, finite precision means that the cancellations
will not be exact. As a result, the net approximation will be contaminated by
high frequency polynomial oscillations, or rather, polynomial noise.

Clearly, in this case it would be wise to represent $|{1}\rangle$ with the
order-$0$ confluent Heun polynomial. This is consistent with the completeness
of each confluent Heun polynomial fixed order basis: an order-$p$ confluent
Heun polynomial basis can exactly represent polynomials of the same order;
while lower order polynomials may also be represented, their representation in
polynomial bases of much larger order increases the likelihood of numerical
noise.

Further, for non-polynomial functions, it will not clear a priori which, if
any, fixed order confluent Heun polynomial basis should be used. For example,
since the QNMs radial functions, $|{f}\rangle$, are generally non-polynomial,
one might attempt to apply a very large order confluent Heun polynomial
representation,

$\displaystyle|{f}\rangle\;=\;\sum_{k=0}^{P}\,|{y_{Pk}}\rangle{\langle{y_{Pk}}\,|\,{f}\rangle}\;,$
(39)

where $P\gg 1$. In practice, the result would be contaminated by polynomial
noise: all ${\langle{y_{Pk}}\,|\,{f}\rangle}$ would have to be given to
arbitrarily high precision such that the sum’s $(P+1)^{2}$ polynomial
coefficients (i.e. there are $P+1$ polynomials, and each with $P+1$
coefficients) add without significant contamination from round-off error.

In contrast, a basis of $P+1$ classical polynomials (i.e. with unique orders
ranging from $0$ to $P$) would only have $\frac{1}{2}(P+1)^{2}$ coefficients
to total. More importantly, unlike the order-$P$ confluent Heun polynomials,
such a basis would have expansion coefficients (analogs of
${\langle{y_{Pk}}\,|\,{f}\rangle}$) that necessarily decrease as $k$
increases, meaning that fewer than $\frac{1}{2}(P+1)^{2}$ coefficients may
ultimately be of practical importance.

### Towards canonical polynomials

– Together, these ideas may be summarized as follows. Since an order-$p$
confluent Heun polynomial basis is complete, it may be used to approximate the
radial problem in a way that explicitly reduced to a polynomial solution in
the limit that $\alpha_{p}=0$. Further, it may be possible to treat general
solutions of the radial problem perturbatively, where
$\alpha_{p}=(p-{p_{\star}^{+}})(p-{p_{\star}^{-}})$ is the perturbative
parameter. However, since $\hat{L}{{}_{p}}^{(\text{cH})}$ will generally be
non-band diagonal, the order-$p$ representation of QNM radial functions,
$|{f}\rangle$, may not quickly converge to a useful approximation. Unlike
approximation with classical polynomials, it is not always possible to reduce
error by increasing the order of the confluent Heun polynomial basis used.

These deficits are underpinned by the fact that each order-$p$ confluent Heun
polynomial basis is a multiplex, i.e. all members of each confluent Heun
polynomial basis are polynomials of the same order. There is, therefore, good
reason to wonder if there exists a basis of polynomials, each with a distinct
order, that is orthogonal under the scalar product (Eq. 14).

##  V Natural & Canonical Polynomials

The previous discussion does not rule out the construction of polynomials that
are (i) complete in the infinite dimensional sense, (ii) orthogonal with
respect to the scalar product (Eq. 14), and (iii) can represent
$\mathcal{D}_{\xi}$ as a symmetric matrix. While such polynomials cannot be
exactly of the Heun type for the reasons stated in Sec. III, they may yet
possess the key features of classical polynomials, namely completeness and
orthogonality. In this sense, polynomials with properties (i) and (ii) will be
referred to as canonical. Polynomials satisfying (iii) will be referred to as
natural to $\mathcal{D}_{\xi}$.

For example, the Chebyshev polynomials defined on $\xi\in[0,1]$ are canonical,
but not natural because $\mathcal{D}_{\xi}$ is not self-adjoint with respect
to the Chebyshev inner-product. In contrast, the confluent Heun polynomials
discussed in Sec. III and Sec. IV are natural to $\mathcal{D}_{\xi}$, but due
to their overcompleteness and atypical orthogonality, they are not canonical.

The possibility that the confluent Heun polynomials may be used to construct
polynomials that are both natural to $\mathcal{D}_{\xi}$, and canonical, is
our present concern. Henceforth, such polynomials will be referred to as
canonical confluent Heun polynomials. When it is clear that confluent Heun
polynomials are being discussed, the canonical confluent Heun polynomials may
be simply referred to as the “canonical polynomials”.

In this section, two equivalent methods to construct canonical confluent Heun
polynomials are described. The first method pertains to the construction of
canonical polynomials directly from confluent Heun polynomials defined in Sec.
III. This method will be referred to as “direct construction”. The second
method is based upon the standard Gram-Schmidt process. This method will be
referred to as “Gram-Schmidt construction”.

This section concludes with comments on the equivalence between direct and
Gram-Schmidt construction, their completeness, and finally two relevant
representations of the canonical confluent Heun polynomials in terms of the
confluent Heun polynomials.

### Direct Construction

– Let a canonical confluent Heun polynomial of order $p$ be denoted
$|{u_{p}}\rangle$, where it is represented in the basis of order $p$ confluent
Heun polynomials,

$\displaystyle|{u_{p}}\rangle\;$
$\displaystyle=\;\sum_{k=0}^{p}\,\text{c}_{pk}\,|{{y}_{pk}}\rangle\;.$ (40)

The direct construction of $|{u_{p}}\rangle$ from $|{y_{pk}}\rangle$ requires
the determination of constants, $\text{c}_{pk}$, such that the resulting
polynomials are canonical (properties i and ii), and natural to
$\mathcal{D}_{\xi}$ (property iii). We begin our construction with property
(ii).

Explicitly, property (ii) means that

$\displaystyle{\langle{u_{p^{\prime}}}\,|\,{u_{p}}\rangle}\;\propto\;\delta_{p^{\prime}p}\;.$
(41)

We will now show that Eqs. (40) and (41) suffice to determine all expansion
coefficients, $c_{pk}$, via an iterative linear process.

Without loss of generality, let $c_{pp}=1$. It follows that, for $p=0$,

$\displaystyle|{u_{0}}\rangle=|{y_{00}}\rangle\;.$ (42)

For $p=1$, applying $c_{11}=1$ to Eq. (40) means that
$|{u_{1}}\rangle=c_{10}|{y_{10}}\rangle+|{y_{11}}\rangle$. Orthogonality (Eq.
41) then requires that

$\displaystyle{\langle{u_{0}}\,|\,{u_{1}}\rangle}\;=\;c_{10}{\langle{y_{00}}\,|\,{y_{10}}\rangle}+{\langle{y_{00}}\,|\,{y_{11}}\rangle}\;=\;0\;.$
(43)

Thus, the single unknown constant $c_{10}$ may be solved for, yielding

$\displaystyle
c_{10}\;=\;-\frac{{\langle{y_{00}}\,|\,{y_{11}}\rangle}}{{\langle{y_{00}}\,|\,{y_{10}}\rangle}}\;.$
(44)

For $p=2$, we may again consider $|{u_{2}}\rangle$ as defined in Eq. (40),
namely
$|{u_{2}}\rangle=c_{20}|{y_{20}}\rangle+c_{21}|{y_{21}}\rangle+|{y_{22}}\rangle$.
Here, $c_{20}$ and $c_{21}$ are the unknowns of interest. Requiring that
${\langle{u_{0}}\,|\,{u_{2}}\rangle}=0$ and
${\langle{u_{1}}\,|\,{u_{2}}\rangle}=0$ results in two coupled equations,

$\displaystyle
c_{20}{\langle{u_{0}}\,|\,{y_{20}}\rangle}+c_{21}{\langle{u_{0}}\,|\,{y_{21}}\rangle}+{\langle{u_{0}}\,|\,{y_{22}}\rangle}\;$
$\displaystyle=\;0\;,$ (45a) $\displaystyle
c_{20}{\langle{u_{1}}\,|\,{y_{20}}\rangle}+c_{21}{\langle{u_{1}}\,|\,{y_{21}}\rangle}+{\langle{u_{1}}\,|\,{y_{22}}\rangle}\;$
$\displaystyle=\;0\;.$ (45b)

Moving the last terms in Eqs. (45a) and (45b) to their respective right hand
sides yields that

$\displaystyle
c_{20}{\langle{u_{0}}\,|\,{y_{20}}\rangle}+c_{21}{\langle{u_{0}}\,|\,{y_{21}}\rangle}\;$
$\displaystyle=\;-{\langle{u_{0}}\,|\,{y_{22}}\rangle}\;,$ (46a)
$\displaystyle
c_{20}{\langle{u_{1}}\,|\,{y_{20}}\rangle}+c_{21}{\langle{u_{1}}\,|\,{y_{21}}\rangle}\;$
$\displaystyle=\;-{\langle{u_{1}}\,|\,{y_{22}}\rangle}\;.$ (46b)

Equation (46) is equivalent to the linear matrix equation,

$\displaystyle\hat{X}_{2}\;\vec{c}_{2}\;=\;\vec{b}_{2}\;,$ (47)

where $\hat{X}_{2}$ has elements
${X}_{2,jk}={\langle{u_{j}}\,|\,{y_{2k}}\rangle}$, $\vec{c}_{2}$ has elements
$c_{2j}$, and $\vec{b}_{2}$ has elements
$b_{2,j}=-{\langle{u_{j}}\,|\,{y_{22}}\rangle}$.

In Eqs. (46) and (47), the essential features of subsequent cases (i.e. $p>2$)
are present: For any $|{u_{p}}\rangle$, there will be $p$ equations of the
form of Eq. (46), with a total of $p$ unknowns, and this system of linear
equations will generally be equivalent to a linear matrix equation,

$\displaystyle\hat{X}_{p}\;\vec{c}_{p}\;=\;\vec{b}_{p}\;.$ (48)

As in Eq. (48), $\hat{X}_{p}$ and $\vec{b}_{p}$ will generally depend on all
polynomials $|{u_{p^{\prime}}}\rangle$ with $p^{\prime}<p$, i.e. all previous
iterations. Since both $\\{y_{pk}\\}_{k=0}^{p}$ and
$\\{u_{p^{\prime}}\\}_{p^{\prime}=0}^{p}$ are complete finite dimensional
spaces, $\hat{X}_{p}$ plays the role of a change-of-basis matrix, allowing for
the determination of $c_{pk}$ via

$\displaystyle\vec{c}_{p}\;=\;\hat{X}_{p}^{-1}\,\vec{b}_{p}\;.$ (49)

In Eq. (49), $\hat{X}_{p}$ and $\vec{b}_{p}$ are implicitly functions of all
$\vec{c}_{p^{\prime}}$, where $p^{\prime}<p$. Thus, given a solution for
$\vec{c}_{p}$, one may then evaluate the right hand side of Eq. (49) to
compute $\vec{c}_{p+1}$ which in turn allows for the evaluation of the
canonical polynomial, $|{u_{p+1}}\rangle$ via Eq. (40). In this way, the
general determination of canonical confluent Heun polynomials follows
inductively from Eqs. (40) and (41).

Having constructed the polynomials, $u_{p}$, to be orthogonal, two scaling
freedoms remain: signature and normalization.

Since the scalar product (Eq. 14) is unchanged with respect to negation of
both of its inputs (e.g.
${\langle{u_{p}}\,|\,{u_{p}}\rangle}={\langle{-u_{p}}\,|\,{-u_{p}}\rangle}$),
there is an overall $\pm 1$ scaling, or rather “signature” freedom. Here, this
is fixed by dividing each $u_{p}$ by its value at $\xi=1$,
$\displaystyle{u_{p}}$ $\displaystyle\leftarrow{u_{p}}/u_{p}(1)\;.$ (50a) The
final freedom is that of normalization. Here, this is fixed in the usual way,
$\displaystyle{u_{p}}$
$\displaystyle\leftarrow{u_{p}}/\sqrt{{\langle{u_{p}}\,|\,{u_{p}}\rangle}}\;.$
(50b)

### Gram-Schmidt Construction

– The well-known Gram-Schmidt process acts on a linearly independent sequence,
producing a new sequence that is orthogonal with respect to a chosen scalar
product Axler (2015). Given a predefined order for the original sequence, the
resulting orthogonal sequence is known to be unique, up to normalization
choices. In this sense, the familiar classical polynomials (e.g. those of
Hermite, Jacobi and Laguerre) are uniquely defined by applying Gram-Schmidt to
the monomials, e.g. $\xi^{p}$, such that the resulting polynomial sequence is
orthogonal with respect to the relevant inner product Abramowitz and Stegun
(1964).

Similarly, the canonical confluent Heun polynomials may be defined by applying
Gram-Schmidt to the monomials using the weight function, Eq. (16).
Orthogonality is guaranteed because the algorithm iteratively removes the non-
orthogonal components between different elements: In particular, given a
sequence of monomials $\\{\xi^{p}\\}_{p=0}^{\infty}$, the algorithm
iteratively defines canonical polynomials, $u_{p}(\xi)$, to be $\xi^{p}$ plus
a lower order polynomial that is constrained by requiring that
${\langle{u_{p}}\,|\,{u_{p^{\prime}}}\rangle}=0$ for all $p^{\prime}<p$. The
results is a sequence of orthogonal polynomials that can then be normalized
according to Eq. (50).

It happens that that the weight function used to define the scalar product in
Eq. (14) is of the Pollaczek-Jacobi type, but with complex valued parameters
(Eqs. 12 and 17). Since our weight function may generally be evaluated using
analytic continuation (see Eq. III), the Gram-Schmidt polynomials described
above may be thought of as analytically continued Pollaczek-Jacobi polynomials
Chen and Dai (2010). Further connections between the polynomials described
above, and the Pollaczek-Jacobi polynomials will be explored in future work
Gurevich and London (2024).

### Equivalence and Completeness

– Both direct and Gram-Schmidt constructions of the canonical confluent Heun
polynomials are unique iterative processes with identical initial conditions,
namely that $|{u_{0}}\rangle$ is a constant. In each construction, the
requirement of orthogonality, i.e.
${\langle{u_{p^{\prime}<p}}\,|\,{u_{p}}\rangle}=0$, uniquely defines
$|{u_{1}}\rangle$, as well as subsequent polynomials. Since, upon each
iteration, there is only one non-trivial result for each $|{u_{p}}\rangle$,
both processes result in normalized canonical confluent Heun polynomials that
are equivalent up to a root-of-unity scaling factor.

Equivalence is easily checked by e.g. writing $|{u_{0}}\rangle$ and
$|{u_{1}}\rangle$ explicitly as polynomials in $\xi$, and then verifying that
both constructions yield the same result, which is

$\displaystyle u_{0}(\xi)\;$
$\displaystyle=\;{\langle\xi^{0}\rangle}^{-1/2}\;,$ (51a) $\displaystyle
u_{1}(\xi)\;$
$\displaystyle=\;\left[{\frac{{\langle\xi^{0}\rangle}}{{\langle\xi^{0}\rangle\langle\xi^{2}\rangle-\langle\xi^{1}\rangle^{2}}}}\right]^{1/2}\;\left(\frac{\langle\xi^{1}\rangle}{\langle\xi^{0}\rangle}+\xi\right)\;.$
(51b)

In Eq. (51), we recall that $\langle\xi^{k}\rangle$ are monomial moments
defined by Eq. (III).

---
Figure 2:  The first five canonical confluent Heun polynomials, i.e.
$u_{p}(\xi)$ with $p\in\\{0,1,2,3,4\\}$. Polynomials are constructed using
field spin weight $s=-2$, BH spin $a/M=0.7$, and Kerr QNM frequency
$M\tilde{\omega}_{220}=0.5326-0.0808i$. For ease of presentation, polynomials
have been divided by their value at $\xi=1$ (i.e. spatial infinity). Top
panel: Real part of each polynomial. Bottom panel: imaginary part of each
polynomial. Related discussion is provided in Sec. VII.

On arrives at Eqs. (51a) and (51b) regardless of whether direct or Gram-
Schmidt construction is used.

Henceforth, we will make no distinction between the canonical confluent Heun
polynomials generated by direct or Gram-Schmidt construction.

Regarding completeness, since the results of direct and Gram-Schmidt
construction are equivalent, if the result of one construction is complete,
then so is the other. The completeness of the Gram-Schmidt constructed
polynomials is perhaps the most straightforward. Any sequence of Gram-Schmidt
polynomials, if constructed from a complete set, will also be complete. The
monomials are complete when all monomial moments are well-defined444See Sec.
VI of Paper I for a detailed discussion. according to Eq. (III) London
(2023a); Adkins and Weintraub (2012); Axler (2015). For example, in Paper I,
it was found that the monomial moments are well-defined for all QNM
frequencies with non-zero real parts, which is the case for most QNM
frequencies of astrophysical relevance. Therefore, we expect that the
canonical confluent Heun polynomials will generally be complete in that
setting.

In practice, since the confluent Heun polynomials are eigenfunctions to
$\mathcal{D}{{}_{\xi}}$, it is useful to represent each $|{u_{p}}\rangle$ in
terms of confluent Heun polynomials, $|{y_{p^{\prime}k^{\prime}}}\rangle$.
However, since the confluent Heun polynomials are overcomplete, there is no
unique representation.

### Single and Multi-order representations

– The overcompleteness of the confluent Heun polynomials means that there are
multiple equivalent representations of any integrable function in terms of the
confluent Heun polynomials, $|{y_{p^{\prime}k}}\rangle$. Here, this fact is
applied to the canonical confluent Heun polynomials.

Since each space of order $p$ confluent Heun polynomials is a $p+1$
dimensional basis, the identity operator for all order $p$ polynomials,
${\mathbb{I}}_{p}$, may be written as it was in Eq. (26),

$\displaystyle{\mathbb{I}}_{p}\;=\;\sum_{k=0}^{p}\,|{y_{pk}}\rangle\langle{y_{pk}}|\;.$
(52)

Having defined the canonical confluent Heun polynomials, $|{u_{p}}\rangle$,
two mathematically different but practically equivalent representations result
from Eq. (52).

The first is simply the ansatz used for direct construction, i.e. Eq. (40),

$\displaystyle|{u_{p}}\rangle\;$
$\displaystyle=\;{\mathbb{I}}_{p}|{u_{p}}\rangle\;$ (53a)
$\displaystyle=\;\sum_{k=0}^{p}|{{y}_{pk}}\rangle{\langle{y_{pk}}\,|\,{u_{p}}\rangle}\;$
(53b) $\displaystyle=\;\sum_{k=0}^{p}c_{pk}\,|{{y}_{pk}}\rangle\;.$ (53c)

In Eqs. (53a) and (53b), the order-$p$ identity operator, ${\mathbb{I}}_{p}$,
has been applied to $|{u_{p}}\rangle$. In Eq. (53c), it is noted that
$c_{pk}={\langle{y_{pk}}\,|\,{u_{p}}\rangle}$. The direct construction
algorithm discussed previously is simply a way to determine
${\langle{y_{pk}}\,|\,{u_{p}}\rangle}$ when $|{u_{p}}\rangle$ is unknown.

Within Eq. (53), each canonical confluent Heun polynomial is represented in
the basis of a different polynomial order, meaning that a set of multiple
canonical confluent Heun polynomials would require the generation of confluent
Heun polynomials from multiple orders. For this reason, Eq. (53) will
henceforth be referred to as the multi-order representation of the canonical
confluent Heun polynomials.

For any finite sequences of canonical confluent Heun polynomials,
$\\{|{u_{p}}\rangle\\}_{p=0}^{p^{\prime}}$, the multi-order representation
requires the evaluation of all polynomials within the set
$\bigcup_{p=0}^{p^{\prime}}\\{|{y_{pk}}\rangle\\}_{k=0}^{p}$. Since there will
be $\frac{1}{2}(p^{\prime}+1)(p^{\prime}+2)$ such polynomials, their
evaluation presents a significant, but ultimately avoidable, computational
expense.

A second representation allowed by Eq. (52) entails choosing a single
confluent Heun polynomial basis of order $p^{\prime}$ to represent all
canonical confluent Heun polynomials with $p\leq p^{\prime}$,

$\displaystyle|{u_{p}}\rangle\;$
$\displaystyle=\;{\mathbb{I}}_{p^{\prime}}|{u_{p}}\rangle\;$ (54a)
$\displaystyle=\;\sum_{k=0}^{p^{\prime}}|{{y}_{p^{\prime}k}}\rangle{\langle{y_{p^{\prime}k}}\,|\,{u_{p}}\rangle}\;$
(54b)
$\displaystyle=\;\sum_{k=0}^{p^{\prime}}d_{p^{\prime}pk}\,|{{y}_{p^{\prime}k}}\rangle\;.$
(54c)

In Eqs. (54a) and (54b), the identity operator for the order-$p^{\prime}$
confluent Heun polynomials has been applied to $|{u_{p}}\rangle$, and in Eq.
(54c), $d_{p^{\prime}pk}={\langle{y_{p^{\prime}k}}\,|\,{u_{p}}\rangle}$.

$\displaystyle
d_{p^{\prime}pk}\;=\;{\langle{y_{p^{\prime}k}}\,|\,{u_{p}}\rangle}$ (55)

In practice, one may first generate $|{u_{p}}\rangle$ via Gram-Schmidt
construction, and then explicitly compute
${\langle{y_{p^{\prime}k}}\,|\,{u_{p}}\rangle}$ using Eq. (19). Equation (54)
will henceforth be referred to as the single-order representation of the
canonical confluent Heun polynomials. If the largest polynomial order of
interest is $p^{\prime}$, then the single-order representation has the
practical advantage of only requiring $p^{\prime}+1$ confluent Heun
polynomials to be computed after Gram-Schmidt construction.

##  VI Tridiagonalization

We may now return to Teukolsky’s radial equation for QNMs,

$\displaystyle\mathcal{L}{{}_{\xi}}\,|{f}\rangle\;$
$\displaystyle=\;A\,|{f}\rangle\;.$ (56)

Having determined canonical polynomials that are natural to
$\mathcal{L}{{}_{\xi}}$, their completeness and orthogonality may be used to
represent Eq. (56) as a simple linear matrix equation.

Since the canonical confluent Heun polynomials are complete, we may use them
to represent the identity operator,

$\displaystyle{\mathbb{I}}\;$
$\displaystyle=\;\sum_{p=0}^{\infty}\,|{u_{p}}\rangle\langle{u_{p}}|\;.$ (57)

This allows each operator and vector within Eq. (56) to be projected into the
basis of canonical confluent Heun polynomials,

$\displaystyle{\mathbb{I}}\,\mathcal{L}{{}_{\xi}}\,{\mathbb{I}}\,|{f}\rangle\;$
$\displaystyle=\;A\,{\mathbb{I}}\,|{f}\rangle\;.$ (58)

Using Eq. (57) to expand Eq. (58) yields

$\displaystyle\sum_{p,q}\,|{u_{p}}\rangle$
$\displaystyle\langle{u_{p}}|\,\mathcal{L}{{}_{\xi}}\,|{u_{q}}\rangle{\langle{u_{q}}\,|\,{f}\rangle}\;\;$
(59a)
$\displaystyle=\;A\,\sum_{p}\,|{u_{p}}\rangle{\langle{u_{p}}\,|\,{f}\rangle}\;,$
$\displaystyle\hat{L}\,\vec{f}\;$ $\displaystyle=\;A\,\vec{f}\;.$ (59b)

In going from Eq. (59a) to Eq. (59b), we recognize $\hat{L}$ to be the matrix
whose elements are

$\displaystyle
L_{pq}\;=\;{\langle{u_{p}}\,|\,{\mathcal{L}{{}_{\xi}}|\,u_{q}}\rangle}\;,$
(60)

and $\vec{f}$ to be the vector whose elements are

$\displaystyle f_{p}\;=\;{\langle{u_{p}}\,|\,{f}\rangle}\;.$ (61)

Recall that since $\mathcal{L}{{}_{\xi}}$ is self-adjoint with respect to the
scalar product, the following equivalences hold:
${\langle{u_{p}}\,|\,{\mathcal{L}{{}_{\xi}}|\,u_{q}}\rangle}\;=\;{\langle{u_{p}}\,|\,{\mathcal{L}{{}_{\xi}}\,u_{q}}\rangle}\;=\;{\langle{\mathcal{L}{{}_{\xi}}\,u_{p}}\,|\,{u_{q}}\rangle}$.

In this section, two mathematically different but practically equivalent
routes to tridiagonalizing $\hat{L}$ are presented. The first route is based
upon the single-order representation of the canonical confluent Heun
polynomials, and it will be discussed in detail. The second route is based
upon their multi-order representation, and it will be described summarily as
the step-by-step reasoning is both similar to that of the single-order
approach, and already detailed in Paper I.

Regardless of the representation used, the result is the same: the confluent
Heun operator $\mathcal{L}_{\xi}$ is approximated by a matrix, $\hat{L}$, that
is of finite size, symmetric and tridiagonal. In effect, the complex Jacobi
matrix of $\mathcal{L}_{\xi}$ is determined Teschl (2000); Kuijlaars  _et al._
(2003).

Our present concern is the structure of $\hat{L}$. Since its matrix elements
are ${\langle{u_{p}}\,|\,{\mathcal{L}{{}_{\xi}}\,u_{q}}\rangle}$, the action
of $\mathcal{L}{{}_{\xi}}$ on $|{u_{q}}\rangle$ must first be understood,

$\displaystyle\mathcal{L}{{}_{\xi}}\,|{u_{q}}\rangle\;=\;(\text{C}_{0}+\text{C}_{1}[1-\xi])\,$
$\displaystyle|{u_{q}}\rangle$ (62a)
$\displaystyle\;+\quad\quad\quad\;\mathcal{D}_{\xi}\,$
$\displaystyle|{u_{q}}\rangle\;.$ (62b)

In Eq. (62), as in Eq. (11), $\mathcal{L}_{\xi}$ has been written in terms of
the physical constants $\text{C}_{0}$ and $\text{C}_{1}$ (Eq. 62a), as well as
the confluent Heun operator $\mathcal{D}_{\xi}$ (Eq. 62b). Plainly, Eq. (62a)
only contains constant rescalings of $|{u_{p}}\rangle$ and terms dependent on
$\xi|{u_{q}}\rangle$, and Eq. (62b) is simply
$\mathcal{D}_{\xi}\,|{u_{q}}\rangle$. Since
${\langle{u_{p}}\,|\,{u_{q}}\rangle}=\delta_{pq}$, constant rescalings of
$|{u_{p}}\rangle$ will only contribute to the diagonal of $\hat{L}$.
Therefore, to understand the structure of $\hat{L}$, we must first refine our
understanding of $\xi|{u_{q}}\rangle$ and
$\mathcal{D}_{\xi}\,|{u_{q}}\rangle$.

It happens that $\xi|{u_{q}}\rangle$ can be represented in terms of only three
canonical polynomials,

$\displaystyle\xi|{u_{q}}\rangle=\sigma_{q0}|{u_{q-1}}\rangle+\sigma_{q1}|{u_{q}}\rangle+\sigma_{q2}|{u_{q+1}}\rangle\;.$
(63)

Equation (63) is a three-term recursion relationship between canonical
confluent Heun polynomials. It is a standard result in polynomial theory, and
can be understood as resulting from completeness, orthogonality, and symmetry
of the scalar product, namely ${\langle{u_{p}}\,|\,{\xi
u_{q}}\rangle}={\langle{\xi u_{p}}\,|\,{u_{q}}\rangle}$.

In Eq. (63), the recursion coefficients, $\sigma_{q0}$, $\sigma_{q1}$ and
$\sigma_{q2}$ are $p$-dependent complex numbers. Unlike many instances in
polynomial theory where the recursion coefficients are known analytically,
here they are to be treated numerically, as determined by evaluation of the
scalar product:

$\displaystyle\sigma_{q0}\;$
$\displaystyle=\;{\langle{u_{q-1}}\,|\,{\xi\,u_{q}}\rangle}\;,$ (64a)
$\displaystyle\sigma_{q1}\;$
$\displaystyle=\;{\langle{u_{q}}\,|\,{\xi\,u_{q}}\rangle}\;,$ (64b)
$\displaystyle\sigma_{q2}\;$
$\displaystyle=\;{\langle{u_{q+1}}\,|\,{\xi\,u_{q}}\rangle}\;.$ (64c)

In the context of $\mathcal{L}{{}_{\xi}}\,|{u_{q}}\rangle$ (i.e. Eq. 62), the
existence of a three-term recursion relationship between polynomials means
that terms within Eq. (62a) only contribute to the central, lower and upper
diagonal of $\hat{L}$. Therefore, Eq. (62a) will generally contribute to the
tridiagonal part of $\hat{L}$.

This should be expected of Eq. (62a), since all canonical type polynomials
(e.g. Chebyshev and Laguerre) have three-term recursions. However, since not
all canonical polynomials are natural to $\mathcal{D}_{\xi}$, it happens that
$\mathcal{D}_{\xi}|{u_{q}}\rangle$ is a problem-specific quantity.

Further, since $\mathcal{D}_{\xi}$ contains derivatives, it is here that the
representation of each $|{u_{q}}\rangle$ has meaningful effect: while the
previously discussed single- and multi-order representations are equivalent
(Eqs. 53 and 54), their behaviors under differentiation are not.

### The single-order approach

– We may first consider the single-order representation of the canonical
confluent Heun polynomials. In that picture, any polynomial of order $q$ is
expanded in terms of confluent Heun polynomials of order $p^{\prime}\geq q$
(Eq. 56). For a canonical confluent Heun polynomial $|{u_{q}}\rangle$, this
expansion is written as

$\displaystyle|{u_{q}}\rangle\;$
$\displaystyle=\;\sum_{k=0}^{p^{\prime}}d_{p^{\prime}qk}\,|{{y}_{p^{\prime}k}}\rangle\;.$
(65)

Therefore the action of $\mathcal{D}{{}_{\xi}}$ on $|{u_{q}}\rangle$ is given
by,

$\displaystyle\mathcal{D}{{}_{\xi}}\,|{u_{q}}\rangle\;$
$\displaystyle=\;\sum_{k=0}^{p^{\prime}}d_{p^{\prime}qk}\,\mathcal{D}{{}_{\xi}}\,|{{y}_{p^{\prime}k}}\rangle\;$
(66a)
$\displaystyle=\;\sum_{k=0}^{p^{\prime}}\,d_{p^{\prime}qk}\,(\lambda_{p^{\prime}k}\,+\,\mu_{p^{\prime}}\,\xi)\,|{{y}_{p^{\prime}k}}\rangle\;,$
(66b)
$\displaystyle=\;|{v_{p^{\prime}q}}\rangle\,+\,\mu_{p^{\prime}}\,\xi\,|{u_{q}}\rangle\;.$
(66c)

In Eq. (66b), the two parameter eigenvalue relation for the confluent Heun
polynomials, Eq. (22a), has been applied. In Eq. (66c), terms that depend on
$\lambda_{p^{\prime}k}$ have been grouped to define a new polynomial,

$\displaystyle|{v_{p^{\prime}q}}\rangle\;=\;\sum_{k=0}^{p^{\prime}}\,{d}_{p^{\prime}qk}\,\lambda_{p^{\prime}k}\,|{{y}_{p^{\prime}k}}\rangle\;,$
(67)

and it has been noticed that remaining terms in Eq. (66b) are simply
proportional to $|{u_{q}}\rangle$. Since $\mathcal{D}_{\xi}$ increases the
order of polynomials by one (see Sec. III), $|{v_{p^{\prime}q}}\rangle$ has
order $q+1$.

It follows (from Eqs. 62 and 66c) that $\mathcal{L}_{\xi}\,|{u_{q}}\rangle$
may be written as

$\displaystyle\mathcal{L}{{}_{\xi}}|{u_{q}}\rangle\;=$
$\displaystyle\;(\text{C}_{0}+\text{C}_{1})|{u_{q}}\rangle$ (68a)
$\displaystyle\;+\;|{v_{p^{\prime}q}}\rangle$ (68b)
$\displaystyle\;+\;\alpha_{p^{\prime}}\,\xi\,|{u_{q}}\rangle.$ (68c)

In Eq. (68c), recall that $\alpha_{p^{\prime}}$ is defined in Eq. (36). The
structure of Eq. (68) enables a few useful observations. Equations (68a-68c)
imply that $\hat{L}$ is the sum of three matrices, $\hat{L}^{(\text{S})}_{0}$,
$\hat{L}^{(\text{S})}_{1}$, and $\hat{L}^{(\text{S})}_{2}$,

$\displaystyle\hat{L}\;=\;\hat{L}^{(\text{S})}_{0}+\hat{L}^{(\text{S})}_{1}+\hat{L}^{(\text{S})}_{2}\;.$
(69)

Their respective matrix elements are

$\displaystyle{L}_{0,pq}^{(\text{S})}\;$
$\displaystyle=\;(\text{C}_{0}+\text{C}_{1})\,\delta_{pq}\;,$ (70a)
$\displaystyle{L}_{1,pq}^{(\text{S})}\;$
$\displaystyle=\;{\langle{u_{p}}\,|\,{v_{p^{\prime}q}}\rangle}\;,$ (70b)
$\displaystyle{L}_{2,pq}^{(\text{S})}\;$
$\displaystyle=\;\alpha_{p^{\prime}}{\langle{u_{p}}\,|\,{\xi\,u_{q}}\rangle}\;.$
(70c)

Since $\mathcal{L}_{\xi}$ is self-adjoint with respect to the scalar product,
$\hat{L}$ will be symmetric, meaning that $L_{pq}=L_{qp}$. It may be observed
that $\hat{L}^{(\text{S})}_{0}$ and $\hat{L}^{(\text{S})}_{2}$ are symmetric,
owing to the symmetry of the Kronecker delta, and of the scalar product. Thus
for
$\hat{L}^{(\text{S})}_{0}+\hat{L}^{(\text{S})}_{1}+\hat{L}^{(\text{S})}_{2}$
to be symmetric, $\hat{L}^{(\text{S})}_{1}$ must be symmetric, meaning that

$\displaystyle{\langle{u_{p}}\,|\,{v_{p^{\prime}q}}\rangle}\;=\;{\langle{v_{p^{\prime}p}}\,|\,{u_{q}}\rangle}\;.$
(71)

This conclusion may also be reached via the consideration of
${\langle{u_{p}}\,|\,{\mathcal{D}{{}_{\xi}}\,u_{q}}\rangle}\;=\;{\langle{\mathcal{D}{{}_{\xi}}\,u_{p}}\,|\,{u_{q}}\rangle}$,
which directly reduces to Eq. (71). Thus, Eq. (71) is only true for scalar
products under which $\mathcal{D}_{\xi}$ is self-adjoint (e.g. Eq. 71 does not
hold for Chebyshev inner-product, meaning that the associated
$\hat{L}^{(\text{S})}_{1}$ matrix would not be symmetric) London (2023a).

Lastly, it may be deduced that $\hat{L}^{(\text{S})}_{1}$ is tridiagonal. This
follows from the fact that each $|{v_{p^{\prime}q}}\rangle$ is of order $q+1$.
In particular, given any order $q+1$ polynomial, e.g.
$|{v_{p^{\prime}q}}\rangle$, its canonical confluent Heun polynomial
representation will not contain contributions from $|{u_{p}}\rangle$ if
$p>q+1$,

$\displaystyle|{v_{p^{\prime}q}}\rangle\;=$
$\displaystyle\;{\mathbb{I}}\,|{v_{p^{\prime}q}}\rangle$ (72a)
$\displaystyle=$
$\displaystyle\;\sum_{p=0}^{\infty}\,|{u_{p}}\rangle{\langle{u_{p}}\,|\,{v_{p^{\prime}q}}\rangle}\;,$
(72b) $\displaystyle=$
$\displaystyle\;\sum_{p=0}^{q+1}\,|{u_{p}}\rangle{\langle{u_{p}}\,|\,{v_{p^{\prime}q}}\rangle}\;.$
(72c)

In Eqs. (72a-72b), a canonical confluent Heun polynomial representation is
applied to $|{v_{p^{\prime}q}}\rangle$. In Eq. (72c), it has been acknowledged
that terms with $p>q+1$ do not contribute; otherwise, the resulting polynomial
would be of order greater than $q+1$. This condition on
${\langle{u_{p}}\,|\,{v_{p^{\prime}q}}\rangle}$, as well as its counterpart on
${\langle{u_{q}}\,|\,{v_{p^{\prime}p}}\rangle}$ may be stated thusly,

$\displaystyle{\langle{u_{p}}\,|\,{v_{p^{\prime}q}}\rangle}\;=\;0\text{, if
}\;\;p>q+1\;,$ (73a)
$\displaystyle{\langle{u_{q}}\,|\,{v_{p^{\prime}p}}\rangle}\;=\;0\text{, if
}\;\;q>p+1\;.$ (73b)

Note that
${\langle{u_{q}}\,|\,{v_{p^{\prime}p}}\rangle}={\langle{v_{p^{\prime}p}}\,|\,{u_{q}}\rangle}$,
not due to Eq. (71), but rather to the multiplicative associative property
applied to the integrand of the scalar product (Eq. 14). Also note that the
conditional in Eq. (73b) is equivalent to $p<q-1$.

Since Eqs. (73a) and (73b) hold simultaneously, it follows that
${\langle{u_{p}}\,|\,{v_{p^{\prime}q}}\rangle}\neq 0$ when $q-1\leq p\leq
q+1$. In other words, each $|{v_{p^{\prime}q}}\rangle$ is exactly equal to a
sum over only three canonical confluent Heun polynomials,

$\displaystyle|{v_{p^{\prime}q}}\rangle\;=\;\nu_{q0}|{u_{q-1}}\rangle+\nu_{q1}|{u_{q}}\rangle+\nu_{q2}|{u_{q+1}}\rangle\;.$
(74)

In Eq. (74), $\nu_{q0}$, $\nu_{q1}$ and $\nu_{q2}$ are complex valued
$q$-dependent numbers defined by the scalar product,

$\displaystyle\nu_{q0}\;$
$\displaystyle=\;{\langle{u_{q-1}}\,|\,{v_{p^{\prime}q}}\rangle}\;,$ (75a)
$\displaystyle\nu_{q1}\;$
$\displaystyle=\;{\langle{u_{q}}\,|\,{v_{p^{\prime}q}}\rangle}\;,$ (75b)
$\displaystyle\nu_{q2}\;$
$\displaystyle=\;{\langle{u_{q+1}}\,|\,{v_{p^{\prime}q}}\rangle}\;.$ (75c)

The elements of $\hat{L}^{(\text{S})}_{1}$, namely
${\langle{u_{p}}\,|\,{v_{p^{\prime}q}}\rangle}$, are only nonzero between its
upper and lower diagonal, so $\hat{L}^{(\text{S})}_{1}$ is a tridiagonal
matrix.

Since $\hat{L}^{(\text{S})}_{1}$ and $\hat{L}^{(\text{S})}_{2}$ are symmetric
and tridiagonal, and $\hat{L}^{(\text{S})}_{0}$ is diagonal, if follows that
$\hat{L}=\hat{L}^{(\text{S})}_{0}+\hat{L}^{(\text{S})}_{1}+\hat{L}^{(\text{S})}_{2}$
is symmetric and tridiagonal. Equivalently,

$\displaystyle\mathcal{L}_{\xi}|{u_{q}}\rangle\;=\;\nu_{q3}\,|{u_{q-1}}\rangle+\nu_{q4}\,|{u_{q}}\rangle+\nu_{q5}\,|{u_{q+1}}\rangle\;,$
(76)

where,

$\displaystyle\nu_{q3}\;$
$\displaystyle=\;\alpha_{p^{\prime}}\,\sigma_{q0}+\nu_{q0}\;,$ (77a)
$\displaystyle\nu_{q4}\;$
$\displaystyle=\;\alpha_{p^{\prime}}\,\sigma_{q1}+\nu_{q1}+(\text{C}_{0}+\text{C}_{1})\;,$
(77b) $\displaystyle\nu_{q5}\;$
$\displaystyle=\;\alpha_{p^{\prime}}\,\sigma_{q2}+\nu_{q2}\;.$ (77c)

In Eq. (77), we recall from Eq. (36) that
$\alpha_{p^{\prime}}=\mu_{p^{\prime}}-\text{C}_{1}$, and that $\sigma_{q0}$,
$\sigma_{q1}$ and $\sigma_{q2}$ are defined by Eq. (64).

### The multi-order approach

– Here, an alternative approach to tridiagonalization is summarized, following
and building upon Sec. VIII of Paper I. This approach uses the multi-order
representation of the canonical confluent Heun polynomials,

$\displaystyle|{u_{q}}\rangle\;$
$\displaystyle=\;\sum_{k=0}^{q}c_{qk}\,|{{y}_{qk}}\rangle\;.$ (78)

While the net result, namely $\hat{L}$, is the same as that found previously,
the multi-order approach enables the derivation of relationships between
$\hat{L}$ and the canonical confluent Heun polynomials’ three-term recursion
that are ostensibly hidden within the single-order approach. This point is
illustrated by first applying the multi-order representation of the canonical
confluent Heun polynomials to the radial problem, and then using the symmetry
of the Jacobi matrix $\hat{L}$ to relate elements of $\hat{L}$ to the
canonical confluent Heun polynomials’ three-term recursion. This proceeds as
follows.

Under the multi-order approach, $\hat{L}$ may be written as the sum of three
matrices,

$\displaystyle\hat{L}\;=\;\hat{L}^{(\text{M})}_{0}+\hat{L}^{(\text{M})}_{1}+\hat{L}^{(\text{M})}_{2}\;.$
(79)

Their respective matrix elements are

$\displaystyle{L}_{0,pq}^{(\text{M})}\;$
$\displaystyle=\;(\text{C}_{0}+\text{C}_{1})\,\delta_{pq}\;,$ (80a)
$\displaystyle{L}_{1,pq}^{(\text{M})}\;$
$\displaystyle=\;{\langle{u_{p}}\,|\,{v_{q}}\rangle}\;,$ (80b)
$\displaystyle{L}_{2,pq}^{(\text{M})}\;$
$\displaystyle=\;\alpha_{q}\,{\langle{u_{p}}\,|\,{\xi\,u_{q}}\rangle}\;.$
(80c)

In Eq. (80b), $|{v_{q}}\rangle$ is an order-$q$ “auxiliary” polynomial defined
as

$\displaystyle|{v_{q}}\rangle\;=\;\sum_{k=0}^{q}\,\text{c}_{qk}\,\lambda_{qk}\,|{{y}_{qk}}\rangle\;,$
(81)

where it is recalled that $\lambda_{qk}$ is the eigenvalue of the confluent
Heun polynomial, $|{y_{qk}}\rangle$. Consideration of
${\langle{u_{p}}\,|\,{\mathcal{D}{{}_{\xi}}\,u_{q}}\rangle}\;=\;{\langle{\mathcal{D}{{}_{\xi}}\,u_{p}}\,|\,{u_{q}}\rangle}$
yields a that $|{v_{q}}\rangle$ may be equated with a sum over only a few
canonical confluent Heun polynomials,

$\displaystyle|{v_{q}}\rangle=\sigma_{q3}|{u_{q-1}}\rangle+\sigma_{q4}|{u_{q}}\rangle\;.$
(82)

In Eq. (82), $\sigma_{q3}$ and $\sigma_{q4}$ are complex valued $q$-dependent
numbers defined by the scalar product,

$\displaystyle\sigma_{q3}\;$
$\displaystyle=\;{\langle{u_{q-1}}\,|\,{v_{q}}\rangle}\;,$ (83a)
$\displaystyle\sigma_{q4}\;$
$\displaystyle=\;{\langle{u_{q}}\,|\,{v_{q}}\rangle}\;.$ (83b)

It may be observed that $\hat{L}^{(\text{M})}_{0}$ is diagonal,
$\hat{L}^{(\text{M})}_{1}$ upper-tridiagonal (due to Eq. 82), and
$\hat{L}^{(\text{M})}_{2}$ is tridiagonal. Therefore, $\hat{L}$ is tridiagonal
as expected.

It is informative to note that, while $\hat{L}^{(\text{M})}_{0}$ is symmetric,
Eqs. (80b) and (80c) are asymmetric: Eq. (80b) is asymmetric due to Eq. (82)
not having a contribution from $|{u_{q+1}}\rangle$ (i.e.
$\hat{L}^{(\text{M})}_{1}$ is only upper-tridiagonal), and Eq. (80c) is
asymmetric due to the $q$-dependent scale factor, $\alpha_{q}$. Thus, since
$\hat{L}$ is symmetric, $\hat{L}^{(\text{M})}_{1}$ and
$\hat{L}^{(\text{M})}_{2}$ must combine to form a quantity that is also
symmetric. In turn, this means that $\hat{L}^{(\text{M})}_{1}$ and
$\hat{L}^{(\text{M})}_{2}$ are not independent objects.

The link between $\hat{L}^{(\text{M})}_{1}$ and $\hat{L}^{(\text{M})}_{2}$ is
perhaps easiest to understand by noting that the multi-order representation of
the canonical confluent Heun polynomials allows
$\mathcal{L}_{\xi}|{u_{q}}\rangle$ to be written as,

$\displaystyle\mathcal{L}_{\xi}|{u_{q}}\rangle\;=\;\sigma_{q5}\,|{u_{q-1}}\rangle+\sigma_{q6}\,|{u_{q}}\rangle+\sigma_{q7}\,|{u_{q+1}}\rangle\;,$
(84)

where

$\displaystyle\sigma_{q5}\;$
$\displaystyle=\;\alpha_{q}\,\sigma_{q0}+\sigma_{q3}\;,$ (85a)
$\displaystyle\sigma_{q6}\;$
$\displaystyle=\;\alpha_{q}\,\sigma_{q1}+\sigma_{q4}+(\text{C}_{0}+\text{C}_{1})\;,$
(85b) $\displaystyle\sigma_{q7}\;$
$\displaystyle=\;\alpha_{q}\,\sigma_{q2}\;.$ (85c)

Plainly stated, since $\hat{L}$ is a symmetric matrix its matrix elements,
$L_{pq}$, are such that $L_{pq}=L_{qp}$. Equation (85) means that $L_{pq}$
will be within the set $\\{\sigma_{q5},\sigma_{q6},\sigma_{q7}\\}$ for
$p\in\\{q-1,q,q+1\\}$, respectively. Since $\hat{L}$ is tridiagonal, the only
nontrivial symmetry statement is $L_{q-1,q}=L_{q,q-1}$, which is equivalent to
$L_{q,q+1}=L_{q+1,q}$ under $q\rightarrow q+1$. Using equation Eqs. (84) and
(85), $L_{q-1,q}=L_{q,q-1}$ becomes,

$\displaystyle\sigma_{q5}\;$ $\displaystyle=\;\sigma_{q-1,7}\;,$ (86a)
$\displaystyle\alpha_{q}\,\sigma_{q0}+\sigma_{q3}\;$
$\displaystyle=\;\;\alpha_{q-1}\,\sigma_{q-1,2}\;.$ (86b)

Figure 3:  Matrices related to canonical confluent Heun polynomials, and
representations of Teukolsky’s radial operator for $(s,\ell,m,n)=(-2,2,2,6)$,
$a/M=0.7$, and $M\tilde{\omega}_{{\ell mn}}=0.4239-1.0954i$. All color bars
correspond to the $\log_{10}$ of the absolute values of quantities noted in
annotation, e.g. for panel (a),
$\log_{10}|{\langle{u_{p}}\,|\,{u_{q}}\rangle}|$ is visualized. Tiles (in
purple) away from each panel’s tridiagonal band are exactly zero. See Sec. VII
for related discussion.

In Eq. (86a), $L_{q-1,q}=L_{q,q-1}$ has been rewritten using Eq. (84). In
going from Eq. (86a) to Eq. (86b), Eq. (85) has been used to expand both
sides. Equation (86b) may be further reduced by noting that symmetry of the
canonical polynomials’ three-term recursion (See Eq. 63) means that
$\sigma_{q-1,2}=\sigma_{q0}$. Therefore, Eq. (86b) allows $\sigma_{q3}$ (i.e.
${\langle{u_{q-1}}\,|\,{v_{q}}\rangle}$) to be written as,

$\displaystyle\sigma_{q3}\;$
$\displaystyle=\;(\alpha_{q-1}-\alpha_{q})\,\sigma_{q0}\;,$ (87a)
$\displaystyle{\langle{u_{q-1}}\,|\,{v_{q}}\rangle}\;$
$\displaystyle=\;(2(q-1)+\text{C}_{4})\,{\langle{u_{q-1}}\,|\,{\xi\,u_{q}}\rangle}\;.$
(87b)

In Eq. (87a), Eq. (86b) as been solved for $\sigma_{q3}$. In Eq. (87b),
$\sigma_{q3}$ and $\sigma_{q0}$ have been replaced with their scalar-product
definitions via Eqs. (62a) and (83a), and Eqs. (22b) and (27) have been used
to rewrite instances of $\alpha_{q}$.

Equation (87b) communicates that the non-zero off-diagonal elements of
$\hat{L}^{(\text{M})}_{1}$, namely ${\langle{u_{q-1}}\,|\,{v_{q}}\rangle}$,
are directly proportional to the three-term-recursion coefficient,
${\langle{u_{q-1}}\,|\,{\xi\,u_{q}}\rangle}$. Recalling Eq. (80c), the same is
clearly true for off-diagonal elements of $\hat{L}^{(\text{M})}_{2}$.
Together, those two observations mean that the off-diagonal elements of
$\hat{L}$ are generally proportional to the three-term-recursion coefficients
$\sigma_{q0}$ and $\sigma_{q2}$.

##  VII Select numerical results

Before being put into practice, the preceding theoretical development benefits
from a few points of guidance London (2023a). For all computations, physical
parameters, namely $\\{a/M,s,\tilde{\omega},m\\}$ must be chosen such that the
derived parameters, ${\text{C}_{1}}$ through ${\text{C}_{4}}$, may then be
computed according to Eq. (12). For nearly all subsequent manipulations
involving polynomials or series expansions, the fundamental operation is the
scalar product (See Eq. 14).

While it is possible to evaluate the scalar product by directly integrating
it, success in that endeavor requires the formulation of a complex path that
avoids singularities in the integrand London (2023a). A much simpler method of
evaluating the scalar product is given by Eqs. (19) and (III), which reduces
the effort to the evaluation of standard special functions via their well-
known analytic continuations.

Given the scalar product, it is subsequently useful to determine the
eigenvectors and values of matrices, such as $\hat{L}$ which encodes the QNM
problem (Eq. 59b). This may be accomplished with any of several standard
linear algebra packages Harris  _et al._ (2020); mpmath development team
(2023). Henceforth, Eq. (19) is used to compute all scalar product, and Eq.
(III) is implemented with $128$ decimal points of general precision arithmetic
via the Python package, mpmath mpmath development team (2023). This package is
also used to solve various matrix eigenvalue problems.

In this section, example applications of the canonical confluent Heun
polynomials are provided along with related practical guidance. While many
cases are considered, most examples are connected by a single fiducial case: a
Kerr BH of dimensionless spin $a/M=0.7$, and its QNM with
$(s,\ell,m,n)=(-2,2,2,6)$. Broadly, examples range from the brief
consideration of the canonical confluent Heun polynomials, to their
application to the confluent Heun polynomial and QNM problems. In the latter
example, the QNMs’ radial functions are explicitly computed and shown to
satisfy Teukolsky’s radial equation (Eq. 13a). This section concludes with a
brief discussion of the transition from Schwarzschild to Kerr, from the
perspective of the QNM radial functions.

### Canonical confluent Heun polynomials

– Examples of the first five canonical confluent Heun polynomials are shown in
Fig. 2.

As seen in the top panel of Fig. 2, the polynomials may be formatted such that
their real parts have zero roots along the real line, meaning that the net
complex valued polynomials will have no roots at real values of the fractional
radial coordinate, $\xi$. This is simply a result of the fact that the
canonical confluent Heun polynomials have complex valued coefficients (Eq.
51), and so will generally have roots situated away from the real line.

### Tridiagonalization

– The canonical confluent Heun polynomials may be used to tridiagonalize
confluent Heun equations such as the QNM’s radial equation (Eq. 13), or the
confluent Heun polynomial equation (Eq. 22).

In Fig. 3 a step-by-step application of the canonical confluent Heun
polynomials to the QNM radial problem is shown. There, the top three panels
(i.e. $(a)$, $(c)$ and $(e)$) show quantities derived directly from the
canonical confluent Heun polynomials. The remaining panels show how the multi-
order representation (Eq. 78) of each quantity contributes to the
representation of the radial differential operator according to Eq. (79).

When written in terms of matrix elements, Eq. (79) is

$\displaystyle
L_{pq}\;=\;{L}^{(\text{M})}_{0,pq}+{L}^{(\text{M})}_{1,pq}+{L}^{(\text{M})}_{2,pq}\;.$
(88)

Under this framing, the qualitative features of panels $(b)$, $(d)$ and $(f)$
are determined by those of panels $(a)$, $(c)$ and $(e)$, respectively.

Conceptually, canonical confluent Heun polynomial orthogonality (Eq. 41) is
demonstrated in panel $(a)$, upper-tridiagonality of auxiliary polynomials
(Eq. 82) is demonstrated in panel $(c)$, and three-term recursion (Eq. 64) is
dempstrated in panel $(e)$. Then end result, panel $(g)$, may be considered to
represent the three-term recursion of orthogonal polynomials which approximate
solutions to the QNM radial problem Teschl (2000); Kuijlaars  _et al._ (2003).

### Polynomial/non-polynomial duality

– The canonical confluent Heun polynomials, $|{u_{p}}\rangle$, may be applied
to the confluent Heun polynomial problem. In doing so, a quantitative example
of polynomial/non-polynomial duality is encountered, along with related
insights into the QNM problem.

The confluent Heun polynomial problem was stated in Eq. (22a) as,

$\displaystyle\mathcal{D}{{}_{\xi}}\;|{y_{pk}}\rangle\;$
$\displaystyle=\;(\lambda_{pk}\;+\;\mu_{p}\,\xi)\;|{y_{pk}}\rangle\;,$ (89a)
To apply the canonical confluent Heun polynomials, Eq. (89a) may be restated
as $\displaystyle\mathcal{L}_{\xi}^{(p)}\;|{y_{pk}}\rangle\;$
$\displaystyle=\;\lambda_{pk}\;|{y_{pk}}\rangle\;,$ (89b) where
$\displaystyle\mathcal{L}_{\xi}^{(p)}\;=\;-\;\mu_{p}\,\xi\;+\;\mathcal{D}{{}_{\xi}}\;.$
(89c)

---
Figure 4:  Example of confluent Heun polynomial/non-polynomial duality for
polynomial order $p=6$, and physical parameters $(s,\ell,m,n)=(-2,2,2,6)$,
$a/M=0.7$, and $M\tilde{\omega}_{{\ell mn}}=0.4239-1.0954i$. All color bars
correspond to the $\log_{10}$ of the absolute values of quantities noted in
annotation. Matrix elements of zeo value are colored purple. Panel $(a)$:
Matrix representation of the radial operator (Eq. 69). Panel $(b)$: Matrix
representation of the radial operator’s eigenvectors (Eq. 90). In panels $(a)$
and $(b)$, vertical and horizontal white lines delineate polynomial and non-
polynomial sectors. Panel $(c)$: The radial operator’s polynomial eigenvalues
(filled circles) and the first $13$ non-polynomial eigenvalues (open circles).

By using the canonical confluent Heun polynomials, Eq. (89b) may be rewritten
as

$\displaystyle\hat{L}^{(p)}\,\vec{y}_{pk}\;=\;\lambda_{pk}\,\vec{y}_{pk}\;,$
(90)

where $\hat{L}^{(p)}$ and $\vec{y}_{pk}$, respectively, have elements given by

$\displaystyle{L}^{(p)}_{\alpha\beta}\;$
$\displaystyle=\;{\langle{u_{\alpha}}\,|\,{\mathcal{L}_{\xi}^{(p)}|\,u_{\beta}}\rangle}\;,$
(91a) $\displaystyle{y}_{pk,\beta}\;$
$\displaystyle=\;{\langle{u_{\beta}}\,|\,{y_{pk}}\rangle}\;.$ (91b)

Since $k$ and $\alpha$ are both integers between $0$ and $p$, the sequence of
all eigenvectors comprise the rows of a matrix whose elements are
${y}_{p\alpha,\beta}$. Similarly, the eigenvalues may be labeled as
$\lambda_{p\alpha}$ rather than $\lambda_{pk}$.

At focus in Fig. 4 are the confluent Heun polynomial problem’s matrix
representation, eigenvectors, and eigenvalues for order-$6$ polynomials.
Physical parameters used are identical to those in Table (1) for
$(\ell,m,n)=(2,2,6)$.

|  |  |
---|---|---|---
Figure 5:  Radial eigenvalue distributions for select quasinormal mode cases,
all with BH dimensionless spin $a/M=0.7$, and angular indices
$(\ell,m)=(2,2)$. In all panels, radial eigenvalues, $A_{n^{\prime}}$, are
noted with open circles, and the known QNM eigenvalue, $A_{\ell mn}$, is shown
as a red filled circle (c.f. Fig. 1). Left to right: Radial eigenvalues for
the respective $n\in\\{0,3,6,12\\}$ QNM overtones. Respective absolute
differences between known QNM eigenvalues (Refs. London  _et al._ (2020);
Stein (2019)) and those computed here (i.e. $|A_{\ell mn}-A_{n}|$) are
$\\{2.72\,\times 10^{-9},4.72\,\times 10^{-9},1.24\,\times
10^{-9},6.28\,\times 10^{-9}\\}$. Center right ($n=6$): Eigenvalues for the
closest confluent Heun polynomial problem are shown as black triangles. See
panel $(c)$ of Fig. 4 for comparison.

In panel $(a)$ of Fig. 4, ${L}^{(p)}_{\alpha\beta}$ are shown. There it may be
observed that

$\displaystyle{L}^{(p)}_{6,7}\;=\;{L}^{(p)}_{7,6}\;=\;0\;.$ (92)

This is consistent with the so-called “$\Delta_{p+1}$ condition” which is used
to define the confluent Heun polynomials via a finite dimensional matrix
eigenvalue problem (see Eq. 79 of Paper I as well as Refs. DLMF ; Cook and
Zalutskiy (2014)). There the basic concept is that a power series solution to
Eq. (89b) will terminate after $p+1$ terms if the $\Delta_{p+1}$ condition on
a related three-term recursion is satisfied.

The $\Delta_{p+1}$ condition is unable to provide insight on the non-
polynomial solutions, simply due to fact that these solutions have series
expansions of infinite length. In other words, the non-polynomial solutions
satisfy an infinite dimensional matrix eigenvalue problem, rather than a
finite dimensional one555 In such cases, standard routes to finding solutions
include the use of continued fractions (see e.g. Ref. Leaver (1985)) or the
representation of the problem with an orthonormal basis, such as the canonical
confluent Heun polynomial basis introduced here. . Unlike the power-series
perspective of the $\Delta_{p+1}$ condition, the canonical confluent Heun
polynomial representation of the problem naturally applies to finite and
infinite dimensional problems: i.e., the monomials are the basis for power
series expansions, but unlike the canonical confluent Heun polynomial,
monomials are not an orthonormal basis, meaning that standard matrix
interpretations do not apply for non-terminating series. Since the canonical
confluent Heun polynomial are orthonormal, the polynomial and non-polynomial
sectors are clearly visible in Fig. 4a.

At focus in Fig. 4b are the eigenvectors; rows labeled in $\alpha$ correspond
to vectors $\vec{y}_{p\alpha}$. There, the polynomial solutions are shown as
rows that have zero values after $\beta=6$, as is equivalent to each confluent
Heun polynomial only being comprised of canonical confluent Heun polynomials
of order $6$ and below. The non-polynomial solutions are seen as the lower
right-hand block, for which values are non-zero only for $\beta>6$, i.e. the
non-polynomial solutions do not contain canonical confluent Heun polynomials
with order $p\leq 6$.

Polynomial and non-polynomial eigenvalues are shown in Fig. 4c as points in
the complex plane. There, polynomial eigenvalues show a distinctly different
distribution than their non-polynomial counterparts. If a QNM is approximately
polynomial as described in Sec. IV, then its eigenvalue distribution should be
expected to also reflect polynomial/non-polynomial duality.

### QNM radial eigenvalue distributions

– It was noted in Sec. IV that the $(\ell,m,n)=(2,2,6)$ QNM may be
approximately polynomial for a BH spin of $a=0.7$ (Table 1). In particular,
this QNM’s upper pseudo-polynomial order (Eq. 33) is

$\displaystyle{p_{\star}^{+}}\;$ $\displaystyle=\;6.2586+0.0743i\;.$ (93)

Since the imaginary part of ${p_{\star}^{+}}$ is small, and the real part is
close to $6$ (c.f. Eq. 37), it is reasonable to expect rough agreement between
that QNM’s radial eigenvalues, and the order-$6$ confluent Heun polynomial’s
polynomial and non-polynomial eigenvalues (i.e. ${p_{\star}^{+}}\approx 6$).

In Fig. 5, the radial eigenvalue distributions for four select QNMs are shown.
There, as $n$ increases, so does the approximately polynomial sector of the
eigenvalue distribution. For the $n=6$ overtone, the order-$6$ confluent Heun
polynomial eigenvalues have been added for comparison. In that case, it is
clear that the qualitative features of the QNM eigenvalue distribution are
precisely captured by that of the confluent Heun polynomial eigenvalues.
Further, it is clear that, for this case, the QNM separation constant,
$A_{226}$, is indeed well approximated by the dominant order-$6$ confluent
Heun polynomial eigenvalue,

$\displaystyle
A_{226}\;\approx\;{\text{C}_{0}}+{\text{C}_{1}}+\lambda_{66}\;.$ (94)

In Eq. (94), Eq. (31) has been restated in the current context, and the
parameters, ${\text{C}_{0}}$ and ${\text{C}_{1}}$, have been computed with the
relevant QNM frequency $\tilde{\omega}_{226}$. The fractional error between
the left and right hand sides of Eq. (94) is $4.6\%$.

### Validation of radial functions

– When represented with canonical confluent Heun polynomials, the QNMs’ radial
problem may be written as

$\displaystyle\mathcal{L}{{}_{\xi}}\,|{f_{n}}\rangle\;$
$\displaystyle=\;A_{n}\,|{f_{n}}\rangle\;.$ (95)

In Eq. (95), we have simply rewritten Eq. (13a) with the overtone label, $n$,
explicitly added. When represented with canonical confluent Heun polynomials,
the radial functions, $f_{n}(\xi)$, are

$\displaystyle|{f_{n}}\rangle\;=\;\sum_{k=0}^{\infty}\,|{u_{k}}\rangle{\langle{u_{k}}\,|\,{f_{n}}\rangle}\;.$
(96)

In coordinate notation, this is simply

$\displaystyle{f_{n}}(\xi)\;=\;\sum_{k=0}^{\infty}\,{u_{k}}(\xi)\,{\langle{u_{k}}\,|\,{f_{n}}\rangle}\;.$
(97)

Upon solving the matrix eigenvalue problem (Eq. 59b), the accuracy of the each
$f_{n}(\xi)$ may be quantified by

$\displaystyle|(\mathcal{L}_{\xi}-A_{n})f_{n}(\xi)|/|f_{n}(\xi)|\;=\;\Delta_{n}(\xi)\;.$
(98)

In Eq. (98), Eq. (95) has been rearranged such that all quantities appear on
the left hand side. The absolute value has then been taken to yield a single,
real-valued error measure.

|  |  |
---|---|---|---
Figure 6:  Radial functions, $f_{n}(\xi)$, and related floating point error
estimates, $\Delta_{n}$, for the same cases shown in Fig. 5. Top row: The real
part of each radial function, $\text{Re}\,f_{n}(\xi)$, is denoted with a solid
orange line. The imaginary part, $\text{Im}\,f_{n}(\xi)$, is denoted with a
dashed purple line, and the function’s positive and negative magnitude,
$\pm|f_{n}(\xi)|$ are denoted with dotted black lines. Bottom row: errors for
each radial function, $\Delta_{n}$ (Eq. 98). See Sec. VII for related
discussion.

If $f_{n}(\xi)$ are eigenfunctions, and the derivatives within
$\mathcal{L}_{\xi}$ can be evaluated with infinite precision, then
$\Delta_{n}(\xi)$ would equal zero for all values of $\xi$. In practice,
$\Delta_{n}(\xi)$ will be very small relative to one, and dominated by
numerical noise.

In Fig. 6, example radial functions are shown along with their errors for the
same QNM cases discussed in Fig. 5. There, errors have been computed using
floating-point precision and spline interpolation for the derivatives Harris
_et al._ (2020); London  _et al._ (2020). Despite the result being
significantly larger than the calculation’s original precision, all cases are
observed to have dimensionless errors on the order of $10^{-7}$. In the top
row of panels, radial functions have been scaled such that their largest
absolute value is $1$. There, it may be observed that as the overtone number
increases, the relative value of the radial function at the BH event horizon
appears to decrease, possibly implying tension with the physical (asymptotic)
boundary conditions, which require $f_{n}(\xi)$ to be non-zero when
$\xi\in\\{0,1\\}$.

### Schwarzschild and Kerr

– Together, the previous examples enable the following observations. When
considered for a single frequency parameter, $\tilde{\omega}$, the radial
functions typically have distinct eigenvalues (Fig. 5), and so are linearly
independent Axler (2015). Although the space of eigenfunctions is countably
infinite, finite dimensional truncations of this space may be solved for using
canonical confluent Heun polynomials. In practice, the radial functions will
be finitely many and linearly independent, and will therefore constitute a
basis.

Further, due to the self-adjoint nature of the radial operator with respect to
the scalar product (Eq. 13), the basis of radial functions is orthogonal
London (2023a). Since the canonical confluent Heun polynomial are an
orthonormal basis, scalar products between different radial functions,
$|{f_{p}}\rangle$ and $|{f_{q}}\rangle$, are equivalent to the simple vector
dot product between their vector representations,

$\displaystyle{\langle{f_{p}}\,|\,{f_{q}}\rangle}\;$
$\displaystyle=\;\sum_{jk}\,{\langle{f_{p}}\,|\,{u_{j}}\rangle}{\langle{u_{j}}\,|\,{u_{k}}\rangle}{\langle{u_{k}}\,|\,{f_{q}}\rangle}\;$
(99a)
$\displaystyle=\;\sum_{k}\,{\langle{f_{p}}\,|\,{u_{k}}\rangle}{\langle{u_{k}}\,|\,{f_{q}}\rangle}\;$
(99b) $\displaystyle=\;\vec{f}_{p}\cdot\vec{f}_{q}\;.$ (99c)

Consequently, the radial functions may be normalized in the standard way,
$|{f_{p}}\rangle\leftarrow|{f_{p}}\rangle/\sqrt{{\langle{f_{p}}\,|\,{f_{p}}\rangle}}$.
In other words, when computed for a single frequency parameter, the radial
functions define an orthonormal basis.

A central question in single BH perturbation theory is whether QNMs for one
physical case may be used to efficiently represent QNMs in another. For
example, it is extremely common to represent the Kerr QNMs’ angular functions
in the basis of those for Schwarzschild London (2023a, b); Berti and Klein
(2014); O’Sullivan and Hughes (2014); Cook and Zalutskiy (2014); Press and
Teukolsky (1973). Here, we are concerned with whether the same is true for the
QNMs’ radial functions. In effect, we ask whether there exists a mixing matrix
between sets of radial functions derived from different BH spins.

To answer this question, it cannot be avoided that, even for a single QNM,
different values of BH spin correspond to different QNM frequencies London and
Fauchon-Jones (2019); Leaver (1985); Berti and Klein (2014); Berti  _et al._
(2016). When combined with the fact that the scalar product’s weight function
(Eq. 16) is frequency dependent, this means that the QNM radial functions for
different BH spins are members of different scalar product spaces. Therefore,
a complete and selfconsistent treatment of the question requires the use of a
frequency-independent weight function. While such a weight function was
defined in Sec. IV-B of Paper I, we consider a complete treatment of its use
to be outside of the current scope.

In lieu of a more rigorous treatment, Eq. (99) may be used to construct a
simplified one. Let the $p^{\text{th}}$ radial function for a given BH spin,
$a$, be denoted $\vec{f}_{p}^{\,(a)}$. As mentioned previously, if
$\vec{f}_{p}^{\,(a)}$ are defined with QNM frequency $\tilde{\omega}_{{\ell
mn}}$, then $\vec{f}_{n}^{\,(a)}$ is the QNM’s radial function. Since each
$a$-dependent canonical confluent Heun polynomial basis is orthonormal, Eq.
(99) may be interpreted to mean that the scalar product for all
$\vec{f}_{p}^{\,(a)}$ is given by the right hand side of Eq. (99c). Since the
right hand side of Eq. (99c) has no weight function, we may superficially
conclude that the standard vector dot product is a frequency-independent
scalar product of potential use. This conclusion is only superficial because
it ignores the fact that a frequency dependent weight is still required to
arrive at Eq. (99c). Nevertheless, use of Eq. (99c) for radial functions of
different BH spins is sufficient for the present example.

Let the scalar product between a Schwarzschild (i.e. $a=0$) and Kerr radial
function be

|  |  |
---|---|---|---
Figure 7:  Mixing matrices between Schwarzschild and Kerr radial functions for
four values of BH spin: $a/M\in\\{0,0.01,0.70,0.94\\}$. For each case, the
first $24$ radial functions have been computed using canonical confluent Heun
polynomials derived from the $s=-2$, $(\ell,m,n)=(2,2,0)$ QNM frequencies.
Mixing matrices have been computed according to Eq. (100).

$\displaystyle\vec{f}_{p}^{(0)}\cdot\vec{f}_{q}^{(a)}\;=\;\sum_{k}\,{\langle{f_{p}^{(0)}}\,|\,{u_{k}^{(0)}}\rangle}{\langle{u_{k}^{(a)}}\,|\,{f_{q}^{(a)}}\rangle}\;.$
(100)

In Eq. (100), $u_{k}^{(a)}$ are canonical confluent Heun polynomials derived
from a BH spin of $a$, where $(s,l,m,n)$ are implicitly fixed. Vectors
$\vec{f}_{p}^{(0)}$ and $\vec{f}_{q}^{(a)}$ will be considered normalized,
making e.g. the set of all $\vec{f}_{p}^{(0)}$ orthonormal. Therefore, the
decomposition of a Kerr radial function onto the basis of Schwarzschild radial
functions is

$\displaystyle\vec{f}_{q}^{(a)}\;=\;\sum_{q=0}\;\vec{f}_{p}^{(0)}\;(\vec{f}_{p}^{(0)}\cdot\vec{f}_{q}^{(a)})\;.$
(101)

In Eq. (101), $(\vec{f}_{p}^{(0)}\cdot\vec{f}_{q}^{(a)})$ are elements of the
mixing matrix between the Schwarzschild and Kerr radial functions. In the case
of $a=0$, Eq. (101) simply reduces to
$\vec{f}_{q}^{(a)}=\sum_{q=0}\vec{f}_{p}^{(0)}\delta_{pq}$, signaling a
diagonal mixing matrix.

In Fig. 7, the $(s,\ell,m,n)=(-2,2,2,0)$ QNM frequencies have been used to
compute Schwarzschild-Kerr mixing matrices for four values of BH spin666 The
successful reproduction of Fig. 7 requires that one keep in mind the overall
sign ambiguity discussed in Sec. V: Given canonical polynomials, $u_{p}^{(a)}$
and $u_{p}^{(a^{\prime})}$, one can avoid the presence of an incongruous minus
sign by using, e.g., the scalar product space of $u_{p}^{(a)}$ to compute the
matrix whose diagonal elements are
${\langle{u_{p}^{(a)}}\,|\,{u_{p}^{(a^{\prime})}}\rangle}$. Diagonal elements
with a negative real part signal that $u_{p}^{(a^{\prime})}$ should be scaled
by $-1$. . As seen from left to right in Fig. 7, as BH spin increases, so too
does the magnitude of off-diagonal matrix elements. Nevertheless, the mixing
matrix appears to be approximately diagonal throughout the range of spins
shown.

This result implies that there exists a one-to-one relationship between
Schwarzschild and Kerr radial functions across a wide range of BH spins.

Such a one-to-one relationship is known to exist between the Schwarzschild and
Kerr angular functions London (2023b); Berti and Klein (2014). These functions
are the spin weighted spherical and spheroidal harmonics, respectively Newman
and Penrose (1966); Goldberg  _et al._ (1967); London (2023b); O’Sullivan and
Hughes (2014); Cook and Zalutskiy (2014). It was found in Ref. London (2023b)
that the existence of a spherical-to-spheroidal mapping underpins the
existence of another set of angular functions called adjoint-spheroidal
harmonics, the significance of which may be understood thusly: Due to their
different QNM frequencies, the QNMs’ spheroidal harmonics are not orthogonal
with each other; however, they are biorthogonal with the adjoint-spheroidal
harmonics. Therefore the adjoint-spheroidal harmonics signal QNM orthogonality
between modes of different $\ell$. This perspective also defines the sense in
which the QNMs are complete in $\ell$ London (2023b).

In this context, Fig. 7 implies the potential bi-orthogonality and
completeness of the QNMs’ radial functions, i.e. the set of all QNM radial
functions, where each radial function has frequency $\tilde{\omega}_{{\ell
mn}}$, and where only $n$ is allowed to vary. However, given the close
relationship between the QNMs’ radial functions and the confluent Heun
polynomials, it is possible that the radial functions are not complete, but
rather overcomplete (Sec. IV). Indeed, overcompleteness is perhaps most
likely, given the known existence of QNM multiplets Cook and Zalutskiy (2014);
Maassen van den Brink (2004).

##  VIII Discussion and conclusions

The present work is based upon a simple premise: to solve Teukolsky’s radial
equation, rather than looking to special functions from other problems, we
need only look to the equation itself London (2023a). The result is a method
for computing global solutions to the confluent Heun equation via a single
expansion in polynomials that are natural to the problem. We have introduced
such polynomials in Sec. V. In Sec. VI, we have shown that these polynomials
enable the radial equation to be represented as a matrix that is both
tridiagonal and symmetric. Lastly, in Sec. VII, we have demonstrated the use
of these polynomials in various numerical examples, the last of which may be
of particular relevance to the time-dependent perturbations of BHs Green  _et
al._ (2023); Cannizzaro  _et al._ (2023); Redondo-Yuste  _et al._ (2023);
Sberna  _et al._ (2022); Lagos and Hui (2022). We now briefly discuss
connections and implications that have not been fully articulated thus far.

As noted previously, one of our goals has been to gain a deeper mathematical
understanding of the QNM’s radial problem. Consequently, our method for
solving Teukolsky’s radial equation is particular to the QNMs. If our method
can be generalized to non-QNM boundary conditions, such as those of a point
particle orbiting a BH, then new techniques of analysis and methods of
evaluation may result. Concurrently, it may also be useful to reproduce our
results for other formulations of BH perturbation theory, and in more exotic
spacetimes Mino  _et al._ (1997); Berti  _et al._ (2009a); Hughes (2000);
Ripley (2022).

Of the techniques presented here, a number are open to refinement, further
exploration, and connection with other recent results. For example, we have
shown that, in some cases, the QNM’s separation constant is well approximated
by the eigenvalue of a confluent Heun polynomial. We expect that more analytic
understanding is tractable when either representing the problem in confluent
Heun polynomials, or their canonical counterparts. Further, there is good
reason to think that more may be learned by a thorough comparison between our
results, and those known for Pollaczek-Jacobi polynomials Chen and Dai (2010);
Chen  _et al._ (2019); Min and Chen (2023); Yu  _et al._ (2022). This effort
is ongoing and near completion Gurevich and London (2024). There is also the
possibility of combining the results presented here with recent iterative
generalized spectral methods for computing QNM frequencies Chung  _et al._
(2023a); Blázquez-Salcedo  _et al._ (2023); Chung  _et al._ (2023b); Ghojogh
_et al._ (2023); Zhu  _et al._ (2023).

Lastly, one clear next step is to expand upon the results of Sec. VII to
better understand whether the QNM’s radial functions are complete. This would
likely require a greater practical understanding of the complex-valued radial
coordinate introduced in Paper I. Furthermore, this may require an analysis
similar to that of Ref. London (2023b) for the spheroidal harmonics.
Intriguingly, if the QNM’s radial functions are complete (in some practical
way), then it would follow that any bounded spatial deviation from Kerr might
be represented in terms of QNMs.

##  IX Acknowledgements

The authors thank Katy Clough, Stephen Green, Scott Hughes, Laura Sberna and
Bernard Whiting for supportive discussions. Lionel London was funded at King’s
College London by the Royal Society URF\R1\211451. Michelle Gurevich was
funded at King’s College London by Royal Society enhancement grant
RF\ERE\210040 and by the King’s College London department of physics.

## References

  * Vishveshwara (1970) C. V. Vishveshwara, Phys. Rev. D 1, 2870 (1970).
  * Abbott  _et al._ (2016a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016a), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ (2016b) B. P. Abbott _et al._ (Virgo, LIGO Scientific), Phys. Rev. Lett. 116, 221101 (2016b), arXiv:1602.03841 [gr-qc] .
  * Abbott  _et al._ (2021a) R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 11, 021053 (2021a), arXiv:2010.14527 [gr-qc] .
  * Abbott  _et al._ (2021b) R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), (2021b), arXiv:2111.03606 [gr-qc] .
  * Abbott  _et al._ (2019) B. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 9, 031040 (2019), arXiv:1811.12907 [astro-ph.HE] .
  * Hughes  _et al._ (2019) S. A. Hughes, A. Apte, G. Khanna, and H. Lim, Phys. Rev. Lett. 123, 161101 (2019), arXiv:1901.05900 [gr-qc] .
  * Khan  _et al._ (2016) S. Khan, S. Husa, M. Hannam, F. Ohme, M. Pürrer, X. Jiménez Forteza, and A. Bohé, Phys. Rev. D93, 044007 (2016), arXiv:1508.07253 [gr-qc] .
  * Husa  _et al._ (2016) S. Husa, S. Khan, M. Hannam, M. Pürrer, F. Ohme, X. Jiménez Forteza, and A. Bohé, Phys. Rev. D93, 044006 (2016), arXiv:1508.07250 [gr-qc] .
  * Pompili  _et al._ (2023) L. Pompili _et al._ , (2023), arXiv:2303.18039 [gr-qc] .
  * Press (1971) W. H. Press, Astrophys. J. Lett. 170, L105 (1971).
  * Press and Teukolsky (1973) W. H. Press and S. A. Teukolsky, Astrophys. J. 185, 649 (1973).
  * Teukolsky (1973) S. A. Teukolsky, Astrophys. J. 185, 635 (1973).
  * Teukolsky and Press (1974) S. A. Teukolsky and W. H. Press, Astrophysical Journal 193, 443 (1974).
  * Leaver (1985) E. Leaver, Proc. Roy. Soc. Lond. A A402, 285 (1985).
  * Berti  _et al._ (2006a) E. Berti, V. Cardoso, and C. M. Will, Phys. Rev. D73, 064030 (2006a), arXiv:gr-qc/0512160 [gr-qc] .
  * London  _et al._ (2014) L. London, D. Shoemaker, and J. Healy, Phys. Rev. D90, 124032 (2014), arXiv:1404.3197 [gr-qc] .
  * Hamilton  _et al._ (2023) E. Hamilton, L. London, and M. Hannam, Phys. Rev. D 107, 104035 (2023), arXiv:2301.06558 [gr-qc] .
  * Carullo  _et al._ (2018) G. Carullo _et al._ , Phys. Rev. D 98, 104020 (2018), arXiv:1805.04760 [gr-qc] .
  * Baibhav  _et al._ (2023) V. Baibhav, M. H.-Y. Cheung, E. Berti, V. Cardoso, G. Carullo, R. Cotesta, W. Del Pozzo, and F. Duque, (2023), arXiv:2302.03050 [gr-qc] .
  * Carullo  _et al._ (2019) G. Carullo, W. Del Pozzo, and J. Veitch, Phys. Rev. D 99, 123029 (2019), [Erratum: Phys.Rev.D 100, 089903 (2019)], arXiv:1902.07527 [gr-qc] .
  * Cotesta  _et al._ (2022) R. Cotesta, G. Carullo, E. Berti, and V. Cardoso, Phys. Rev. Lett. 129, 111102 (2022), arXiv:2201.00822 [gr-qc] .
  * Mishra  _et al._ (2010) C. K. Mishra, K. G. Arun, B. R. Iyer, and B. S. Sathyaprakash, Phys. Rev. D 82, 064010 (2010), arXiv:1005.0304 [gr-qc] .
  * Bhagwat  _et al._ (2020) S. Bhagwat, X. J. Forteza, P. Pani, and V. Ferrari, Phys. Rev. D 101, 044033 (2020), arXiv:1910.08708 [gr-qc] .
  * Lagos and Hui (2022) M. Lagos and L. Hui, (2022), arXiv:2208.07379 [gr-qc] .
  * Green  _et al._ (2023) S. R. Green, S. Hollands, L. Sberna, V. Toomani, and P. Zimmerman, Phys. Rev. D 107, 064030 (2023), arXiv:2210.15935 [gr-qc] .
  * Cannizzaro  _et al._ (2023) E. Cannizzaro, L. Sberna, S. R. Green, and S. Hollands, (2023), arXiv:2309.10021 [gr-qc] .
  * Sberna  _et al._ (2022) L. Sberna, P. Bosch, W. E. East, S. R. Green, and L. Lehner, Phys. Rev. D 105, 064046 (2022), arXiv:2112.11168 [gr-qc] .
  * Redondo-Yuste  _et al._ (2023) J. Redondo-Yuste, D. Pereñiguez, and V. Cardoso, (2023), arXiv:2312.04633 [gr-qc] .
  * Arun  _et al._ (2022) K. G. Arun _et al._ (LISA), Living Rev. Rel. 25, 4 (2022), arXiv:2205.01597 [gr-qc] .
  * Maselli  _et al._ (2017) A. Maselli, K. Kokkotas, and P. Laguna, Phys. Rev. D95, 104026 (2017), arXiv:1702.01110 [gr-qc] .
  * Berti  _et al._ (2016) E. Berti, A. Sesana, E. Barausse, V. Cardoso, and K. Belczynski, Phys. Rev. Lett. 117, 101102 (2016), arXiv:1605.09286 [gr-qc] .
  * Barausse  _et al._ (2015) E. Barausse, V. Cardoso, and P. Pani, J. Phys. Conf. Ser. 610, 012044 (2015), arXiv:1404.7140 [astro-ph.CO] .
  * Klein  _et al._ (2016) A. Klein _et al._ , Phys. Rev. D 93, 024003 (2016), arXiv:1511.05581 [gr-qc] .
  * Bailes  _et al._ (2021) M. Bailes, B. K. Berger, P. R. Brady, M. Branchesi, K. Danzmann, M. Evans, K. Holley-Bockelmann, B. R. Iyer, T. Kajita, S. Katsanevas, M. Kramer, A. Lazzarini, L. Lehner, G. Losurdo, H. Lück, D. E. McClelland, M. A. McLaughlin, M. Punturo, S. Ransom, S. Raychaudhury, D. H. Reitze, F. Ricci, S. Rowan, Y. Saito, G. H. Sanders, B. S. Sathyaprakash, B. F. Schutz, A. Sesana, H. Shinkai, X. Siemens, D. H. Shoemaker, J. Thorpe, J. F. J. van den Brand, and S. Vitale, Nature Reviews Physics 3, 344 (2021).
  * Vallisneri  _et al._ (2015) M. Vallisneri, J. Kanner, R. Williams, A. Weinstein, and B. Stephens,  _Proceedings, 10th International LISA Symposium: Gainesville, Florida, USA, May 18-23, 2014_ , J. Phys. Conf. Ser. 610, 012021 (2015), arXiv:1410.4839 [gr-qc] .
  * Karnesis  _et al._ (2022) N. Karnesis _et al._ , (2022), arXiv:2209.04358 [gr-qc] .
  * Kokkotas and Schmidt (1999) K. D. Kokkotas and B. G. Schmidt, Living Review in Relativity 2 (1999).
  * Nollert (1999) H.-P. Nollert, Class. Quant. Grav. 16, R159 (1999).
  * Berti  _et al._ (2009a) E. Berti, V. Cardoso, and A. O. Starinets, Class.Quant.Grav. 26, 163001 (2009a), arXiv:0905.2975 [gr-qc] .
  * Cook and Zalutskiy (2014) G. B. Cook and M. Zalutskiy, Phys. Rev. D 90, 124021 (2014), arXiv:1410.7698 [gr-qc] .
  * Leaver (1986) E. W. Leaver, Journal of Mathematical Physics 27, 1238 (1986).
  * Leaver (1986) E. W. Leaver, Phys. Rev. D 34, 384 (1986).
  * Berti  _et al._ (2003) E. Berti, V. Cardoso, K. D. Kokkotas, and H. Onozawa, Phys. Rev. D 68, 124018 (2003), arXiv:hep-th/0307013 .
  * Berti  _et al._ (2006b) E. Berti, V. Cardoso, and M. Casals, Phys. Rev. D73, 024013 (2006b), [Erratum: Phys. Rev.D73,109902(2006)], arXiv:gr-qc/0511111 [gr-qc] .
  * Berti and Cardoso (2006) E. Berti and V. Cardoso, Physical Review D 74 (2006), 10.1103/physrevd.74.104020.
  * Yang  _et al._ (2013) H. Yang, F. Zhang, A. Zimmerman, D. A. Nichols, E. Berti, and Y. Chen, Phys. Rev. D87, 041502 (2013), arXiv:1212.3271 [gr-qc] .
  * Yang  _et al._ (2012) H. Yang, D. A. Nichols, F. Zhang, A. Zimmerman, Z. Zhang, and Y. Chen, Phys. Rev. D 86, 104006 (2012), arXiv:1207.4253 [gr-qc] .
  * Zimmerman  _et al._ (2015) A. Zimmerman, H. Yang, F. Zhang, D. A. Nichols, E. Berti, and Y. Chen, (2015), arXiv:1510.08159 [gr-qc] .
  * Zimmerman and Mark (2016) A. Zimmerman and Z. Mark, Phys. Rev. D93, 044033 (2016), [Erratum: Phys. Rev.D93,no.8,089905(2016)], arXiv:1512.02247 [gr-qc] .
  * Maggiore (2008) M. Maggiore, Phys. Rev. Lett. 100, 141301 (2008), arXiv:0711.3145 [gr-qc] .
  * Andersson (1995) N. Andersson, Phys. Rev. D 51, 353 (1995).
  * Andersson (1997) N. Andersson, Phys. Rev. D 55, 468 (1997), arXiv:gr-qc/9607064 .
  * London (2023a) L. London, “A radial scalar product for kerr quasinormal modes,” (2023a), arXiv:2312.17678 [gr-qc] .
  * Fackerell and Crossman (2008) E. D. Fackerell and R. G. Crossman, Journal of Mathematical Physics 18, 1849 (2008), https://pubs.aip.org/aip/jmp/article-pdf/18/9/1849/11184786/1849_1_online.pdf .
  * London (2023b) L. T. London, Phys. Rev. D 107, 044056 (2023b), arXiv:2006.11449 [gr-qc] .
  * Ronveaux and Arscott (1995) A. Ronveaux and F. Arscott,  _Heun’s Differential Equations_, Oxford science publications (Oxford University Press, 1995).
  * (58) DLMF, “Solutions Analytic at Three Singularities: Heun Polynomials,” https://dlmf.nist.gov/, Release 1.1.10 of 2023-06-15, f. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.
  * Fiziev (2009a) P. P. Fiziev, Phys. Rev. D 80, 124001 (2009a), arXiv:0906.5108 [gr-qc] .
  * Hortacsu (2012) M. Hortacsu,  , 23 (2012), arXiv:1101.0471 [math-ph] .
  * Dariescu  _et al._ (2021) C. Dariescu, M.-A. Dariescu, and C. Stelea, Adv. High Energy Phys. 2021, 5512735 (2021), arXiv:2102.03850 [hep-th] .
  * Magnus  _et al._ (2021) A. P. Magnus, F. Ndayiragije, and A. Ronveaux, Journal of Approximation Theory 263, 105522 (2021).
  * Pinchover and Rubinstein (2005) Y. Pinchover and J. Rubinstein, “Equations in high dimensions,” in  _An Introduction to Partial Differential Equations_ (Cambridge University Press, 2005) p. 226–281.
  * Abramowitz and Stegun (1964) M. Abramowitz and I. A. Stegun,  _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_ , ninth dover printing, tenth gpo printing ed. (Dover, New York, 1964).
  * Arfken  _et al._ (2013) G. B. Arfken, H. J. Weber, and F. E. Harris, in  _Mathematical Methods for Physicists (Seventh Edition)_, edited by G. B. Arfken, H. J. Weber, and F. E. Harris (Academic Press, Boston, 2013) seventh edition ed., pp. 401–445.
  * Courant and Hilbert. (1954) R. Courant and D. Hilbert., Quarterly Journal of the Royal Meteorological Society 80, 485 (1954).
  * Chen and Dai (2010) Y. Chen and D. Dai, Journal of Approximation Theory 162, 2149 (2010).
  * Axler (2015) S. Axler,  _Linear algebra done right (eBook)_, 3rd ed. (Springer, Cham, 2015).
  * Chung  _et al._ (2023a) A. K.-W. Chung, P. Wagle, and N. Yunes, Phys. Rev. D 107, 124032 (2023a), arXiv:2302.11624 [gr-qc] .
  * Blázquez-Salcedo  _et al._ (2023) J. L. Blázquez-Salcedo, F. S. Khoo, J. Kunz, and L. M. González-Romero, (2023), arXiv:2312.10754 [gr-qc] .
  * Chung  _et al._ (2023b) A. K.-W. Chung, P. Wagle, and N. Yunes, (2023b), arXiv:2312.08435 [gr-qc] .
  * Ghojogh  _et al._ (2023) B. Ghojogh, F. Karray, and M. Crowley, “Eigenvalue and generalized eigenvalue problems: Tutorial,” (2023), arXiv:1903.11240 [stat.ML] .
  * Zhu  _et al._ (2023) H. Zhu, J. L. Ripley, A. Cárdenas-Avendaño, and F. Pretorius, (2023), arXiv:2309.13204 [gr-qc] .
  * Teschl (2000) G. Teschl,  _Jacobi Operators and Completely Integrable Nonlinear Lattices_, Mathematical surveys and monographs (American Mathematical Society, 2000).
  * Kuijlaars  _et al._ (2003) A. B. J. Kuijlaars, A. Martinez-Finkelshtein, and R. Orive, “Orthogonality of jacobi polynomials with general parameters,” (2003), arXiv:math/0301037 [math.CA] .
  * Magaña Zertuche  _et al._ (2021) L. Magaña Zertuche _et al._ , (2021), arXiv:2110.15922 [gr-qc] .
  * Berti  _et al._ (2009b) E. Berti, V. Cardoso, and A. O. Starinets, Classical and Quantum Gravity 26, 163001 (2009b).
  * Jaramillo  _et al._ (2021) J. L. Jaramillo, R. Panosso Macedo, and L. Al Sheikh, Phys. Rev. X 11, 031003 (2021), arXiv:2004.06434 [gr-qc] .
  * Boyanov  _et al._ (2023) V. Boyanov, V. Cardoso, K. Destounis, J. L. Jaramillo, and R. Panosso Macedo, (2023), arXiv:2312.11998 [gr-qc] .
  * Destounis and Duque (2023) K. Destounis and F. Duque (2023) arXiv:2308.16227 [gr-qc] .
  * Courty  _et al._ (2023) A. Courty, K. Destounis, and P. Pani, Phys. Rev. D 108, 104027 (2023), arXiv:2307.11155 [gr-qc] .
  * Torres (2023) T. Torres, Phys. Rev. Lett. 131, 111401 (2023), arXiv:2304.10252 [gr-qc] .
  * Leung  _et al._ (1994) P. T. Leung, S. Y. Liu, and K. Young, Phys. Rev. A 49, 3057 (1994).
  * Wald (2008) R. M. Wald, Journal of Mathematical Physics 20, 1056 (2008), https://pubs.aip.org/aip/jmp/article-pdf/20/6/1056/8150384/1056_1_online.pdf .
  * Beyer (1995) H. R. Beyer, Journal of Mathematical Physics 36, 4815 (1995), https://pubs.aip.org/aip/jmp/article-pdf/36/9/4815/7505756/4815_1_online.pdf .
  * O’Sullivan and Hughes (2014) S. O’Sullivan and S. A. Hughes, Phys. Rev. D 90, 124039 (2014), [Erratum: Phys.Rev.D 91, 109901 (2015)], arXiv:1407.6983 [gr-qc] .
  * Berti and Klein (2014) E. Berti and A. Klein, Phys. Rev. D 90, 064012 (2014), arXiv:1408.1860 [gr-qc] .
  * Blanchet (2014) L. Blanchet, Living Rev. Rel. 17, 2 (2014), arXiv:1310.1528 [gr-qc] .
  * Boyer and Lindquist (2004) R. H. Boyer and R. W. Lindquist, Journal of Mathematical Physics 8, 265 (2004), https://pubs.aip.org/aip/jmp/article-pdf/8/2/265/11041897/265_1_online.pdf .
  * Newman and Penrose (1962) E. Newman and R. Penrose, Journal of Mathematical Physics 3, 566 (1962).
  * Newman and Penrose (1966) E. T. Newman and R. Penrose, Journal of Mathematical Physics 7, 863 (1966), https://doi.org/10.1063/1.1931221 .
  * Ruiz  _et al._ (2008) M. Ruiz, R. Takahashi, M. Alcubierre, and D. Nunez, Gen. Rel. Grav. 40, 2467 (2008), arXiv:0707.4654 [gr-qc] .
  * Merlin  _et al._ (2016) C. Merlin, A. Ori, L. Barack, A. Pound, and M. van de Meent, Phys. Rev. D 94, 104066 (2016), arXiv:1609.01227 [gr-qc] .
  * Whiting (1989) B. F. Whiting, Journal of Mathematical Physics 30, 1301 (1989), https://pubs.aip.org/aip/jmp/article-pdf/30/6/1301/8158708/1301_1_online.pdf .
  * Fiziev (2009b) P. P. Fiziev, Journal of Physics A: Mathematical and Theoretical 43, 035203 (2009b).
  * Berti  _et al._ (2007) E. Berti, J. Cardoso, V. Cardoso, and M. Cavaglia, Phys.Rev. D76, 104044 (2007), arXiv:0707.1202 [gr-qc] .
  * Lax (2002) P. Lax,  _Functional analysis_ (Wiley, New York, 2002).
  * García-Quirós  _et al._ (2020) C. García-Quirós, M. Colleoni, S. Husa, H. Estellés, G. Pratten, A. Ramos-Buades, M. Mateu-Lucena, and R. Jaume, (2020), arXiv:2001.10914 [gr-qc] .
  * Kristensson (2010) G. Kristensson,  _Second order differential equations - special functions and their classification_ (Springer, Germany, 2010).
  * Chihara (2011) T. Chihara,  _An Introduction to Orthogonal Polynomials_, Dover Books on Mathematics (Dover Publications, 2011).
  * Cho  _et al._ (2009) H. T. Cho, A. S. Cornell, J. Doukas, and W. Naylor, Phys. Rev. D 80, 064022 (2009), arXiv:0904.1867 [gr-qc] .
  * Sakurai (1993) J. J. Sakurai,  _Modern Quantum Mechanics (Revised Edition)_, 1st ed. (Addison Wesley, 1993).
  * Christensen (2003) O. Christensen,  _An Introduction to Frames and Riesz Bases_ (Birkhäuser Boston, 2003).
  * Gurevich and London (2024) M. Gurevich and L. London, “Properties of natural polynomials for Kerr quasi-normal modes,” (2024), in Prep.
  * Adkins and Weintraub (2012) W. Adkins and S. Weintraub,  _Algebra: An Approach via Module Theory_, Graduate Texts in Mathematics (Springer New York, 2012).
  * Harris  _et al._ (2020) C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. Fernández del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Nature 585, 357–362 (2020).
  * mpmath development team (2023) T. mpmath development team,  _mpmath: a Python library for arbitrary-precision floating-point arithmetic (version 1.3.0)_ (2023), transverse–tracelessh (TTh)ttp://mpmath.org/.
  * London  _et al._ (2020) L. London, E. Fauchon, and E. Hamilton, “llondon6/positive: map,”  (2020).
  * Stein (2019) L. C. Stein, J. Open Source Softw. 4, 1683 (2019), arXiv:1908.10377 [gr-qc] .
  * London and Fauchon-Jones (2019) L. London and E. Fauchon-Jones, Class. Quant. Grav. 36, 235015 (2019), arXiv:1810.03550 [gr-qc] .
  * Goldberg  _et al._ (1967) J. N. Goldberg, A. J. MacFarlane, E. T. Newman, F. Rohrlich, and E. C. G. Sudarshan, J. Math. Phys. 8, 2155 (1967).
  * Maassen van den Brink (2004) A. Maassen van den Brink, J. Math. Phys. 45, 327 (2004), arXiv:gr-qc/0303095 .
  * Mino  _et al._ (1997) Y. Mino, M. Sasaki, M. Shibata, H. Tagoshi, and T. Tanaka, Progress of Theoretical Physics Supplement 128, 1 (1997), http://oup.prod.sis.lan/ptps/article-pdf/doi/10.1143/PTPS.128.1/5438984/128-1.pdf .
  * Hughes (2000) S. A. Hughes, Phys. Rev. D 62, 044029 (2000), [Erratum: Phys.Rev.D 67, 089902 (2003)], arXiv:gr-qc/0002043 .
  * Ripley (2022) J. L. Ripley, Class. Quant. Grav. 39, 145009 (2022), arXiv:2202.03837 [gr-qc] .
  * Chen  _et al._ (2019) Y. Chen, G. Filipuk, and L. Zhan, Journal of Mathematical Physics 60 (2019), 10.1063/1.5102102.
  * Min and Chen (2023) C. Min and Y. Chen, Mathematical Methods in the Applied Sciences 46, 15270–15284 (2023).
  * Yu  _et al._ (2022) J. Yu, C. Li, M. Zhu, and Y. Chen, Journal of Mathematical Physics 63, 063504 (2022), https://pubs.aip.org/aip/jmp/article-pdf/doi/10.1063/5.0062949/16564109/063504_1_online.pdf .