# A radial scalar product for Kerr quasinormal modes

Lionel T. London  King’s College London, Strand, London WC2R 2LS, United
Kingdom  Institute for High-Energy Physics, University of Amsterdam, Science
Park 904, 1098 XH Amsterdam, Netherlands

###### Abstract

A scalar product for quasinormal mode solutions to Teukolsky’s homogeneous
radial equation is presented. Evaluation of this scalar product can be
performed either by direct integration, or by evaluation of a confluent
hypergeometric functions. For direct integration, it is explicitly shown that
the quasinormal modes’ radial functions are regular on a family of physically
bounded complex paths. The related scalar product will be useful for better
understanding analytic solutions to Teukolsky’s radial equation, particularly
the quasi-normal modes, their potential spatial completeness, and whether the
quasi-normal mode overtone excitations may be estimated by spectral
decomposition rather than fitting. With that motivation, the scalar product is
applied to confluent Heun polynomials where it is used to derive their
peculiar orthogonality and eigenvalue properties. A potentially new
relationship is derived between the confluent Heun polynomials’ scalar
products and eigenvalues. Using these results, it is shown for the first time
that Teukolsky’s radial equation (and perhaps similar confluent Heun
equations) are, in principle, exactly tri-diagonalizable. To this end,
“canonical” confluent Heun polynomials are conjectured.

##  I Introduction

The aftermath of binary black hole (BBH) merger is thought to be a spinning
and perturbed black hole (BH) that rings down like a struck, and shimmering
gong Vishveshwara (1970); Abbott  _et al._ (2016a, b, 2021a, 2021b, 2019a).
The corresponding gravitational radiation encodes information about the
progenitor BBH system (i.e. masses, spins, momenta and environment), and about
the true nature of gravity Kamaretsos  _et al._ (2012); London  _et al._
(2014); Hughes  _et al._ (2019); Berti  _et al._ (2006a, 2009); Maselli  _et
al._ (2017); Kokkotas and Schmidt (1999); Evstafyeva  _et al._ (2023); Arun
_et al._ (2022); Pacilio and Bhagwat (2023). Ongoing and planned experiments
endeavor to detect and learn from that radiation with increasing precision
Bailes  _et al._ (2021); Sathyaprakash  _et al._ (2012); Berti  _et al._
(2006a); Vallisneri  _et al._ (2015); Karnesis  _et al._ (2022); Arun  _et
al._ (2022). At all stages, those efforts are guided by our mathematical
understanding of BHs, their gravitational waves, and the special functions
best suited to them Abbott  _et al._ (2017, 2016c); Misner  _et al._ (1973);
Maggiore (2007); Newman and Penrose (1966, 1962); Teukolsky (1972); Whiting
(1989). Recently, it has become clear that there are limitations in our
mathematical understanding of ringdown radiation, and that these limitations
may actively confound the interpretation of numerical and experimental data
Giesler  _et al._ (2019); McWilliams (2019); Carullo  _et al._ (2019); Ota and
Chirenti (2020); Bhagwat  _et al._ (2020); Forteza  _et al._ (2020); Leong
_et al._ (2023); Khera  _et al._ (2023); Franchini and Völkel (2023); Nee  _et
al._ (2023); Cotesta  _et al._ (2022); Baibhav  _et al._ (2023).

The purpose of this article is to introduce a framework (via a radial scalar
product) that may enable a deeper practical understanding of ringdown’s
“overtones”, and perhaps gravitational waves from compact binaries broadly.

To date, the LIGO-Virgo-Kagra collaboration (LVK) has confidently observed
gravitational radiation from $83$ BBH mergers; among them, the ringdown
radiation of $31$ have been compared against gravitational wave signal models
(theory predictions) to test for physics beyond General Relativity (GR) Abbott
_et al._ (2021a, b, 2019a, c, 2019b, d, e). Thus far, all BBH remnants have
been found to be consistent with spinning (i.e. Kerr) BHs Abbott  _et al._
(2019b, 2021d, 2021e).

Planned gravitational wave detectors such as the Einstein Telescope (ET),
Cosmic Explorer (CE) and the Laser Interferometer Space Antenna (LISA) are,
for example, expected to detect signal to noise ratios of $\sim 10^{3}$ (ET &
CE) and $\sim 10^{5}$ (LISA) from GW150914-like111This example refers to BBHs
that have similar configurations to that of the first gravitational wave
detection: approximately non-spinning progenitors and having a luminosity
distance of $\sim 410~{}\mathrm{Mpc}$ Abbott  _et al._ (2016c, a). BBHs of
$10^{2}~{}\mathrm{M}_{\odot}$ and $10^{6}~{}\mathrm{M}_{\odot}$, respectively
Berti  _et al._ (2016); Mishra  _et al._ (2010); Robson  _et al._ (2019). Such
detections are expected to dramatically increase the potential for ringdown
data analyses to verify proposed environmental, dark matter and quantum
effects Klein  _et al._ (2016); Franchini and Völkel (2023); Macedo  _et al._
(2013); Berti  _et al._ (2006a); Barausse  _et al._ (2015); Zhao  _et al._
(2023). However, the correct interpretation of analyses involving ringdown is
contingent upon our understanding of ringdown’s underlying mathematical
structure Berti  _et al._ (2003); Berti and Klein (2014); Berti  _et al._
(2006b); London (2020); London  _et al._ (2014); London (2023); London and
Fauchon-Jones (2019); Husa  _et al._ (2016); Hamilton  _et al._ (2021, 2023);
Cotesta  _et al._ (2020); Pompili  _et al._ (2023); Ossokine  _et al._ (2020);
Hughes  _et al._ (2019); Lagos and Hui (2022).

The term “ringdown” specifically refers to linear combinations of BH Quasi-
Normal Modes (QNMs) Press (1971); Leaver (1985); Press and Teukolsky (1973).
Thus the mathematical structure of ringdown is that of the QNMs. Due to their
astrophysical relevance, the QNMs of Kerr BHs are perhaps the most studied
Berti  _et al._ (2003); Berti and Klein (2014); Berti  _et al._ (2006b);
London (2020); London  _et al._ (2014); London (2023); London and Fauchon-
Jones (2019); Husa  _et al._ (2016); Hamilton  _et al._ (2021, 2023); Cotesta
_et al._ (2020); Pompili  _et al._ (2023); Ossokine  _et al._ (2020); Hughes
_et al._ (2019); Ghosh  _et al._ (2018); Yang  _et al._ (2012, 2013);
Zimmerman and Mark (2016); Lagos and Hui (2022).

The Kerr QNMs emerge within the study of Einstein’s equations linearized
around the Kerr metric Vishveshwara (1970); Press and Teukolsky (1973); Leaver
(1985). Within that topic, a central result is that each single QNM solution
is a function of time (a damped sinusoid), times a function of spherical polar
angles (a spheroidal harmonic), and times a function of radial distance (a
confluent Heun function) Press and Teukolsky (1973); Dariescu  _et al._
(2021); Fiziev (2010); Hatsuda (2021).

The time domain’s sinusoidal damping is determined by discrete QNM (or
ringdown) frequencies Leaver (1985). For gravitational perturbations of Kerr,
QNM frequencies are labeled by three integers: $\ell$, $m$ and $n$, where
$\ell\geq 2$, $|m|\leq\ell$ and $n\geq 0$ Berti  _et al._ (2006a). Typically,
$\ell$ is referred to as a polar index, $m$ an azimuthal index, and $n$ an
“overtone” index. Respectively, $\ell$, $m$ and $n$ originate from the polar,
azimuthal and radial dimensions, as defined in a BH centered frame Leaver
(1985).

For physically relevant scenarios, i.e. where the QNM frequencies depend on
$\ell$, the spheroidal harmonics have recently been shown to not only be
biorthogonal222Two sequences, $\\{a_{i}\\}_{i=0}^{\infty}$ and
$\\{b_{j}\\}_{j=0}^{\infty}$, are biorthogonal if
${\langle{a_{i}}\,|\,{b_{j}}\rangle}=\delta_{ij}$ Brauer (1964); Christensen
(2003). with the so-called adjoint-spheroidal harmonics, but also complete
London (2023). This conclusion is based on the fact that $\ell$ is directly
related to the eigenvalues of e.g. Jacobi polynomials Fackerell and Crossman
(2008). As a result, the spheroidal harmonics may be thought of as a basis of
special functions, capable of representing gravitational waves from arbitrary
sources.

This is well known to be true of the spherical harmonics, which are a special
case of the spheroidal harmonics. Consequently, spherical harmonics play a
variety of key roles in gravitational wave theory and data analysis, from the
development of perturbative approximations (e.g. Post-Newtonian theory and
gravitational wave self force), and the storage of simulated gravitational
waves (Numerical Relativity), to the evaluation of gravitational wave signal
models for use in parameter inference Newman and Penrose (1966); Goldberg  _et
al._ (1967); Blanchet (2014); Abbott  _et al._ (2016d); O’Sullivan and Hughes
(2014); Ruiz  _et al._ (2008); Breuer  _et al._ (1977). Analogous applications
for the spheroidal and adjoint-spheroidal harmonics have been proposed London
and Hughes (2022). Perhaps surprisingly, there has been comparatively little
research into whether the QNM’s radial functions are special functions of
potentially broad use in gravitational wave theory Green  _et al._ (2023);
Lagos and Hui (2022); Yang  _et al._ (2012); Leaver (1986).

In $1973$, Teukolsky derived the differential equation that the radial
functions must satisfy, and at the time noted that his radial equation has two
regular singular points and one irregular one, meaning that it cannot be
transformed into a hypergeometric-type equation for which there are many well
known analytic methods of solutions Teukolsky (1973). This appears to have
informed a branch of research which seeks to solve Teukolsky’s radial equation
using numerical methods, e.g. O’Sullivan and Hughes (2014); Hughes (2000a);
Taracchini  _et al._ (2014); Hughes  _et al._ (2019); Islam  _et al._ (2022);
Rifat  _et al._ (2019); Campanelli  _et al._ (2001); Krivan  _et al._ (1996);
Vishal  _et al._ (2023).

However, in 1985 Leaver presented what was perhaps the first Frobenius-type
series solution for the QNM’s radial functions Leaver (1985). Leaver later
pointed out that, although such series solutions may be locally convergent,
they may diverge at either the event horizon or at spatial infinity Leaver
(1986); Ronveaux and Arscott (1995). It wasn’t until 1996 that Mano, Suzuki
and Takasugi (MST) presented a kind of “trick” that allows for the analytic
matching between one series expansion valid at the BH horizon, and another
valid near spatial infinity Mano  _et al._ (1996); Sasaki and Tagoshi (2003).
Mostly used in the context of gravitational wave self-force and particle
perturbation theory, various developments upon the MST method have been
presented Isoyama  _et al._ (2022); Casals and Zimmerman (2019); Sasaki and
Tagoshi (2003); Fujita and Tagoshi (2004).

In parallel, there has been broad and growing interest in a generalization of
the hypergeometric equation called the Heun equation Ronveaux and Arscott
(1995); Chen and Jing (2023); Fiziev (2009a); Hortacsu (2012); Dariescu  _et
al._ (2021); Magnus  _et al._ (2021). Unlike the hypergeometric equation,
which has three regular singular points, the Heun equation has four Ronveaux
and Arscott (1995). It was perhaps first pointed out by Leaver and then Fiziev
that Teukolsky’s radial equation is an instance of the confluent Heun
equation, also known as a generalized spheroidal equation Fiziev (2009a).
Fiziev’s work appears to coincide with what is perhaps a recent and multi-
disciplinary interest in Heun equationsHortacsu (2012). For example, it has
very recently been observed that a careful consideration of the confluent Heun
function’s asymptotic behavior allow them to be applied directly to
astrophysically relevant scenarios in particle perturbation theory (e.g.
gravitational waves from extreme mass ratio inspirals) Chen and Jing (2023).

Together, methods for evaluating the QNM’s radial functions that are derived
from Teukolsky, MST and Heun underlay a vast (and growing) body of research.
However, thus far, this work appears to have not fully addressed the following
questions:

  1. (_Q_ 1)

The properties of many special functions in physics may be framed using ideas
from Sturm-Liouville theory and functional analysis Arfken  _et al._ (2013);
Kristensson (2010); Courant and Hilbert. (1954); Adkins and Weintraub (2012).
There, solutions to differential equations are conceptualized as eigenvectors,
where the vector scalar products (i.e. symmetric nondegenerate bilinear forms
Adkins and Weintraub (2012)) are defined by an integral whose kernel is
weighted such that the differential equations is formally self-adjoint. For
example, the scalar product for the QNM’s angular functions can be framed as
that of the Jacobi polynomials (i.e. with the Jacobi weight function), and the
Jacobi differential equation is self-adjoint with respect to that scalar
product Fackerell and Crossman (2008). In analogy:

Is there an appropriate scalar product for the QNM’s radial functions?

  2. (_Q_ 2)

Classical polynomials (i.e. those of Hermite, Laguerre and Jacobi) are known
to be solutions to special cases of the hypergeometric equation Abramowitz and
Stegun (1964). It is well known that the QNM angular functions are underpinned
by spin-weighted spherical harmonics, which are themselves underpinned by
Jacobi polynomials Fackerell and Crossman (2008). Therefore, there is good
reason to wonder whether an analogous conclusion be developed for the QNM’s
radial functions:

Are the QNM’s radial functions underpinned by classical polynomials, or some
generalization thereof? Equivalently, is the overtone label, i.e. $n$,
directly related to the order of a complete and orthonormal polynomial
sequence?

Regarding (_Q_ 1), many studies have encountered the need to perform integrals
over solutions to Teukolsky’s radial equation, e.g. Green  _et al._ (2023);
Motl and Neitzke (2003); Hughes (2000a); Leaver (1986); Andersson (1997);
Jaramillo  _et al._ (2021). Perhaps most relevant to the present article,
Leaver (Ref. Leaver (1986)) studied Schwarzschild (i.e. non-spinning) BHs, and
used a Green’s function approach to derive various spectral-type expansions,
including those for QNMs. Similar approaches have recently been used to study
non-linear (particularly quadratic) QNMs of Schwarzschild, e.g. Lagos and Hui
(2022). This branch of research intersects with (_Q_ 1), however it does not
directly draw from Sturm-Liouville theory due to apparent divergences in QNM-
type problems that are not present in classical resonance
problems333Similarly, Ref. Nollert and Price (1999) implies that the QNM
problem may not be compatible with Sturm-Liouville theory because QNMs damp in
the time domain. In Sec. IV of this article, this is shown to not be the case.
Leaver (1986). Interestingly, Ref. Leaver (1986) observes that exactly such
divergence may be overcome by allowing for a radial integration path that is
complex valued. In this article, the same idea is encountered for Kerr, and is
motivated by coordinate dependence.

In an adjacent branch of research, many studies have made use of what’s called
the “formal adjoint” of operators related to QNMs where, in effect, a weight
function of $1$ is considered, e.g. Teukolsky and Press (1974); Wald (1979);
Green  _et al._ (2023); Dias  _et al._ (2022); Fiziev (2009a); Jaramillo  _et
al._ (2021). From that perspective, Teukolsky’s radial equation is not self-
adjoint due to the weight function or scalar product used. Consequently, it
appears that an investigation of Teukolsky’s radial equation in the context of
Sturm-Liouville theory with a non-trivial weight has not yet been pursued.

Regarding (_Q_ 2), a number studies have attempted to use classical
polynomials to represent solutions to the QNM’s radial functions, e.g. Chung
_et al._ (2023a, b); Ghojogh  _et al._ (2023). These efforts encounter the
following mathematical difficulty: the asymptotic (i.e. large distance from
the BH) behavior of the classical polynomials is incompatible with the QNM’s
radial functions, specifically the constraint that gravitational radiation
from QNMs must be purely outgoing towards spatial infinity Teukolsky (1973);
Press and Teukolsky (1973); Teukolsky and Press (1974). As a result, one
encounters a “generalized spectral problem”, which may result in the
appearance of non-physical solutions Chung  _et al._ (2023a). In contrast, a
purely spectral formulation would result in a matrix representation of the
problem that is e.g. tridiagonal, or generally band-diagonal, and has no
artificial eigensolutions Cook and Zalutskiy (2014); Fackerell and Crossman
(2008). Although it has been recognized that Teukolsky’s radial equation is a
confluent Heun equation, it appears that confluent Heun polynomials have not
been used to frame the QNM problem as purely spectral (until this work).

In 2014, Cook and Zalutskiy studied the connection between “algebraically
special” QNMs of nonspinning BH’s, and “total transmission” QNMs in the
presence of BH spin Cook and Zalutskiy (2014). There, it was found that for
special QNM frequencies, the associated radial functions can be directly
identified with select confluent Heun polynomials. With (_Q_ 2) in mind, there
is good reason to ask whether the confluent Heun polynomials are only relevant
to the total transmission QNMs, or whether there exists a single formulation
that connects confluent Heun polynomials to all QNM solutions?

This article provides a first-principles development of a radial scalar
product for QNMs, thereby affirmatively answering (_Q_ 1). Towards (_Q_ 2),
this article applies the scalar product to confluent Heun polynomials, where
it is shown that the polynomials may be defined for any wave frequency,
including the complex valued QNM frequencies. While a full treatment of (_Q_
2) is determined to be outside of this article’s scope, significant progress
is provided towards an affirmative answer. In particular, this article
concludes by showing that Teukolsky’s radial equation is, at least in
principle, exactly tridiagonalizable. To this end, the existence of
potentially new “canonical confluent Heun polynomials” is conjectured.

##  II Scope and organization

###  II.1 Scope

This article discusses QNM solutions to Teukolsky radial equation in the
context of isolated Kerr BHs, absent any effects beyond GR. Some aspects of
the discussion make minimal physical assumptions, while others are specific to
gravitational wave ringdown. Broadly, the reader is cautioned that the author
does not aspire to the language of formal mathematical rigor (e.g. the text
will not include theorems and lemmas).

Most mathematical arguments make no assumption about spin weight, $s$.
However, $s\in\\{-2,+2\\}$ are given particular attention due to their
correspondence with gravitational radiation. Further, Teukolsky’s radial
equation, even when formatted for the QNMs, may be considered for any
frequency. Therefore, much of this article’s content is not specific to the
QNM’s frequencies, but rather the structure of Teukolsky’s radial equation
when formatted with QNM boundary conditions.

This article’s central result is a frequency dependent scalar product;
however, it is also noted in Sec. IV.2 that this is but one of a number of
equivalent formulations. For example, a frequency independent scalar product
is shown in Eq. (31), and its relation to this article’s frequency dependent
scalar product is described in adjacent text. That scalar product may be of
future interest to an investigation analogous to Ref. London (2023), where
different QNM frequencies are studied simultaneously to learn about the
completeness and biorthogonality properties of the QNM’s spheroidal harmonics.

Throughout this article, QNM frequencies are used to evaluate numerical
examples due to their relevance to gravitational wave theory and experiment.
QNM frequencies for this article have been generated using Leaver’s method via
the publicly available positive and qnm Python packages London  _et al._
(2020); Stein (2019). All numerical results presented have been computed with
$64$ significant figures of precision via the mpmath package mpmath
development team (2023).

###  II.2 Organization

Section III provides an overview of Kerr QNMs and Teukolsky’s radial equation.
Section IV discusses a simplified version of Teukolsky’s radial equation in
the context of Sturm-Liouville theory, and uses ideas therefrom to derive a
scalar product for the Kerr QNMs. Section V focuses on two complementary
numerical methods (direct integration and analytic continuation) for
evaluating the scalar product. Section VI presents the limitations of direct
integration and analytic continuation, and then compares and contrasts them.

Section VII presents two pedagogical applications of the scalar product.
There, a pedagogical construction of the confluent Heun polynomials is
provided, along with a derivation their most relevant properties. Section VIII
discusses the implications of the scalar product and confluent Heun
polynomials on Teukolsky’s radial equation. There, a class of “canonical
confluent Heun polynomials” is proposed, and used to argue that Teukolsky’s
radial equation may be exactly tridiagonalizable.

Finally, in Sec. IX, further implications of the scalar product and related
results are discussed in the context of future work.

##  III Teukolsky’s equation with QNM boundary conditions

In this section we connect the observable gravitational wave polarizations
with this article’s topical starting point, Teukolsky’s radial equation with
QNM boundary conditions. Along the way, we will encounter spin weighted fields
Newman and Penrose (1966, 1962), the definition of QNMs Nollert (1999); Berti
_et al._ (2016), the importance of a scalar product and polynomials to
Teukolsky’s angular equation, and how QNM solutions to Teukolsky’s radial
equation may be framed in close analogy with the angular case Leaver (1985).
These ideas converge on Teukolsky’s radial equation in a simplified form Cook
(2020); Fiziev (2010); Chen and Jing (2023). That equation sets the stage for
our application of Sturm-Liouville theory Pinchover and Rubinstein (2005);
Abramowitz and Stegun (1964); Arfken  _et al._ (2013); Courant and Hilbert.
(1954), and the development of a radial scalar product in Sec. IV.

###  III.1 From gravitational wave strain to the perturbative wave quantity

Any gravitational wave detector will be sensitive to the two strain
polarizations predicted by GR, $h_{+}$ and $h_{\times}$ Ruiz  _et al._ (2008);
Thorne (1980). The study of gravitational wave plane waves motivates an
extremely useful complex valued strain given by

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

In Eq. (1a), $h$ is the complex valued gravitational wave strain Ruiz  _et
al._ (2008). We will consider it to depend on Boyer–Lindquist coordinates,
where time, radius, polar angle and azimuthal angle are denotes by $t$, $r$,
$\theta$ and $\phi$ respectively Boyer and Lindquist (2004). In Eq. (1b),
$\psi_{4}$ is the $5$’th Weyl scalar in Newman-Penrose notation Ruiz  _et al._
(2008); Newman and Penrose (1962). In that formalism, $\psi_{4}$ is a complex
valued scalar field of spin weight $s=-2$ Newman and Penrose (1966); Teukolsky
(1973). Henceforth, complex valued qualities encode or modify the relative
importance of $h_{+}$ and $h_{\times}$ within $h$. Within the study of
gravitationally perturbed isolated BHs, two Weyl scalars are of particular
interest: $\psi_{4}$ ($s=-2$), and the closely related $\psi_{0}$ ($s=+2$).
Both encode information about gravitational waves Teukolsky (1973).

For linear perturbations of Kerr, $\psi_{4}$ is treated in terms of a generic
wave quantity, $\psi$:

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

In Eq. (2), the quantity scaling $\psi$ is a spin coefficient Teukolsky
(1973), where $a=|J|/M$ is the BH spin magnitude. We will hold that $G=c=1$,
leaving $a$ to have units of the black hole mass, $M$.

The use of $\psi$ is that it satisfies exactly the same partial differential
equation as other physically relevant wave quantities. In the case of the
first Weyl scalar, $\psi_{0}$, we simply have that $\psi=\psi_{0}$ Newman and
Penrose (1962). Whether $\psi$ refers to $\psi_{0}$ or $\psi_{4}$ is toggled
by the spin weight of the field, $s$. Thus for $\psi_{0}$ ($s=+2$) or
$\psi_{4}$ (s=-2), $\psi$ satisfies Teukolsky’s equation, which is one way of
representing Einstein’s equations linearized about the Kerr solution Hughes
(2000b); Teukolsky (1973, 1972); Teukolsky and Press (1974); Press and
Teukolsky (1973); Mino  _et al._ (1997).

###  III.2 Teukolsky’s equation and its QNM solutions

If we denote the partial differential operator for Teukolsky’s equation as
$\mathcal{L}_{tr\theta\phi}$, with

$\displaystyle\mathcal{L}_{tr\theta\phi}\;=\;\mathcal{L}_{tr\theta\phi}(s,r,\theta,\partial_{t},\partial_{r},\partial_{\theta},\partial_{\phi}),$
(3)

then the schematic form of Teukolsky’s equation is

$\displaystyle\mathcal{L}_{tr\theta\phi}\;\;\psi(t,r,\theta,\phi)\;=\;4\pi\Sigma
T\;.$ (4)

In Eq. (3),
$\\{\partial_{t},\partial_{r},\partial_{\theta},\partial_{\phi}\\}$ denote
differentiation with respect to the the Boyer-Lindquist coordinates Boyer and
Lindquist (2004). An explicit definition of $\mathcal{L}_{tr\theta\phi}$ is
provided in Appx. (A). In Eq. (4), $\Sigma$ is the usual quantity that appears
in the Kerr metric,

$\displaystyle\Sigma\equiv r^{2}+a^{2}\cos^{2}\theta\;,$ (5)

and $T$ is a source term derived from the stress-energy tensor. Since we are
interested in isolated BHs, we have that

$\displaystyle T(t,r,\theta,\phi)\;=\;0\;.$ (6)

The fact that $\mathcal{L}_{tr\theta\phi}$ (see Eq. 3) does not explicitly
depend on $t$ and $\phi$ means that $\psi$ of the form

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

may separate444The ansatz provided in Eq. (7) obscures the fact that our net
solution will depend on the tortoise coordinate via the retarded time Press
and Teukolsky (1973); Teukolsky (1973). The radial part of the retarded time
is ultimately encoded within $R(r)$. $\mathcal{L}_{tr\theta\phi}$ into two
coupled equations for $R(r)$ and $S(\theta)$. In Eq. (7), $\tilde{\omega}$ is
generally a complex valued frequency:

$\displaystyle\tilde{\omega}\;=\;\omega\,-\,i/\tau\;.$ (8)

In Eq. (8), one should think of any physical situation in which an isolated BH
experiences a momentary and small perturbation. The post-merger of two
sufficiently compact objects is an example of particular importance to
astrophysics Helfer  _et al._ (2017); Abbott  _et al._ (2021b, a); Baker  _et
al._ (2006). For the simplest scenarios, the resulting gravitational radiation
will, qualitatively, ring down Nollert (1999); Press (1971); Berti  _et al._
(2007); Kokkotas and Schmidt (1999); Nollert and Price (1999). The central
frequency of that ringing will be $\omega$, and the exponential decay rate of
that ringing will be $1/\tau$, where $\tau$ is strictly positive; otherwise,
Eq. (7) would diverge at late times, and the BH would be linearly unstable
Vishveshwara (1970); Press and Teukolsky (1973); Press (1971); Whiting (1989).

In general, gravitational radiation will be a sum of such mode-like terms (Eq.
7) weighted by constant factors that depend on spacetime geometry and the
perturber Berti and Cardoso (2006); Kamaretsos  _et al._ (2012); London  _et
al._ (2014); London (2020). The net spatial dependence of each term in the sum
defines a single QNM Berti  _et al._ (2007); Kokkotas and Schmidt (1999);
Nollert and Price (1999). There will always be a period in time for which QNM
radiation dominates over other solutions to the homogeneous radial equation,
such as power-law tails, prompt emission, and nonlinear QNMs Leaver (1986);
Lagos and Hui (2022); Andersson (1997). This inherently poorly defined region
of QNM dominance is colloquially referred to as ringdown Press (1971); Nollert
(1999).

Properties of QNM solutions can be determined by applying the solution ansatz
defined in Eq. (7), and then studying the resulting coupled eigenvalue
problems. These are separate radial and angular differential equations,
connected by the separation constant, $A$:

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

Only the schematic structure of Eqs. (9a) and (9b) are important here, but
reader may find full expressions for $\mathcal{L}{{}_{u}}$ and
$\mathcal{L}{{}_{r}}$ in Appx. (A). Henceforth, we will refer to Eq. (9a) as
the angular equation, and Eq. (9b) as the radial equation. Respectively,
angular problem and radial problem will refer to the solving of each
eigenvalue equation.

In Eq. (9), $u=\cos(\theta)$, $\mathcal{L}{{}_{u}}$ will be referred to as the
angular differential operator, and $\mathcal{L}{{}_{r}}$ will be referred to
as the radial differential operator. The angular functions, $S(u)$, are called
spheroidal harmonics, and we will refer to $R(r)$ as Teukolsky’s radial
functions. Both are known to be closely related to the confluent Heun
functions555As noted by Leaver, they are also called generalized spheroidal
wave functions Leaver (1986). Fiziev (2010); Cook and Zalutskiy (2014); Cho
_et al._ (2009); Chen and Jing (2023).

###  III.3 Lessons from Teukolsky’s angular problem: the roles of the scalar
product and Jacobi polynomials

Although both angular and radial equations are known to be similar in form
(e.g. Ref. Leaver (1986)), there are good reasons to first study the angular
equation, and then apply lessons therefrom to the radial equation. Relative to
the radial equation, the angular equation is well known to be simpler; as a
result, there are many well known and highly accurate methods for computing
its global666i.e. series solutions that apply across the physical domain
analytic solutions Leaver (1985); Cook and Zalutskiy (2014); O’Sullivan and
Hughes (2014); Berti  _et al._ (2006b); Yang  _et al._ (2012). Pertinent to
this article is the spectral method presented in Ref. Fackerell and Crossman
(2008), and whether two key concepts therein may be applied to the radial
problem.

For the convenience of the reader, a modern and self-contained presentation
this method is provided in Appx. (B). There, the two key concepts are:

  * (i)

The structure of the differential equation allows the development of a scalar
product that is natural for the solution space. This concept motivates (_Q_
1).

  * (ii)

The structure of the differential equation points the way to problem specific
special functions, namely the Jacobi polynomials. In turn, the Jacobi
polynomials may be used to exactly represent analytic solutions. This concept
motivates (_Q_ 2).

If applicable to the radial problem, concept (i) may enable greater
understanding of QNM (spatial) orthogonality, bi-orthogonality, and
completeness London (2023). Similarly, concept (ii) may yield a representation
of the radial functions that may not only be computationally efficient, but
may also lend clear insight into the mathematical meaning of the so-called
overtone label, $n$. For these reasons, concepts (i) and (ii) will guide
forthcoming discussion of the radial equation.

###  III.4 QNM solutions of Teukolsky’s radial equation

The study of BH QNMs is traditionally contained within the study of spacetime
exterior to the event horizon Regge and Wheeler (1957); Vishveshwara (1970);
Press and Teukolsky (1973). In this context, QNM boundary conditions are
defined on two asymptotic surfaces: the event horizon and spatial infinity
Teukolsky (1973); Press and Teukolsky (1973); Teukolsky and Press (1974). This
makes it computationally and algebraically useful to consider a compactified
radial coordinate Leaver (1985, 1986),

$\displaystyle\xi\;$ $\displaystyle=\;(r-r_{+})/(r-r_{-})\;.$ (10)

In Eq. (10), $r_{+}=M+\sqrt{M^{2}-a^{2}}$ is the outer Kerr horizon (i.e. the
event horizon), and $r_{-}=M-\sqrt{M^{2}-a^{2}}$ is the inner Kerr horizon.
The compactified coordinate, $\xi$, has the effect of mapping $r=r_{+}$ to
$\xi=0$, $r=r_{-}$ to $\xi\rightarrow-\infty$, and $r\rightarrow\infty$ to
$\xi=1$.

The QNM boundary conditions are that the gravitational radiation is purely
ingoing (towards the BH) at $\xi=0$ and purely outgoing near $\xi=1$. A
derivation of these conditions is provided in e.g. Refs.Press and Teukolsky
(1973); Cook and Zalutskiy (2014). Since these conditions are imposed between
$r=r_{+}$ and $r\rightarrow\infty$, we will only consider
$r\in[r_{+},\infty)$, and so $\xi\in[0,1)$.

The QNM boundary conditions are asymptotic in the sense that they are
restrictions on the functional form of $\psi$ near the event horizon and
spatial infinity. As used by Leaver, these boundary conditions amount to a
similarity transformation on the radial problem Ronveaux and Arscott (1995);
Leaver (1985, 1986). This transformation proceeds by seeking solutions to the
radial problem of the form

$\displaystyle R(r(\xi))\;$ $\displaystyle=\;\mu(\xi)f(\xi)\;.$ (11)

In Eq. (11), $\mu(\xi)$ is a product of functional forms required by the
asymptotic boundary conditions,

$\displaystyle\mu(\xi)\;$
$\displaystyle=\;e^{\frac{2i\delta\tilde{\omega}}{1-\xi}}(1-\xi)^{1+2(s-iM\tilde{\omega})}$
(12)
$\displaystyle\quad\quad\times\xi^{-(iM\tilde{\omega}+s)+\frac{i\left(am-2M^{2}\tilde{\omega}\right)}{2\delta}}\;.$

The factors that make up $\mu$ are determined by structure of Frobenius
(series) solutions to Eq. (9b). In Eq. (12), the argument of the exponential,
as well as the powers of $\xi$ and $1-\xi$ are selected by the physical
boundary conditions Cook and Zalutskiy (2014); Fiziev (2010); Leaver (1986).
In turn, Eq. (11) strictly requires that $f(\xi)$ be non-zero at $\xi=0$ and
$\xi=1$; otherwise the scalings imparted by $\mu$ (i.e. the QNM boundary
conditions) would be spoiled. This requirement makes $f(\xi)$ the natural
quantity for which to seek a power series solution.

The application of Eq. (11) to radial eigenproblem (Eq. 9b) has the effect of
transforming $\mathcal{L}_{r}$ into a new differential operator,
$\mathcal{L}_{\xi}$,

$\displaystyle\left[\mu(\xi)^{-1}\,\mathcal{L}{{}_{r}}\,\mu(\xi(r))\right]f(\xi)\;$
$\displaystyle=\;\left[\mu(\xi)^{-1}A\,\mu(\xi)\right]f(\xi)\;$ (13a)
$\displaystyle\mathcal{L}{{}_{\xi}}\,f(\xi)\;$ $\displaystyle=\;A\,f(\xi)\;.$
(13b)

In Eq. (13a), we have applied the solution ansatz given by Eq. (11) to the
radial eigenproblem, Eq. (9b). We have then divided by $\mu$. Since $\mu$ is
closely related to the singular exponents used in Frobenius’ method, the
division is exact Leaver (1985). The result, Eq. (13b), is a transformed
eigenproblem that has the same eigenvalues, $A$. It is in this sense that
$\mathcal{L}_{\xi}$ is formally similar to $\mathcal{L}_{r}$ Axler (2015); Lax
(2002); Christensen (2003),

$\displaystyle\mathcal{L}{{}_{\xi}}\;$
$\displaystyle=\;\mu(\xi)^{-1}\,\mathcal{L}{{}_{r}}\,\mu(\xi(r))\;.$ (14)

While we refer to the effect of Eqs. (13) and (14) as a similarity
transformation, it is sometimes referred to as an elementary power
transformation, or a homotopic transformation Ronveaux and Arscott (1995);
Chen and Jing (2023); Cho (2020). The reader should note that Eqs. (13a) and
(14) contain a coordinate transformation ($r$ to $\xi$), as well as a
similarity transformation.

In Eq. (14), the new differential operator, $\mathcal{L}_{\xi}$, has the form
given by Eq. (15).

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

There, the constants, $\text{C}_{0}$ to $\text{C}_{4}$, are defined in Eq.
(16) using the physical parameters discussed previously.

$\displaystyle\delta\;$ $\displaystyle=\;\sqrt{M^{2}-a^{2}}$ (16a)
$\displaystyle\text{C}_{0}\;$
$\displaystyle=\;-2am\tilde{\omega}-2i\tilde{\omega}(-\delta+M(2s+1))$ (16b)
$\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$
(16c)
$\displaystyle\quad\;\;-(4M\tilde{\omega}+i)\frac{\left(am-2M^{2}\tilde{\omega}\right)}{\delta}\;$
$\displaystyle\text{C}_{2}\;$ $\displaystyle=\;4i\delta\tilde{\omega}$ (16d)
$\displaystyle\text{C}_{3}\;$
$\displaystyle=\;-2(s+1)+4i\tilde{\omega}(M-\delta)$ (16e)
$\displaystyle\text{C}_{4}\;$
$\displaystyle=\;s+3-6iM\tilde{\omega}+\frac{i\left(am-2M^{2}\tilde{\omega}\right)}{\delta}$
(16f)

Equation (15) has a few special properties. We have started with Boyer-
Lindquist coordinates, but one also arrives at Eq. (15) if one instead begins
with Kerr Ingoing or Kerr Outgoing coordinates Teukolsky and Press (1974). One
also arrives at Eq. (15) (albeit with $s\rightarrow-s$ and $A\rightarrow
A-2s$) if one uses the alternative tetrad choice discussed in Ref. Teukolsky
and Press (1974), which amounts to $t\rightarrow-t$ and
$\phi\rightarrow-\phi$. As noted in e.g. Refs. Teukolsky and Press (1974);
Whiting (1989), the use of either $s=2$ or $s=-2$ results in mathematically
different descriptions of the same gravitational radiation.

Equation (15) has regular singular points at $\xi\in\\{-\infty,0\\}$ and one
irregular singular point at $\xi=1$. Thus Eq. (15) is an instance of the
confluent Heun differential operator Ronveaux and Arscott (1995). This is also
true in the Schwarzschild limit. While it is sometimes said that Ref. Leaver
(1986) treats Eq. (15) as what is called a generalized spheroidal operator, we
note here that this is a slight misnomer Ronveaux and Arscott (1995); Hortacsu
(2012). Equation (15) may only be transformed into a generalized spheroidal
operator if one applies an asymptotic boundary condition at $r=r_{-}$ that is
inconsistent with QNMs Leaver (1985); Cook and Zalutskiy (2014). Appendix (C)
provides an expanded disambiguation between different forms of the confluent
Heun equation and the radial equation.

Equations (13b-16) are commonly used to define series solutions for $f(\xi)$.
The requirement that these series solutions converge on $\xi\in[0,1)$ happens
to be equivalent to finding eigenvalues, $A$, such that certain continued
fraction are zero Leaver (1985, 1986); Mano  _et al._ (1996). In particular,
see Ref.Leaver (1986) for a proof of convergence for series solutions, example
applications of continued fractions to the radial and angular eigen-problems,
and how the eigenvalues may be determined by searching for roots of continued
fractions.

From the use of continued fractions, values of $A$ are well known to be
discrete, and by convention are labeled by the size of their imaginary part:
$n=0$ corresponds to the eigenvalue with the smallest value of $|1/\tau|$, and
$n=1$ has the next-smallest value Berti  _et al._ (2016); Nollert (1999);
Kokkotas and Schmidt (1999). However, the use of continued fractions has
arguably obscured deeper understanding of these eigenvalues, and their related
solutions. In search of that understanding, it appears that there has been
relatively little application of Sturm-Liouville theory to Eq. (15) in the BH
context Leaver (1985); Nollert (1999); Berti and Cardoso (2006); Ripley
(2022); Lin and Qian (2023); Press (1971); Motl and Neitzke (2003); Ronveaux
and Arscott (1995). Nevertheless, since the QNM eigenvalues are well known to
be unique for most physical scenarios 777There is some ambiguity in how one
defines the eigenvalues. In particular, the fact that $\text{C}_{0}$ is non-
zero means that the eigenvalue that one starts with, $A$, need not strictly
define the eigenvalue that is used in Eq. (9a). Thus in e.g. the limit that
$a=0$, where $A=\ell(\ell+1)-s(s+1)$ for every radial eigenfunction, the
effective eigenvalue, $A-\text{C}_{0}$, may be a unique function of
$\tilde{\omega}$. (e.g. London and Fauchon-Jones (2019); Cook (2020); Seidel
(1989); Yang  _et al._ (2012)), the QNM’s themselves are necessarily linearly
independent888More precisely, finite sets of QNMs will be linearly independent
when their eigenvalues are unique. Infinite sets of QNMs will be minimal if
their eigenvalues remain distinct as some order label on the set goes to
infinity Christensen (2003); London (2023)., and thus well suited to Sturm-
Liouville theory Courant and Hilbert. (1954).

##  IV Definition of a radial scalar product

The aim of Sturm-Liouville theory is to understand the existence and
asymptotic behavior of eigenvalues and eigenfunctions Pinchover and Rubinstein
(2005); Arfken  _et al._ (2013); Courant and Hilbert. (1954). Here, we will
focus on a facet of the theory that pertains to whether there exists a scalar
product and an associated weight function, for which $\mathcal{L}_{\xi}$ is
self-adjoint.

###  IV.1 The scalar product and the weight function

By scalar product, we mean a symmetric bilinear form that takes two non-
singular functions999The reader should not confuse the non-singular but
otherwise arbitrary function, $\mathrm{a}(\xi)$, with the BH spin, $a$.,
$\mathrm{a}(\xi)$ and $\mathrm{b}(\xi)$, and performs an operation on them
that results in a number Adkins and Weintraub (2012); Gohberg  _et al._
(2006). In less technical parlance, we simply mean a bra-ket operation that is
similar to what one encounters in quantum mechanics:

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

In Eq. (17), $\text{W}(\xi)$ is a called a measure, or weight function. It
will soon be the focus of our discussion.

Since there is no conjugation involved in Eq. (17), the scalar-product shown
is not conjugate-symmetric. This feature distinguishes it from an inner-
product, and allows analysis of $\mathcal{L}_{\xi}$, irrespective of its
complex valued coefficients. The analysis of immediate interest pertains
whether $\text{W}(\xi)$ may be defined such that $\mathcal{L}_{\xi}$ is
formally self-adjoint Pinchover and Rubinstein (2005); Courant and Hilbert.
(1954).

The adjoint of $\mathcal{L}_{\xi}$ is defined to be an operator,
${\mathcal{L}_{\xi}}^{\dagger}$, such that

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

Chapter $5$ Section 3 of Ref. Courant and Hilbert. (1954) provides a detailed
overview of how $\mathcal{L}_{\xi}$ may be explicitly defined using
integration by parts. In that derivation, the fact that the coefficient of the
second derivative term in Eq. (15) vanishes at $\xi=0$ and $\xi=1$ is key to
Sturm-Liouville boundary conditions being satisfied Pinchover and Rubinstein
(2005). There, it is also shown that $\mathcal{L}_{\xi}$ is formally self-
adjoint, i.e.

$\displaystyle\mathcal{L}{{}_{\xi}}\;$
$\displaystyle=\;{\mathcal{L}{{}_{\xi}}}^{\dagger}\;,$ (19)

if it also holds that

$\displaystyle\mathcal{L}_{\xi}\;$
$\displaystyle=\;p_{0}(\xi)\;+\;\text{W}(\xi)^{-1}\,\partial_{\xi}\left[\,(\text{W}(\xi)\,p_{2}(\xi))\,\times\,\partial_{\xi}\,\right]$
(20a)
$\displaystyle=\;p_{0}(\xi)\;+\;\frac{\partial_{\xi}(\text{W}(\xi)\,p_{2}(\xi))}{\text{W}(\xi)}\,\partial_{\xi}\;+\;p_{2}(\xi)\partial^{2}_{\xi}\;.$
(20b)

In Eq. (20), $p_{0}(\xi)$ and $p_{2}(\xi)$ are the coefficient functions seen
in zero’th and second derivative terms of Eq. (15), respectively. If we define
$p_{1}$ to be the coefficient function of $\partial_{\xi}$ in Eq. (15), then

$\displaystyle p_{0}\;$ $\displaystyle=\;\text{C}_{0}+\text{C}_{1}(1-\xi)\;,$
(21a) $\displaystyle p_{1}\;$
$\displaystyle=\;\text{C}_{2}+\text{C}_{3}(1-\xi)+\text{C}_{4}(1-\xi)^{2}$
(21b) $\displaystyle p_{2}\;$ $\displaystyle=\;\xi(\xi-1)^{2}\;.$ (21c)

In Eq. (20a), the appearance of the weight function is a direct result of Eqs.
(18) and (19). To go from Eq. (20a) to Eq. (20b), we have applied the product
rule (a.k.a. Leibniz’s rule), to
$\partial_{\xi}\,\left[(\text{W}p_{2})\times\partial_{\xi}\right]$. For
$\mathcal{L}_{\xi}$ in Eq. (15) to be equivalent to that in Eq. (20b), it must
be the case that coefficients of $\partial_{\xi}$ match,

$\displaystyle p_{1}\;$
$\displaystyle=\;\text{W}^{-1}\,\partial_{\xi}(\text{W}\,p_{2})$ (22a)
$\displaystyle\text{W}^{-1}\,\partial_{\xi}\text{W}\;$
$\displaystyle=\;p_{2}^{-1}(p_{1}-\partial_{\xi}p_{2})\;.$ (22b)

In Eq. (22a) we have simply equated the first derivative coefficient of Eq.
(20b) with that of Eq. (15). In Eq. (22b), we have used Leibniz’s rule to
expand the right-hand-sid of Eq. (22a), and then we have rearranged terms to
emphasize that Eq. (22b) is a linear differential equation whose left-hand-
side is equal to $\partial_{\xi}\ln\text{W}$.

Equation (22b) can be solved for W using the coefficient functions defined in
Eq. (21),

$\displaystyle\text{W}(\xi)\;$
$\displaystyle=\;\,\frac{1}{p_{2}(\xi)}\,\exp\left(\int^{\xi}\frac{p_{1}(\xi^{\prime})}{p_{2}(\xi^{\prime})}\,d\xi^{\prime}\,\right)$
(23a)
$\displaystyle=\;e^{\frac{\text{C}_{2}}{1-\xi}}\;(1-\xi)^{-\text{C}_{2}-\text{C}_{3}-2}\;\xi^{\text{C}_{2}+\text{C}_{3}+\text{C}_{4}-1}\;.$
(23b)

In Eq. (23a), the integration constant has been chosen to be zero, and in Eq.
(23b) we explicitly evaluate the integral using Eqs. (21b) and (21c).

It will soon be useful to simplify our notation for $\mathrm{W}(\xi)$ by
defining a shorthand for combinations of parameters,

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

Consequently, the weight function takes a less cumbersome form,

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

and the scalar product between two any two arbitrary functions,
$\text{a}(\xi)$ and $\text{b}(\xi)$, may now be written as,

$\displaystyle{\langle{\mathrm{a}}\,|\,{\mathrm{b}}\rangle}\;=\;\int_{0}^{1}\,\mathrm{a}(\xi)\,\mathrm{b}(\xi)\;\;\xi^{\text{B}_{0}}(1-\xi)^{\text{B}_{1}}e^{\frac{\text{B}_{2}}{1-\xi}}\;d\xi\;.$
(26)

###  IV.2 Properties and non-uniqueness of the weight function

With Eq. (26), we arrive at one of this article’s key results. Equation (26)
along with Eqs. (16) and (24) define a scalar product for which Teukolsky’s
radial equation, Eq. (15), is self-adjoint. In turn, self-adjointness, in
principle, implies that as long as the eigenvalues may be defined, and that
they are unique, then the eigenfunctions are orthogonal, and very likely
complete Ronveaux and Arscott (1995). However, there are good reasons why the
explicit demonstration of orthogonality and completeness are beyond the scope
of this article. The most pressing reason has to do with the fact that
$\text{W}(\xi)$, naively, appears to diverge on the event horizon, and at
spatial infinity. Understanding that apparent divergence and its consequences
is preliminary to further investigation.

Visual inspection of Eq. (23) points to the origin of these divergences. In
Eq. (23a), the factor

$\displaystyle p_{2}(\xi)^{-1}\;=\;\xi^{-1}(1-\xi)^{-2}$ (27)

clearly creates a potential for the weight function to become singular at
$\xi=0$ and $\xi=1$. For this reason, the radial eigen-problem may be said to
be of the singular101010Note that e.g. the Jacobi and Legendre polynomial
problems are also singular in this way Arfken  _et al._ (2013). Sturm-
Liouville type Pinchover and Rubinstein (2005). In Eq. (23b), we see that the
net exponents of $\xi$ and $1-\xi$ may still prevent divergence if their total
real part is positive. Similarly, the exponential factor may diverge if the
real part of $\text{C}_{2}=4i\delta\tilde{\omega}$ is greater than zero.

Thus, to determine the physical conditions in which $\text{W}(\xi)$ diverges,
one might embark on an analysis of the sign of parameters within Eq. (23b) or,
equivalently, Eq. (25). While this is a somewhat cumbersome exercise, it is
easy to see that the exponential factor,

$\displaystyle
e^{\text{C}_{2}/(1-\xi)}\;=\;e^{4i\delta\tilde{\omega}/(1-\xi)}\;,$ (28)

will generally diverge at spatial infinity. For that we need to recall the
discussion below Eq. (8): Since the imaginary part of $\tilde{\omega}$ must be
negative for time-domain stability, the real part of
$4i\delta\tilde{\omega}/(1-\xi)$ is generally positive for the QNMs,

$\displaystyle\text{Re}\,\frac{4i\delta\tilde{\omega}}{1-\xi}\;=\;\frac{4\delta}{\tau(1-\xi)}\;>\;0\;.$
(29)

The inequality in Eq. (29) holds because $\delta=\sqrt{M^{2}-a^{2}}>0$ and
$\tau>0$. Thus it would appear that $\text{W}(\xi)$ will always diverge on at
least one boundary. This appearance happens to be just that.

Above, we have implicitly assumed that $\xi=0$ and $\xi=1$ are approached from
the real line. While that is a physically meaningful assumption, previous
studies have found it to be an unnecessary practical impediment Leaver (1986);
Green  _et al._ (2023); Motl and Neitzke (2003); Mano  _et al._ (1996). In
particular, when evaluating the scalar product, $\xi$ need not be a real
valued quantity. This point is explored in the two sections that follow. A
kind of “physicist’s analytic continuation” is detailed in Sec. V.1, and then
a complementary analytic continuation technique is presented in Sec. V.2.

We conclude this section by pointing out two alternative perspectives on the
weight function. The first hinges on the fact that Teukolsky’s equation allows
for the definition of a weight function prior to applying the QNM boundary
conditions. The second pertains to whether the weight function is unique, once
the QNM boundary conditions have been applied.

The first perspective is that the form of the weight function is closely
linked to the scaling function, $\mu(\xi)$, defined in Eq. (12). In
particular, it is easy show (see Appx. B) that

$\displaystyle\text{W}(\xi)\;$
$\displaystyle=\;\mu(\xi)^{2}\,\text{W}_{0}(\xi)\;.$ (30)

In Eq. (31), $\text{W}_{0}$ is the radial equation’s weight function before
applying the similarity transformation shown in Eq. (15),

$\displaystyle\text{W}_{0}(\xi)\;$
$\displaystyle=\;\xi^{s}\,(1-\xi)^{-2(s+1)}\;.$ (31)

This weight function is independent of $\tilde{\omega}$, $m$ and $a$, and it
is directly relevant for the original radial functions, $R(\xi)$ (see Eq. 9b).
In essence, Eq. (30) defines how the weight function changes under similarity
transformations on the main problem. One way to arrive at Eq. (30) is to note
that such transformations preserve the value of scalar products. While
$\text{W}_{0}$ may be of interest for future work, it is not directly
considered here because the corresponding differential operator does not
support polynomial solutions.

The second perspective is that $\text{W}(\xi)$ is not the only weight function
for which $\mathcal{L}_{\xi}$ is self-adjoint. It is left as an exercise for
the reader to show that one can actually enforce a strictly positive and real
valued weight function,

$\displaystyle\mathcal{W}(\xi)\;=\;|\xi^{\text{B}_{0}}(1-\xi)^{\text{B}_{1}}e^{\frac{\text{B}_{2}}{1-\xi}}|\;.$
(32)

This exercise follows exactly the steps discussed in Eqs. (20-23), but with
the requirement that the weight function is wrapped within an absolute value.
While $\mathcal{W}(\xi)$ is amenable to the direct integration approach that
we will discuss shortly, it is not at all clear to the author how it may be
used with the analytic continuation approach that will be discussed in Sec.
V.2. Further, there do not appear to be significant gains in using
$\mathcal{W}(\xi)$ rather than $\text{W}(\xi)$. Since complexification of the
radial coordinate is also needed for $\mathcal{W}(\xi)$, scalar products will
in general still be complex valued. $\mathcal{W}(\xi)$ is noted for potential
future interest. Henceforth, this article will focus on the evaluation of
scalar products with the comparatively simpler weight function,
$\text{W}(\xi)$.

|  |
---|---|---
Figure 1: Examples integration paths (white curves) against weight function
contours (colored contours) for the $(s,\ell,m,n)=(-2,2,2,3)$ QNM at three BH
spins: (left) $a/M=0.7$, $M\tilde{\omega}=0.4713-0.5843i$, (center)
$a/M=0.99$, $M\tilde{\omega}=0.8695-0.2058i$, and (right) $a/M=0.999$,
$M\tilde{\omega}=0.9558-0.07299i$. Contours for the real and imaginary parts
of the weight function,
$\text{W}(\xi)=\mathrm{Re}(\text{W})+i\,\mathrm{Im}(\text{W})$, are shown. For
$\mathrm{Re}(\text{W})$, colors between cyan and magenta represent
dimensionless values between $-50$ and $50$ respectively. For
$\mathrm{Im}(\text{W})$, colors between magenta and yellow represent values
between $-50$ and $50$ respectively. Convergence of contours at spatial
infinity ($\mathrm{Re}(\xi)=1$) corresponds to the divergence of the weight
function. Horizontal contours along $\mathrm{Re}(\xi)<0$ and
$\mathrm{Re}(\xi)>0$ are branch cuts. For left, center and right panels,
coordinate parameters, $(K_{0},K_{2})$, are $(0.8530,23.4621)$,
$(0.7105,1.4010\times 10^{-4})$ and $(0.1585,2.5088\times 10^{-4})$,
respectively. Dotted horizontal and vertical lines mark the respective
locations of the real and imaginary axes.

##  V Evaluation of scalar products

The scalar product developed in the previous section has been copied below,

$\displaystyle{\langle{\text{a}}\,|\,{\text{b}}\rangle}\;=\;\int_{0}^{1}\,\text{a}(\xi)\,\text{b}(\xi)\;\;\xi^{\text{B}_{0}}(1-\xi)^{\text{B}_{1}}e^{\frac{\text{B}_{2}}{1-\xi}}\;d\xi\;.$
(33)

One of the simplest ways to evaluate Eq. (33) is direct integration (i.e.
explicit numerical integration via discretization of the domain). However,
divergence of the weight function at spatial infinity or the BH horizon means
that naive integration along real values of $\xi$ is bound to fail. So other
means must be considered.

In Sec. V.1, we make use of coordinate dependence to develop a method to
directly integrate Eq. (33). In Sec. V.2, we note that Eq. (33) coincides with
the integral representation of the Tricomi confluent hypergeometric function,
also known as the Kummer function of the second kind Arfken  _et al._ (2013);
Chen and Dai (2010). That fact, along with the well known analytic
continuation of e.g. the Gamma function to complex inputs, is used to
analytically evaluate Eq. (33). This section concludes with a brief comparison
of the pros and cons of direct integration and analytic continuation.

###  V.1 Direct integration

An extremely useful lesson from calculus is that, given an integral, we may
change the coordinate representation of its integrand without changing the
value of the integral. So while the kernel of Eq. (33) is coordinate
dependent, its value, ${\langle{\text{a}}\,|\,{\text{b}}\rangle}$, is
coordinate independent. This concept leads us to seek a new coordinate, $z$,
such that the kernel of Eq. (33) does not diverge. Henceforth, we will use
$\xi(z)$ to refer to a complex valued integration path, and $\xi$ to refer to
the quantity defined by Eq. (10).

The following definition of $\xi(z)$ has the desired effect of augmenting
singular behavior at the horizon and spatial infinity,

$\displaystyle\xi(z)\;=\;\frac{{x_{2}}z^{{x_{0}}}}{({x_{2}}+1)z^{{x_{0}}}-1}\;.$
(34)

In Eq. (34), $x_{0}$ and $x_{2}$ are coordinate parameters that modify
$\text{B}_{0}$ and $\text{B}_{2}$, respectively. We will consider
$\text{Re}(x0)>0$, such that $\xi(0)=0$. The dependence of Eq. (34) on $x_{2}$
is motivated by the complex radial coordinate used in Ref. Mano  _et al._
(1996). It may also be developed by placing the Laguerre polynomial equation
into a $\xi$-like compactified coordinate, and then studying that system’s
weight function. The dependence of Eq. (34) on $x_{0}$ is motivated by the
observation that $z^{x_{0}}$ is the simplest apparent coordinate modification
that impacts singular behavior arising from $\text{B}_{0}$.

The effect of this coordinate choice is found by applying Eq. (34) to the
scalar product’s integrand,

$\displaystyle{\langle{\text{a}}\,|\,{\text{b}}\rangle}=\int_{0}^{1}\text{a}(\xi(z))\,\text{b}(\xi(z))\;\left[\text{W}(\xi(z))\;\partial_{z}\xi(z)\right]\;dz\;.$
(35)

In Eq. (35), we have used the fact that $d\xi=\partial_{z}\xi(z)\,dz$ to write
the integral explicitly in $z$. According to Eq. (34), $\xi(z=0)=0$, and
$\xi(z=1)=1$; thus, in Eq. (35), the limits of integration do not change.

Within Eq. (34), the quantity in square brackets is of particular practical
use. It is an effective weight function,

$\displaystyle\text{W}_{\text{eff.}}(z)\;$
$\displaystyle=\;\text{W}(\xi(z))\;\partial_{z}\xi(z)$ (36a)
$\displaystyle=\;e^{\frac{\text{B}_{2}\text{x}_{2}}{z^{\text{x}_{0}}-1}}\;z^{\text{x}_{0}(\text{B}_{0}+1)-1}\;\left(1-z^{\text{x}_{0}}\right)^{\text{B}_{1}}\;$
(36b)
$\displaystyle\quad\quad\times\;\left(1-(\text{x}_{2}+1)z^{\text{x}_{0}}\right)^{-\text{B}_{0}-\text{B}_{1}-2}\;.$

While Eq. (34) provides a coordinate that modifies the boundary behavior of
$\text{W}(\xi)$, it is Eq. (36) that must be analyzed to determine the special
values of $x_{0}$ and $x_{2}$ that prevent the effective weight function from
diverging at the domain boundaries. In essence, we have two independent
conditions: (i) $\text{W}_{\text{eff.}}(0)=0$, and (ii)
$\text{W}_{\text{eff.}}(1)=0$. We also have two unknowns: $x_{0}$ and $x_{1}$.
Thus, $x_{0}$ and $x_{1}$ may be determined (up to positive constants) by
careful analysis of $\text{W}_{\text{eff.}}$.

This analysis is somewhat tedious, but tractable. First, one should focus on
the factor of $z^{\text{x}_{0}(\text{B}_{0}+1)-1}$ within
$\text{W}_{\text{eff.}}$. If $x_{0}$ is chosen appropriately, then the real
part of $\text{x}_{0}(\text{B}_{0}+1)-1$ should equal some positive constant,
$\text{K}_{0}$, thereby allowing $\text{W}_{\text{eff.}}$ to go to zero as $z$
goes to zero. Since, $\text{B}_{0}$ is complex valued, $x_{0}$ must be also.
These ideas result in the following constraints,

$\displaystyle\text{K}_{0}$
$\displaystyle=-1+\left[1+\text{Re}(\text{B}_{0})\right]\text{Re}(x_{0})-\text{Im}(\text{B}_{0})\text{Im}(x_{0}),$
(37a) $\displaystyle\text{K}_{0}$ $\displaystyle>0\;.$ (37b)

If we assume that $\text{Re}(x_{0})$ only depends on
$\text{Re}(\text{B}_{0})$, then then the result for $x_{0}$ follows,

$\displaystyle\text{x}_{0}\;$
$\displaystyle=\;\text{K}_{0}\,\frac{\text{sgn}(\mathrm{Re}(\text{B}_{0})+1)}{\mathrm{Re}(\text{B}_{0})+1}$
(38)
$\displaystyle\quad\;+\;i\,\frac{-1\,+\,\text{K}_{0}\,[\text{sgn}(\;\mathrm{Re}(\text{B}_{0})+1\;)-1]}{\mathrm{Im}(\text{B}_{0})}\;.$

Note that, in Eq. (38), $\text{Re}(\text{x}_{0})>0$, as is needed for
$\xi|_{z=0}$ to be zero.

Given, $x_{0}$, one may then determine $x_{2}$. Returning to
$\text{W}_{\text{eff.}}$, if its exponential factor has a negative argument as
$z\rightarrow 1$, then $\text{W}_{\text{eff.}}$ is assured to be regular in
that limit. To that end, one should consider
$e^{\frac{\text{B}_{2}\text{x}_{2}}{z^{\text{x}_{0}}-1}}$, and seek $x_{2}$
such that

$\displaystyle\lim_{z\rightarrow
1}\frac{\text{B}_{2}\text{x}_{2}}{z^{\text{x}_{0}}-1}\;$
$\displaystyle\sim\;-\frac{\text{K}_{2}}{1-z}\;,$ (39a)
$\displaystyle\text{K}_{2}\;$ $\displaystyle>0\;.$ (39b)

To evaluate the left-hand-side of Eq. (39a), one can use the binomial
theorem111111For this it is useful to work with $z^{\prime}=1-z$.; this
results in a quantity with the same form as the right-hand-side of Eq. (39a).
Upgrading the asymptotic equivalence in Eq. (39a) to an equality yields the
following solution for $\text{x}_{2}$,

$\displaystyle\text{x}_{2}\;$
$\displaystyle=\;\text{K}_{2}\;\text{x}_{0}\,/\,\text{B}_{2}\;.$ (40)

Note that, since $\text{x}_{0}$ is complex valued, the factor of
$\left(1-(\text{x}_{2}+1)z^{\text{x}_{0}}\right)^{-\text{B}_{0}-\text{B}_{1}-2}$
within $\text{W}_{\text{eff.}}$ will never encounter a root. Thus it is of
little practical concern.

In Eqs. (38) and (40) we define the values of the $\text{x}_{0}$ and
$\text{x}_{2}$ are allowed to take such that $\text{W}_{\text{eff.}}$ is
regular (actually zero) at the event horizon and spatial infinity. In doing
so, we have transferred the coordinate dependence to two new, strictly
positive, parameters, $\text{K}_{0}$ and $\text{K}_{2}$. What remains is the
possibility that $|\xi(z)|\gg 1$.

If $\text{K}_{0}$ and $\text{K}_{2}$ are not chosen such that
$|\xi(z)|\lessapprox 1$, then numerical evaluation of Eq. (35) can encounter
extremely large numbers. This numerically problematic possibility can be
avoided if one determines $\text{K}_{0}$ and $\text{K}_{2}$ by requiring that
$\xi(z)$ is contained as much as possible within a unit circle centered about
$\xi=1/2$. This may be accomplished via the following minimization,

$\displaystyle\\{\text{K}_{0},\text{K}_{2}\\}=\text{argmin}_{\text{K}_{0},\text{K}_{2}}S(\text{K}_{0},\text{K}_{2},z)$
(41a) where $\displaystyle S(\text{K}_{0},\text{K}_{2},z)\;$
$\displaystyle=\;|1/2-\text{mean}_{z}(|\xi(z)|)|$ (41b)
$\displaystyle\quad\;+\;\text{max}_{z}|\text{Im}\,\xi(z)|\;+\;|\text{min}_{z}\;\text{Re}\,\xi(z)|\;.$

In Eq. (41), argmin outputs the values for which $S$ is minimized in
$\text{K}_{0}$ and $\text{K}_{2}$, $\text{mean}_{z}$ computes the unweighted
mean over values of $z$, $\text{min}_{z}$ returns the smallest value of a
quantity in $z$, and $\text{max}_{z}$ returns the largest value of a quantity
in $z$.

In practice, one has a number of remaining options for actually evaluating Eq.
(35). For the evaluation of Eq. (41), it is found to be useful to densely
sample, on log-scale, $z$ near $0$ and $1$. The coordinate freedom associated
with $x_{0}$ introduces oscillations into $\xi(z)$ which must be well resolved
for the accurate determination of $\text{K}_{0}$ and $\text{K}_{2}$.
Similarly, for the evaluation of Eq. (35), one will want to densely sample
near the boundaries; it is found that $z=\frac{1}{2}(1+\sin(t\pi/2))$, with
$t\in(-1,1)$ is sufficient for this purpose. For the determination of
$\text{K}_{0}$ and $\text{K}_{2}$, and for the evaluation of Eq. (35), it is
found that at least $10^{4}$ points in $z$ are needed for robust results. With
these details in hand, it is now possible to compute the integration path,
$\xi(z)$, and related scalar-products, Eq. (35), using e.g. $5^{\text{th}}$
order spline integration Virtanen  _et al._ (2020).

Figure 1 shows three example integration paths, $\xi(z)$, set against contours
of the weight function, $\text{W}(\xi)$, in the complex plane. In focus is the
effect of BH spin $a$. The result for $a=0$ is qualitatively similar to that
for $a=0.7$. Branch cuts are visible as overlapping horizontal contours along
the real line for $\xi<0$ and $\xi>1$. The shape of integration paths are seen
to vary significantly with increasing BH spin, changing most drastically near
the extremal limit, $a=M$.

For $a=0.99$ (the central panel), the integration path appears to cross the
branch cut many times. The top panel of Fig. 2 shows this in detail. There,
the apparent crossing of branch cuts is an artifact of presentation: $\xi(z)$
actually runs along all branches the Riemann surface defined by $\text{W}(z)$,
but in Fig. 1, only a single branch of $\text{W}(\xi)$ is shown Kuijlaars  _et
al._ (2003). This can be understood by noting that the transformed scaler-
product, Eq. (35), uses an integration path that follows the real line, where
there are no branch cuts. This is shown explicitly in the bottom panel of Fig.
2.

###  V.2 Analytic continuation

---
Figure 2: Examples integration paths (white curves) against weight function
contours (colored contours) for the $(s,\ell,m,n)=(-2,2,2,3)$ QNM at BH spin
$a/M=0.99$. (top) A zoom-in of the central panel of Fig. 1 with axes centered
around the BH horizon at $\xi=0$. (bottom) Visualization of the coordinate
transformed integration path and effective weight function for the same case.
Contour formatting and coordinate parameters are identical to those in Fig. 1.

The scalar product developed in Sec. IV is,

$\displaystyle{\langle{\text{a}}\,|\,{\text{b}}\rangle}\;=\;\int_{0}^{1}\,\text{a}(\xi)\,\text{b}(\xi)\;\;\xi^{\text{B}_{0}}(1-\xi)^{\text{B}_{1}}e^{\frac{\text{B}_{2}}{1-\xi}}\;d\xi\;.$
(42)

The method described in the previous section uses coordinate freedom to extend
the domain upon which Eq. (42) can be evaluated. This is the basic purpose of
the set of techniques known as analytic continuation Ahlfors ; Ronveaux and
Arscott (1995). It happens that there is a more direct, and ultimately more
robust, route to evaluating Eq. (42) with analytic continuation.

This alternative route begins with the observation that if
$\text{a}(\xi)=\text{b}(\xi)=1$, then the left-hand-side of Eq. (42) coincides
exactly with the integral representation of the Tricomi121212The nomenclature
for these functions appears to be particularly inconsistent across sources.
Some readers may find Refs. DLMF (a) and Chen and Dai (2010) of particular
use. confluent hypergeometric function Arfken  _et al._ (2013); Chen and Dai
(2010); Abramowitz and Stegun (1964); DLMF (a); Leaver (1986),

$\displaystyle\int_{0}^{1}\,$
$\displaystyle\xi^{\text{B}_{0}}\,(1-\xi)^{\text{B}_{1}}\,e^{\frac{\text{B}_{2}}{1-\xi}}\;d\xi$
(43)
$\displaystyle=\;e^{\text{B}_{2}}\,\Gamma(\text{B}_{0}+1)\,U(1+\text{B}_{0},-\text{B}_{1},-\text{B}_{2})\;.$

In Eq. (43), $\Gamma(x)$ is the usual Euler-Bernoulli Gamma function, and
$U(\text{x},\text{y},\text{z})$ is the Tricomi confluent hypergeometric
function,

$\displaystyle U(\text{x},\text{y},\text{z})\;=\;$
$\displaystyle\frac{\Gamma(1-\text{y})}{\Gamma(\text{x}-\text{y}+1)}M(\text{x},\text{y},\text{z})$
(44)
$\displaystyle\;+\;\frac{\Gamma(\text{y}-1)}{\Gamma(\text{x})}\text{z}^{1-\text{y}}M(\text{x}-\text{y}+1,2-\text{y},\text{z})\;.$

Both may be evaluated for complex inputs by considering the appropriate
integration contour Whittaker and Watson (1996). Finally, in Eq. (44),
$M(\text{x},\text{y},\text{z})$ is the confluent hypergeometric function of
the first kind (See e.g. Eq. $13.2.42$ of Ref. DLMF (a)).

The Gamma functions is constructed to extend the factorial function to the
complex numbers. The confluent hypergeometric function is a product of Gamma
functions and a power series in z, and so it is straightforward to evaluate.

This parade of special functions has the following use: Every $\text{a}(\xi)$
and $\text{b}(\xi)$ of interest to Eq. (42) will be smooth, and so posses a
series expansion,

$\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}\;.$ (45)

This allows every scalar product to be expressed as a sum over scalar products
of monomials,

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

In Eq. (46), $\langle\xi^{p}\rangle$ is the $p^{\text{th}}$ monomial moment
defined as,

$\displaystyle\left<\xi^{p}\right>\;$
$\displaystyle=\;\int_{0}^{1}\,\xi^{p}\,\text{W}(\xi)\;d\xi$ (47a)
$\displaystyle=\;\int_{0}^{1}\,\xi^{\text{B}_{0}+p}\,(1-\xi)^{\text{B}_{1}}\,e^{\frac{\text{B}_{2}}{1-\xi}}\;d\xi$
(47b)
$\displaystyle=\;e^{\text{B}_{2}}\,\Gamma(\text{B}_{0}+p+1)\,U(1+\text{B}_{0}+p,-\text{B}_{1},-\text{B}_{2})\;.$
(47c)

In Eq. (47a), we have effectively defined $\left<\xi^{p}\right>$ to be
${\langle{1}\,|\,{\xi^{p}}\rangle}$. In Eq. (47b), we have applied the
definition of the scalar product, and in Eq. (47c), we have equated the
integral over $\xi$ with the product of special functions found in Eq. (43),
but with $\text{B}_{0}\rightarrow\text{B}_{0}+p$. The result is that the
scalar product of any two smooth functions may be equated with a sum over
confluent hypergeometric terms,

$\displaystyle{\langle{\text{a}}\,|\,{\text{b}}\rangle}\;$
$\displaystyle=\;e^{\text{B}_{2}}\,\sum_{j,k}\;\text{a}_{j}\,\text{b}_{k}\,\Gamma(\text{B}_{0}+j+k+1)\,$
(48) $\displaystyle\hskip
56.9055pt\times\,U(1+\text{B}_{0}+j+k,-\text{B}_{1},-\text{B}_{2})\;.$

In practice, it is useful to precompute monomial moments using Eq. (47), and
then compute scalar products using Eq. (46). The monomial moments are defined
up to an overall constant. Henceforth, we chose to normalize all
$\left<\xi^{p}\right>$ such that $\left<\xi^{1}\right>=1$,

$\displaystyle\left<\xi^{p}\right>\;\leftarrow\;\left<\xi^{p}\right>/\left<\xi^{1}\right>\;.$
(49)

##  VI Limitations of direct integration and analytic continuation

Direct integration, analytic continuation, and their application to monomial
moments and confluent Heun polynomials are all subject to various caveats and
limitations. Here, these facets are reviewed. At times, technical details are
omitted for brevity, but particular attention is given to practical
workarounds. While analytic continuation and direct integration are found to
be complementary, analytic continuation is determined to be preferred.

###  VI.1 Direct Integration: disallowed BH spins, and the zero frequency
limit

The effective weight function introduced in Sec. V.1 was

$\displaystyle\text{W}_{\text{eff.}}(z)\;$
$\displaystyle=\;e^{\frac{\text{B}_{2}\text{x}_{2}}{z^{\text{x}_{0}}-1}}\;z^{\text{x}_{0}(\text{B}_{0}+1)-1}\;\left(1-z^{\text{x}_{0}}\right)^{\text{B}_{1}}\;$
(50)
$\displaystyle\quad\quad\times\;\left(1-(\text{x}_{2}+1)z^{\text{x}_{0}}\right)^{-\text{B}_{0}-\text{B}_{1}-2}\;.$

In Sec. V.1, a direct integration approach was introduced to evaluate scalar
product using $\text{W}_{\text{eff.}}(z)$. This approach has two limiting
features: (i) for a given value of $\tilde{\omega}$, there are values of the
BH spin, $a$, such that an integration path cannot be constructed, and (ii)
when $\tilde{\omega}=0$, the weight function may diverge at $\xi=1$.

For any value of $\tilde{\omega}=\omega-i/\tau$, there may exist values of the
BH spin,

$\displaystyle a=a_{*}\;,$ (51)

such that an integration path cannot be constructed as described in Sec. V.1.
This limitation results from the requirement that $\text{Re}(\text{x}_{0})>0$.

In Sec. V.1, the coordinate parameter $\text{x}_{0}$ was introduced so that
$\text{W}_{\text{eff.}}(z)$ is zero at $z=0$. The value found for
$\text{x}_{0}$ to have that effect was,

$\displaystyle\text{x}_{0}\;$
$\displaystyle=\;\text{K}_{0}\,\frac{\text{sgn}(\mathrm{Re}(\text{B}_{0})+1)}{\mathrm{Re}(\text{B}_{0})+1}$
(52)
$\displaystyle\quad\;+\;i\,\frac{-1\,+\,\text{K}_{0}\,[\text{sgn}(\;\mathrm{Re}(\text{B}_{0})+1\;)-1]}{\mathrm{Im}(\text{B}_{0})}$

where $\text{sgn}(x)$ is defined by $x=\text{sgn}(x)\,|x|$, and

$\displaystyle\text{B}_{0}\;=$
$\displaystyle\;-\left(\frac{2M(\delta+M)}{\tau\delta}+s\right)$ (53)
$\displaystyle\hskip
14.22636pt\;+\;i\left(\frac{am-2M^{2}\omega}{\delta}-2M\omega\right)\;.$

and $\text{K}_{0}>0$. As seen in Eq. (52), the real part of $\text{x}_{0}$ is
generally regular, since $\text{sgn}(\mathrm{Re}(\text{B}_{0})+1)$ is defined
to be zero when $\text{Re}(\text{B}_{0})+1$ is zero. However, the imaginary
part of $\text{x}_{0}$ may diverge if $\text{Im}(\text{B}_{0})=0$. If we had
constructed $\text{x}_{0}$ differently, e.g. by allowing
$\text{Re}(\text{x}_{0})$ to depend on both $\text{Re}(\text{B}_{0})$ and
$\text{Im}(\text{B}_{0})$, then the requirement that
$\text{Re}(\text{x}_{0})>0$ would still have resulted in a
$1/\text{Im}(\text{B}_{0})$ divergence.

Taking $\omega$ as known, there may exist values of the BH spin, namely
$a_{*}$, such that

$\displaystyle\text{Im}{(\text{B}_{0})}|_{a=a_{*}}\;$ $\displaystyle=\;0$
(54a) Applying the definition of $\text{B}_{0}$ (see Eq. 53), and recalling
that $\delta=\sqrt{M^{2}-a^{2}}$ allows Eq. (54a) to be solved for $a_{*}$,
$\displaystyle a_{*}\;$
$\displaystyle=\;\frac{4mM^{2}w}{m^{2}+4M^{2}w^{2}}\;.$ (54b)

For physical applications, it happens that for any QNM $(\ell,m,n)$, having
frequency

$\displaystyle\tilde{\omega}(a)\;=\;\omega(a)-i/\tau(a)\;,$ (55)

there is a discrete set BH spins where $\omega$ is such that $a=a_{*}$. For a
given QNM $\ell$, $m$ and $n$, this set may be defined by first evaluating
$\omega(a)$ using e.g. positive.physics.qnmobjLondon  _et al._ (2020), and
then searching for values of $a$ such that

$\displaystyle|a-a_{*}(\tilde{\omega}{(a)})|\;=\;0\;.$ (56)

For $(\ell,m,n)=(2,2,0)$, there is only one such BH spin,

$\displaystyle a_{*}\;\approx\;0.996632\;.$ (57)

It is instructive to contemplate the behavior of $\text{x}_{0}$ when $a$ is
above this critical value.

Towards extremal BH spin, $a\rightarrow 1$, $\text{x}_{0}$ typically remains
regular, meaning that an integration path, $\xi(z)$, may still be constructed.
This is a result of how $\text{x}_{0}$ depends on $\text{B}_{0}$. Since
$\lim_{a\rightarrow 1}|\text{Re}(\text{B}_{0})|\rightarrow\infty$,
$\text{Re}(\text{x}_{0})$ tends to zero in that same limit. Similarly, it may
be deduced that $\text{Im}(\text{x}_{0})$ tends to zero in that limit. Recall
from Eq. (34) that the integration path, $\xi(z)$, is not defined when
$\text{x}_{0}=0$. Despite this, the reasoning above implies that $\xi(z)$ may
be constructed in the near extremal limit, within the limits of numerical
precision; i.e. when $\text{x}_{0}$ is small, but strictly non-zero.

Lastly, a second coordinate parameter, $\text{x}_{2}$, was introduced so that
exponential divergence of the weight function was avoided (see Sec. V.1). The
full expression for $\text{x}_{2}$ was

$\displaystyle\text{x}_{2}\;$
$\displaystyle=\;\frac{\text{K}_{2}\,\text{x}_{0}}{\text{B}_{2}}\;$ (58a)
$\displaystyle=\;\frac{\text{K}_{2}\,\text{x}_{0}}{4i\delta\tilde{\omega}}\;.$
(58b)

In Eq. (58b), we have used the fact that
$\text{B}_{2}=4i\delta\tilde{\omega}$. Since $\text{x}_{2}$ is proportional to
$\tilde{\omega}^{-1}$, it generally diverges in the zero frequency limit.
However, since only
$\text{x}_{2}\,\text{B}_{2}\,=\,\text{K}_{2}\,\text{x}_{0}$ appears in the
weight function (see e.g. Eq. 50), $\text{W}_{\text{eff.}}(z)$ may still be
evaluated.

Further, in the zero frequency limit there is no exponential divergence of the
weight function at $\xi=1$, and so there is no formal need for $\text{x}_{2}$
at all in that limit. Seeing as

$\displaystyle\text{W}_{\text{eff.}}(z)|_{\tilde{\omega}=0}\;\propto\;\left(1-z^{\text{x}_{0}}\right)^{\text{B}_{1}}\;,$
(59)

$\text{W}_{\text{eff.}}$ will have a pole at $z=1$ if
$\text{Re}(\text{B}_{1})<0$. This is of no concern when $\tilde{\omega}\neq
0$, since $\text{x}_{2}$ is designed to enforce exponential damping at $z=1$;
otherwise, when $\tilde{\omega}=0$, it is possible to construct a different
coordinate, $z^{\prime}$, such that this pole is avoided. This work-around may
only be viable when $\tilde{\omega}=0$; otherwise, it is found that when
$\tilde{\omega}\neq 0$, the required coordinate change necessarily introduces
oscillations in $\xi(z^{\prime})$ that grow exponentially as $z$ increases.
Thus, when $\tilde{\omega}=0$, one can only use coordinate freedom to either
avoid the essential singularity at $\xi(z)=1$ (i.e. that resulting from the
exponential), or the pole at $\xi(z)=1$, but not both.

In this, direct integration’s underlying limitation is one of mathematical
perspective: it happens to be more robust to use analytic continuation via
integration over closed contours, rather than singularity avoiding integration
over finite paths. For example, Ref. Ronveaux and Arscott (1995) details how
single and double Pochhammer contours may be used in similar mathematical
settings. This approach makes use of standard theorems in complex analysis,
and is ultimately equivalent to analytic continuation as outlined in Sec. V.2.

###  VI.2 Analytic continuation and limitations of the Gamma function

It was found in Sec. V.2 that scalar products may be evaluated analytically
using monomial moments, $\left<\xi^{p}\right>$,

$\displaystyle\left<\xi^{p}\right>\;$
$\displaystyle=\;\Gamma(1+\text{B}_{0}+p)\,U(1+\text{B}_{0}+p,-\text{B}_{1},-\text{B}_{2})\;.$
(60a) The qualitative behavior of the monomial moments can be described by
$\displaystyle\left<\xi^{p}\right>\;$
$\displaystyle\approx\;\frac{\Gamma(1+\text{B}_{0}+p)\,\Gamma(1+\text{B}_{1})}{\Gamma(2+\text{B}_{0}+p+\text{B}_{1})}\;,$
(60b)

where the leading term in the power series definition for
$U(1+\text{B}_{0}+p,-\text{B}_{1},-\text{B}_{2})$ has been used to approximate
it with ${\Gamma(1+\text{B}_{1})}/{\Gamma(2+\text{B}_{0}+p+\text{B}_{1})}$
DLMF (a). Since $\text{B}_{2}=4i\delta\tilde{\omega}{}$, Eq. (60b) is exact in
the zero frequency limit, $\tilde{\omega}=0$, or in the extremal spin limit,
$\delta=\sqrt{M^{2}-a^{2}}=0$.

For either Eq. (60a) or Eq. (60b), $\left<\xi^{p}\right>$ will be regular when
instances of the Gamma function are regular. Since $\Gamma(z)$ is well known
to have simple poles when $z$ is a negative integer, it follows that
$\left<\xi^{p}\right>$ can diverge when,

$\displaystyle 1+\text{B}_{0}+p$
$\displaystyle\in\\{\;z\in\mathbb{Z}\;|\;z<0\;\\}\;,$ (61a) and $\displaystyle
1+\text{B}_{1}$ $\displaystyle\in\\{\;z\in\mathbb{Z}\;|\;z<0\;\\}\;.$ (61b)

In Eq. (61), $\\{z\in\mathbb{Z}\,|\,z<0\\}$ denotes the set of all negative
integers. The practical implications of Eq. (61) may be investigated by
rewriting $\text{B}_{0}$ and $\text{B}_{1}$ in terms of physical parameters,

$\displaystyle
p-s+1-\frac{2M(\delta+M)}{\tau\,\delta}\quad\quad\quad\quad\;\;$ (62a)
$\displaystyle+\;\frac{i(am-2M\omega(\delta+M))}{\delta}$
$\displaystyle\;=\;-z_{0}\;,$ $\displaystyle 2s+1-\frac{4M}{\tau}-4iM\omega$
$\displaystyle\;=\;-z_{1}\;.$ (62b)

In Eq. (62), $z_{0}$ and $z_{1}$ are natural numbers, and we have used Eq.
(8), Eqs. (16) and (24) to expand in terms of basic physical quantities.
Equation (62a) is equivalent to Eq. (61a), and Eq. (62b) is equivalent to Eq.
(61b).

Inspection of Eqs. (62a) and (62b) reveals that they can only hold when
$M\tilde{\omega}{}$ is purely imaginary (i.e. $\omega=0$). This limits the
physical relevance of these equations to the total transmission QNMs Cook and
Zalutskiy (2014) (LL: Check me). Furthermore, Eq. (62a) can only be true if
$a=0$, thus it is only relevant to the algebraically special modes of
Schwarzschild Cook and Zalutskiy (2014); Berti  _et al._ (2003).

Both Eqs. (62a) and (62b) have non-trivial implications for evaluation of the
scalar product with analytic continuation. If Eq. (62a) holds, then there
exist some monomials, $\xi^{p}$, such that the respective monomial moments are
infinite when defined by Eq. (60a). If Eq. (62b) holds, then no monomial
moments are finite according to Eq. (60a).

These limitations are, again, a result of mathematical perspective: evaluation
of the scalar product with Gamma and confluent hypergeometric functions may
fail, depending on $\text{Re}(\text{B}_{0})$ or $\text{Re}(\text{B}_{1})$;
however, we saw in the previous section that the viability of direct
integration only depends on $\text{Im}(\text{B}_{0})$. Thus the two
perspectives have complementary limitations.

###  VI.3 Direct Integration vs. Analytic continuation

In practice, analytic continuation has the following relative attributes that
may make it preferable in most circumstances. Since analytic continuation does
not require the determination of optimal coordinate parameters (see Sec. V.1),
its implementation requires fewer steps. In addition, analytic continuation
requires monomial moments to compute scalar products, ans do does not require
the storage of arrays (discretizations of $\xi(z)$) of length $10^{5}$ or
greater (see Sec. V.1). Instead, monomial moments with degrees up to twice
that of the largest polynomial order of interest may be precomputed, ideally
using arbitrary precision arithmetic mpmath development team (2023), and then
used for all scalar products (Sec. V.2). Thus, while direct integration may
serve as a valuable alternative, particularly in cases where the Gamma
function cannot be evaluated, analytic continuation is expected to be
structurally simpler and computationally more efficient.

All forthcoming numerical results will make use of analytic continuation with
$64$ digits of precision via the public mpmath package mpmath development team
(2023).

|  |
---|---|---
Figure 3:  Example distributions of absolute monomial moments,
$|\langle\xi^{k}\rangle|$ v.s. monomial order, $k$, for the $(\ell,m)=(2,2)$
QNM at BH dimensionless spin of $a/M=0.7$. Left, central and right panels show
results for $n$ of $0$, $2$ and $6$, respectively. The corresponding QNM
frequencies are $M\tilde{\omega}=0.5326-0.0808i$,
$M\tilde{\omega}=0.4999-0.4123i$, and $M\tilde{\omega}=0.4239-1.0954i$
respectively. At focus is the effect of negating the spin weight, $s$, for
increasing overtone label. For all cases shown, $|\langle\xi\rangle|=1$ has
been enforced. Closed circles show points for $s=+2$, and open circles show
points for $s=-2$.

##  VII Pedagogical Applications

Previous sections provide a foundation for applying the scalar product, Eq.
(63), to physical scenarios in which the series expansion for radial functions
is either known (analytic continuation), or unknown (direct integration).
Direct integration is applicable in both scenarios.

$\displaystyle{\langle{\text{a}}\,|\,{\text{b}}\rangle}\;=\;\int_{0}^{1}\,\text{a}(\xi)\,\text{b}(\xi)\;\;\xi^{\text{B}_{0}}(1-\xi)^{\text{B}_{1}}e^{\frac{\text{B}_{2}}{1-\xi}}\;d\xi\;.$
(63)

Before leaping into applications to Teukolsky’s radial equation, there are
good reasons to first look closely at two pedagogical cases.

In this section, we will first apply the scalar product to the monomial
moments for different QNM values of the frequency $\tilde{\omega}$, paying
particular attention to the effect of spin weight, $s$. For a broad range of
scenarios, monomial moments with $s=2$ decrease monotonically with monomial
order, while this is less so the case for $s=-2$. We will then apply the
scalar product to the confluent Heun polynomials.

Previously, Cook and Zalutskiy considered confluent Heun polynomials for BHs
in Ref. Cook and Zalutskiy (2014); Fiziev (2009a). There, the polynomials were
defined such that a special choice $\tilde{\omega}{}$ resulted in
$\mathcal{L}_{\xi}$ (Eq. 15) having polynomial eigenfunctions. These are the
well known algebraically special and total transmission QNMs Cook and
Zalutskiy (2014); Andersson (1994); Dias and Reall (2013); Berti  _et al._
(2003); Whiting (1989); Fiziev (2009a).

Here, we take a different perspective: We will a priori define
$\tilde{\omega}{}$ to be any QNM frequency 131313In the language of Ref.
Ronveaux and Arscott (1995), we choose to focus on Class I of confluent Heun
polynomials.. This requires that, in effect, we redefine the potential of
$\mathcal{L}_{\xi}$ such that its eigenfunctions are polynomials. The
consequences of this choice will have a significant impact on forthcoming
results and discussion.

We will use the scalar product to derive the orthogonality between confluent
Heun polynomials of like order. Lastly, we will use the scalar product to
derive a relationship between the confluent Heun polynomials’ eigenvalues, and
a specific average of $\xi$. This relationship provides some insight into the
structure of the eigenvalues.

###  VII.1 Monomial moments

The scalar-product of any two smooth functions may be evaluated using a linear
combination of monomial moments,

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

For many common scalar products, monomial moments $\left<\xi^{p}\right>$
decrease monotonically with $p$. This is easily demonstrated by considering a
hypothetical constant weight function:

$\displaystyle\text{W}(\xi)_{\text{hypothetical}}\;$ $\displaystyle=\;1\;,$
(65a)
$\displaystyle{\langle{\text{a}}\,|\,{\text{b}}\rangle}_{\text{hypothetical}}\;$
$\displaystyle=\;\int_{0}^{1}\,\text{a}(\xi)\,\text{b}(\xi)d\xi\;,$ (65b)
$\displaystyle\left<\xi^{p}\right>_{\text{hypothetical}}\;$
$\displaystyle=\;(p+1)^{-1}\sim p^{-1}\;.$ (65c)

In Eq. (65a), a weight function of $1$ results in a scalar product equivalent
to that used for the shifted Legendre polynomials, Eq. (65b). In turn, this
results in monomial moments that are asymptotically equivalent to one over the
degree of the monomial, Eq. (65c). In practice, this means that if the
sequences comprised of $\text{a}_{j}$ and $\text{b}_{k}$ are convergent, then
so too is $\text{a}_{j}\,\text{b}_{k}\,\left<\xi^{j+k}\right>$. We will now
examine this argument in the context of the QNM’s scalar product. Our goal is
to gain a basic and qualitative understanding of how monomial moments depend
on two physical parameters: $s$ and $M\tilde{\omega}$.

We previously observed that monomial moments may be computed via analytic
continuation of the confluent hypergeometric function, $U$, and the gamma
function, $\Gamma$,

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

The confluent hypergeometric function, $U$, is defined as a series expansion
in its third input DLMF (a). Thus some understanding of $\left<\xi^{p}\right>$
may be found by examining $U$ expanded in powers of $\text{B}_{2}$.
Examination of the zeroth order term happens to be sufficient for qualitative
understanding of how $\left<\xi^{p}\right>$ behaves when $p$ is large. For
precise and predictive understanding, one would have to extend the forthcoming
analysis to higher orders. This will not be done here so that the following
discussion remains fairly simple.

To zero’th order in its third input, $U$ takes the approximate form,

$\displaystyle
U(\text{x},\text{y},\text{z})\;\approx\;\frac{\Gamma(1-\text{y})}{\Gamma(\text{x}-\text{y}+1)}+\mathcal{O}(z)\;.$
(67)

When applied to Eq. (66), this yields

$\displaystyle\left<\xi^{p}\right>\;$
$\displaystyle\approx\;\frac{\Gamma(1+\text{B}_{0}+p)\,\Gamma(1+\text{B}_{1})}{\Gamma(2+\text{B}_{0}+p+\text{B}_{1})}$
(68a) $\displaystyle\sim\;p^{-1-\text{B}_{1}}\;\Gamma(1+\text{B}_{1})\;.$
(68b)

To go from Eq. (68a) to Eq. (68b), we have used the fact that, as
$p\rightarrow\infty$, $\Gamma(p+\alpha)\,\sim\,p^{\alpha}\,\Gamma(p)$
Abramowitz and Stegun (1964). The physics of the situation is revealed by
using Eqs. (24) and (16) to rewrite $\text{B}_{1}$ in terms of $M$,
$\tilde{\omega}$ and $s$,

$\displaystyle\left<\xi^{p}\right>\;$
$\displaystyle\sim\;p^{-1-2s+4iM\tilde{\omega}}\;$ (69a)
$\displaystyle\sim\;p^{-1-2s+4M/\tau}e^{4iM\omega\,\ln(p)}\;.$ (69b)

In Eq. (69a), we see that the previous result of $p^{-1}$ is effectively
generalized to depend on the spin weight, $s$, and the dimensionless
frequency, $M\tilde{\omega}$. Since the net exponent of $p$ is now complex
valued, it is the real part of the exponent that determines the large-$p$
behavior of $|\left<\xi^{p}\right>|$. In Eq. (69b), this is shown explicitly
by using $\tilde{\omega}=w-i/\tau$. Recall that $\tau$ is strictly positive as
is needed for time domain stability.

The presence of $e^{4iM\omega\,\ln(p)}$ in Eq. (69b) means that, in general,
the real and imaginary parts of $\left<\xi^{p}\right>$ will trace out a
logarithmic spiral.

The presence of $p^{-1-2s+4M/\tau}$ within Eq. (69b) allows us to discern
three additional features:

  1. (i)

There exist some combination of $s$ and $M/\tau$ such that monomial moments
either decrease, asymptote to a constant, or increase as $p$ tends towards
infinity.

  2. (ii)

Since Eq. (69b) depends on $p^{-2s}$, monomial moments for $s=+2$ will be much
more inclined to decrease as $p\rightarrow\infty$:

$\displaystyle\left<\xi^{p}\right>|_{s=+2}\;\sim\;\
p^{-8}\,\left<\xi^{p}\right>|_{s=-2}$ (70)

In Eq. (70), we have used the fact that if $\tilde{\omega}{}$ is a QNM
frequency, then $\tilde{\omega}|_{s=+2}=\tilde{\omega}|_{s=-2}$.

  3. (iii)

Lastly, the presence of $4M/|\tau|$ in Eq. (69b) means that, as $\tau$
decreases, the rate of change of $\left<\xi^{p}\right>$ will typically
increase. This is perhaps most relevant for the QNM overtones, for which
$|\tau|$ decreases as the overtone label, $n$, increases Zimmerman and Mark
(2016); Berti  _et al._ (2003).

These three qualitative features are visible in Fig. 3:

  1. (i)

In the left panel ($n=0$), and for $s=-2$, we see that the monomial moments
$|\left<\xi^{p}\right>|$ decrease monotonically with $p$. However, for $s=-2$
and $n\in\\{2,6\\}$, moments are non-monotonic, and at $n=6$, they increase as
$p$ increases.

  2. (ii)

In each panel, we see that monomial moments for spin weight $s=+2$ are
significantly depressed relative to those for $s=-2$.

  3. (iii)

The scale of each panel is identical so that they effect of increasing
overtone label (i.e. decreasing $\tau$) is easily visible for both spin
weights.

The fact that absolute monomial moments, $|\left<\xi^{p}\right>|$, are not
strictly decreasing with $p$ has implications for the numerical behavior of
scalar products. Naively, it might seem that, if $|\left<\xi^{p}\right>|$
increases with $p$, then scalar products might diverge. This happens to not be
the case because the monomial moments are generally complex valued.

|
---|---
|
Figure 4:  The first four confluent Heun polynomials for $p=k$, spin weight
$s=-2$, black hole dimensionless spin $a/M=0.86$, and $\tilde{\omega}$ defined
by the $(\ell,m,n)=(3,3,1)$ QNM. The corresponding QNM frequency is
$M\tilde{\omega}=0.9840-0.2140i$. Polynomial orders $p\in\\{0,1,2,3\\}$ are
shown as a function of the physically valued fractional radius $\xi$. The top
left and right panels shows the real and imaginary parts of each polynomial,
respectively. Bottom left and right panels show the corresponding polynomial
amplitudes and phases, respectively.

We have previously noted that the real and imaginary parts of
$\left<\xi^{p}\right>$ will trace out a logarithmic spiral in $p$. This
implies that linear combinations of monomial moments will, in general, have a
telescoping quality: as $p$ increases, successive terms will interfere
destructively, meaning that although $|\left<\xi^{p}\right>|$ may increase,
its net contribution must be considered along with adjacent terms. This
telescoping quality behoves the use of arbitrary precision (e.g. mpmath mpmath
development team (2023)) when evaluating moments by analytic
continuation141414Leaver came to the same conclusion during his analysis of
series expansions in $\xi$, which he refers to as Jaffé-type expansion Leaver
(1986). Leaver (1986); Ronveaux and Arscott (1995). However, this quality also
prompts a question.

Since the weight function, $\text{W}(\xi)$, allows the computation of
arbitrary monomial moments, the monomials are a basis for the scalar product
space defined by Eq. (64) Adkins and Weintraub (2012); Gohberg  _et al._
(2006). Thus, it is fair to wonder whether there exists a related basis that
is more suitable for numerical approximation? This is the basic question that
links the study of monomial moments to applied polynomial theory Chihara
(2011). A complete investigation of this question is beyond the current scope.
Nevertheless, a first step towards an answer may be found by applying the
scalar product to confluent Heun polynomials Whiting (1989); Cook and
Zalutskiy (2014); Fiziev (2009b); DLMF (b); Magnus  _et al._ (2021); Fiziev
(2009a).

###  VII.2 Confluent Heun polynomials

The differential operator that results from applying QNM boundary conditions
to Teukolsky’s radial equation was

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

Since the operator given in Eq. (71) is of the confluent Heun type, its eigen-
solutions are the aptly named confluent Heun functions Ronveaux and Arscott
(1995); Magnus  _et al._ (2021).

The subset of these with terminating series solutions are, unsurprisingly,
named confluent Heun polynomials DLMF (a); Cook and Zalutskiy (2014); Fiziev
(2009b); Ronveaux and Arscott (1995).

A confluent Heun polynomial of order $p$ is typically considered to result
from special values of $\text{C}_{0}$ and $\text{C}_{1}$ Cook and Zalutskiy
(2014); Fiziev (2009b). In the context of Teukolsky’s radial equation, these
special values have been interpreted as a constraint on $\tilde{\omega}$. This
is the case for the algebraically special modes of Schwarzschild and the total
transmission modes of Kerr Cook and Zalutskiy (2014); Berti  _et al._ (2003);
Fiziev (2009a). Here, a different perspective on the confluent Heun polynomial
conditions for BHs will be used. This perspective is based on two
observations:

  1. (i)

Classical polynomials are typically defined to be eigenfunctions of a purely
differential operator (i.e. one with no $\partial_{\xi}^{0}$ term). For the
confluent Heun operator, Eq. (71), this purely differential part is denoted by
Eqs. (71b) and (71c), and repeated below in Eq. (72),

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

  2. (ii)

Any attempt to construct order $p$ polynomials that are eigenfunctions of
$\mathcal{D}_{\xi}$ requires that the “eigenvalue” be a linear function of the
domain variable,

$\displaystyle\mathcal{D}{{}_{\xi}}\;\mathbf{y}_{pk}(\xi)\;$
$\displaystyle=\;\sigma_{pk}(\xi)\;\mathbf{y}_{pk}(\xi)\;,$ (73a)
$\displaystyle\sigma_{pk}(\xi)\;$
$\displaystyle=\;\lambda_{pk}\;+\;\mu_{p}\,\xi\;.$ (73b)

Plainly, $\sigma_{pk}(\xi)$ is not an eigenvalue in the traditional sense.
Instead, it is an instance of a two parameter eigenvalue problem Ronveaux and
Arscott (1995); Atkinson and Mingarelli (2011). We will see that $\mu_{p}$ is
uniquely determined by the order of the polynomial, meaning that it need not
be conceptually separated from the traditional eigenvalue, $\lambda_{pk}$.

In what follows, we will apply these ideas to construct confluent Heun
polynomials that are valid for any value of $\tilde{\omega}$ (i.e. any QNM
frequency). We will see that the operator, $\mathcal{D}{{}_{\xi}}$, is self-
adjoint under the same scalar product that we have discussed previously. We
will then use the self-adjointness of this operator to prove and numerically
demonstrate the orthogonality of confluent Heun polynomials for fixed
polynomial order. Lastly, we will use self-adjointness to derive an expression
for the confluent Heun polynomials’ eigenvalues in terms of ratios of scalar
products.

####  VII.2.1 Construction

Let $\mathbf{y}_{pk}(\xi)$ be a confluent Heun polynomial of order $p$,

$\displaystyle\mathbf{y}_{pk}(\xi)\;=\;\sum_{q=0}^{p}\,\text{a}_{pkq}\,\xi^{q}\;.$
(74)

We will soon see that the label, $k$, is useful because there are will
generally be $p+1$ order $p$ polynomials.

The polynomial coefficients, $\text{a}_{pkq}$, are determined by applying the
form above for $\mathbf{y}_{pk}(\xi)$ to the eigen-relationship given by Eq.
(73a). One then asserts that the coefficients of each power in $\xi^{q}$ must
be equal on both sides of Eq. (73a). The result is a three term recursion
relationship or, equivalently, a discrete eigenvalue relationship. To simplify
our notation, we will henceforth drop two indices, $p$ and $k$, from
$\text{a}_{pkq}$,

$\displaystyle\text{a}_{q}\;\leftarrow\;\text{a}_{pkq}\;.$ (75)

Thus the discrete eigenvalue relationship takes the schematic form,

$\displaystyle\text{a}_{q-1}\,\alpha_{q}\;+\;\text{a}_{q}\,\beta_{q}\;+\;\text{a}_{q+1}\,\gamma_{q}\;=\;\lambda_{pk}\,\text{a}_{q}\;,$
(76)

where

$\displaystyle\alpha_{q}\;$ $\displaystyle=\;-(q-1)(\text{C}_{4}+q-2)-\mu_{p}$
(77a) $\displaystyle\beta_{q}\;$
$\displaystyle=\;q\,(\text{C}_{3}+2(\text{C}_{4}+q-1))$ (77b)
$\displaystyle\gamma_{q}\;$
$\displaystyle=\;-(q+1)(\text{C}_{2}+\text{C}_{3}+\text{C}_{4}+q)\;.$ (77c)

An initial boundary condition,

$\displaystyle\text{a}_{0}\,\beta_{0}\;+\;\text{a}_{1}\,\gamma_{0}\;=\;\lambda_{pk}\,\text{a}_{0}\;,$
(78)

is imposed for consistency with Eq. (74). In Eq. (78), $\text{a}_{0}=1$ is a
common initial normalization choice.

A necessary but insufficient condition for the three-term recursion to
terminate at the $p^{\text{th}}$ term, is simply that $\alpha_{p+1}=0$. This
is the so-called $\Delta_{p+1}=0$ condition DLMF (b); Cook and Zalutskiy
(2014). In our notation, this takes the form,

$\displaystyle\alpha_{p+1}\;$
$\displaystyle=\;p\,(\,p+\text{C}_{4}-1)-\mu_{p}\;=\;0\;.$ (79)

In turn, this can be interpreted as a constraint that allows the determination
of $\mu_{p}$,

$\displaystyle\mu_{p}\;$ $\displaystyle=\;p\,(\,p+\text{C}_{4}-1\,)\;.$ (80)

Thus $\mu_{p}$ must be non-zero for the existence of non-trivial polynomial
solutions to Eq. (73a).

Furthermore, if $\mu_{p}$ is given by Eq. (80), then the discrete eigenvalue
relationship given by Eq. (76) defines a finite dimensional matrix eigenvalue
problem151515 Note that, along with the series’ boundary conditions (Eqs.
78-79), it is the finite nature of the system of equations that justifies the
eigenproblem interpretation. In contrast, non-polynomial solutions to the
radial problem will not satisfy Eq. (79), meaning that e.g. a continued
fraction based method, and a related boundary condition at the infinite end of
the series are required Leaver (1986, 1985).  Ronveaux and Arscott (1995);
DLMF (b); Cook and Zalutskiy (2014). The matrix of interest,
$\hat{\mathbf{A}}_{p}$, is asymmetric and tridiagonal,

$\hat{\mathbf{A}}_{p}\;=\;\left(\begin{array}[]{ccccc}\beta_{0}&\gamma_{0}&0&\cdots&0\\\
\alpha_{1}&\beta_{1}&\gamma_{1}&\ddots&0\\\
0&\alpha_{2}&\beta_{2}&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\gamma_{p-1}\\\
0&0&\cdots&\alpha_{p}&\beta_{p}\end{array}\right)\;.$ (81)

By construction, the eigenvectors of $\hat{\mathbf{A}}_{p}$ are the polynomial
coefficients, and the eigenvalues of $\hat{\mathbf{A}}_{p}$ are the polynomial
eigenvalues,

$\displaystyle\\{\text{a}_{pkq}\\}_{q=0}^{p}\;$
$\displaystyle=\;{k}^{\text{th}}\text{ Eigenvector of
}\hat{\mathbf{A}}_{p}\;,$ (82a) $\displaystyle\lambda_{pk}\;$
$\displaystyle=\;{k}^{\text{th}}\text{ Eigenvalue of }\hat{\mathbf{A}}_{p}\;.$
(82b)

Since $\hat{\mathbf{A}}_{p}$ is a $(p+1)\times(p+1)$ matrix, there will be
$p+1$ eigenvalues, $\lambda_{pk}$. If the eigenvalues are unique, then we may
impose an order on $\mathbf{y}_{pk}$, whereby

$\displaystyle k\in\\{0,1,2,...,p\\}\;.$ (83)

In general, the $p+1$ values of $\lambda_{pk}$ will correspond to the roots of
$\hat{\mathbf{A}}_{p}$’s characteristic polynomial. Lacking a deeper
understanding how the roots of all such polynomials are distributed, the
following ordering will henceforth be used

$\displaystyle\text{Im}\,\lambda_{p0}\,<\,\text{Im}\,\lambda_{p1}\,<\,\text{Im}\,\lambda_{p2}\,<\,...\,<\,\text{Im}\,\lambda_{pp}\;.$
(84)

In this sense, for a given order $p$, there is a multiplex comprised of $p+1$
solutions, where each solution is labeled by the multiplex index, $k$.

####  VII.2.2 Orthogonality between polynomials of like order

If each $\lambda_{pk}$ is unique, then each corresponding polynomial,
$\mathbf{y}_{pk}$, will be linearly independent Axler (2015). Since
$\mathcal{L}_{\xi}=\mathcal{D}_{\xi}+(\text{C}_{0}+(1-\xi)\text{C}_{1})$ is
self-adjoint with respect to the scalar product, Eq. (64), $\mathcal{D}_{\xi}$
is also self-adjoint,

$\displaystyle\mathcal{D}{{}_{\xi}}\;$
$\displaystyle=\;{\mathcal{D}{{}_{\xi}}}^{\dagger}\;.$ (85)

Here, we will use the self-adjointness of $\mathcal{D}_{\xi}$ to show that
when $\lambda_{pk}$ are unique, $\mathbf{y}_{pk}$ must be orthogonal. To that
end we will use a standard argument that is based on the behavior of
${\langle{\mathbf{y}_{p^{\prime}k^{\prime}}}\,|\,{\mathcal{D}_{\xi}\mathbf{y}_{pk}}\rangle}$
when $(p^{\prime},k^{\prime})$ swaps places with $(p,k)$. Along the way, we
will encounter a twist that results from the fact that $\mathbf{y}_{pk}$
satisfy a non-standard eigenvalue equation, Eq. (73a).

Let us begin by analyzing the matrix whose elements are
${\langle{\mathbf{y}_{p^{\prime}k^{\prime}}}\,|\,{\mathcal{D}_{\xi}\mathbf{y}_{pk}}\rangle}$.
Since $\mathcal{D}_{\xi}$ is self-adjoint with respect to the scalar product
(Eq. 26), we have that

$\displaystyle{\langle{\mathbf{y}_{p^{\prime}k^{\prime}}}\,|\,{\mathcal{D}{{}_{\xi}}\,\mathbf{y}_{pk}}\rangle}\;$
$\displaystyle=\;{\langle{\mathcal{D}{{}_{\xi}}\,\mathbf{y}_{p^{\prime}k^{\prime}}}\,|\,{\mathbf{y}_{pk}}\rangle}\;.$
(86)

Applying the eigenvalue relationship,

$\displaystyle\mathcal{D}{{}_{\xi}}\;\mathbf{y}_{pk}(\xi)\;$
$\displaystyle=\;(\lambda_{pk}\;+\;\mu_{p}\,\xi)\;\mathbf{y}_{pk}(\xi)$ (87)

to both sides of Eq. (86) yields

$\displaystyle\lambda_{pk}$
$\displaystyle{\langle{\mathbf{y}_{p^{\prime}k^{\prime}}}\,|\,{\mathbf{y}_{pk}}\rangle}\;+\;\mu_{p}\,\langle{\mathbf{y}_{p^{\prime}k^{\prime}}}|\xi\;|{\mathbf{y}_{pk}}\rangle$
(88)
$\displaystyle=\;\lambda_{p^{\prime}k^{\prime}}{\langle{\mathbf{y}_{p^{\prime}k^{\prime}}}\,|\,{\mathbf{y}_{pk}}\rangle}\;+\;\mu_{p^{\prime}}\,\langle{\mathbf{y}_{p^{\prime}k^{\prime}}}|\xi\;|{\mathbf{y}_{pk}}\rangle\;.$

A useful rearrangement of Eq. (88) is given by Eq. (89),

$\displaystyle(\lambda_{pk}-\lambda_{p^{\prime}k^{\prime}})\,$
$\displaystyle{\langle{\mathbf{y}_{p^{\prime}k^{\prime}}}\,|\,{\mathbf{y}_{pk}}\rangle}\;=\;(\mu_{p^{\prime}}-\mu_{p})\,\langle{\mathbf{y}_{p^{\prime}k^{\prime}}}|\xi\;|{\mathbf{y}_{pk}}\rangle\;.$
(89)

The aforementioned twist is that, unlike in single parameter eigenvalue
problems, the right-hand-side of Eq. (89) is non-zero. Nevertheless, the
standard orthogonality argument is recovered by restricting Eq. (89) to
polynomials of like order,

$\displaystyle p\;$ $\displaystyle=\;p^{\prime}\;.$ (90)

The resulting version of Eq. (89) is,

$\displaystyle(\lambda_{pk}-\lambda_{pk^{\prime}})\,{\langle{\mathbf{y}_{pk^{\prime}}}\,|\,{\mathbf{y}_{pk}}\rangle}\;$
$\displaystyle=\;0\;.$ (91)

If $\lambda_{pk}$ are distinct for different $k$, then Eq. (91) can only hold
if

$\displaystyle{\langle{\mathbf{y}_{pk^{\prime}}}\,|\,{\mathbf{y}_{pk}}\rangle}\;$
$\displaystyle=\;\eta_{pk}^{2}\,\delta_{k^{\prime}k}\;.$ (92)

In Eq. (92), $\eta_{pk}$ is an inverse normalization constant. Henceforth, we
will work with the normalized confluent Heun polynomials,

$\displaystyle{y}_{pk}\;$ $\displaystyle=\;\eta_{pk}^{-1}\;\mathbf{y}_{pk}\;,$
(93)

where,

$\displaystyle{\langle{{y}_{pk^{\prime}}}\,|\,{{y}_{pk}}\rangle}\;$
$\displaystyle=\;\delta_{k^{\prime}k}\;.$ (94)

|  |
---|---|---
Figure 5:  Distribution of roots for confluent Heun polynomials,
$y_{pk}(\xi)$, for three cases, all with $s=-2$, $a/M=0.7$, and $\ell=m=2$.
The three cases differ in polynomial order, $p$, and the frequency,
$\tilde{\omega}$, used: (left) $y_{33}$, the $n=0$ QNM frequency,
$M\tilde{\omega}=0.5326-0.0808i$, was used, (center) $y_{77}$, the $n=0$ QNM
frequency, $M\tilde{\omega}=0.4713-0.5843i$, was used, and (right)
$y_{20,20}$, the $n=12$ QNM frequency, $M\tilde{\omega}=0.4155-2.5050i$, was
used. Polynomial roots are marked by white open circles. Roots were calculated
using numpy.roots Horn and Johnson (1990); Harris  _et al._ (2020). Values of
polynomial phase, $\arg(y_{pk})$, are shaded from blue to yellow, indicating
values between $-\pi$ and $\pi$, respectively. Dotted horizontal and vertical
lines mark the respective locations of the real and imaginary axes.

####  VII.2.3 A scalar product expression for eigenvalues

For polynomials of unlike order, Eq. (89) may be recast as a relationship
between polynomial eigenvalues. This is done by solving for $\lambda_{pk}$,

$\displaystyle\lambda_{pk}\;=\;\lambda_{p^{\prime}k^{\prime}}\,+\,(\mu_{p^{\prime}}-\mu_{p})\,\frac{\langle{{y}_{p^{\prime}k^{\prime}}}|\xi\;|{{y}_{pk}}\rangle}{{\langle{{y}_{p^{\prime}k^{\prime}}}\,|\,{{y}_{pk}}\rangle}}\;.$
(95)

Equation (95) applies to all $(p,k)$ and $(p^{\prime},k^{\prime})$, where
$p\neq p^{\prime}$. An instructive special case is that of

$\displaystyle p^{\prime}=0\;.$ (96)

This corresponds to ${y}_{00}$, which must be constant. If ${y}_{00}$ is a
constant, then for the eigenvalue relation (Eq. 87) to hold, it must be true
that

$\displaystyle\lambda_{00}\;$ $\displaystyle=\;\mu_{0}\;=\;0\;.$ (97)

Applying Eq. (97) to Eq. (95) gives

$\displaystyle\lambda_{pk}\,\;$
$\displaystyle=\;\mu_{p}\,\frac{\langle{y_{00}}|\xi\;|{y_{pk}}\rangle}{{\langle{y_{00}}\,|\,{y_{pk}}\rangle}}\;.$
(98)

Recalling that $\mu_{p}=p\,(\,p+\text{C}_{4}-1\,)$, Eq. (98) is equivalent to

$\displaystyle\lambda_{pk}\,\;$
$\displaystyle=\;p\,(\,p+\text{C}_{4}-1\,)\,\frac{\langle{y_{00}}|\xi\;|{y_{pk}}\rangle}{{\langle{y_{00}}\,|\,{y_{pk}}\rangle}}\;.$
(99)

If $y_{pk}$ are known, then Eq. (99) allows the computation of the related
eigenvalues via the scalar product.

####  VII.2.4 Numerical examples

Numerical computation of the confluent Heun polynomials may proceed by first
computing $\text{C}_{2}$, $\text{C}_{3}$ and $\text{C}_{4}$. Then, for a
chosen polynomial order, one would compute $\hat{\mathbf{A}}_{p}$ (Eq. 81), so
that its eigenvalues and vectors may be determined. Since scalar products
between polynomials is a weighted sum over monomial moments (Eq. 64), it is
useful to precompute and store all of the monomial moments relevant to a given
polynomial order. This allows all polynomial scalar products to be evaluated
without explicit integration (Eq. 47). From this perspective, numerical
evaluation of confluent Heun polynomials may be accomplished on-the-fly via
the appropriate sum over monomials (Eq. 74). Here, this procedure is used to
provide select numerical examples.

Figure 4 shows example confluent Heun polynomials of orders
$p\in\\{0,1,2,3\\}$. Only the polynomials for which $p=k$ are shown; others
have qualitatively similar behavior. As expected, the real and imaginary parts
of each case behave as polynomials in $\xi$. However, since each polynomial
has complex value coefficients, the polynomials’ amplitudes and phases are
non-polynomial functions of $\xi$. Similarly, since the polynomials are
complex valued, we should not expect the roots of their respective real and
imaginary parts to be simultaneous. For example, the top left panel of Fig. 4
shows $\text{Re}\,y_{11}$ having a root near $\xi=0.5$, but
$\text{Im}\,y_{11}$ having no root at all on the real line.

Figure 5 shows the distribution of polynomial roots for three cases. In focus
is the distribution of roots in the complex plane as polynomial order and
overtone label increase. As expected, roots are typically not real valued. For
large $n$ and large polynomial orders, roots are found to form two families,
one situated around the origin, and another dispersed across the upper-right
quadrant. This is seen in the right panel mof Fig. 5. In all panels of Fig. 5,
each root is seen along with its associated phase boundary. Any smooth curve,
e.g. $\xi(z)$, may be defined to have a continuous phase by adding the
appropriate integer factor of $2\pi$ at each phase boundary.

Figure 6 shows examples of orthogonality at fixed polynomial order,

$\displaystyle p\;=\;20\;.$ (100)

for three cases. At focus, is the effect of increasing QNM overtone label or,
equivalent, decreasing the QNM decay time $\tau$, such that $1/\tau$ increases
within $\tilde{\omega}=\omega-i/\tau$. The $p=20$ polynomials are considered
to be somewhat extreme cases: for this and all lower polynomial orders, no
appreciable numerical artifacts are encountered for the numerical and physical
parameters used here.

The left columns of Fig. 6 show the absolute values of Gramian matrix
elements, $\langle y_{pk^{\prime}}|y_{pk}\rangle$, and the right columns show
absolute values of $\langle y_{pk^{\prime}}|\xi|y_{pk}\rangle$. Scalar
products have been evaluated using analytic continuation (Eq. 64) with $64$
significant figures of precision via the dps setting of mpmath mpmath
development team (2023).

Orthogonality of the polynomials at fixed order is seen through the diagonal
structure of each Gramian. The top left panel, having $n=0$, displays a slight
deviation away from orthogonality (i.e. values near $10^{-40}$ rather than
$10^{-50}$ and below) as $k$ nears it maximum value of $20$.

This effect has two causes, one technical and the other mathematical. The
technical cause of this feature is simply the value of the precision
parameter, dps$=64$, which means that numerical results near $10^{-64}$ are
more prone to round-off error. The mathematical cause of this feature is that,
like the legendre polynomials, the confluent Heun polynomials have
normalization constants that generally decrease as polynomial order
increases161616For the Legendre polynomials, the analog of $\eta_{pk}^{2}$
(see Eq. 92) is $(2p+1)^{-1}$.. Consequently the action of normalization
effectively scales each polynomial by a very large number, exacerbating finite
precision errors. At worst, square normalization constants might be $10^{-64}$
or smaller and would thereby be corrupted by finite precision error. In turn,
this would correspond to normalization constants $10^{-32}$ or smaller, and
the act of normalization would be to scale the non-normalized polynomials by a
corrupted number of order $10^{32}$. This situation can be avoided by
increasing numerical precision.

The confluent Heun polynomials, while orthogonal at fixed order, do not
display a feature common among classical polynomials. In particular, they do
not possess three-term-recursion relations between polynomials of different
values of $k$ Chihara (2011); Arfken  _et al._ (2013). This may be understood
thusly: Since $\xi\times y_{pk}$ is a polynomial or order $p+1$, it is to be
expected that it cannot be written in terms of only polynomials of order $p$.
Consequently, the matrix whose elements are $\langle
y_{pk^{\prime}}|\xi|y_{pk}\rangle$ will in general have no non-zero values,
and therefore will not be band-diagonal.

This feature of the confluent Heun polynomials is shown in the right column of
Fig. 6. There, the structure of $|\langle y_{pk^{\prime}}|\xi|y_{pk}\rangle|$
is seen to increase non-trivially with overtone label, $n$. This may be
partially understood in the context of the polynomial eigenvalues.

It was previously noted that the small deviations from orthogonality are
linked to small normalization constants, and that increasing structure in the
distribution of $|\langle y_{pk^{\prime}}|\xi|y_{pk}\rangle|$ are linked to
the polynomial eigenvalues. Polynomial normalization constants and eigenvalues
are shown, respectively, in the left and right panels of Fig. 7. There, the
same three cases seen in Fig. 6 are considered.

At focus in the left panel of Fig. 7 is the distribution of normalization
constants. While normalization constants are found to generally decrease with
polynomial order (not shown in figure), they are also found to generally
decrease as $k$ increases, and they are found to generally increase as the
imaginary part $\tilde{\omega}$ increases. Branching in the distribution of
normalization constants is linked to branching in the respective polynomial
eigenvalues. Similarly, Eq. (95) links eigenvalue branching to the
distribution of $|\langle y_{pk^{\prime}}|\xi|y_{pk}\rangle|$ seen in the
right column of Fig. 6.

At focus in the right panel of Fig. 7 is how the distribution of eigenvalues
changes as the imaginary part of $\tilde{\omega}$ becomes more negative (i.e.
as $n$ increases). For $n=0$, either $\text{Re}\,\tilde{\omega}$ or
$\text{Im}\,\tilde{\omega}$ may be used to define equivalent orders (i.e. maps
between polynomials and the index $k$). This is not the case for
$n\in\\{6,12\\}$, where branches appear in the eigenvalue distribution. In
these cases, the eigenvalues are still unique, so a well defined order can
still be constructed.

|
---|---
|
|
Figure 6:  Gramian matrices for three QNM overtone labels: $n=0$ (top row), $n=6$ (middle row), and $n=12$ (bottom row). Considered here are the order $p=20$ normalized confluent Heun polynomials with $s=2$, $a/M=0.70$ and QNM frequency labels $(\ell,m)=(2,2)$. The QNM frequencies for top middle and bottom rows are $M\tilde{\omega}=0.5326-0.0808i$, $M\tilde{\omega}=0.4239-1.0954i$ and $M\tilde{\omega}=0.4098-3.1270i$, respectively. The left column shows the absolute values of $\langle y_{pk^{\prime}}|y_{pk}\rangle$; there, note that the color bar is log scaled between $10^{-50}$ and $10^{0}$. The right column shows absolute values of $\langle y_{pk^{\prime}}|\xi|y_{pk}\rangle$; there, the color bar is logs scaled between $10^{-12}$ and $10^{0}$ to accentuate graph structure.  |
---|---
Figure 7:  Distribution of normalization constants, $\eta_{pk}$, and
eigenvalues,
$\lambda_{pk}=\mathrm{Re}\lambda_{pk}+i\,\mathrm{Im}\lambda_{pk}$, for three
different QNM overtone labels: $n\in\\{0,6,12\\}$. For all values of $n$, the
spin weight is $s=2$, the black hole dimensionless spin is $a/M=0.70$, and the
QNM $(\ell,m)$ are $(2,2)$. The left panel shows normalization constants, and
the right panel shows eigenvalues. Black dots mark $n=0$, orange circles mark
$n=6$, and magenta squares mark $n=12$.

##  VIII Implications for analytic solutions to Teukolsky’s radial equation

We previously set out to apply two concepts to Teukolsky’s radial problem.
When stated in the context of the radial equation, these two concepts are:

  * (i)

That a scalar product may be developed from the structure of the radial
equation (i.e. (_Q_ 1)).

  * (ii)

That special functions may be developed from the structure of the radial
equation, and that these functions may be used to exactly represent analytic
solutions (i.e. (_Q_ 2)).

Previous sections show that concept (i) is indeed applicable to the radial
equation, thus affirmatively answering (_Q_ 1). A radial scalar product was
developed in Sec. IV. Methods for evaluation of the scalar product were
developed in Sec. V. Connections between the scalar product and concept (ii)
were discussed in Sec. VII. In particular, the scalar products connection with
problem specific special functions (i.e. the confluent Heun polynomials) was
detailed in Sec. VII.2. Here, we take a first step towards applying confluent
Heun polynomials to the representation of analytic solutions to Teukolsky’s
radial equation for QNMs. A non-classical feature of the confluent Heun
polynomials means that, for now, only a preliminary and conceptual step is
taken.

Recall that Teukolsky’s radial equation is a confluent Heun equation of the
form,

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

In Eq. (101), the transformed radial operator, $\mathcal{L}{{}_{\xi}}$, is
simply rewritten such the differential part is explicitly identified as
$\mathcal{D}_{\xi}$ (see Eq. 72). In Sec. VII.2, it was found that
$\mathcal{D}_{\xi}$ generates confluent Heun polynomials, $y_{pk}(\xi)$, such
that each $y_{pk}(\xi)$ satisfies a two parameter eigenvalue problem,

$\displaystyle\mathcal{D}{{}_{\xi}}\;|{{y}_{pk}}\rangle\;$
$\displaystyle=\;(\lambda_{pk}\;+\;\mu_{p}\,\xi)\;|{{y}_{pk}}\rangle\;,$ (102)

where, for a given polynomial order $p$, there are $p+1$ solutions, $y_{pk}$.
This situation differs starkly from that of classical (e.g. Jacobi)
polynomials, for which every single polynomial order maps to a single
polynomial solution. In other words, for each polynomial order, classical
polynomials have a simplex of solutions, while the confluent Heun polynomials
have a multiplex.

The existence of a multiplex for each polynomial order prompts a nontrivial
conceptual tension: The Jacobi polynomials’ recursion relationships play a key
role in the angular problem (see Appx. B); however, there cannot exist similar
recursion relationships between the confluent Heun polynomials, $y_{pk}$. This
is because such recursion relationships require an orthogonal polynomial
sequence, wherein each polynomial has a unique order (i.e. is a simplex)
Chihara (2011). Further, since each confluent Heun polynomial multiplex is
composed of polynomials of like order, the collection of all such multiplexes
would appear to contain too many polynomials, meaning that the they allow for
multiple equivalent representations (i.e. decompositions) of radial functions.

This does not, however, preclude the possibility that there are special
combinations of confluent Heun polynomials allow for unique representations,
and classical recursion relationships. While this possibility is perhaps too
open ended to be fully pursued here, we may assume (somewhat adventurously)
that such polynomials exist, and then contemplate implications for Teukolsky’s
radial equation:

Suppose that solutions to the radial equation, $|{f}\rangle$, can be exactly
represented in terms of orthonormal polynomials, $u_{n}$,

$\displaystyle|{f}\rangle\;=\;\sum_{p=0}^{\infty}a_{p}\,|{u_{p}}\rangle\;,$
(103a) where it is additionally supposed that $\displaystyle|{u_{p}}\rangle\;$
$\displaystyle=\;\sum_{k=0}^{p}\,\text{c}_{pk}\,|{{y}_{pk}}\rangle\;,$ (103b)
$\displaystyle{\mathbb{I}}\;$
$\displaystyle=\;\sum_{p=0}^{\infty}\,|{u_{p}}\rangle\langle{u_{p}}|\;,$
(103c) and
$\displaystyle\xi|{u_{p}}\rangle=\sigma_{p0}|{u_{p-1}}\rangle+\sigma_{p1}|{u_{p}}\rangle+\sigma_{p2}|{u_{p+1}}\rangle\;.$
(103d)

In Eqs. (103a-103d), despite not being classical, $u_{n}$ are constructed to
have the essential properties of classical polynomials, namely completeness,
orthonormality and three-term recursion Chihara (2011); Chen and Dai (2010).
In this sense, we may think of $u_{n}$ as having canonical properties. This
motivates their being referred to as canonical confluent Heun polynomials.

Further, in Eqs. (103a-103d), $a_{p}$, $c_{pk}$, $\sigma_{p0}$, $\sigma_{p1}$
and $\sigma_{p2}$ are assumed to be generally non-zero complex valued
constants. Henceforth, $\sigma_{pj}$ will refer to complex valued constants
that are assumed to be known if $u_{p}$ are known. In Eq. (103c),
${\mathbb{I}}$ is the identity operator for the space of all confluent Heun
functions of the fractional radial coordinate, $\xi$. Equation (103c) requires
that the expansion coefficients present in Eq. (103a) are
$a_{p}={\langle{u_{p}}\,|\,{f}\rangle}$. The three-term-recursion relation
given by Eq. (103d) is a universal feature of all (simplexical) orthogonal
polynomial sequences, and it is a direct consequence of Eq. (103c) (i.e.
completeness) along with the fact that
${\langle{u_{j}}\,|\,{\xi|u_{k}}\rangle}={\langle{\xi
u_{j}}\,|\,{u_{k}}\rangle}={\langle{u_{j}}\,|\,{\xi u_{k}}\rangle}$ Chihara
(2011).

Naively, if each order $p$ polynomial multiplex is a $p+1$ dimensional
orthonormal basis, then it should be capable of exactly representing any order
$p$ polynomial on the same domain Axler (2015). Thus it may be possible to
construct $u_{p}$ such that Eqs. (103b-103d) hold.

In effect, Eq. (103) fully defines a solution ansatz. This solution ansatz may
be used to investigate whether the scalar product and confluent Heun
polynomials imply that Teukolsky’s radial equation is, like the angular
equation, exactly tridiagonalizable. To this end, we may attempt to represent
the radial equation in the basis defined by $\\{u_{p}\\}_{p=0}^{\infty}$.

In bra-ket notation, the transformed radial equation is

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

Applying Eq. (103c) such that $\mathcal{L}_{\xi}$ is represented as a matrix,
$\hat{L}$, and $|{f}\rangle$ a vector, $\vec{f}$, gives

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

In Eq. (105b), the elements of $\hat{L}$ are
$\langle{u_{p^{\prime}}}|\,\mathcal{L}{{}_{\xi}}\,|{u_{p}}\rangle$, and the
elements of $\vec{f}$ are ${\langle{u_{p}}\,|\,{f}\rangle}$. Henceforth,
$p^{\prime}$ will label rows, and $p$ will label columns, meaning that $p\geq
p^{\prime}$ denotes the upper triangle of $\hat{L}$, and $p\leq p^{\prime}$
its lower. To compute the matrix elements,
$\langle{u_{p^{\prime}}}|\,\mathcal{L}{{}_{\xi}}\,|{u_{p}}\rangle$, the action
of $\mathcal{L}{{}_{\xi}}$ on $|{u_{p}}\rangle$ must be considered in more
detail.

Rewriting $\mathcal{L}{{}_{\xi}}$ in terms of $\mathcal{D}_{\xi}$ (see Eq.
101) allows $\mathcal{L}{{}_{\xi}}\,|{u_{p}}\rangle$ to be expanded as

$\displaystyle\mathcal{L}{{}_{\xi}}\,|{u_{p}}\rangle\;=\;(\text{C}_{0}+\text{C}_{1}[1-\xi])\,|{u_{p}}\rangle+\mathcal{D}_{\xi}\,|{u_{p}}\rangle\;.$
(106)

In Eq. (106), $\mathcal{D}_{\xi}\,|{u_{p}}\rangle$ can be rewritten using the
definition of $u_{p}$ in terms of confluent Heun polynomials (Eq. 103b), and
the two parameter eigenvalue relationship (Eq. 102). Doing so yields

$\displaystyle\mathcal{D}_{\xi}\,|{u_{p}}\rangle\;$
$\displaystyle=\;\sum_{k=0}^{p}\,\text{c}_{pk}\,\mathcal{D}_{\xi}\,|{{y}_{pk}}\rangle$
(107a)
$\displaystyle=\;\sum_{k=0}^{p}\,\text{c}_{pk}\,(\lambda_{pk}\,+\,\mu_{p}\,\xi)\,|{{y}_{pk}}\rangle\;,$
(107b) $\displaystyle=\;|{v_{p}}\rangle\,+\,\mu_{p}\,\xi\,|{u_{p}}\rangle\;.$
(107c)

In going from Eq. (107b) to Eq. (107c), the term dependent on the confluent
Heun eigenvalues, $\lambda_{pk}$, has been used to define a new polynomial,

$\displaystyle|{v_{p}}\rangle\;=\;\sum_{k=0}^{p}\,\text{c}_{pk}\,\lambda_{pk}\,|{{y}_{pk}}\rangle\;,$
(108)

and it has been noticed that the remaining term in Eq. (107b) is simply
proportional to $|{u_{n}}\rangle$.

It follows from Eq. (107c) that $\mathcal{L}{{}_{\xi}}\,|{u_{p}}\rangle$ may
be rewritten as

$\displaystyle\mathcal{L}{{}_{\xi}}|{u_{p}}\rangle=(\text{C}_{0}+\text{C}_{1})|{u_{p}}\rangle+|{v_{p}}\rangle+(\mu_{p}-\text{C}_{1})\xi|{u_{p}}\rangle\;,$
(109)

and that the matrix element,
$\langle{u_{p^{\prime}}}|\mathcal{L}{{}_{\xi}}|{u_{p}}\rangle$, is

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

Since $|{u_{p}}\rangle$ are orthonormal (i.e.
${\langle{u_{p^{\prime}}}\,|\,{u_{p}}\rangle}=\delta_{p^{\prime}p}$), Eq.
(110a) only contributes to diagonal matrix elements. Since $|{u_{p}}\rangle$
possess a three term recursion relation given by Eq. (103d), Eq. (110b) will
only contribute when $|p^{\prime}-p|\leq 1$.

It happens that ${\langle{u_{p^{\prime}}}\,|\,{v_{p}}\rangle}$ is not only
non-zero when $|p^{\prime}-p|\leq 1$, but also restricted to the upper-
triangle (i.e. $p\geq p^{\prime}$). This follows from the fact that
$\mathcal{D}_{\xi}$ is self-adjoint w.r.t. the scalar product, meaning that

$\displaystyle{\langle{u_{p^{\prime}}}\,|\,{\mathcal{D}_{\xi}\,u_{p}}\rangle}\;=\;{\langle{\mathcal{D}_{\xi}\,u_{p^{\prime}}}\,|\,{u_{p}}\rangle}\;.$
(111)

Using Eq. (107c) to expand both sides of Eq. (111) yields,

$\displaystyle{\langle{u_{p^{\prime}}}\,|\,{v_{p}}\rangle}-{\langle{v_{p^{\prime}}}\,|\,{u_{p}}\rangle}\;=\;(\mu_{p}-\mu_{p^{\prime}}){\langle{u_{p^{\prime}}}\,|\,{\xi|u_{p}}\rangle}\;.$
(112)

Equation (112) may be understood thusly: Since $u_{p^{\prime}}$ is a member of
a complete orthonormal sequence, its scalar product with any order $p$
polynomial will be non-zero if $p\geq p^{\prime}$ and zero otherwise Chihara
(2011). In other words, any order $p$ polynomial may be expressed as a linear
combination of $u_{p^{\prime}}$, where $p^{\prime}\leq p$ is required since
the order of interest is $p$. Thus
${\langle{u_{p^{\prime}}}\,|\,{v_{p}}\rangle}$ and
${\langle{v_{p^{\prime}}}\,|\,{u_{p}}\rangle}$ only have non-zero values in
the respective upper ($p\geq p^{\prime}$) and lower ($p\leq p^{\prime}$)
triangle of $\hat{L}$. Further, $\xi\,u_{p}$ is a linear combination of
$\\{u_{p-1},u_{p},u_{p+1}\\}$ due to Eq. (103d), so the matrix whose elements
are ${\langle{u_{p^{\prime}}}\,|\,{\xi|u_{p}}\rangle}$ will be tridiagonal.
Together, these points mean that the only way that Eq. (112) can hold, i.e.
the only way that the upper and lower triangular structure of
${\langle{u_{p^{\prime}}}\,|\,{v_{p}}\rangle}$ and
${\langle{v_{p^{\prime}}}\,|\,{u_{p}}\rangle}$ can be equal to the tridiagonal
structure of ${\langle{u_{p^{\prime}}}\,|\,{\xi|u_{p}}\rangle}$, is if
$\\{{\langle{u_{p^{\prime}}}\,|\,{v_{p}}\rangle}\\}_{p^{\prime}p}$ is upper-
triangular and tridiagonal. This means that $|{v_{p}}\rangle$ must be a linear
combination of $u_{p-1}$ and $u_{p}$,

$\displaystyle|{v_{p}}\rangle=\sigma_{p3}|{u_{p-1}}\rangle+\sigma_{p4}|{u_{p}}\rangle\;,$
(113)

where $\sigma_{p3}={\langle{u_{p-1}}\,|\,{v_{p}}\rangle}$ and
$\sigma_{p4}={\langle{u_{p}}\,|\,{v_{p}}\rangle}$ are complex valued constants
which are determined by $c_{pk}$ (i.e. $u_{p}$) and the confluent Heun
eigenvalues $\lambda_{pk}$.

Together, Eq. (103d), Eq. (110) and Eq. (113) allow
$\mathcal{L}_{\xi}|{u_{p}}\rangle$ to be written as

$\displaystyle\mathcal{L}_{\xi}|{u_{p}}\rangle\;=\;\sigma_{p5}\,|{u_{p-1}}\rangle+\sigma_{p6}\,|{u_{p}}\rangle+\sigma_{p7}\,|{u_{p+1}}\rangle\;,$
(114)

where

$\displaystyle\sigma_{p5}\;$
$\displaystyle=\;(\mu_{p}-\text{C}_{1})\sigma_{p0}+\sigma_{p3}\;,$ (115a)
$\displaystyle\sigma_{p6}\;$
$\displaystyle=\;(\mu_{p}-\text{C}_{1})\sigma_{p1}+\sigma_{p4}+(\text{C}_{0}+\text{C}_{1})\;,$
(115b) $\displaystyle\sigma_{p7}\;$
$\displaystyle=\;(\mu_{p}-\text{C}_{1})\sigma_{p2}\;.$ (115c)

The implication of Eqs. (114-115) is that, if the canonical polynomials
$u_{n}$ can be determined, then Teukolsky’s radial operator is exactly
tridiagonalizable.

##  IX Discussion

The primary results of this article are (i) the introduction of a radial
scalar product for QNMs (Sec. IV), (ii) the development of QNM specific
confluent Heun polynomials that may be defined for arbitrary frequency
parameter (Sec. VII.2), and (iii) the demonstration that, in principle, the
scalar product and proposed canonical confluent Heun polynomials enable the
exact tridiagonalization of Teukolsky’s radial equation for QNMs (Sec. VIII).

It is also of note that the scalar product may be evaluated by direct
integration with a complex valued radius, or by more standard analytic
extension using confluent hypergeometric functions (See Sec. VII.1). While
analytic extension was found to be the more practical route, the ability to
construct a complex radial coordinate, $\xi(z)$, for evaluation of the scalar
product also implies (via Eq. (30)) that the QNMs themselves are spatially
bounded on $\xi(z)$. Since $\xi(z)$ must be constructed to coincide with
Sturm-Liouville boundary conditions, a QNM on $\xi(z)$ would have physically
relevant values at both the event horizon and spatial infinity.

Regarding the canonical polynomials, a clear next step is to investigate their
computation. Like the Gram-Schmidt polynomials, which are unique once certain
processing choices are made, one might expect that there are multiple
constructions of the canonical polynomials. This a topic of ongoing research
that is near conclusion London and Gurevich (2023a).

Similarly, the tridiagonalization of Teukolsky’s radial equation directly
implies that its solutions (for fixed frequency parameter) are orthogonal.
This would be directly analogous to the orthogonality of the spheroidal
harmonics at fixed frequency (i.e. fixed “oblateness”) London (2023). It would
also be compatible with the conclusions of complementary studies of linear
perturbations of Kerr that have encountered a kind of QNM orthogonality Green
_et al._ (2023).