# Radiation fluxes of gravitational, electromagnetic, and scalar perturbations
in type-D black holes: an exact approach

Changkai Chen  0000-0002-4023-0682 Department of Physics, Key Laboratory of
Low Dimensional Quantum Structures and Quantum Control of Ministry of
Education, Synergetic Innovation Center for Quantum Effects and Applications,
Hunan Normal University, Changsha, 410081, Hunan, China  Jiliang Jing
0000-0002-2803-7900<EMAIL_ADDRESS>Department of Physics, Key Laboratory
of Low Dimensional Quantum Structures and Quantum Control of Ministry of
Education, Synergetic Innovation Center for Quantum Effects and Applications,
Hunan Normal University, Changsha, 410081, Hunan, China

###### Abstract

We present a novel method that solves Teukolsky equations with the source to
calculate radiation fluxes at infinity and event horizon for any perturbation
fields of type-D black holes. For the first time, we use the confluent Heun
function to obtain the exact solutions of ingoing and outgoing waves for the
Teukolsky equation. This benefits from our derivation of the asymptotic
analytic expression of the confluent Heun function at infinity. It is
interesting to note that these exact solutions are not subject to any
constraints, such as low-frequency and weak-field. To illustrate the
correctness, we apply these exact solutions to calculate the gravitational,
electromagnetic, and scalar radiations of the Schwarzschild black hole.
Numerical results show that the proposed exact solution appreciably improves
the computational accuracy and efficiency compared with the 23rd post-
Newtonian order expansion and the Mano-Suzuki-Takasugi method.

##  I Introduction

Black hole (BH) perturbation theory [1, 2] is a method used to study various
real relativistic objects, such as massive compact objects, jets, supernova
explosions, binary systems, etc. Originally developed as a metric perturbation
theory, for Schwarzschild BHs, Regge [3] and Zerilli [4] decoupled and
separated a single master equation for the metric perturbation into odd and
even parity parts, respectively. However, no such equation has been
established for rotating Kerr BHs. This led Bardeen and Press [5] to derive a
master equation for the curvature perturbation of a Schwarzschild BH without
the source ($T_{\ell m\omega}=0$) using the Newman-Penrose null-tetrad
formalism, where the tetrad components of the curvature tensor serve as
fundamental variables. Expanding to a Kerr BH with the source ($T_{\ell
m\omega}\neq 0$), Teukolsky [6] derived the curvature perturbation equation.
This resulting equation, known as the Teukolsky equation, represents a wave
equation for the null-tetrad component of the Weyl tensors $\psi_{0}$ and
$\psi_{4}$. The Teukolsky equation describes the dynamics of various fields of
different spins as perturbations (scalar, neutrino, electromagnetic, and
gravitational perturbations) to the Kerr metric. Therefore, the Teukolsky
equation can be used as a mathematical model of gravitational waves (GWs) to
construct GW waveform templates. In recent years there is an increased
interest in GW detection, and we anticipate gaining a deeper understanding of
BHs’ demographics and properties in the coming years [7, 8, 9, 10, 11]. Future
space-based GW detectors, such as the Laser Interferometer Space Antenna
(LISA) [12], TianQin [13, 14], and Taiji [15, 16], will be built with the
specific purpose of detecting GW signals from sources that radiate in the
millihertz ($mHz$) bandwidth, namely extreme mass ratio inspirals (EMRIs) [17,
18, 19]. BH perturbation theory has been utilized in modeling EMRI systems and
their associated GWs [20, 21]. Under the influence of gravitational radiation
reaction, the compact secondary body undergoes a slow inspiral into the
primary supermassive black hole, emitting GWs that propagate to infinity. Due
to potential interference from other sources concurrently emitting GWs in the
$mHz$ bandwidth [12], detection and parameter estimation of EMRI signals will
rely on matched filtering techniques. Accurate calculation of GW waveform
templates is thus of paramount importance, requiring precision on the order of
fractions of radians in GW phases [22, 23]. Upon precise modeling and
interpretation, the obtained data presents an opportunity for unparalleled
experiments examining general relativity and detecting hitherto unknown
astrophysical phenomena [24, 25, 26].

The Teukolsky equation without the source simplifies the perturbation problem
significantly because it does not need to construct Green’s function. This
homogenous equation has been widely adopted to study the quasinormal modes
(QNMs) of BHs and other relativistic celestial bodies, including massive dense
stars, jets, supernova explosions, and more. By exploiting infinite series of
special functions, Leaver derived the analytical expressions for the solutions
of the Regge-Wheeler (RW) and Teukolsky equations [27], which enabled precise
computation of the QNMs spectrum of black holes [28] and discussed their
excitation using Green’s function method [29]. But Green’s function
constructed by Leaver cannot be applied to calculate gravitational radiation.

It is known, but historically somewhat under-appreciated in the physics
literature, that both the RW and Teukolsky equations without the source are
examples of the confluent Heun equation [30, 31]. One possible reason for the
historical lack of attention given to the confluent Heun equation in physics
literature is the difficulty in numerically calculating its associated Heun
class functions in the past, leading to a prevalence of Leaver’s solution.
However, recent advances in numerical algorithms have enabled the computation
of Heun class functions using various mathematical software packages. Maple
software version 7 introduced numerical calculations of Heun class functions
in 2001. And Motygin provided MATLAB codes for general and confluent Heun
functions in 2015 [32] and 2018 [33], respectively. The 12.1 version of
Mathematica software released in 2020 also includes numerical calculation
capabilities for Heun class functions. It should be noted, however, that the
performance and accuracy of Heun class functions can differ across software
packages. Fiziev provided the analytical solutions [34, 35, 36, 37] to the
source-free perturbation equations (RW or Teukolsky equations) in terms of the
confluent Heun function. Furthermore, these solutions have been utilized to
calculate QNMs of the Teukolsky equation describing the Schwarzschild black
hole[38], including its continuous spectrum[39], as well as for the central
engine of Gamma-ray bursts (GRB) and cosmic jet of the Kerr black hole[40,
41]. Following Fiziev’s work, numerous researchers have employed Heun class
functions to study source-free perturbation equations in various spacetimes,
such as QNMs [42] of Dirac field in $2+1$ dimensional GW background [43], QNMs
and the reflection coefficient [44] of massless fields in Kerr-Newman–de
Sitter BH [45], and QNMs of the massive scalar field in Kerr-AdS5 BH [46].
Moreover, Cook et al. [47] converted the modes equations of Kerr BH into the
confluent Heun equation and then solved the radial equation using the
continued fraction method, and the angular equation using the spectral method,
resulting in improved accuracy for the QNMs.

Figure 1: The difference between our method and other methods in black hole
perturbation

In previous research, the construction of Green’s function for the Teukolsky
equation with the source was deemed analytically infeasible. As a result,
numerical integration methods were widely employed by most scholars initially
to investigate the radiation phenomenon associated with black hole
perturbation. These numerical methods involve simulating the propagation of
GWs in a background spacetime distorted by a rotating BH with an external
particle or other perturbation included. Such simulations furnish invaluable
insight into the behavior of gravitational waves and their interactions with
black holes[48, 49, 50, 51, 52]. Subsequently, using the post-Newtonian (PN)
expansion, Sasaki derived a part of Green’s function [53] for the RW equation.
This derivation involved only the ingoing wave solution $X_{\ell
m\omega}^{{\rm{in}}}$ and its asymptotic amplitudes at infinity. Utilizing the
Sasaki-Nakamura (SN) transformation, Sasaki transformed $X_{\ell
m\omega}^{{\rm{in}}}$ into the ingoing wave solution $R_{\ell
m\omega}^{{\rm{in}}}$ for the Teukolsky equation with the source in a
Schwarzschild BH [54] and a Kerr BH [55]. This approach is applied to
determining $R^{\rm in}_{\ell m\omega}$ of a Schwarzschild BH up to 5.5PN
order [56], as well as $R^{\rm in}_{\ell m\omega}$ of a Kerr BH up to 4PN
order [1, 57]. For the absorption and tail correction of a Schwarzschild BH,
it is necessary to construct another part of Green’s function containing
$R_{\ell m\omega}^{{\rm{up}}}$, which represents the homogeneous solution of
the Teukolsky equation of pure outgoing waves. Therefore, Poisson and Sasaki
[58] utilized the spherical Hankel function to construct only 1PN $X_{\ell
m\omega}^{{\rm{up}}}$, but this solution is controversial and does not satisfy
the conservation of the Wronskian for Kerr BHs 111It follows from the
conservation of the Wronskian that $X_{\ell m\omega}^{{\rm{up}}}(r\to
2M)\sim-\bar{A}_{\ell m\omega}^{{\rm{up}}}{e^{-i\omega{r^{*}}}}+A_{\ell
m\omega}^{{\rm{in}}}{e^{i\omega{r^{*}}}},A_{\ell m\omega}^{{\rm{in}}}\neq 0$.
But $X_{\ell m\omega}^{{\rm{up}}}$ constructed using spherical Hankel
functions has the amplitude $A_{\ell m\omega}^{{\rm{in}}}=0$, which fails to
satisfy this conservation relation..

The Mano-Suzuki-Takasugi (MST) method is a second analytical approach for
constructing Green’s function of the Teukolsky equation with the source.
Japanese researchers ‘revamped’ Leaver’s solutions and developed new series
solutions for the radial solutions [59, 60, 2]. The primary distinction
between the MST solution and Leaver’s solutions is that the former can obtain
analytical expressions for asymptotic amplitudes, which the five solutions
produced by Leaver do not satisfy. The series solutions employed by the MST
method, known as the MST expansions, are naturally adapted to carrying out
low-frequency expansions. Casals believed that the MST series converge
theoretically for any frequency value, albeit their convergence speed
diminishes as the frequency magnitude increases [61]. Initially, Mano et al.
[59, 60] only presented the MST method as the second-order post-Minkowskian
expansion result. Subsequently, Fujita numerically computed the renormalized
angular momentum of the MST method and achieved remarkably precise solutions
[62, 63]. Additionally, Fujita utilized the MST method for computing
gravitational radiation at arbitrary high PN order. For instance, Fujita
provided 22PN results in a Schwarzschild BH [64] and 11PN results in a Kerr
BH[65]. Other than that, Fujita presented 5.5PN GW polarizations and
associated factorized resumed waveforms[66]. The MST method can also be
employed to compute the self-force (SF) acting on point particles [20, 67, 68,
69, 70].

Apart from the MST method, there exist alternative approaches for computing
the gravitational wave energy flux. Fully relativistic GW fluxes from orbits
of non-spinning particles were initially computed for eccentric orbits around
a Schwarzschild BH [71] and circular equatorial orbits around a Kerr BH [72].
Glampedakis et al. [73] calculated Fluxes from eccentric orbits in the Kerr
spacetime. Fully generic GW fluxes from non-spinning particles were derived by
Drasco[74], and Hughes et al. computed adiabatic waveforms of EMRIs [75].
Circular orbits in a BH spacetime with the spin of the secondary taken into
account were investigated via GW fluxes in Refs. [76, 77, 78, 79, 80], while
quasi-circular adiabatic evolution of such orbits that includes spin was
presented in Refs. [81, 82, 83, 84]. The first-order SF for circular orbits in
the Schwarzschild spacetime was calculated in Ref. [85]. Moreover, Skoupy et
al. [86] computed the fluxes from spinning bodies on eccentric equatorial
orbits around a Kerr BH, and provided the linear in spin approximation for the
adiabatic evolution [87]. Referring to the frequency-domain method developed
in Refs. [88, 89] for generic orbits of spinning bodies around a Kerr black
hole, Hughes et al. [90] computed the asymptotic GW fluxes from a spinning
body moving on such orbits up to linear order in the secondary spin in a Kerr
BH.

In this work, we solve the Teukolsky equation with the source for the first
time using the confluent Heun function to obtain an exact method for
calculating gravitational, electromagnetic and scalar radiation fluxes for any
type-D BHs. Our method can get the complete Green’s functions including
$R_{\ell m\omega}^{{\rm{in,up}}}$. Different from the PN expansion and MST
method in the BH perturbation approach, our strategy involves expressing the
general solution of the homogeneous Teukolsky equation as a linear combination
of two linearly independent particular solutions (constructed in the form of
confluent Heun functions), and then using the boundary conditions to obtain
$R_{\ell m\omega}^{{\rm{in,up}}}$. The construction difference between our
method and other methods is shown in Fig. 1.

The remainder of the paper is arranged as follows. In Sec. II, the Teukolsky
equation is rewritten into a more general form, which can cover many curvature
perturbation equations of Type-D spacetimes. In Sec. III, we revamped Green’s
function method to obtain radiation fluxes for any perturbation fields. For
Teukolsky equations without the source, we proposed new analytical solution
constructed by confluent Heun functions in Sec. IV. We present the analytical
expression for the confluent Heun function at infinity for the first time in
Sec. V. To verify the correctness of our method, our method is applied to
calculate the fluxes of scalar, electromagnetic, and gravitational
perturbations for Schwarzschild BHs in Sec. VI. Sec. VII is devoted to
assessing the accuracy and effectiveness of the method proposed in this work.
To this end, we compare our results with those of post-Newtonian expansion and
MST methods to validate the applicability of our approach. The conclusions are
presented in Sec. VIII. In this paper, we use geometrized units: $c=G=1$.

##  II Teukolsky Equations of Type-D BHs

We now reformulate radial Teukolsky equations (RTE) into a more comprehensive
form that can encompass all previously identified Teukolsky equations of
Type-D BHs, which, names the general form of radial Teukolsky equation
(GFRTE), is expressed as

$\left[{\Delta_{n}^{-s+1}\frac{d}{{dr}}\Delta_{n}^{s+1}\frac{d}{{dr}}+V(r)}\right]{R_{\ell
m\omega}}={\Delta_{n}}{{}_{s}T_{\ell m\omega}},$ (1)

where

$\displaystyle\Delta_{n}=\sum\limits_{i=0}^{n}{{b_{i}}{r^{2-i}}}=\prod\limits_{i=1}^{n}{{{\left({r-{r_{i}}}\right)}}},$
(2) $\displaystyle V(r)=\sum\limits_{i=0}^{\infty}{{{\bf v}_{i}}{r^{i}}}.$ (3)

The Newman-Penrose formalism allows for the determination of the explicit form
(3) of $V(r)$, and it can be observed that the simplest expansion of $V(r)$
does not contain any terms with fractional powers of $r$. The common approach
is to approximate $\Delta_{n}$ by $\Delta_{2}$ and $\Delta_{4}$ 222The
Teukolsky equation of Kerr-Newman (anti-)de Sitter BHs [91, 92, 93] without
the source corresponds to $\Delta_{4}$-type., resulting in various types of
the GFRTEs with different forms in the potential term $V(r)$. This paper
specifically focuses on the Teukolsky equation of $\Delta_{2}$-type, that is
$\Delta_{2}=\left({r-{r_{-}}}\right)\left({r-{r_{+}}}\right)$. Here, $r_{-}$
is the inner horizon, and $r_{+}$ is the outer (event) horizon.

For different black holes, the values of $\Delta_{2}$ and $V(r)$ in the
Teukolsky equation will change. Now, we collect information on seven black
holes, and the corresponding $\Delta_{2}$ and $V(r)$ are organized as follows:

I. Schwarzschild BHs

The potential $V(r)$ of Schwarzschild BHs for all perturbation fields [5, 54]
is

${V_{{\rm{Sch}}}}={\omega^{2}}{r^{4}}+2is\omega{r^{2}}(r-3M)-{\Delta_{2}}{\lambda},$
(4)

where ${r_{-}}=0$, ${r_{+}}=r_{\rm H}=2M$ and ${\lambda}$ is given in Eq.
(112) of Ref. [2].

II. Reissner-Nordström BHs

The potential $V(r)$ of Reissner-Nordström BHs for the massive charged scalar
perturbation field [94] is

${V_{{\rm{RN}}}}={{(\omega
r^{2}-eQr)}^{2}}-\Delta_{2}\left({{\lambda}+{\mu^{2}}{r^{2}}}\right),$ (5)

with $r_{\pm}=M\pm\sqrt{{M^{2}}-{Q^{2}}}$. Here, $Q$ is the charge of the
black hole and $e$ is the elementary charge.

III. high-dimensional Schwarzschild BHs

The potential $V(r)$ of the background of a ($4+1$)-dimensional, non-rotating,
neutral BHs projected onto a 3-brane with the perturbation fields of
$s=0,\frac{1}{2},1$ [95], that is

$\begin{array}[]{*{20}{l}}{{V_{(4+1)}}={\omega^{2}}{r^{4}}-is\omega\hat{r}_{\rm{H}}^{2}r}\\\
{\quad\quad+{\Delta_{2}}\left({2i\omega
sr+s{\Delta^{\prime\prime}}-2s-\lambda}\right)},\end{array}$ (6)

with ${r_{\pm}}=\pm{\hat{r}}_{\rm H}$. The horizon radius ${\hat{r}}_{\rm H}$
is shown in Eq. (2.4) of Ref. [95].

IV. Kerr BHs

The potential $V(r)$ of Kerr BHs for all perturbation fields [6] is

$V_{\rm
K}={K^{2}}-isK\Delta_{2}^{\prime}+\Delta_{2}\left({2isK^{\prime}-\lambda}\right),$
(7)

where ${r_{\pm}}=M\pm\sqrt{{M^{2}}-{a^{2}}}$ and
$K=\left({{r^{2}}+{a^{2}}}\right)\omega-ma$.

V. Sixth-derivative correction to Kerr BHs

$V(r)$ of sixth-derivative correction to Kerr BHs for all perturbation fields
[96, 97] is

$\displaystyle{V_{{\rm{K6}}}}$
$\displaystyle={K^{2}}+\Delta_{2}\left({isK^{\prime}-\hat{\lambda}}\right)$
$\displaystyle-
is\left(2M\omega({{r^{2}}-{a^{2}}})-am{{\Delta^{\prime}_{2}}}\right)$ (8)

where $\hat{\lambda}={a^{2}}{\omega^{2}}-s+{B_{lm}}$ and ${B_{lm}}$ is the
angular separation constant in Refs. [96, 97].

VI. Kerr-Sen BHs

The potential $V(r)$ of Kerr-Sen BHs [98] for the massive charged scalar
perturbation field [99, 100, 101] is

$\begin{array}[]{l}V_{\rm KS}={(\omega(\Delta_{2}+2Mr)-eQr-am)^{2}}\\\
\quad\quad-{\Delta_{2}}\left({{\mu^{2}}(\Delta_{2}+2Mr)+{\lambda}}\right),\end{array}$
(9)

where ${r_{\pm}}=M-b\pm\sqrt{{(M-b)^{2}}-{a^{2}}}$ and $b=Q^{2}/2M$.

Table 1:  General Form of Teukolsky Equations  Potential Term $V(r)$ | $\\{{r_{\pm}}\\}$ | Perturbation Fields
---|---|---
${V_{{\rm{Sch}}}}={\omega^{2}}{r^{4}}+2is\omega{r^{2}}(r-3M)-{\Delta_{2}}{\lambda}$ | $\left\\{M\pm M\right\\}$ | all spins
${V_{{\rm{RN}}}}={{(\omega r^{2}-eQr)}^{2}}-\Delta_{2}\left({{\lambda}+{\mu^{2}}{r^{2}}}\right)$ | $\left\\{M\pm\sqrt{{M^{2}}-{Q^{2}}}\right\\}$ | $s=0$
$\begin{array}[]{*{20}{l}}{{V_{(4+1)}}={\omega^{2}}{r^{4}}-is\omega\hat{r}_{\rm{H}}^{2}r}\\\ {\quad\quad+{\Delta_{2}}\left({2i\omega sr+s{\Delta^{\prime\prime}}-2s-\lambda}\right)}\end{array}$ | $\\{{\pm{\hat{r}}_{\rm H}}\\}$ | $|s|=0,\frac{1}{2},1\quad$
$V_{\rm K}={K^{2}}-isK\Delta_{2}^{\prime}+\Delta_{2}\left({2isK^{\prime}-\lambda}\right)$ | $\left\\{{M\pm\sqrt{{M^{2}}-{a^{2}}}}\right\\}$ | all spins
$\begin{aligned} {V_{{\rm{K6}}}}&={K^{2}}+\Delta_{2}\left({isK^{\prime}-\hat{\lambda}}\right)\\\ &-is\left(2M\omega({{r^{2}}-{a^{2}}})-am{{\Delta_{2}^{\prime}}}\right)\end{aligned}$ | $\left\\{{M\pm\sqrt{{M^{2}}-{a^{2}}}}\right\\}$ | all spins
$\begin{array}[]{l}V_{\rm KS}={(\omega(\Delta_{2}+2Mr)-eQr-am)^{2}}\\\ \quad\quad-{\Delta_{2}}\left({{\mu^{2}}(\Delta_{2}+2Mr)+{\lambda}}\right)\end{array}$ | $\left\\{{M-b\pm\sqrt{{{(M-b)}^{2}}-{a^{2}}}}\right\\}$ | $s=0$
${V_{{\rm{KN1}}}}={K^{2}}-{\rm{i}}sK\Delta_{2}^{\prime}+\Delta_{2}\left({2{\rm{i}}sK^{\prime}-\lambda}\right)$ | $\left\\{{M\pm\sqrt{{M^{2}}-{a^{2}}-{Q^{2}}}}\right\\}$ | all spins
${V_{{\rm{KN2}}}}=K^{2}-\lambda{\Delta_{2}}$ | $\left\\{{M\pm\sqrt{{M^{2}}-{a^{2}}-{Q^{2}}}}\right\\}$ | all spins
$V_{\rm KN3}={\left({K-eQr}\right)^{2}}-\Delta_{2}\left({{\mu^{2}}\left({{r^{2}}+{a^{2}}}\right)+\lambda}\right)\quad$ | $\left\\{{M\pm\sqrt{{M^{2}}-{a^{2}}-{Q^{2}}}}\right\\}$ | $s=0$
$\begin{array}[]{l}{V_{{\rm{KN4}}}}={\left({K-eQr}\right)^{2}}-is\left({K-eQr}\right){\Delta^{\prime}_{2}}\\\ \quad\quad\quad+{\Delta_{2}}\left({2{\rm{i}}sK^{\prime}+\lambda}\right)\end{array}$ | $\left\\{{M\pm\sqrt{{M^{2}}-{a^{2}}-{Q^{2}}}}\right\\}$ | $|s|=1,2$

VII. Kerr-Newman BHs

For the curvature perturbation equation, due to the coupling between different
types of perturbation fields, it seems impossible to transform the general
perturbation of KN BHs into a single equation except for some limit cases or
the scalar perturbation field. We collected four examples with
${r_{\pm}}=M\pm\sqrt{{M^{2}}-{a^{2}}-{Q^{2}}}$ as follow:

VII-1. Using weakly charged approximation, Dudley and Finley [102, 103]
derived a Teukolsky-like equation (named as Dudley-Finley (DF) equation) for
all spin fields. And $V(r)$ of DF equation for all perturbation fields is

${V_{{\rm{KN1}}}}={K^{2}}-{\rm{i}}sK\Delta_{2}^{\prime}+\Delta_{2}\left({2{\rm{i}}sK^{\prime}-\lambda}\right),$
(10)

While the QNMs of DF equation yield exact solutions for scalar perturbations,
they are considered a conceptually questionable approximation for
gravitational and electromagnetic modes[104, 105].

VII-2. $V(r)$ of the eikonal limit $\ell\gg 1$ for all perturbation fields
[106] is

${V_{{\rm{KN2}}}}=K^{2}-\lambda{\Delta_{2}},$ (11)

VII-3. $V(r)$ of Kerr-Newman BH for charged massive scalar perturbation field
[107, 108] is

$V_{\rm
KN3}={\left({K-eQr}\right)^{2}}-\Delta_{2}\left({{\mu^{2}}\left({{r^{2}}+{a^{2}}}\right)+\lambda}\right),$
(12)

VII-4. $V(r)$ of Kerr-Newman BH for charged photons and gravitons perturbation
fields [109] is

$\begin{array}[]{l}{V_{{\rm{KN4}}}}={\left({K-eQr}\right)^{2}}-is\left({K-eQr}\right){\Delta^{\prime}_{2}}\\\
\quad\quad\quad+{\Delta_{2}}\left({2{\rm{i}}sK^{\prime}+\lambda}\right),\end{array}$
(13)

All the black holes given above satisfy Type-D black holes, that is, the four
Weyl tensors (${\Psi_{0}},{\Psi_{1}},{\Psi_{3}},{\Psi_{4}}$) are 0. The
potential terms of the above ten Teukolsky equations are summarized in Table
1.

##  III Fluxes and Green’s Function Method

Teukolsky equations with the source have found widespread application in the
study of Radiation fluxes and GW waveforms [2, 20, 21, 17, 18, 19, 110],
absorption [111] of gravitational waves and tail correction [61] of linear
field perturbations, as well as scalar self-force [112] and electromagnetic
self-force [113] acting on a charged particle in Kerr spacetime. In this
section, the Green’s function method is employed to solve the Teukolsky
equation with the source, in order to derive radiation fluxes. Using the Green
function method, it is necessary to construct two homogeneous solutions of Eq.
1, denoted as $R^{\rm in}_{\ell\omega}(r)$ and $R^{\rm up}_{\ell\omega}(r)$,
that satisfy the following boundary conditions:

$\displaystyle R_{\ell
m\omega}^{{\rm{in}}}\to\left\\{{\begin{array}[]{*{20}{l}}{B_{\ell
m\omega}^{{\rm{trans}}}{\Delta^{-s}}{{\rm{e}}^{-iP{r^{*}}}},}&{r\to{r_{+}}}\\\
{{B_{\ell
m\omega}^{{\rm{ref}}}}{r^{1-2s}}{{\rm{e}}^{i\omega{r^{*}}}}+{{B_{\ell
m\omega}^{{\rm{inc}}}}}r^{-1}{{\rm{e}}^{-i\omega{r^{*}}}},}&{r\to+\infty,}\end{array}}\right.$
(16) $\displaystyle R_{\ell
m\omega}^{{\rm{up}}}\to\left\\{{\begin{array}[]{*{20}{l}}{C_{\ell
m\omega}^{{\rm{up}}}{{\rm{e}}^{iP{r^{*}}}}+{C_{\ell
m\omega}^{{\rm{ref}}}}\Delta^{-s}{{\rm{e}}^{-iP{r^{*}}}},}&{r\to{r_{+}},}\\\
{{C_{\ell
m\omega}^{{\rm{trans}}}}{r^{1-2s}}{{\rm{e}}^{i\omega{r^{*}}}},}&{r\to+\infty,}\end{array}}\right.$
(19)

where $P=P(\omega,a,m,r_{\pm},Q,\mu)$ and $\Delta_{2}$ is abbreviated as
$\Delta$. The tortoise coordinate $r^{*}$ can be defined by

$\displaystyle{r^{*}}$
$\displaystyle=r+\frac{{{r_{+}}+{r_{-}}}}{{{r_{+}}-{r_{-}}}}{r_{+}}\ln\frac{{r-{r_{+}}}}{{{r_{+}}+{r_{-}}}}$
$\displaystyle-\frac{{{r_{+}}+{r_{-}}}}{{{r_{+}}-{r_{-}}}}{r_{-}}\ln\frac{{r-{r_{-}}}}{{2M}}.$
(20)

With the two homogeneous solutions $R_{\ell m\omega}^{{\rm{in,up}}}$, we can
easily obtain the solution of the inhomogeneous Teukolsky equation (1) with
purely ingoing behavior at the horizon and purely outgoing behavior at
infinity, which is described by

${R_{\ell m\omega}}(r)={}_{s}\tilde{Z}_{\ell m\omega}^{\infty}R_{\ell
m\omega}^{{\rm{up}}}(r)+{}_{s}\tilde{Z}_{\ell m\omega}^{{\rm{H}}}R_{\ell
m\omega}^{{\rm{in}}}(r),$ (21)

The asymptotic amplitudes $\tilde{Z}_{\ell m\omega}^{{\infty,\rm{H}},s}$ are
radial integrals over the source term:

${}_{s}\tilde{Z}_{\ell m\omega}^{{\infty,\rm{H}}}=\frac{1}{{W_{\rm
C}}}\int_{{r_{+}}}^{\infty}dr^{\prime}\frac{{{{}_{s}T_{\ell
m\omega}}(r^{\prime})R_{\ell
m\omega}^{{\rm{in,up}}}(r^{\prime})}}{{{\Delta^{2}}(r^{\prime})}},$ (22)

where ${{\rm{W}}_{C}}$ is the conserved Wronskian, that is

$\displaystyle{W_{\rm C}}$ $\displaystyle=R_{\ell
m\omega}^{{\rm{up}}}\frac{d}{{d{r^{*}}}}R_{\ell m\omega}^{{\rm{in}}}-R_{\ell
m\omega}^{{\rm{in}}}\frac{d}{{d{r^{*}}}}R_{\ell m\omega}^{{\rm{up}}}$
$\displaystyle=2i\omega B_{\ell m\omega}^{{\rm{inc}}}C_{\ell
m\omega}^{{\rm{trans}}}.$ (23)

The amplitudes ${}_{s}\tilde{Z}_{\ell m\omega}^{{\rm H},\infty}$ of
inhomogeneous equations (1) can be employed to precisely obtain energy and
angular momentum fluxes for scalar, electromagnetic, and gravitational
perturbations. In the case of circular orbits, ${}_{s}\tilde{Z}_{\ell
m\omega}^{{\rm H},\infty}$ in Eq. 22 takes the form

${}_{s}\tilde{Z}_{\ell m\omega}^{{\rm H},\infty}={}_{s}\tilde{Z}_{\ell
m\omega}^{{\rm H},\infty}\,\delta(\omega-m\,\Omega).$ (24)

The time-averaged luminosity (energy flux) at infinity is provided by [111]

${{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\biggl{|}}}_{s}{}\left\langle{\frac{{dE}}{{dt}}}\right\rangle_{\infty}=\sum\limits_{\ell=2}^{\infty}{\sum\limits_{m=-\ell}^{\ell}{{}_{s}{\beta_{\ell
m\omega}}}}\frac{{|{}_{s}\tilde{Z}_{\ell
m\omega}^{\infty}{|^{2}}}}{{4\pi{\omega^{2(|s|-1)}}}},$ (25)

where $\left<\cdots\right>$ represents the time average and $\omega=m\Omega$.

Similarly, the time-averaged luminosity at the horizon becomes[111, 114]

${{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\biggl{|}}}_{s}{}\left\langle{\frac{{dE}}{{dt}}}\right\rangle_{\rm{H}}=\sum\limits_{\ell=2}^{\infty}{\sum\limits_{m=-\ell}^{\ell}{{}_{s}{\alpha_{\ell
m\omega}}}}\frac{{|{}_{s}Z_{\ell
m\omega}^{{\rm{H}}}{|^{2}}}}{{4\pi{\omega^{2}}}},$ (26)

where $\alpha^{s}_{\ell m\omega}$ and $\beta^{s}_{\ell m\omega}$ are the
coefficients that contain the Starobinsky-Teukolsky constant [48, 67].

The angular momentum flux at infinity and event horizon for our case are given
by

${{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\biggl{|}}}_{s}{}{\left\langle{\frac{{dJ}}{{dt}}}\right\rangle_{{\rm{H,}}\infty}}=\frac{1}{\Omega}{{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\biggl{|}}}_{s}\left\langle{\frac{{dE}}{{dt}}}\right\rangle_{{\rm{H,}}\infty}}$
(27)

Moreover, the paper presents several functions with symmetry, which we
summarize as follows:

$\displaystyle{}_{s}Z_{\ell,-m\omega}^{\infty,{\rm{H}}}={(-1)^{\ell}}{}_{s}\bar{Z}_{\ell
m\omega}^{\infty,{\rm{H}}},$ (28a) $\displaystyle{{}_{s}}{S_{\ell
m}}(\theta)={(-1)^{(s+\ell)}}{{\kern 1.0pt}_{s}}{S_{\ell,-m}}(\pi-\theta),$
(28b) $\displaystyle{{\bar{R}}_{\ell,-m,-\omega}}={R_{\ell m\omega}},$ (28c)

where the bar denotes complex conjugation.

##  IV Analytical Solution of Teukolsky Equations without the source

Teukolsky equations without the source have been extensively utilized to study
various physical phenomena related to black holes, such as quasinormal modes
[28, 47, 104, 105, 106], Hawking radiation [95], near-superradiant scattering
[109], scalar clouds [100, 101, 107, 108], the central engine [40] of Gamma-
ray bursts, and cosmic jets [41]. The homogenous radial Teukolsky equations
(HRTEs) can be expressed as

$\left[{\Delta^{-s+1}\frac{d}{{dr}}\Delta^{s+1}\frac{d}{{dr}}+V(r)}\right]{R_{\ell
m\omega}}=0,$ (29)

###  IV.1 General Solution of HRTEs

Eq. 29 is a second-order ordinary differential equation (ODE) whose general
solution can be expressed as:

${R_{\ell m\omega}}(r)={C_{1}}R_{0}^{\beta}(r)+{C_{2}}R_{0}^{-\beta}(r),$ (30)

where $R_{0}^{\pm\beta}(r)$ are two linear independent particular solutions of
Eq. 29, $C_{1}$ and $C_{2}$ are constants that should be determined based on
different boundary conditions.

The HRTEs (29) is not the standard ODE that matches the known special
function. By utilizing the so-called S-homotopic transformation [31], the
HRTEs (29) can be transformed into an ODE that corresponds to a special
function known as the Heun class equation. Thus, two linear independent
particular solutions of Eq. 29 over the entire range
$r\in{\cal{R}}=[r_{H},\infty)$ can be obtained as

$R_{0}^{\pm\beta}(r)=S_{\rm{0}}^{\pm\beta}(x)\mathbb{H}^{\pm\beta}(x),$ (31)

where $x$ is defined as a new coordinate that is obtained by applying a Möbius
(isomorphic) transformation, which is

$x=-\frac{{r-r_{+}}}{{r_{+}}-{r_{-}}}.$ (32)

Introducing the unnormalized S-homotopic transformation

$S_{\rm{0}}^{\pm\beta}(x)={\left({-x}\right)^{\frac{1}{2}(\pm\beta-s)}}{\left({1-x}\right)^{\frac{1}{2}(\gamma-s)}}{{\rm{e}}^{\frac{1}{2}\alpha
x}},$ (33)

which can be regarded as the asymptotic behaviors of the function
$R_{0}^{\pm\beta}(r)$ at regular singularities, then substituting Eqs. 31, 32
and 33 into Eq. 29, we can obtain standard confluent Heun equation

$\displaystyle\mathbb{H}^{\prime\prime}-\frac{{\left({-{x^{2}}\alpha+\left({-2-\beta-\gamma+\alpha}\right)x+1+\beta}\right)}}{{x\left({x-1}\right)}}\mathbb{H}^{\prime}$
$\displaystyle-\left(\begin{array}[]{l}\left({\left({-2-\beta-\gamma}\right)\alpha-2\delta}\right)x\\\
+\left({\beta+1}\right)\alpha+\left({-\gamma-1}\right)\beta-\gamma-2\eta\end{array}\right)\frac{\mathbb{H}}{{2x\left({x-1}\right)}}=0.$
(36)

Two linear independent particular solutions
$\mathbb{H}_{\rm{0}}^{\pm\beta}(x)$ of sec. IV.1 can be expressed as

$\displaystyle\mathbb{H}_{\rm{0}}^{\beta}(x)={\rm{HeunC}}(\alpha,\beta,\gamma,\delta,\eta;x),$
(37a)
$\displaystyle\mathbb{H}_{\rm{0}}^{-\beta}(x)={(-x)^{-\beta}}{\rm{HeunC}}(\alpha,-\beta,\gamma,\delta,\eta;x),$
(37b)

where HeunC is the confluent Heun333In this paper, the HeunC function are
implemented using the corresponding functions and symbolic notations provided
by the computational software Maple. function [30, 31, 115], and
$\alpha,\beta,\gamma,\delta$ and $\eta$ are parameters that should be
determined for given black holes. In Table 2, we calculate the parameters
corresponding to the ten Teukolsky equations proposed in Sec. II .

Table 2:  parameters that should be determined for given black holes

$V(r)$ | $\alpha$ | $\beta$ | $\gamma$ | $\delta$ | $\eta$
---|---|---|---|---|---
${V_{{\rm{Sch}}}}$ | $2i\omega{r_{+}}$ | $-s-2i\omega{r_{+}}$ | $s$ | $-2is\omega{r_{+}}-2{\omega^{2}}r_{+}^{2}$ | $\begin{array}[]{l}2is\omega{r_{+}}+2{\omega^{2}}r_{+}^{2}\\\ -\frac{1}{2}{s^{2}}-s-\lambda\end{array}$
$V_{{\rm{RN}}}$ | $-2\sqrt{{\mu^{2}}-{\omega^{2}}}{r_{\rm{x}}}$ | $\frac{{2i{r_{+}}}}{{{r_{\rm{x}}}}}(\omega{r_{+}}-eQ)$ | $\frac{{2i{r_{-}}}}{{{r_{\rm{x}}}}}(eQ-\omega{r_{-}})$ | $\begin{array}[]{l}2r_{\rm{x}}^{2}{\omega^{2}}-2eQ\omega{r_{\rm{x}}}\\\ -{\mu^{2}}(r_{-}^{2}-r_{+}^{2}){r_{\rm{x}}}\end{array}$ | $\begin{array}[]{l}-\lambda-\frac{{{r_{+}}}}{{r_{\rm{x}}^{2}}}(2{Q^{2}}{e^{2}}{r_{-}}+\\\ 2Qe\omega r_{+}^{2}-2{\omega^{2}}r_{+}^{3}\\\ +{\mu^{2}}r_{-}^{2}{r_{+}}-2{\mu^{2}}{r_{-}}r_{+}^{2}\\\ +{\mu^{2}}r_{+}^{3}+4{\omega^{2}}{r_{-}}r_{+}^{2}\\\ -6Qe\omega{r_{-}}{r_{+}})\end{array}$
$V_{(4+1)}$ | $-2i\omega{r_{\rm{x}}}$ | $\begin{array}[]{l}-\frac{1}{{{r_{\rm{x}}}}}(r_{\rm{x}}^{2}{s^{2}}+\\\ 4is\omega r_{+}^{3}-4{\omega^{2}}r_{+}^{4}{)^{\frac{1}{2}}}\end{array}$ | $\begin{array}[]{l}\frac{1}{{{r_{\rm{x}}}}}(r_{\rm{x}}^{2}{s^{2}}+\\\ 4is\omega{r_{-}}r_{+}^{2}-4{\omega^{2}}r_{-}^{4}{)^{\frac{1}{2}}}\end{array}$ | $2\omega{r_{\rm{x}}}\left({\omega{r_{\rm{x}}}+is}\right)$ | $\begin{array}[]{l}-\frac{{{s^{2}}}}{2}+3s-\lambda+\\\ \frac{{\omega{r_{+}}}}{{r_{\rm{x}}^{2}}}(2isr_{-}^{2}-3is{r_{-}}{r_{+}}+\\\ 3isr_{+}^{2}-4\omega{r_{-}}r_{+}^{2}+2\omega r_{+}^{3})\end{array}$
$V_{\rm K}$ | $-2i\omega r_{\rm x}$ | $-s+\frac{{2i\omega(r_{+}^{2}+{a^{2}})-2iam}}{{{r_{\rm{x}}}}}$ | $s+\frac{{2i\omega(r_{-}^{2}+{a^{2}})-2iam}}{{{r_{\rm{x}}}}}$ | $2\omega r_{\rm x}({\omega r_{\rm{w}}+is})$ | $\begin{array}[]{l}2is\omega{r_{+}}-\frac{1}{2}{s^{2}}-s-\lambda\\\ -\frac{2}{{r_{\rm{x}}^{2}}}({\omega r_{\rm m}-am})({\omega r_{\rm o}-am})\end{array}$
$V_{\rm K6}$ | $-2i\omega r_{\rm x}$ | $\begin{array}[]{l}-\frac{1}{{{r_{\rm{x}}}}}\\{r_{\rm{x}}^{2}{s^{2}}-4r_{\rm{o}}^{2}{\omega^{2}}\\\ +8is[am(M-{r_{+}})\\\ -M\omega({a^{2}}-r_{+}^{2})]\\\ +16Mam\omega{r_{+}}\\\ -4{m^{2}}{a^{2}}{\\}^{\frac{1}{2}}}\end{array}$ | $\begin{array}[]{l}-\frac{1}{{{r_{\rm{x}}}}}\\{r_{\rm{x}}^{2}{s^{2}}\\\ -4{\left({{a^{2}}+r_{-}^{2}}\right)^{2}}{\omega^{2}}\\\ +8is[am(M-{r_{-}})\\\ -M\omega({a^{2}}-r_{-}^{2})]\\\ +16Mam\omega{r_{-}}\\\ -4{m^{2}}{a^{2}}{\\}^{\frac{1}{2}}}\end{array}$ | $2\omega{r_{\rm{x}}}\left({\omega r_{\rm{w}}+is}\right)$ | $\begin{array}[]{l}2is\omega{r_{+}}-s-\frac{1}{2}{s^{2}}-\hat{\lambda}\\\ +\frac{1}{{r_{\rm{x}}^{2}}}\\{-2{a^{4}}{\omega^{2}}-2(2iMs\omega\\\ +2{r_{-}}{r_{+}}{\omega^{2}}+{m^{2}}){a^{2}}\\\ +ma[4\omega Mr_{\rm{w}}+2is(2M-r_{\rm{w}})]\\\ +2\omega{r_{+}}[2iMs{r_{-}}\\\ -\omega r_{+}^{2}(2{r_{-}}-{r_{+}})]\\}\end{array}$
$V_{\rm KS}$ | $-2\sqrt{{\mu^{2}}-{\omega^{2}}}{r_{\rm{x}}}$ | $\begin{array}[]{l}-\frac{{2i}}{{{r_{\rm{x}}}}}(2M\omega{r_{+}}\\\ -ma-eQ{r_{+}})\end{array}$ | $\begin{array}[]{l}\frac{{2i}}{{{r_{\rm{x}}}}}(2M\omega{r_{-}}\\\ -ma-eQ{r_{-}})\end{array}$ | $\begin{array}[]{l}-2{r_{\rm{x}}}(eQ\omega\\\ +({\mu^{2}}-2{\omega^{2}})M)\end{array}$ | $\begin{array}[]{l}-2{r_{+}}({\mu^{2}}M+2{\omega^{2}}M+eQ\omega)\\\ -\lambda-2am\omega+\frac{1}{{r_{\rm{x}}^{2}}}\left\\{{[-2{a^{2}}{m^{2}}}\right.\\\ +2{r_{-}}{r_{+}}(4M\omega-eQ)\\\ -2ma({r_{-}}+{r_{+}})]eQ\\\ -8M\omega(-\omega r_{+}^{3}+2\omega{r_{-}}r_{+}^{2}\\\ +M\omega{r_{-}}{r_{+}}-{r_{+}}r_{-}^{2}\omega\\\ \left.{-\frac{{ma}}{2}r_{\rm{w}})}\right\\}\end{array}$
$V_{\rm{KN1}}$ | $-2i\omega r_{\rm x}$ | $-s+\frac{{2i\omega(r_{+}^{2}+{a^{2}})-2iam}}{{{r_{\rm{x}}}}}$ | $s+\frac{{2i\omega(r_{-}^{2}+{a^{2}})-2iam}}{{{r_{\rm{x}}}}}$ | $2\omega r_{\rm x}({\omega(r_{-}+r_{+})+is})$ | $\begin{array}[]{l}2is\omega{r_{+}}-\frac{1}{2}{s^{2}}-s-\lambda\\\ -\frac{2}{{r_{\rm{x}}^{2}}}\left[{\omega r_{\rm m}-am}\right]\times\\\ \left[{\omega r_{\rm o}-am}\right]\end{array}$
$V_{\rm{KN2}}$ | $-2i\omega r_{\rm x}$ | $\begin{array}[]{l}-\frac{2}{{{r_{\rm{x}}}}}[r_{\rm{x}}^{2}{s^{2}}-\\\ {(-ma+{r_{\rm{m}}}\omega)^{2}}{]^{\frac{1}{2}}}\end{array}$ | $\begin{array}[]{l}\frac{2}{{{r_{\rm{x}}}}}[r_{\rm{x}}^{2}{s^{2}}-\\\ {(-ma+{r_{\rm{m}}}\omega)^{2}}{]^{\frac{1}{2}}}\end{array}$ | $2{\omega^{2}}r_{-}^{2}-2{\omega^{2}}r_{+}^{2}$ | $\begin{array}[]{l}-\frac{1}{2}{s^{2}}-s-\lambda\\\ -\frac{2}{{r_{\rm{x}}^{2}}}\left[{\omega r_{\rm m}-am}\right]\times\\\ \left[{\omega r_{\rm o}-am}\right]\end{array}$
$V_{\rm{KN3}}$ | $-2\sqrt{{\mu^{2}}-{\omega^{2}}}{r_{\rm{x}}}$ | $\begin{array}[]{l}-\frac{{2{\kern 1.0pt}i}}{{{r_{\rm{x}}}}}({a^{2}}\omega+\omega r_{+}^{2}\\\ -ma-eQ{r_{+}})\end{array}$ | $\begin{array}[]{l}\frac{{2{\kern 1.0pt}i}}{{{r_{\rm{x}}}}}({a^{2}}\omega+\omega r_{-}^{2}\\\ -ma-eQ{r_{-}})\end{array}$ | $\begin{array}[]{l}-{r_{\rm{x}}}[2eQ\omega+\\\ ({r_{-}}+{r_{+}})({\mu^{2}}-2{\omega^{2}})]\end{array}$ | $\begin{array}[]{l}-{\mu^{2}}r_{\rm m}-\lambda\\\ -\frac{2}{{r_{\rm{x}}^{2}}}(\omega r_{\rm m}-ma-{r_{+}}eQ)\\\ \times(\omega r_{\rm o}-ma-eQ{r_{-}})\end{array}$
$V_{\rm{KN4}}$ | $-2i\omega r_{\rm x}$ | $\begin{array}[]{l}-\frac{1}{{r_{\rm{x}}^{2}}}[-4{a^{4}}{\omega^{2}}\\\ +8{a^{3}}m\omega\\\ +4{a^{2}}(2eQ\omega{r_{+}}\\\ -is\omega{r_{\rm{x}}}\\\ -2{\omega^{2}}r_{+}^{2}-{m^{2}})\\\ +4ma(2i\omega r_{+}^{2}\\\ -2eQ{r_{+}}+s{r_{\rm{x}}})\\\ +r_{\rm{x}}^{2}{s^{2}}-4is\omega r_{+}^{2}r_{\rm{x}}^{2}\\\ -4r_{+}^{2}{\left({eQ-\omega{r_{+}}}\right)^{2}}{]^{\frac{1}{2}}}\end{array}$ | $\begin{array}[]{l}\frac{1}{{{r_{\rm{x}}}}}[-4{a^{4}}{\omega^{2}}\\\ +8{a^{3}}m\omega\\\ +4{a^{2}}(2eQ\omega{r_{-}}\\\ +is\omega{r_{\rm{x}}}\\\ -2{\omega^{2}}r_{-}^{2}-{m^{2}})\\\ -4ma(is{r_{\rm{x}}}\\\ +2eQ{r_{-}}-2\omega r_{-}^{2})\\\ +r_{\rm{x}}^{2}{s^{2}}+4is\omega r_{-}^{2}{r_{\rm{x}}}\\\ -4r_{-}^{2}{\left({eQ-\omega{r_{-}}}\right)^{2}}{]^{\frac{1}{2}}}\end{array}$ | $\begin{array}[]{l}2\omega{r_{\rm{x}}}[({r_{-}}+{r_{+}})\omega\\\ +is-eQ]\end{array}$ | $\begin{array}[]{*{20}{l}}{2is\omega{r_{+}}-\frac{1}{2}{s^{2}}-s-\lambda}\\\ {-\frac{2}{{r_{\rm{x}}^{2}}}(\omega{r_{\rm{m}}}-ma-{r_{+}}eQ)}\\\ {\times(\omega{r_{\rm{o}}}-ma-eQ{r_{-}})}\end{array}$

  * $\clubsuit$

In Table 2, some symbols are defined as: ${r_{\rm{x}}}={r_{-}}-{r_{+}}$,
${r_{\rm{w}}}={r_{-}}+{r_{+}}$, $r_{\rm o}={a^{2}}+2{r_{-}}{r_{+}}-r_{+}^{2}$,
$r_{\rm m}=r_{+}^{2}+{a^{2}}$.

Using the $\rm HeunC$ function (37) and the S-homotopic transformation (33),
the general solution of Eq. 29 for $\Delta_{2}$-type can be expressed as:

$\displaystyle{R_{\ell m\omega}}$
$\displaystyle={C_{1}}S_{0}^{\beta}(x){\rm{HeunC}}(\alpha,\beta,\gamma,\delta,\eta;x)$
$\displaystyle+{C_{2}}S_{0}^{-\beta}(x){\rm{HeunC}}(\alpha,-\beta,\gamma,\delta,\eta;x).$
(38)

The general solution (IV.1) can be effectively applied to a specific physical
model by fixing the parameters ($\alpha,\beta,\gamma,\delta$, $\eta$) for a
given black hole. This allows us to tailor the solution to accurately describe
the characteristics and properties of the particular black hole under
consideration. $C_{1}$ and $C_{2}$ are the combination coefficients that are
explicitly solved by three significant boundary conditions (Dirichlet,
Neumann, and Robbin) in this paper. The application of these boundary
conditions requires the asymptotic behavior of the confluent Heun function at
the outer horizon and infinity, respectively. Expanding the confluent Heun
function in power series for the independent variable $x$ around the regular
singular point $x=0$ [30], yields the following asymptotic behavior at the
outer horizon:

$\mathop{\lim}\limits_{x\to
0}{\rm{HeunC}}(\alpha,\beta,\gamma,\delta,\eta;x)=1,\quad r\rightarrow r_{+},$
(39)

Expanding the confluent Heun function in a sector around the irregular
singular point at infinity [30], the asymptotic behavior at infinity can be
expressed as:

$\displaystyle\mathop{\lim}\limits_{\left|x\right|\to\infty}{\rm{HeunC}}(\alpha,\beta,\gamma,\delta,\eta;x)\to$
$\displaystyle
D_{\odot}^{\beta}\;{x^{-\frac{{\beta+\gamma+2}}{2}-\frac{\delta}{\alpha}}}+D_{\otimes}^{\beta}{{\rm{e}}^{-\alpha
x}}{x^{-\frac{{\beta+\gamma+2}}{2}+\frac{\delta}{\alpha}}},\ \ \ \
r\to\infty,$ (40)

where ${{D}_{\otimes}^{\beta}}$ and ${{D}_{\odot}^{\beta}}$ are undetermined
constants. Only when the constants ${{D}_{\otimes}^{\beta}}$ and
${{D}_{\odot}^{\beta}}$ are known, the value of the combination coefficient
$C_{1}$ and $C_{2}$ can be determined according to the boundary conditions.
Because the calculation process of the constants ${{D}_{\otimes}^{\beta}}$ and
${{D}_{\odot}^{\beta}}$ is very complicated and lengthy, we will introduce
their calculation process in detail in Sec. V.

While Fiziev proposed a general solutions similar to sec. IV.1 for solving
QNMs [38], he was not acquainted with the constants ${{D}_{\otimes}^{\beta}}$
and ${{D}_{\odot}^{\beta}}$. Therefore, he believes that providing explicit
analytical expressions for $C_{1}$ and $C_{2}$ is an unsolved difficult
problem in mathematics[34, 116]. In the past decade, Bezerra and Vieira have
dedicated their efforts to the exploration of Hawking radiation in scalar
fields[117, 118, 119, 120, 121]. They also constructed a general solution
similar to sec. IV.1 for the Klein-Gordon equation. Since the decay rate
$\Gamma_{+}$ does not involve the calculation of ${{D}_{\otimes}^{\beta}}$ and
${{D}_{\odot}^{\beta}}$, they can simply use asymptotic behavior (39) to
obtain the Hawking radiation spectrum.

Solving ${{D}_{\otimes}^{\beta}}$ and ${{D}_{\odot}^{\beta}}$ analytically is
one of the key advantages of our general solutions, which is different from
the general solutions of these literatures [34, 37, 38, 117, 118, 119, 120,
121]. It is also beneficial to extend the general solutions of homogeneous
equations to inhomogeneous equations (1). In contrast, the general solutions
presented by Fiziev and Vieira are specific cases within our broader
framework. sec. IV.1 covers the general solution of QNM, Hawking radiation and
other physical problems. Furthermore, their general solution is unable to
construct the Green’s function required by the inhomogeneous equation and the
outgoing wave solution $R_{\ell m\omega}^{{\rm{up}}}$. Consequently, we apply
the general solution to construct the solution of inhomogeneous equations and
$R_{\ell m\omega}^{{\rm{up}}}$, thereby broadening its application scope.

###  IV.2 Ingoing and Outgoing Wave Solutions

To get the general solution of the Teukolsky equation with the source under
given boundary conditions, we should first find the ingoing wave and outgoing
wave solutions $R_{\ell m\omega}^{\rm in,up}$ of the homogenous Teukolsky
equation, then obtain the general solution of the inhomogeneous Teukolsky
equation utilizing the Green’s function method [122]. The most common methods
for solving the homogenous solutions $R^{\rm in,up}_{\ell m\omega}$ are the PN
expansion [54] of the SN equation and the MST method[2, 59, 60]. The former
method involves utilizing the Chandrasekhar-Sasaki-Nakamura transformation
[123, 124, 125, 54] to convert the HRTEs (29) into the SN equation, and then
derives the PN expansion of $R^{\rm in}_{\ell m\omega}$. On the other hand, in
the MST method, a series solution of the hypergeometric function that
converges within a finite region is first constructed, and then a series
solution of the Coulomb wave function that converges at infinity is generated.
The two solutions are subsequently matched to yield a convergent solution that
extends from the horizon to infinity.

Now, we commence with the construction of the ingoing and outgoing wave
solutions from the solution (IV.1). Using the asymptotic properties (39) and
(IV.1) of the confluent Heun function and the boundary condition of the
ingoing wave, we can construct the ingoing wave solution. Noting the ingoing
wave solution at the horizon is purely ingoing wave, we have $C_{2}^{\rm
in}=0$. Thus, the ingoing wave solution $R_{\ell m\omega}^{{\rm{in}}}$ is
given by

$\displaystyle R_{\ell m\omega}^{{\rm{in}}}$
$\displaystyle=C_{1}^{{\rm{in}}}R_{\rm{0}}^{\beta}+C_{2}^{{\rm{in}}}R_{\rm{0}}^{-\beta}=C_{1}^{{\rm{in}}}R_{\rm{0}}^{\beta}$
$\displaystyle=C_{1}^{{\rm{in}}}S_{0}^{\beta}(x){\rm{HeunC}}(\alpha,\beta,\gamma,\delta,\eta;x).$
(41)

Similarly, we can construct the outgoing wave solution as

$\displaystyle R_{\ell m\omega}^{{\rm{up}}}$
$\displaystyle=C_{1}^{{\rm{up}}}R_{\rm{0}}^{\beta}+C_{2}^{{\rm{up}}}R_{\rm{0}}^{-\beta}$
$\displaystyle=C_{1}^{{\rm{up}}}S_{0}^{\beta}(x){\rm{HeunC}}(\alpha,\beta,\gamma,\delta,\eta;x)$
$\displaystyle+C_{2}^{{\rm{up}}}S_{0}^{-\beta}(x){\rm{HeunC}}(\alpha,-\beta,\gamma,\delta,\eta;x).$
(42)

According to the asymptotic properties (39) and (IV.1) of the confluent Heun
function, and noting that the outgoing wave solution at infinity is purely
outgoing wave, we have

$\displaystyle
C_{1}^{{\rm{up}}}={\left({-1}\right)^{\beta+1}}\Big{(}{\frac{{D_{\odot}^{-\beta}}}{{D_{\odot}^{\beta}}}}\Big{)}C_{2}^{{\rm{up}}}.$
(43)

Thus, the outgoing wave solution is described by

$\displaystyle R_{\ell m\omega}^{{\rm{up}}}$
$\displaystyle=C_{2}^{{\rm{up}}}\Big{[}{\left({-1}\right)^{\beta+1}}\Big{(}{\frac{{D_{\odot}^{-\beta}}}{{D_{\odot}^{\beta}}}}\Big{)}S_{0}^{\beta}(x){\rm{HeunC}}(\alpha,\beta,\gamma,\delta,\eta;x)$
$\displaystyle+S_{0}^{-\beta}(x){\rm{HeunC}}(\alpha,-\beta,\gamma,\delta,\eta;x)\Big{]}.$
(44)

Compared with PN expansion results, our ingoing wave solution (IV.2) is
considered a complete solution without the need for series expansion. And our
outgoing wave solution (IV.2) is undisputedly accurate, satisfying the
conservation of the Wronskian determinant for any BHs. Additionally, the
solutions (IV.2) and (IV.2) are not constrained by limitations stemming from
slow motion and weak-field approximations, rendering its findings superior to
those of the PN expansion results near the horizon. Compared with the MST
method, our method employs a special function to construct the solutions
$R_{\ell m\omega}^{{\rm{in,up}}}$ which does not involve computing two-sided
infinite series or renormalized angular momentum solutions $\nu$ ( MST method
faces solving the transcendental equation of renormalized angular momentum
$\nu$ ). Therefore, the solutions impose no limitations, while the MST method
is constrained by low-frequency approximations.

##  V Asymptotic Formula of HeunC Function at Infinity

Up to now, an analytic asymptotic expression for the confluent Heun function
$\mathbb{H}(x)$ at infinity $|x|$ has not been reported in the literature.
While most literature provides a linear combination of the two asymptotic
solutions of the confluent Heun function at infinity, as in sec. IV.1, the
coefficients ${{D}_{\otimes}^{\beta}}$ and ${{D}_{\odot}^{\beta}}$ remain
undetermined. The solution $Y(x)$ to the generalized spherical wave equation
(GSWE444Some researchers[126, 127] have named GSWE as confluent Heun
equations, which is not used in this paper. ) is related to $\mathbb{H}(x)$ by
$Y(x)={{\rm{e}}^{i\omega x}}\mathbb{H}(x)$. Additionally, $Y(x)$ can be
expressed as a series in terms of Coulomb wave functions ${F_{n+\nu}}(x)$ [27]
that converges for $x>0$. Therefore, an analytic asymptotic expression for
$\mathbb{H}(x)$ at infinity can also be constructed in terms of a series
solution in terms of ${F_{n+\nu}}(x)$. However, the series in terms of
${F_{n+\nu}}(x)$ does not converge at $x=0$, which prohibits normalization
with the asymptotic expression of $\mathbb{H}(x)$ at $x=0$. Instead, a
proportionality relation between the two can be established, but the
proportionality coefficient $\Xi$ is undetermined. To determine this
coefficient $\Xi$, we can represent $\mathbb{H}(x)$ as a series solution in
terms of hypergeometric functions ${}_{2}{F_{1}}(x)$ [115], which converge at
$x=0$. Therefore, $\mathbb{H}(x)$ can be expressed as a series in terms of
${}_{2}{F_{1}}(x)$ and as a proportionality coefficient multiplied by a series
solution in terms of ${F_{n+\nu}}(x)$. By expanding both series in the
interval $0<x<\infty$, the proportionality coefficient $\Xi$ can be
determined. This mathematical technique is similar to the approach used in
Refs. [59, 60] to determine the analytic asymptotic amplitudes of $R_{\ell
m\omega}^{{\rm{in}}}$ at infinity.

###  V.1 Expansion in Series of Hypergeometric Function

The series expansion of the hypergeometric function can be utilized to
represent the confluent Heun function555This paper assumes
${\rm{Im}}(\alpha)>0$. .

$\displaystyle\mathbb{H}(x)$
$\displaystyle={\rm{HeunC}}(\alpha,\beta,\gamma,\delta,\eta;x)$
$\displaystyle=\mathbb{F}\sum\limits_{n=-\infty}^{+\infty}{f{{{}_{n}^{\nu}}}\,{}_{2}{F_{1}}\left({a,{b},{c};x}\right)}$
(45)

with
$\mathbb{F}={\left(\sum\nolimits_{n=-\infty}^{\infty}{f_{n}^{\nu}}\right)^{-1}}$
is the normalized function at $x=0$. Here,
$a=n+\nu+1+\frac{{\beta+\gamma}}{2},{b}=-n-\nu+\frac{{\beta+\gamma}}{2}$ and
${c}=\beta+1$. Substituting the series solution, as given by sec. V.1, into
the confluent Heun equation (IV.1), we can derive the following three-term
recurrence relation for the expansion coefficients $f_{n}^{\nu}$.

${{\hat{\alpha}}_{n}}f_{n+1}^{\nu}+{{\hat{\beta}}_{n}}f_{n}^{\nu}+{{\hat{\gamma}}_{n}}f_{n-1}^{\nu}=0$
(46)

where

$\displaystyle{{\hat{\alpha}}_{n}}$
$\displaystyle=\frac{{\left({2n+2\nu+2-\beta+\gamma}\right)}}{{8\left({n+\nu+1}\right)\left({2n+2\nu+3}\right)}}$
$\displaystyle\times\left({\alpha
n+\alpha\nu+\alpha-\delta}\right)\left({2n+2\nu+2-\gamma-\beta}\right),$ (47a)
$\displaystyle{{\hat{\beta}}_{n}}$
$\displaystyle=\eta+\frac{\delta}{2}-\frac{{{\beta^{2}}}}{4}-\frac{{{\gamma^{2}}}}{4}+\left({n+\nu}\right)\left({n+\nu+1}\right)$
$\displaystyle+\frac{{\delta\left({\gamma+\beta}\right)\left({\beta-\gamma}\right)}}{{8\left({n+\nu}\right)\left({n+\nu+1}\right)}},$
(47b) $\displaystyle{{\hat{\gamma}}_{n}}$
$\displaystyle=-\frac{{\left({2n+2\nu+\beta-\gamma}\right)}}{{8\left({n+\nu}\right)\left({2n+2\nu-1}\right)}}$
$\displaystyle\times\left({\alpha
n+\alpha\nu+\delta}\right)\left({2n+2\nu+\gamma+\beta}\right).$ (47c)

The phase parameter $\nu$, also known as the renormalized angular momentum,
may be obtained by solving a characteristic equation expressed as the sum of
two infinite continued fractions [27, 128, 126, 127].

${{\hat{\beta}}_{0}}=\frac{{{{\hat{\alpha}}_{-1}}{{\hat{\gamma}}_{0}}}}{{{{\hat{\beta}}_{-1}}-}}\frac{{{{\hat{\alpha}}_{-2}}{{\hat{\gamma}}_{-1}}}}{{{{\hat{\beta}}_{-2}}-}}\frac{{{{\hat{\alpha}}_{-3}}{{\hat{\gamma}}_{-2}}}}{{{{\hat{\beta}}_{-3}}-}}+\cdots+\frac{{{{\hat{\alpha}}_{0}}{{\hat{\gamma}}_{1}}}}{{{{\hat{\beta}}_{1}}-}}\frac{{{{\hat{\alpha}}_{1}}{{\hat{\gamma}}_{2}}}}{{{{\hat{\beta}}_{2}}-}}\frac{{{{\hat{\alpha}}_{2}}{{\hat{\gamma}}_{3}}}}{{{{\hat{\beta}}_{3}}-}}\cdots.$
(48)

Solving the solution $\nu$ of Eq. 48 is a complex task, because Eq. 48 is a
transcendental equation. Nevertheless, two approaches have been developed to
determine $\nu$. The first method involves presenting a series expansion of
$\nu$, without having to solve the transcendental equation directly [59, 60,
61]. However, this approach requires enforcing low-frequency constraints. The
second method, originally introduced by Fujita and Tagoshi[62, 63], utilizes
the Steed algorithm for continued fractions to numerically solve the
transcendental equation (48) and obtain $\nu$, without requiring any
constraints. This paper selects the unconstrained second method over the first
due to the limitations of the series expansion approach.

The series representation of the hypergeometric function is given by [115]

${}_{2}{F_{1}}({{a}},{{b}},{{c}},\tilde{x})=\sum\limits_{j=0}^{\infty}{\frac{{{{(a)}_{j}}{{(b)}_{j}}}}{{{{(c)}_{j}}}}{{\frac{{\tilde{x}}}{{j!}}}^{j}}},\quad
0<\left|{\tilde{x}}\right|<1$ (49)

where $({{a}})_{n}$ denotes the Pochhammer symbol defined as
$({{a}})_{j}={{a}}({{a}}+1)({{a}}+2)\cdots({{a}}+j-1)$ with $({{a}})_{0}=1$.
And $\tilde{x}$ is new variable, defined as $\tilde{x}=1/(1-x)$.

Applying the linear transformation of hypergeometric functions ( Eq. (15.3.8)
in Ref. [129] ), we can derive a relation between ${}_{2}{F_{1}}(x)$ and
${}_{2}{F_{1}}(\tilde{x})$.

$\displaystyle{}_{2}{F_{1}}$
$\displaystyle\left({{{a}},{{b}},{{c}};x}\right)=$
$\displaystyle{{\tilde{x}}^{{a}}}\frac{{\Gamma({{c}})\Gamma({{b}}-{{a}})}}{{\Gamma({{b}})\Gamma({{c}}-{{a}})}}{{\kern
1.0pt}_{2}}{F_{1}}\left({{{a}},{{c}}-{{b}},{{a}}-{{b}}+1;\tilde{x}}\right)$
$\displaystyle+$
$\displaystyle{{\tilde{x}}^{{b}}}\frac{{\Gamma({{c}})\Gamma({{a}}-{{b}})}}{{\Gamma({{a}})\Gamma({{c}}-{{b}})}}{{\kern
1.0pt}_{2}}{F_{1}}\left({{{b}},{{c}}-{{a}},{{b}}-{{a}}+1;\tilde{x}}\right).$
(50)

The confluent Heun function can be expressed in an alternate form by utilizing
Sec. V.1 and V.1, as shown below:

$\mathbb{H}(x)={\mathbb{H}_{n,\nu}}\left({\tilde{x}}\right)+{\mathbb{H}_{-n,-\nu-1}}\left({\tilde{x}}\right),$
(51)

where

$\displaystyle\begin{array}[]{l}{\mathbb{H}_{n,\nu}}\left({\tilde{x}}\right)=\mathbb{F}{{\tilde{x}}^{-\nu+\frac{{\beta+\gamma}}{2}}}\sum\limits_{n=-\infty}^{\infty}{\frac{{\Gamma({{c}})\Gamma({{b}}-{{a}})}}{{\Gamma({{b}})\Gamma({{c}}-{{a}})}}}\\\
\quad\quad\times f_{n}^{\nu}{{\kern
1.0pt}_{2}}{F_{1}}\left({{{a}},{{c}}-{{b}};{{a}}-{{b}}+1;\tilde{x}}\right),\end{array}$
(54)
$\displaystyle\begin{array}[]{l}{\mathbb{H}_{-n,-\nu-1}}\left({\tilde{x}}\right)=\mathbb{F}{{\tilde{x}}^{\nu+1+\frac{{\beta+\gamma}}{2}}}\sum\limits_{n=-\infty}^{\infty}{\frac{{\Gamma({{c}})\Gamma({{a}}-{{b}})}}{{\Gamma({{a}})\Gamma({{c}}-{{b}})}}}\\\
\quad\quad\times f_{n}^{-\nu-1}{{\kern
1.0pt}_{2}}{F_{1}}\left({{{b}},{{c}}-{{a}};{{b}}-{{a}}+1;\tilde{x}}\right).\end{array}$
(57)

Accordingly, the recurrence relation (46) possesses a structure such that
${f_{n}^{-\nu-1}}$ satisfies an identical recurrence relation to that of
${f_{n}^{\nu}}$. Introducing a new coordinate variable
$\hat{z}=2i/(\alpha\tilde{x})$ and applying Eq. 49 into Eq. 54 , we can expand
${\mathbb{H}_{n,\nu}}(\tilde{x})$ as a series in terms of $z$.

${\mathbb{H}_{n,\nu}}=\mathbb{F}{\left({\frac{{2i\hat{z}}}{\alpha}}\right)^{\nu-\frac{{\beta+\gamma}}{2}}}\sum\limits_{k=-\infty}^{\infty}{\sum\limits_{n=k}^{\infty}{{C_{n,n-k}}{{\hat{z}}^{k}}}},$
(58)

where

$\displaystyle{C_{n,n-k}}$
$\displaystyle=\frac{{{{\left({-n-\nu+\frac{{\beta+\gamma}}{2}}\right)}_{n-k}}{{\left({-n-\nu+\frac{\delta}{\alpha}}\right)}_{n-k}}}}{{\Gamma\left({n+\nu+1+\frac{{\beta+\gamma}}{2}}\right)\Gamma\left({n+\nu+1+\frac{{\beta-\gamma}}{2}}\right)}}$
$\displaystyle\times\frac{{\Gamma\left({\beta+1}\right)\Gamma\left({2n+2\nu+1}\right)}}{{{{\left({-2n-2\nu}\right)}_{n-k}}\left({n-k}\right)!}}{\left({-\frac{{i\alpha}}{2}}\right)^{-k}}f_{n}^{\nu},$
(59)

###  V.2 Expansion in Series of Coulomb Wave Function

Assuming that ${\cal F}_{n,\nu}^{C}$ is a non-trivial solution of the
confluent Heun equation (IV.1), which can be written as a series expansion of
Coulomb wave functions. Here, $\hat{z}=\frac{i}{2}\alpha\left({x-1}\right)$.
And ${\cal F}_{n,\nu}^{C}$ is proportional to
$\mathbb{H}={\rm{HeunC}}(\alpha,\beta,\gamma,\delta,\eta;x)$, such that
$\mathbb{H}\sim{\cal F}_{n,\nu}^{C}$. The expression of ${\cal F}_{n,\nu}^{C}$
can written as

${\cal
F}_{n,\nu}^{C}(\hat{z})={\left({\frac{{\alpha}}{2}}\right)^{\tau}}{{\rm{e}}^{-\frac{{i\pi\tau}}{2}-\frac{\alpha}{2}+i\hat{z}}}{{\hat{z}}^{-\frac{{\gamma+\beta+2}}{2}}}f_{n,\nu}^{C}(\hat{z}),$
(60)

where $\tau=\frac{1}{4}\left({3\beta+\gamma+\frac{{2\delta}}{\alpha}}\right)$,
and $f_{n,\nu}^{C}$ can be defined as

$f_{n,\nu}^{C}=\mathbb{F}\sum\limits_{n=-\infty}^{\infty}{{{(-i)}^{n}}\frac{{{{(\nu+1+i\hat{\eta})}_{n}}}}{{{{(\nu+1-i\hat{\eta})}_{n}}}}f_{n}^{\nu}}{F_{n+\nu}}(\hat{\eta},\hat{z}),$
(61)

in which the Coulomb wave function ${F_{n+\nu}}(\hat{\eta},\hat{z})$ can be
defined as

${F_{n+\nu}}(\hat{\eta},\hat{z})={2^{n+\nu}}{{\hat{z}}^{n+\nu+1}}{e^{-i\hat{z}}}\frac{{\Gamma(\hat{a})}}{{\Gamma(\hat{c})}}\Phi\left({\hat{a},\hat{c};2i\hat{z}}\right),$
(62)

with $\hat{a}=n+\nu+1+\frac{{\delta}}{\alpha}$ , $\hat{c}=2n+2\nu+2,$ and
$i{\hat{\eta}}=-\frac{{\delta}}{\alpha}$, $\Phi$ is the regular confluent
hypergeometric function which can be represented by the following series
expansion [115]:

$\Phi\left({\hat{a},\hat{c};2i\hat{z}}\right)=\sum\limits_{j=0}^{\infty}{\frac{{{{(\hat{a})}_{j}}}}{{{{(\hat{c})}_{j}}}}{{\frac{{\left({2i\hat{z}}\right)}}{{j!}}}^{j}}}.$
(63)

Expanding $f_{n,\nu}^{C}$ to series of $\hat{z}$

$f_{n,\nu}^{C}=\mathbb{F}{e^{-i\hat{z}}}{2^{\nu}}{{\hat{z}}^{\nu+1}}\sum\limits_{k=-\infty}^{\infty}{\sum\limits_{n=-\infty}^{k}{{D_{n,n-k}}}}{{\hat{z}}^{k}},$
(64)

where

$\displaystyle{D_{n,k-n}}$
$\displaystyle=\frac{{\Gamma\left({n+\nu+1+\frac{\delta}{\alpha}}\right){{\left({\nu+1-\frac{\delta}{\alpha}}\right)}_{n}}}}{{\Gamma\left({2n+2\nu+2}\right){{\left({\nu+1+\frac{\delta}{\alpha}}\right)}_{n}}}}$
$\displaystyle\times\frac{{{{\left({n+\nu+1+\frac{\delta}{\alpha}}\right)}_{k-n}}}}{{{{(2n+2\nu+2)}_{k-n}}(k-n)!}}{(-1)^{n}}{(2i)^{k}}f_{n}^{\nu}.$
(65)

Because $f_{n,\nu}^{C}({\hat{z}})$ is convergent at infinity, then we discuss
its analytic asymptotic formula at infinity. There is the analytic property
[126, 127, 130] of the confluent hypergeometric function

$\displaystyle\Phi(\hat{a},\hat{c};2i\hat{z})$
$\displaystyle=\frac{{\Gamma(\hat{c})}}{{\Gamma(\hat{c}-\hat{a})}}{e^{i\epsilon\hat{a}\pi}}\Psi(\hat{a},\hat{c};2i\hat{z})$
(66)
$\displaystyle+\frac{{\Gamma(\hat{c})}}{{\Gamma(\hat{a})}}{e^{i\pi(\hat{a}-\hat{c})\epsilon}}{e^{2i\hat{z}}}\Psi(\hat{c}-\hat{a},\hat{c};-2i\hat{z}),$

with

${\epsilon={\rm{sgn}}\left({{\rm{Im}}(2i\hat{z})}\right)=\left\\{{\begin{array}[]{*{20}{l}}{1,\quad\quad{\rm{if}}\quad{\rm{Im}}(2i\hat{z})>0,}\\\
{-1,\quad{\rm{if}}\quad{\rm{Im}}(2i\hat{z})<0,}\end{array}}\right.}$ (67)

where $\Psi$ is the irregular confluent hypergeometric function. Then, Eq.
(61) can be expressed as

$\displaystyle f_{n,\nu}^{C}=f_{n,\nu}^{\odot}+f_{n,\nu}^{\otimes}$ (68a)
$\displaystyle\begin{array}[]{l}f_{n,\nu}^{\odot}={2^{\nu}}\mathbb{F}{e^{i\pi\left({\nu+1+\frac{\delta}{\alpha}}\right)}}{{\hat{z}}^{\nu+1}}{e^{-i\hat{z}}}\frac{{\Gamma\left({\nu+1+\frac{\delta}{\alpha}}\right)}}{{\Gamma\left({\nu+1-\frac{\delta}{\alpha}}\right)}}\\\
\quad\quad\times\sum\limits_{n=-\infty}^{\infty}{f_{n}^{\nu}{{(2i\hat{z})}^{n}}}\Psi\left({\hat{a},\hat{c},2i\hat{z}}\right),\end{array}$
(68d)
$\displaystyle\begin{array}[]{l}f_{n,\nu}^{\otimes}={2^{\nu}}\mathbb{F}{e^{-i\pi\left({\nu+1-\frac{\delta}{\alpha}}\right)}}{{\hat{z}}^{\nu+1}}{e^{i\hat{z}}}\sum\limits_{n=-\infty}^{\infty}{\frac{{{{\left({\nu+1-\frac{\delta}{\alpha}}\right)}_{n}}}}{{{{\left({\nu+1+\frac{\delta}{\alpha}}\right)}_{n}}}}}\\\
\quad\quad\times{(-1)^{n}}f_{n}^{\nu}{(-2i\hat{z})^{n}}\Psi(\hat{c}-\hat{a},\hat{c};-2i\hat{z}).\end{array}$
(68g)

By taking into account the asymptotic behavior [115] of
$\Psi\left({\hat{a},\hat{c};2i\hat{z}}\right)$ at large $|x|$:

$\mathop{\lim}\limits_{x\to\infty}\Psi\left({\hat{a},\hat{c};2i\hat{z}}\right)\rightarrow{(2i\hat{z})^{-\hat{a}}},$
(69)

we can obtain the asymptotic analytic expression of $f_{n,\nu}^{C}$ at
infinity:

$\displaystyle
f_{n,\nu}^{\odot}=A_{n,\nu}^{\odot}{x^{-\frac{\delta}{\alpha}}}{{\rm{e}}^{\frac{\alpha}{2}x}},\quad
f_{n,\nu}^{\otimes}=A_{n,\nu}^{\otimes}{x^{\frac{\delta}{\alpha}}}{{\rm{e}}^{-\frac{\alpha}{2}x}},$
(70a) $\displaystyle
A_{n,\nu}^{\odot}={{\rm{e}}^{-\frac{\alpha}{2}}}{\left({\frac{{i\alpha}}{2}}\right)^{-\frac{\delta}{\alpha}}}\tilde{A}_{n,\nu}^{\odot},{\kern
1.0pt}A_{n,\nu}^{\otimes}={{\rm{e}}^{\frac{\alpha}{2}}}{\left({\frac{{i\alpha}}{2}}\right)^{\frac{\delta}{\alpha}}}\tilde{A}_{n,\nu}^{\otimes}.$
(70b)

where

$\displaystyle\tilde{A}_{n,\nu}^{\odot}={2^{-1-\frac{\delta}{\alpha}}}{{\rm{e}}^{\frac{{i\pi}}{2}\left({\nu+1+\frac{\delta}{\alpha}}\right)}}\frac{{\Gamma\left({\nu+1+\frac{\delta}{\alpha}}\right)}}{{\Gamma\left({\nu+1-\frac{\delta}{\alpha}}\right)}},$
(71a)
$\displaystyle\begin{array}[]{l}\tilde{A}_{n,\nu}^{\otimes}={2^{-1+\frac{\delta}{\alpha}}}{{\rm{e}}^{-\frac{{i\pi}}{2}\left({\nu+1-\frac{\delta}{\alpha}}\right)}}{\left({\sum\limits_{n=-\infty}^{+\infty}{f_{n}^{\nu}}}\right)^{-1}}\\\
\quad\quad\times\left({\sum\limits_{n=-\infty}^{+\infty}{{{(-1)}^{n}}}\frac{{{{\left({\nu+1-\frac{\delta}{\alpha}}\right)}_{n}}}}{{{{\left({\nu+1+\frac{\delta}{\alpha}}\right)}_{n}}}}f_{n}^{\nu}}\right).\end{array}$
(71d)

Analogous to Eq. 68a, we can derive the asymptotic analytical expression of
$f_{-n,-\nu-1}^{C}$ as $|x|\rightarrow\infty$,

${f_{-n,-\nu-1}^{C}={{\rm{e}}^{i\pi\left({\nu+\frac{1}{2}}\right)}}f_{n,\nu}^{\otimes}}+\frac{{\sin\pi\left({\nu+\frac{\delta}{\alpha}}\right)}}{{\sin\pi\left({\nu-\frac{\delta}{\alpha}}\right)}}{{\rm{e}}^{-i\pi\left({\nu+\frac{1}{2}}\right)}}f_{n,\nu}^{\odot}.$
(72)

###  V.3 Proportionality Coefficient

Both solutions, ${\mathbb{H}_{n,\nu}}$ and ${{\cal F}_{n,\nu}^{C}}$, converge
within an extremely wide region of $\hat{z}$. As $k$ is an arbitrary integer,
we set $k=0$, and find that the series representations ( Eqs. 58, 60 and 64 )
of ${\mathbb{H}_{n,\nu}}$ and ${{\cal F}_{n,\nu}^{C}}$ are proportional to the
same single-valued function of $\hat{z}$. Therefore, the analytical properties
of ${\mathbb{H}_{n,\nu}}$ and ${{\cal F}_{n,\nu}^{C}}$ are identical,
indicating that these two solutions are equivalent up to a multiplicative
constant. The proportional coefficient between ${{\mathbb{H}_{n,\nu}}}$ and
${{{\cal F}_{n,\nu}^{C}}}$ can now be determined

$\Xi_{n,\nu}^{\beta}=\frac{{{\mathbb{H}_{n,\nu}}}}{{{\cal
F}_{n,\nu}^{C}}}={2^{-\nu}}{\left({\frac{\alpha}{2}}\right)^{-\hat{\tau}}}{{\rm{e}}^{\frac{{i\pi\hat{\tau}+\alpha}}{2}}}\frac{{\sum\limits_{j=0}^{\infty}{{C_{n,n-k}}}}}{{\sum\limits_{n=-\infty}^{k}{{D_{n,n-k}}}}},$
(73)

where $\hat{\tau}=\frac{{\beta-\gamma}}{4}+\nu+\frac{\delta}{{2\alpha}}.$

For the convenience of calculation, we set $k=0$, so Eq. (73) is simplified as

$\displaystyle\begin{array}[]{*{20}{l}}{\Xi_{n,\nu}^{\beta}=\frac{{{2^{-\nu}}{{\left({\frac{\alpha}{2}}\right)}^{-\hat{\tau}}}{{\rm{e}}^{\frac{{i\pi\hat{\tau}+\alpha}}{2}}}\Gamma\left({\beta+1}\right)\Gamma\left({2\nu+2}\right)}}{{\Gamma\left({\nu+1+\frac{\delta}{\alpha}}\right)\Gamma\left({\nu+1-\frac{{\beta+\gamma}}{2}}\right)\Gamma\left({\nu+1+\frac{{\gamma-\beta}}{2}}\right)}}}\\\
{\times{{\Big{(}{\sum\limits_{n=-\infty}^{0}{\frac{{{{(-1)}^{n}}{{\left({\nu+1-\frac{\delta}{\alpha}}\right)}_{n}}}}{{(-n)!{{\left({2\nu+2}\right)}_{n}}{{\left({\nu+1+\frac{\delta}{\alpha}}\right)}_{n}}}}f_{n}^{\nu}}}\Big{)}}^{-1}}}\end{array}$
(76)
$\displaystyle\begin{array}[]{l}{\times\Big{(}{\sum\limits_{n=0}^{\infty}{{{\left({-1}\right)}^{n}}\frac{{\Gamma\left({n+2\nu+1}\right)\Gamma\left({n+\nu+1+\frac{{\gamma-\beta}}{2}}\right)\Gamma\left({n+\nu+1-\frac{{\beta+\gamma}}{2}}\right)}}{{(n!)\Gamma\left({n+\nu+1+\frac{{\beta-\gamma}}{2}}\right)\Gamma\left({n+\nu+1+\frac{{\beta+\gamma}}{2}}\right)}}f_{n}^{\nu}}}\Big{)}}.\end{array}$
(78)

To obtain the proportional coefficients $\Xi_{-n,-\nu-1}^{\beta}$ for
${\mathbb{H}_{-n,-\nu-1}}$ and ${{\cal F}_{-n,-\nu-1}^{C}}$, we can simply
substitute $n\Rightarrow-n$ and $\nu\Rightarrow-\nu-1$ into Eq. 76.

###  V.4 Infinite Asymptotic Behavior of HeunC Function

Based on the previous sections, we are now able to derive the analytical
asymptotic expression for the confluence Heun function at infinity.
Specifically, we can rewrite sec. IV.1 as follows.

$\displaystyle\mathop{\lim}\limits_{|x|\to\infty}\mathbb{H}(x)$
$\displaystyle=\Xi_{n,\nu}^{\beta}{\cal
F}_{n,\nu}^{C}+\Xi_{-n,-\nu-1}^{\beta}{\cal F}_{-n,-\nu-1}^{C}$
$\displaystyle=D_{\odot}^{\beta}\;{x^{-\frac{{\beta+\gamma+2}}{2}-\frac{\delta}{\alpha}}}+D_{\otimes}^{\beta}{{\rm{e}}^{-\alpha
x}}{x^{-\frac{{\beta+\gamma+2}}{2}+\frac{\delta}{\alpha}}}.$ (79)

Substituting the results from sec. V.2 and V.3 into sec. V.4, the constants
$D_{\odot}^{\beta}$ and $D_{\otimes}^{\beta}$ in the asymptotic behavior at
infinity (IV.1) are given by

$\displaystyle\begin{array}[]{l}D_{\odot}^{\beta}=\Xi_{n,\nu}^{\beta}D_{\odot,n,\nu}^{\beta}\\\
\quad\quad+{{\rm{e}}^{-i\pi\left({\nu+\frac{1}{2}}\right)}}\frac{{\sin\pi\left({\nu+\frac{\delta}{\alpha}}\right)}}{{\sin\pi\left({\nu-\frac{\delta}{\alpha}}\right)}}\Xi_{-n,-\nu-1}^{\beta}D_{\odot,-n,-\nu-1}^{\beta},\end{array}$
(82)
$\displaystyle{D_{\otimes}^{\beta}=\Xi_{n,\nu}^{\beta}D_{\otimes,n,\nu}^{\beta}+{{\rm{e}}^{i\pi\left({\nu+\frac{1}{2}}\right)}}\Xi_{-n,-\nu-1}^{\beta}D_{\otimes,-n,-\nu-1}^{\beta},}$
(83)

with

$\displaystyle D_{\odot,n,\nu}^{\beta}$
$\displaystyle={\left({-1}\right)^{\frac{{\gamma+\beta+2}}{2}+\frac{\delta}{\alpha}}}{\big{(}{\frac{\alpha}{2}}\big{)}^{\tau}}{\big{(}-{\frac{{i\alpha}}{2}}\big{)}^{-\frac{{\gamma+\beta+2}}{2}-\frac{\delta}{\alpha}}}{{\rm{e}}^{-\frac{{i\pi\tau+\alpha}}{2}}}$
$\displaystyle\times{2^{-1-\frac{\delta}{\alpha}}}{{\rm{e}}^{\frac{{i\pi}}{2}\big{(}{\nu+1+\frac{\delta}{\alpha}}\big{)}}}\Xi_{n,\nu}^{\beta}\frac{{\Gamma\big{(}{\nu+1+\frac{\delta}{\alpha}}\big{)}}}{{\Gamma\big{(}{\nu+1-\frac{\delta}{\alpha}}\big{)}}},$
(84) $\displaystyle D_{\otimes,n,\nu}^{\beta}$
$\displaystyle={\big{(}{-1}\big{)}^{\frac{{\gamma+\beta+2}}{2}-\frac{\delta}{\alpha}}}{\left({\frac{\alpha}{2}}\right)^{\tau}}{\left(-{\frac{{i\alpha}}{2}}\right)^{-\frac{{\gamma+\beta+2}}{2}+\frac{\delta}{\alpha}}}$
$\displaystyle\times{{\rm{e}}^{-\frac{{i\pi\tau-\alpha}}{2}}}\Xi_{n,\nu}^{\beta}\frac{{2^{-1+\frac{\delta}{\alpha}}}{{\rm{e}}^{-\frac{{i\pi}}{2}\big{(}{\nu+1-\frac{\delta}{\alpha}}\big{)}}}}{\sum\limits_{n=-\infty}^{+\infty}{f_{n}^{\nu}}}$
$\displaystyle\times\Big{(}{\sum\limits_{n=-\infty}^{+\infty}{{{(-1)}^{n}}}\frac{{{{\big{(}{\nu+1-\frac{\delta}{\alpha}}\big{)}}_{n}}}}{{{{\big{(}{\nu+1+\frac{\delta}{\alpha}}\big{)}}_{n}}}}f_{n}^{\nu}}\Big{)}.$
(85)

##  VI Application to Schwarzschild BHs

For simplicity but without loss of generality, we shall consider the
gravitational, electromagnetic, and scalar fluxes of the Schwarzschild black
hole as an illustrative example, and we shall provide the complete solution
for the purely ingoing wave at the horizon and the purely outgoing wave at
infinity, along with their respective amplitudes. We consider the case when a
test particle of mass $\mu$ travels a circular orbit around a Schwarzschild BH
of mass $M\gg\mu$. To calculate the radiation fluxes, the radial Teukolsky
equation of a Schwarzschild BH with the point source was given by Eq. 29,
reduces to:

$\displaystyle\Delta R_{\ell m\omega}^{\prime\prime}+2(r-M)(s+1)R_{\ell
m\omega}^{\prime}$ (86)
$\displaystyle+\left[{r^{2}}\left({{\omega^{2}}{r^{2}}-4i\omega(r-3M)}\right)-\Delta{\lambda}\right]{R_{\ell
m\omega}}={}_{s}T_{\ell m\omega},$

where $\Delta=r(r-{r_{\rm H}})$, and ${r_{\rm H}}=2M$. The source term
${}_{s}T_{\ell m\omega}$ for the gravitational, electromagnetic, and scalar
perturbations are given in Ref. [6].

Using Green’s function method, it is necessary to construct two linear
independent solutions of Eq. 86, denoted as $R^{\rm in}_{\ell m\omega}(r)$ and
$R^{\rm up}_{\ell m\omega}(r)$, that satisfy the following boundary
conditions:

$\displaystyle R_{\ell
m\omega}^{{\rm{in}}}\to\left\\{{\begin{array}[]{*{20}{l}}{B_{\ell
m\omega}^{{\rm{trans}}}{\Delta^{-s}}{{\rm{e}}^{-i\omega{r^{*}}}},}&{r\to
r_{\rm H}},\\\ {{{B_{\ell
m\omega}^{{\rm{ref}}}}}{r^{1-2s}}{{\rm{e}}^{i\omega{r^{*}}}}+{{B_{\ell
m\omega}^{{\rm{inc}}}}}{r^{-1}}{{\rm{e}}^{-i\omega{r^{*}}}},}&{r\to+\infty,}\end{array}}\right.$
(89) $\displaystyle R_{\ell
m\omega}^{{\rm{up}}}\to\left\\{{\begin{array}[]{*{20}{l}}{C_{\ell
m\omega}^{{\rm{up}}}{{\rm{e}}^{i\omega{r^{*}}}}+{C_{\ell
m\omega}^{{\rm{ref}}}}{\Delta^{-s}}{{\rm{e}}^{-i\omega{r^{*}}}},}&{r\to{r_{\rm
H}},}\\\ {{{C_{\ell
m\omega}^{{\rm{trans}}}}}{r^{1-2s}}{{\rm{e}}^{i\omega{r^{*}}}},}&{r\to+\infty,}\end{array}}\right.$
(92)

where $r^{*}=r+r_{\rm H}\ln(r/r_{\rm H}-1)$ and $s=-2$.

In the case of a circular orbit, the specific energy $\tilde{E}$ and angular
momentum $\tilde{L}_{z}$ of the particle are given by

$\displaystyle\tilde{E}=(r_{0}-2M)/\sqrt{r_{0}(r_{0}-3M)},$ (93)
$\displaystyle\tilde{L}_{z}=\sqrt{Mr_{0}}/\sqrt{1-3M/r_{0}},$ (94)

where $r_{0}$ is the orbital radius. The angular frequency is given by
$\Omega=(M/r_{0}^{3})^{1/2}$ , so the orbital frequency is $\omega=m\Omega$
and the orbital frequency in Mino time is
$\Upsilon_{t}=({\frac{r_{0}^{5}}{r_{0}-3M}})^{1/2}$.

###  VI.1 Normalized solutions

The solutions $R_{\ell m\omega}^{{\rm{in,up}}}$ proposed in Sec. IV.2 are the
analytical solutions satisfying the boundary conditions (89) and (92) of the
Schwarzschild BH. To facilitate comparisons with other methods, we need
normalize $R_{\ell m\omega}^{{\rm{in}}}$ at the event horizon and $R_{\ell
m\omega}^{{\rm{up}}}$ at infinity, obtaining normalized asymptotic amplitudes
of homogeneous equations (29). Using the normalized condition $B_{\ell
m\omega}^{{\rm{trans}}}=1$ at event horizon, the ingoing wave solution (IV.2)
is normalized to obtain

$\tilde{R}_{\ell m\omega}^{{\rm{in}}}=S_{\rm
H}^{\beta}\left(x\right){\rm{HeunC}}\left({\alpha,\beta,\gamma,\delta,\eta;x}\right),$
(95)

where the normalized S-homotopic transformation at horizon is

$S_{\rm
H}^{\beta}\left(x\right)=r^{-2s}_{\rm{H}}{{\rm{e}}^{\frac{1}{2}\alpha\left({x-1}\right)}}{{{\left({-x}\right)}^{\frac{1}{2}\left({\beta-s}\right)}}{{\left({1-x}\right)}^{\frac{1}{2}\left({\gamma-s}\right)}}}.$
(96)

By utilizing the asymptotic behavior (89) of the solution $\tilde{R}_{\ell
m\omega}^{{\rm{in}}}$ as $r\rightarrow\infty$, we derive analytic expressions
of the asymptotic amplitudes $B_{\ell m\omega}^{\rm inc}$ and $B_{\ell
m\omega}^{\rm ref}$ for $\tilde{R}_{\ell m\omega}^{{\rm{in}}}$.

$\displaystyle B_{\ell m\omega}^{{\rm{inc}}}$
$\displaystyle=-r_{\rm{H}}^{5}D_{\odot}^{\beta},$ (97a) $\displaystyle B_{\ell
m\omega}^{{\rm{ref}}}$
$\displaystyle={\left({-1}\right)^{2i\omega{r_{\rm{H}}}-1}}{r_{\rm{H}}}{{\rm{e}}^{-2i\omega{r_{\rm{H}}}}}D_{\otimes}^{\beta}.$
(97b)

Meanwhile, using the normalized condition $C_{\ell m\omega}^{{\rm{trans}}}=1$
at infinity, the outgoing wave solution is normalized to obtain

$\displaystyle\tilde{R}_{\ell m\omega}^{{\rm{up}}}=$
$\displaystyle{\left({-1}\right)^{\beta+1}}\frac{{D_{\odot}^{-\beta}}}{{D_{\odot}^{\beta}}}\tilde{D}S_{\infty}^{\beta}(x){\rm{HeunC}}(\alpha,\beta,\gamma,\delta,\eta;x)$
$\displaystyle+\tilde{D}S_{\infty}^{-\beta}(x){\rm{HeunC}}(\alpha,-\beta,\gamma,\delta,\eta;x),$
(98)

where the normalized S-homotopic transformation at infinity is

$S_{\infty}^{\beta}(x)=-r_{\rm{H}}^{3}{{\rm{e}}^{\frac{1}{2}\alpha\left({x+1}\right)}}{\left({-x}\right)^{\frac{1}{2}(\beta-s)}}{\left({1-x}\right)^{\frac{1}{2}(\gamma-s)}},$
(99)

and

$\tilde{D}={\left({D_{\otimes}^{-\beta}-\frac{{D_{\odot}^{-\beta}}}{{D_{\odot}^{\beta}}}D_{\otimes}^{\beta}}\right)^{-1}}.$

By utilizing the asymptotic behavior (92) of the solution $\tilde{R}_{\ell
m\omega}^{{\rm{up}}}$ as $r\rightarrow r_{\rm{H}}$, we derive analytic
expressions of the asymptotic amplitudes $C_{\ell m\omega}^{{\rm{ref}}}$ and
$C_{\ell m\omega}^{{\rm{up}}}$ for $\tilde{R}_{\ell m\omega}^{{\rm{up}}}$.

$\displaystyle C_{\ell m\omega}^{{\rm{up}}}$
$\displaystyle=-{r}^{3}_{\rm{H}}\tilde{D},$ (100a) $\displaystyle C_{\ell
m\omega}^{{\rm{ref}}}$
$\displaystyle=\frac{{{{\left({-1}\right)}^{-2i\omega{r_{\rm{H}}}}}}}{{{r_{\rm{H}}}}}\frac{{D_{\odot}^{-\beta}}}{{\;D_{\odot}^{\beta}}}{{\rm{e}}^{2i\omega{r_{\rm{H}}}}}\tilde{D}.$
(100b)

For the ingoing wave solution $\tilde{R}_{\ell m\omega}^{{\rm{in}}}$ and
outgoing wave solution $\tilde{R}_{\ell m\omega}^{{\rm{up}}}$, the parameters
($\alpha,\beta,\gamma,\delta$, $\eta$) of these two solutions for the
Schwarzschild BH can be seen in Table 2.

###  VI.2 Energy Fluxes of Schwarzschild BHs

####  VI.2.1 Gravitational field

According to BH perturbation theory, gravitational waves are described by
perturbations in the Weyl scalars $\psi_{0}$ and $\psi_{4}$. For vacuum
solutions corresponding to pure gravitational waves, both $\psi_{0}$ and
$\psi_{4}$ contain the same information about the propagation of gravitational
waves within the Kerr spacetime. However, both for near null infinity and near
the horizon, $\psi_{4}$ exhibits dominance over $\psi_{0}$ for gravitational
waves propagating along the positive radial ($+r$) direction, whereas
$\psi_{0}$ prevails over $\psi_{4}$ for waves propagating along the negative
radial ($-r$) direction. Specifically, $\psi_{0}$ corresponds to gravitational
waves with spin-weight $s=+2$, while $\psi_{4}$ corresponds to gravitational
waves with spin-weight $s=-2$. In studies focusing on gravitational waves
generated by compact binary coalescence that propagate towards future null
infinity, it is natural to compute $\psi_{4}$. Thus, Teukolsky equations have
been extensively applied for accurate evaluations of $\psi_{4}$ in the
pertinent scientific literature.

In order to calculate the radiative energy and angular momentum flues, the
expression of the asymptotic amplitude ${}_{s}\tilde{Z}_{\ell m\omega}^{{\rm
H},\infty}$ in inhomogeneous solution (21) needs to be solved. The expressions
of ${}_{s}\tilde{Z}_{\ell m\omega}^{{\rm H},\infty}$, ${{}_{s}{\beta_{\ell
m\omega}}}$, and ${{}_{s}{\alpha_{\ell m\omega}}}$ can be found in Appendix A.
Finally, theses coefficients and functions are substituted into Eqs. 25 and
26, and the expressions of energy fluxes are obtained. Therefore, we obtain
energy fluxes of the gravitational perturbation field of the spin-weight
$s=-2$.

$\displaystyle{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\biggl{|}}}_{-2}\left\langle{\frac{{dE}}{{dt}}}\right\rangle_{\infty}$
$\displaystyle=\sum\limits_{\ell=2}^{\infty}{\sum\limits_{m=-\ell}^{\ell}{\frac{{16\pi}}{{{\omega^{2}}}}\left|{\frac{1}{{{W_{\rm{C}}}{\Upsilon_{t}}}}\left[{\left({{A_{nn0}}+{A_{\bar{m}n0}}+{A_{\bar{m}\bar{m}0}}}\right)R_{\ell
m\omega}^{{\rm{in}}}}\right.}\right.}}$
$\displaystyle{\left.{-\left({{A_{\bar{m}n1}}+{A_{\bar{m}\bar{m}1}}}\right){{\left({R_{\ell
m\omega}^{{\rm{in}}}}\right)}^{\prime}}\left.{+{A_{\bar{m}\bar{m}2}}{{\left({R_{\ell
m\omega}^{{\rm{in}}}}\right)}^{\prime\prime}}}\right]}\right|^{2}},$ (101a)
$\displaystyle{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\biggl{|}}}_{-2}\left\langle{\frac{{dE}}{{dt}}}\right\rangle_{\rm{H}}$
$\displaystyle=\sum\limits_{\ell=2}^{\infty}{\sum\limits_{m=-\ell}^{\ell}{\frac{{{\rm{4096}}\pi{{(2M{r_{\rm{H}}})}^{5}}({\omega^{2}}+4{{\tilde{\epsilon}}^{2}})({\omega^{2}}+16{{\tilde{\epsilon}}^{2}}){\omega^{2}}}}{{{\lambda^{2}}{{(\lambda+2)}^{2}}+144{\omega^{2}}{M^{2}}.}}}}$
$\displaystyle\times\left|{\frac{1}{{{W_{\rm{C}}}{\Upsilon_{t}}}}\left[{\left({{A_{nn0}}+{A_{\bar{m}n0}}+{A_{\bar{m}\bar{m}0}}}\right)R_{\ell
m\omega}^{{\rm{up}}}}\right.}\right.$
$\displaystyle\left.{-\left({{A_{\bar{m}n1}}+{A_{\bar{m}\bar{m}1}}}\right){{\left({R_{\ell
m\omega}^{{\rm{up}}}}\right)}^{\prime}}\left.{+{A_{\bar{m}\bar{m}2}}{{\left({R_{\ell
m\omega}^{{\rm{up}}}}\right)}^{\prime\prime}}}\right]}\right|^{2}.$ (101b)

where $\tilde{\epsilon}={(4{r_{\rm{H}}})^{-1}}$ and the explicit form of
$A_{nn0}$ and other terms can be found and derived in Ref. [2, 6].

Similarly, substituting ${}_{s}\tilde{Z}_{\ell m\omega}^{{\rm H},\infty}$,
${{}_{s}{\beta_{\ell m\omega}}}$, and ${{}_{s}{\alpha_{\ell m\omega}}}$ into
the radiative energy fluxes (25) and (26), we obtain energy fluxes of the
gravitational perturbation field of the spin-weight $s=+2$.

$\displaystyle{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\biggl{|}}}_{+2}\left\langle{\frac{{dE}}{{dt}}}\right\rangle_{\infty}$
$\displaystyle=\sum\limits_{\ell=2}^{\infty}{\sum\limits_{m=-\ell}^{\ell}{\frac{{{\rm{256}}\pi{\omega^{6}}}}{{{{(\lambda+4)}^{2}}{{(\lambda+6)}^{2}}+144{M^{2}}{\omega^{2}}}}}}$
$\displaystyle\times\left|{\frac{1}{{{W_{\rm{C}}}{\Upsilon_{t}}}}\left[{\left({{A_{ll0}}+{A_{lm0}}+{A_{mm0}}}\right)R_{\ell
m\omega}^{{\rm{in}}}}\right.}\right.$
$\displaystyle{\left.{-\left({{A_{lm1}}+{A_{mm1}}}\right){{\left({R_{\ell
m\omega}^{{\rm{in}}}}\right)}^{\prime}}\left.{+{A_{mm2}}{{\left({R_{\ell
m\omega}^{{\rm{in}}}}\right)}^{\prime\prime}}}\right]}\right|^{2}},$ (102a)
$\displaystyle{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\biggl{|}}}_{+2}\left\langle{\frac{{dE}}{{dt}}}\right\rangle_{\rm{H}}$
$\displaystyle=\sum\limits_{\ell=2}^{\infty}{\sum\limits_{m=-\ell}^{\ell}{\frac{\pi}{{8r_{\rm{H}}^{3}\left({{\omega^{2}}+4{\epsilon^{2}}}\right)}}}}$
$\displaystyle\times\left|{\frac{1}{{{W_{\rm{C}}}{\Upsilon_{t}}}}\left[{\left({{A_{ll0}}+{A_{lm0}}+{A_{mm0}}}\right)R_{\ell
m\omega}^{{\rm{up}}}}\right.}\right.$
$\displaystyle{\left.{-\left({{A_{lm1}}+{A_{mm1}}}\right){{\left({R_{\ell
m\omega}^{{\rm{up}}}}\right)}^{\prime}}\left.{+{A_{mm2}}{{\left({R_{\ell
m\omega}^{{\rm{up}}}}\right)}^{\prime\prime}}}\right]}\right|^{2}}.$ (102b)

####  VI.2.2 Electromagnetic field

Electromagnetic waves are described in terms of perturbations in the
electromagnetic scalars $\phi_{0}$ and $\phi_{2}$. The electromagnetic scalars
$\phi_{0}$ and $\phi_{2}$, similar to Weyl scalars $\psi_{0}$ and $\psi_{4}$,
both contain the same information of wave propagation in Kerr spacetime.
However, when the electromagnetic wave propagates along the positive radial
($+r$) direction, $\phi_{2}$ dominates ; when the electromagnetic wave
propagates along the negative radial ($-r$) direction, $\phi_{0}$ dominates.
Here, $\phi_{0}$ corresponds to the electromagnetic field with $s=+1$, and
$\phi_{2}$ corresponds to the electromagnetic perturbation field with $s=-1$.

Similar to energy fluxes of GW radiation, ${}_{s}\tilde{Z}_{\ell
m\omega}^{{\rm H},\infty}$, ${{}_{s}{\beta_{\ell m\omega}}}$, and
${{}_{s}{\alpha_{\ell m\omega}}}$ are substituted into the radiative energy
fluxes (25) and (26), we obtain energy fluxes of the electromagnetic
perturbation field of the spin-weight $s=-1$.

${{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\biggl{|}}}_{-1}\left\langle{\frac{{dE}}{{dt}}}\right\rangle_{\infty}=\sum\limits_{\ell=2}^{\infty}{\sum\limits_{m=-\ell}^{\ell}{\frac{1}{{2\pi}}{{\left|{\frac{\mathfrak{A}_{-}}{{{W_{\rm{C}}}\Delta}}R_{\ell
m\omega}^{{\rm{in}}}-\frac{\mathfrak{B}_{-}}{{{W_{\rm{C}}}\Delta}}{{\left({R_{\ell
m\omega}^{{\rm{in}}}}\right)}^{\prime}}}\right|}^{2}}}},$ (103a)
$\displaystyle{{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\biggl{|}}}_{-1}\left\langle{\frac{{dE}}{{dt}}}\right\rangle_{\rm{H}}$
$\displaystyle=\sum\limits_{\ell=2}^{\infty}{\sum\limits_{m=-\ell}^{\ell}{\frac{{16{r_{\rm{H}}}{M^{3}}{\omega^{2}}\left({4{r_{\rm{H}}}^{2}{\omega^{2}}+1}\right)}}{{\pi{\lambda^{2}}}}}}$
$\displaystyle\times{\left|{\frac{\mathfrak{A}_{-}}{{{W_{\rm{C}}}\Delta}}R_{\ell
m\omega}^{{\rm{up}}}-\frac{\mathfrak{B}_{-}}{{{W_{\rm{C}}}\Delta}}{{\left({R_{\ell
m\omega}^{{\rm{up}}}}\right)}^{\prime}}}\right|^{2}}.$ (103b)

where the explicit form of ${\mathfrak{A}_{\pm}}$ and ${\mathfrak{B}_{\pm}}$
can be found in Eq. (48) of Ref. [113]

Similarly, substituting ${}_{s}\tilde{Z}_{\ell m\omega}^{{\rm H},\infty}$,
${{}_{s}{\beta_{\ell m\omega}}}$, and ${{}_{s}{\alpha_{\ell m\omega}}}$ into
the radiative energy fluxes (25) and (26), we obtain energy fluxes of the
electromagnetic perturbation field of the spin-weight $s=+1$.

${{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\biggl{|}}}_{+1}\left\langle{\frac{{dE}}{{dt}}}\right\rangle_{\infty}={\sum\limits_{\ell=2}^{\infty}{\sum\limits_{m=-\ell}^{\ell}{\frac{{2{\omega^{4}}}}{{\pi{{(\lambda+2)}^{2}}}}\left|{\frac{\mathfrak{A}_{+}}{{{W_{\rm{C}}}}}R_{\ell
m\omega}^{{\rm{in}}}-\frac{\mathfrak{B}_{+}}{{{W_{\rm{C}}}\Delta}}{{\left({\Delta
R_{\ell m\omega}^{{\rm{in}}}}\right)}^{\prime}}}\right|}}^{2}},$ (104a)
${{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\biggl{|}}}_{+1}\left\langle{\frac{{dE}}{{dt}}}\right\rangle_{\rm{H}}=\sum\limits_{\ell=2}^{\infty}{\sum\limits_{m=-\ell}^{\ell}{\frac{1}{{16\pi
r_{\rm H}}}{{\left|{\frac{\mathfrak{A}_{+}}{{{W_{\rm{C}}}}}R_{\ell
m\omega}^{{\rm{up}}}-\frac{\mathfrak{B}_{+}}{{{W_{\rm{C}}}\Delta}}{{\left({\Delta
R_{\ell m\omega}^{{\rm{up}}}}\right)}^{\prime}}}\right|}^{2}}}}.$ (104b)

####  VI.2.3 Scalar field

Substituting ${}_{s}\tilde{Z}_{\ell m\omega}^{{\rm H},\infty}$,
${{}_{s}{\beta_{\ell m\omega}}}$, and ${{}_{s}{\alpha_{\ell m\omega}}}$ into
the radiative energy fluxes (25) and (26), we obtain energy fluxes of the
scalar perturbation field of the spin-weight $s=0$.

${{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\biggl{|}}}_{0}\left\langle{\frac{{dE}}{{dt}}}\right\rangle_{\infty}=\sum\limits_{\ell=2}^{\infty}{\sum\limits_{m=-\ell}^{\ell}{4\pi{\omega^{2}}{{\left|{\frac{{r_{0}^{2}}}{{{W_{\rm{C}}}{\Upsilon_{t}}}}{{\kern
1.0pt}_{0}}{Y_{\ell m\omega}}R_{\ell m\omega}^{{\rm{in}}}}\right|}^{2}}}},$
(105a)
${{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}\biggl{|}}}_{0}\left\langle{\frac{{dE}}{{dt}}}\right\rangle_{\rm{H}}=\sum\limits_{\ell=2}^{\infty}{\sum\limits_{m=-\ell}^{\ell}{8\pi
r_{\rm
H}{\omega^{2}}{{\left|{\frac{{r_{0}^{2}}}{{{W_{\rm{C}}}{\Upsilon_{t}}}}{{\kern
1.0pt}_{0}}{Y_{\ell m\omega}}R_{\ell m\omega}^{{\rm{up}}}}\right|}^{2}}}}.$
(105b)

where ${}_{s}Y_{\ell m\omega}$ is the spin-weighted spherical harmonics.

The significant application of our method is to calculate radiation fluxes in
which both $R_{\ell m\omega}^{{\rm{in}}}$ and $R_{\ell m\omega}^{{\rm{up}}}$
can be obtained using Eq. (IV.1) directly with the boundary conditions (89)
and (92). Our method differs from traditional results obtained through post-
Newtonian expansion and post-Minkowskian expansion, as it provides complete
results without requiring series expansion. Theoretically, our method should
be more efficient and accurate than the expansion method, and it has a broader
scope of application as it is not limited by physical constraints such as low
frequency, slow-motion or weak-field limits.

##  VII Comparisons with other methods

Figure 2: The 3D figures of logarithmic luminosity at infinity with ${\bf
N}=50$, over the domain $[5M,15M]$.  Figure 3: The 3D figures of logarithmic
luminosity at horizon with ${\bf N}=50$, over the domain $[5M,15M]$.  Figure
4: The 2D figures of logarithmic luminosity at infinity with ${\bf N}=50$ at
$r_{\rm ISCO}$.  Figure 5: The 2D figures of logarithmic luminosity at horizon
with ${\bf N}=50$ at $r_{\rm ISCO}$.  Figure 6: Logarithmic relative errors of
energy fluxes ${}_{s}\left<{dE_{22}/dt}\right>_{\infty,{\rm H}}$ between four
methods and the exact solution with the different floating-point numbers $\bf
N$ at $r_{\rm ISCO}$.  Table 3: The comparison of four methods for energy
fluxes ${}_{s}\left<{dE_{22}/dt}\right>_{\infty,{\rm H}}$ (angular momentum
fluxes ${}_{s}\left<{dJ_{22}/dt}\right>_{\infty,{\rm H}}$) with the different
floating-point numbers $\bf N$ at $r_{\rm ISCO}$. In our tables,
${}_{s}\left<{dE_{22}/dt}\right>_{\infty,{\rm H}}$ is abbreviated as
${\dot{E}}_{\infty,{\rm H}}$.

$s$ | Flux | Method | $\bf N$ | 20 | 30 | 40 | 50 | 60 | 70 | 80 | 90 | 100 | Mean
---|---|---|---|---|---|---|---|---|---|---|---|---|---
-2 | ${\bf{RE}}({\dot{E}}_{\infty})$ | HeunC | $\bf RE$ | 2.006E-14* | 2.168E-24 | 4.097E-34 | 4.845E-46 | 4.589E-55 | 1.882E-64 | 4.296E-74 | 6.043E-84 | 1.422E-94 | -
| ROC | - | 9.966 | 9.724 | 11.927 | 9.024 | 9.387 | 9.642 | 9.852 | 10.628 | 10.019
MST-RTE | $\bf RE$ | 1.892E-11 | 7.584E-17 | 8.281E-22 | 4.790E-27 | 3.400E-32 | 6.227E-37 | 9.808E-43 | 4.476E-47 | 1.422E-52 | -
| ROC | - | 5.397 | 4.962 | 5.238 | 5.149 | 4.737 | 5.803 | 4.341 | 5.498 | 5.141
MST-RWE | $\bf RE$ | 8.894E-12 | 4.043E-16 | 4.448E-21 | 5.755E-26 | 5.062E-30 | 9.023E-36 | 5.219E-41 | 8.896E-46 | 9.094E-51 | -
| ROC | - | 4.342 | 5.178 | 4.787 | 5.009 | 4.969 | 5.036 | 5.221 | 4.777 | 4.915
23PN | $\bf RE$ | 5.417E-07 | 5.417E-07 | 5.417E-07 | 5.417E-07 | 5.417E-07 | 5.417E-07 | 5.417E-07 | 5.417E-07 | 5.417E-07 | -
| ROC | - | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000
${\bf{RE}}({{{\dot{E}}_{\rm H}}})$ | HeunC | $\bf RE$ | 6.543E-11 | 1.300E-20 | 7.178E-31 | 4.771E-42 | 8.066E-52 | 5.255E-61 | 2.221E-70 | 1.545E-82 | 4.779E-92 | -
| ROC | - | 9.702 | 10.258 | 11.177 | 9.772 | 9.186 | 9.374 | 12.158 | 9.510 | 10.142
MST-RTE | $\bf RE$ | 1.047E-10 | 1.440E-15 | 1.051E-20 | 7.309E-26 | 4.858E-31 | 3.720E-36 | 6.107E-41 | 3.496E-46 | 1.820E-51 | -
| ROC | - | 4.862 | 5.137 | 5.158 | 5.177 | 5.116 | 4.785 | 5.242 | 5.284 | 5.095
MST-RWE | $\bf RE$ | 5.510E-11 | 1.042E-15 | 1.877E-21 | 1.417E-25 | 3.886E-32 | 7.783E-36 | 5.495E-41 | 1.024E-45 | 6.537E-51 | -
|  | ROC | - | 4.723 | 4.959 | 4.888 | 5.056 | 3.749 | 6.238 | 4.768 | 4.990 | 4.874
+2 | ${\bf{RE}}({\dot{E}}_{\infty})$ | HeunC | $\bf RE$ | 1.989E-14 | 2.385E-23 | 4.101E-34 | 2.988E-43 | 4.784E-55 | 1.884E-64 | 4.358E-74 | 6.028E-84 | 5.396E-94 | -
| ROC | - | 8.921 | 10.765 | 9.138 | 11.796 | 9.405 | 9.636 | 9.859 | 10.048 | 9.946
MST-RTE | $\bf RE$ | 8.864E-12 | 3.031E-17 | 2.847E-22 | 1.452E-27 | 6.169E-32 | 2.483E-38 | 1.006E-41 | 9.288E-48 | 7.466E-54 | -
| ROC | - | 5.466 | 5.027 | 5.292 | 4.372 | 6.395 | 3.392 | 6.035 | 6.095 | 5.259
MST-RWE | $\bf RE$ | 2.195E-11 | 4.043E-16 | 4.448E-21 | 5.755E-26 | 5.062E-31 | 9.023E-35 | 5.219E-41 | 8.896E-46 | 9.094E-51 | -
|  | ROC | - | 4.735 | 4.959 | 4.888 | 5.056 | 3.749 | 6.238 | 4.768 | 4.990 | 4.923
${\bf{RE}}({{{\dot{E}}_{\rm H}}})$ | HeunC | $\bf RE$ | 1.071E-10 | 2.911E-19 | 1.506E-30 | 2.264E-39 | 1.143E-50 | 2.823E-60 | 3.886E-70 | 2.468E-80 | 5.450E-91 | -
| ROC | - | 8.566 | 11.286 | 8.823 | 11.297 | 9.607 | 9.861 | 10.197 | 10.656 | 10.037
MST-RTE | $\bf RE$ | 8.198E-11 | 1.137E-15 | 4.538E-21 | 8.956E-26 | 1.322E-30 | 1.855E-35 | 3.141E-41 | 6.077E-46 | 1.820E-51 | -
| ROC | - | 4.858 | 5.399 | 4.705 | 4.831 | 4.853 | 5.771 | 4.713 | 5.524 | 5.082
MST-RWE | $\bf RE$ | 1.096E-10 | 1.042E-15 | 1.877E-21 | 1.417E-25 | 3.886E-31 | 7.783E-36 | 5.495E-41 | 1.024E-45 | 6.537E-51 | -
|  | ROC | - | 5.022 | 5.744 | 4.122 | 5.562 | 4.698 | 5.151 | 4.730 | 5.195 | 5.028
$-1$ | ${\bf{RE}}({\dot{E}}_{\infty})$ | HeunC | $\bf RE$ | 2.533E-15 | 2.733E-25 | 1.084E-34 | 2.272E-45 | 1.694E-54 | 3.944E-66 | 6.398E-76 | 2.077E-85 | 1.775E-95 | -
| ROC | - | 9.967 | 9.402 | 10.679 | 9.127 | 11.633 | 9.790 | 9.489 | 10.068 | 10.019
MST-RTE | $\bf RE$ | 7.240E-13 | 2.935E-17 | 1.191E-22 | 3.002E-27 | 1.181E-32 | 6.267E-37 | 4.501E-42 | 8.646E-48 | 3.701E-53 | -
| ROC | - | 4.392 | 5.392 | 4.599 | 5.405 | 4.275 | 5.144 | 5.716 | 5.368 | 5.036
${\bf{RE}}({{{\dot{E}}_{\rm H}}})$ | HeunC | $\bf RE$ | 4.540E-11 | 1.371E-20 | 1.127E-30 | 2.188E-41 | 4.866E-52 | 2.097E-61 | 1.684E-70 | 4.477E-80 | 6.692E-90 | -
| ROC | - | 9.520 | 10.085 | 10.712 | 10.653 | 9.366 | 9.095 | 9.575 | 9.825 | 9.854
MST-RTE | $\bf RE$ | 7.824E-11 | 3.102E-16 | 5.299E-21 | 8.606E-26 | 2.776E-31 | 4.241E-36 | 6.368E-41 | 6.556E-46 | 2.807E-51 | -
|  | ROC | - | 5.402 | 4.767 | 4.789 | 5.491 | 4.816 | 4.823 | 4.987 | 5.368 | 5.056
$+1$ | ${\bf{RE}}({\dot{E}}_{\infty})$ | HeunC | $\bf RE$ | 2.531E-15 | 2.696E-25 | 8.520E-35 | 1.049E-43 | 1.696E-54 | 3.947E-66 | 7.068E-76 | 2.085E-85 | 2.495E-95 | -
| ROC | - | 9.973 | 9.500 | 8.910 | 10.791 | 11.633 | 9.747 | 9.530 | 9.922 | 10.001
MST-RTE | $\bf RE$ | 2.822E-12 | 2.979E-17 | 2.804E-22 | 1.436E-26 | 4.333E-33 | 9.535E-38 | 7.576E-44 | 9.074E-48 | 1.175E-52 | -
| ROC | - | 4.976 | 5.026 | 4.291 | 6.520 | 4.657 | 6.100 | 3.922 | 4.888 | 5.048
${\bf{RE}}({{{\dot{E}}_{\rm H}}})$ | HeunC | $\bf RE$ | 5.493E-11 | 1.285E-20 | 1.166E-30 | 2.105E-42 | 1.125E-50 | 4.089E-60 | 1.812E-70 | 4.657E-80 | 6.835E-90 | -
| ROC | - | 9.631 | 10.042 | 11.743 | 8.272 | 9.440 | 10.353 | 9.590 | 9.833 | 9.863
MST-RTE | $\bf RE$ | 4.628E-11 | 6.097E-17 | 4.777E-21 | 9.611E-26 | 3.575E-31 | 8.646E-36 | 8.646E-41 | 7.472E-46 | 1.068E-50 | -
|  | ROC | - | 5.880 | 4.106 | 4.696 | 5.429 | 4.616 | 5.000 | 5.063 | 4.845 | 4.955
$0$ | ${\bf{RE}}({\dot{E}}_{\infty})$ | HeunC | $\bf RE$ | 1.059E-15 | 1.296E-22 | 2.308E-29 | 7.493E-35 | 2.467E-40 | 8.204E-46 | 2.749E-51 | 3.757E-56 | 5.153E-61 | -
| ROC | - | 6.912 | 6.749 | 5.489 | 5.482 | 5.478 | 5.475 | 4.864 | 4.863 | 5.664
MST-RTE | $\bf RE$ | 1.708E-11 | 1.192E-18 | 4.098E-24 | 7.184E-30 | 9.991E-32 | 2.120E-37 | 9.476E-43 | 2.534E-49 | 2.300E-52 | -
| ROC | - | 7.156 | 5.464 | 5.756 | 1.857 | 5.673 | 5.350 | 6.573 | 3.042 | 5.109
${\bf{RE}}({{{\dot{E}}_{\rm H}}})$ | HeunC | $\bf RE$ | 4.250E-11 | 5.811E-20 | 3.617E-26 | 1.175E-31 | 3.867E-37 | 1.286E-42 | 4.310E-48 | 5.890E-53 | 8.077E-58 | -
| ROC | - | 8.864 | 6.206 | 5.488 | 5.483 | 5.478 | 5.475 | 4.864 | 4.863 | 5.840
MST-RTE | $\bf RE$ | 5.423E-11 | 4.007E-16 | 5.853E-22 | 1.812E-27 | 1.646E-31 | 1.479E-36 | 2.901E-41 | 5.421E-46 | 3.250E-51 | -
|  | ROC | - | 5.131 | 5.835 | 5.509 | 4.042 | 5.046 | 4.707 | 4.728 | 5.222 | 5.028

  * *

The table data in this paper are abbreviated using scientific notation. For
instance, $2.006\times 10^{-14}$ is represented as 2.006E-14.

Table 4: Relative errors between total energy fluxes ${}_{s}\left<{dE/dt}\right>_{\infty,{\rm H}}$ (total angular momentum fluxes ${}_{s}\left<{dJ/dt}\right>_{\infty,{\rm H}}$) of four methods (${\bf N}=100$) and the exact solution at $r_{\rm ISCO}$(NA = not available).  |  | Relative errors of ${}_{s}\left<{dE/dt}\right>_{\infty}$ | Relative errors of ${}_{s}\left<{dE/dt}\right>_{\rm H}$
---|---|---|---
$s$ | Method | HeunC | MST-RTE | MST-RWE | 23PN | HeunC | MST-RTE | MST-RWE | 23PN
$-2$ | $\bf RE$ | 6.354E-73 | 2.769E-52 | 9.458E-51 | 1.496E-4 | 3.135E-76 | 2.308E-51 | 5.370E-51 | 1.750E-3
CPU(s) | 197.9 | 591.7 | 588.1 | NA | 53.0 | 91.2 | 91.9 | NA
$+2$ | $\bf RE$ | 7.421E-72 | 7.993E-53 | 1.303E-50 | NA | 6.222E-74 | 9.458E-51 | 5.370E-51 | NA
CPU(s) | 193.6 | 571.3 | 576.1 | NA | 54.7 | 88.8 | 93.1 | NA
$-1$ | $\bf RE$ | 6.274E-72 | 1.136E-52 | NA | NA | 4.633E-76 | 5.302E-52 | NA | NA
CPU(s) | 192.0 | 574.8 | NA | NA | 51.2 | 94.9 | NA | NA
$+1$ | $\bf RE$ | 3.880E-72 | 2.132E-52 | NA | NA | 9.675E-75 | 9.867E-51 | NA | NA
CPU(s) | 189.5 | 599.2 | NA | NA | 51.4 | 98.0 | NA | NA
$0$ | $\bf RE$ | 5.539E-59 | 1.343E-52 | NA | NA | 4.802E-53 | 3.451E-51 | NA | NA
CPU(s) | 139.2 | 295.5 | NA | NA | 29.2 | 52.5 | NA | NA
Table 5: Relative errors between energy fluxes ${}_{s}\left<{dE_{22}/dt}\right>_{\infty,{\rm H}}$ of the four methods (${\bf N}=100$) and the exact solution in different orbital radius $r_{0}$. In our tables, ${}_{s}\left<{dE_{22}/dt}\right>_{\infty,{\rm H}}$ is abbreviated as ${\dot{E}}_{\infty,{\rm H}}$. $s$ | Flux | $r_{0}$ | $5M$ | $10M$ | $20M$ | $50M$ | $100M$ | $200M$
---|---|---|---|---|---|---|---|---
$-2$ | ${\bf{RE}}({{{\dot{E}}_{\infty}}})$ | HeunC | 1.597E-93 | 1.716E-97 | 8.967E-100 | 9.232E-100 | 1.809E-100 | 4.346E-100
MST-RTE | 2.106E-52 | 2.002E-52 | 4.321E-52 | 4.066E-52 | 9.754E-54 | 1.015E-53
MST-RWE | 2.123E-50 | 5.023E-52 | 2.120E-51 | 7.431E-52 | 1.593E-52 | 1.924E-52
23PN | 2.664E-05 | 7.473E-12 | 1.467E-18 | 1.512E-27 | 1.673E-34 | 1.629E-42
${\bf{RE}}({{{\dot{E}}_{\rm H}}})$ | HeunC | 1.029E-89 | 1.139E-91 | 4.150E-91 | 2.178E-88 | 1.688E-86 | 2.870E-85
MST-RTE | 1.391E-50 | 2.807E-51 | 2.242E-51 | 1.287E-51 | 6.314E-53 | 1.248E-52
MST-RWE | 1.307E-50 | 3.793E-51 | 2.724E-51 | 2.342E-51 | 8.606E-52 | 9.625E-52
$+2$ | ${\bf{RE}}({{{\dot{E}}_{\infty}}})$ | HeunC | 1.413E-91 | 2.554E-97 | 1.601E-99 | 1.141E-100 | 4.792E-101 | 1.883E-100
MST-RTE | 4.323E-53 | 2.095E-52 | 3.307E-52 | 2.149E-53 | 8.198E-54 | 2.720E-54
MST-RWE | 2.123E-50 | 5.023E-52 | 2.120E-51 | 7.431E-52 | 1.593E-52 | 1.925E-52
${\bf{RE}}({{{\dot{E}}_{\rm H}}})$ | HeunC | 8.785E-90 | 9.953E-91 | 6.755E-91 | 5.800E-88 | 3.029E-85 | 9.663E-84
MST-RTE | 1.282E-50 | 2.780E-51 | 7.089E-52 | 2.396E-52 | 1.798E-52 | 4.947E-53
| MST-RWE | 1.307E-50 | 3.793E-51 | 2.724E-51 | 2.342E-51 | 8.606E-52 | 9.625E-52
$-1$ | ${\bf{RE}}({{{\dot{E}}_{\infty}}})$ | HeunC | 1.253E-93 | 1.080E-96 | 3.858E-99 | 3.665E-100 | 2.181E-98 | 1.043E-100
MST-RTE | 5.168E-52 | 5.281E-54 | 1.161E-55 | 2.090E-53 | 2.500E-53 | 3.472E-54
${\bf{RE}}({{{\dot{E}}_{\rm H}}})$ | HeunC | 1.258E-89 | 1.325E-91 | 1.056E-90 | 5.104E-88 | 1.731E-86 | 5.447E-84
| MST-RTE | 6.124E-49 | 2.571E-51 | 1.711E-51 | 1.244E-52 | 3.509E-52 | 1.629E-53
$+1$ | ${\bf{RE}}({{{\dot{E}}_{\infty}}})$ | HeunC | 1.248E-93 | 1.792E-97 | 1.361E-99 | 1.179E-101 | 2.179E-98 | 1.308E-100
MST-RTE | 3.377E-53 | 2.792E-54 | 2.844E-52 | 1.730E-53 | 7.509E-53 | 2.153E-53
${\bf{RE}}({{{\dot{E}}_{\rm H}}})$ | HeunC | 2.158E-90 | 3.051E-91 | 1.517E-90 | 6.192E-88 | 4.198E-86 | 8.736E-84
| MST-RTE | 1.485E-50 | 1.066E-51 | 3.775E-51 | 3.217E-52 | 1.009E-52 | 1.195E-53
$0$ | ${\bf{RE}}({{{\dot{E}}_{\infty}}})$ | HeunC | 2.622E-61 | 3.798E-61 | 6.767E-62 | 5.422E-63 | 7.733E-64 | 1.085E-64
MST-RTE | 1.903E-53 | 1.857E-52 | 3.542E-53 | 4.881E-53 | 2.547E-54 | 3.736E-55
${\bf{RE}}({{{\dot{E}}_{\rm H}}})$ | HeunC | 1.034E-58 | 1.527E-56 | 1.134E-55 | 8.936E-55 | 3.812E-54 | 1.572E-53
| MST-RTE | 1.142E-50 | 1.898E-51 | 8.722E-52 | 2.350E-52 | 3.121E-53 | 1.668E-52

This section presents several numerical results to validate the adaptability
of our method, and compare their accuracy with those which are already
available in the literature for calculating scalar, electromagnetic, and
gravitational radiations. These methods in the literature include numerical
integration (NI) method, high-order post-Newtonian expansion [114] , MST-RTE
method [60, 2] and MST-RWE method [59, 61], whose codes are provided by the
black hole perturbation toolkit (BHPT) [131]. In this subsection, there are
abbreviations of some methods, MST-RTE and MST-RWE 666MST-RWE first utilizes
the MST method to solve the RW equation, and then uses the Chandrasekhar-
Sasaki-Nakamura transformation to convert it into the solution of the
homogeneous Teukolsky equation. are the MST methods for solving the radial
Teukolsky equation (RTE) and Regge-Wheeler equation (RWE), respectively.
Meanwhile, 23PN represents 23th post-Newtonian order expansion, and the HeunC
method is the method proposed in this paper. The BHPT is a collection of
multiple scattered black hole perturbation theory codes, which have been
developed by various individuals or groups over several decades. Although the
BHPToolkit yields a formal solution of the confluent Heun function that solves
the homogenous Teukolsky equation, it solely provides the ingoing wave
solution $R_{\ell m\omega}^{{\rm{in}}}$, but the outgoing wave solution
$R_{\ell m\omega}^{{\rm{up}}}$has not been given so far. Therefore, this
formal solution (This solution is composed of the confluent Heun function, but
it is different from our solution.) of the BHPToolkit is an imperfect method
and cannot calculate radiation fluxes.

In numerical calculations for partial differential equations, the rate of
convergence (ROC) is a measure of how well a numerical solution approaches the
exact solution as the spatial or temporal resolution is increased. It provides
insight into the accuracy and reliability of the method. Therefore, we
introduce the ROC of the floating-point numbers777$\bf N$ is the software
floating-point numbers, which can affect the computational efficiency. The
larger $\bf N$, the more calculation time is required. $\bf N$ to evaluate the
speed of convergence for the computational methods. And the ROC of floating-
point numbers is defined as [132, 133, 134]

${\rm{ROC}}={\log_{10}}\left(\frac{{\bf{RE}}_{\bf N}}{{\bf{RE}}_{{\bf
N}+10}}\right),$ (106)

where ${{\bf{RE}}_{{\bf N}}}$ and ${{\bf{RE}}_{{\bf N}+10}}$ are the relative
error (RE) of ${\bf N}$ and ${\bf N}+10$, respectively. The exact solution in
this paper is the results of the MST-RTE method with ${\bf N}=200$. In our
Tables, CPU is computing time. All the numerical experiments are conducted on
an Intel Core i7-12700H 2.70 GHz processor.

To validate the correctness and applicability of the HeunC method, we test two
kind of radiation fluxes ${}_{s}\left<{dE_{\ell m}/dt}\right>_{\infty,{\rm
H}}$ propagates along the ($+r$) direction. Fig. 2 and 3 show logarithmic
values of ${}_{s}\left<{dE_{\ell m}/dt}\right>_{\infty,{\rm H}}$ of the three
perturbation fields (scalar, electromagnetic, and gravitational perturbations)
with different modes, over the domain $[5M,15M]$. In order to see the value of
the energy fluxes at $r_{\rm{ISCO}}$888The innermost stable circular
orbit(ISCO) is only defined in the equatorial plane, that is $r_{\rm
ISCO}=6M$., we draw a slice diagram of the three perturbation fields as shown
in Fig. 4 and 5. It can be seen from Fig. 2, 3, 4 and 5 that when ${\bf
N}=50$, the energy fluxes obtained by the HeunC method exhibit a commendable
agreement with the results obtained by the MST-RTE method. This numerical
simulation confirms the physical phenomenon that the closer a point particle
approaches a black hole, the greater the amount of radiative energy flux it
generates. Furthermore, when $r_{0}$ is fixed, the maximum energy flux across
all modes are ${}_{s}\left<{dE_{22}/dt}\right>_{\infty,{\rm H}}$. According to
the properties of the spin-weighted spherical harmonic function ${}_{s}Y_{\ell
m\omega}$, when $|\ell|-|m|={\rm odd}$, ${}_{0}Y_{\ell m\omega}=0$. From Eqs.
105a and 105b, when ${}_{0}Y_{\ell m\omega}$ is equal to zero, the energy
fluxes are also zero. Therefore, when $|\ell|-|m|={\rm odd}$, the logarithmic
energy flux of the scalar perturbation field tends to infinity. This property
can be verified by the logarithmic energy fluxes of the scalar perturbation
field in Fig. 4 and 5.

To facilitate a comprehensive comparison of the accuracy and efficiency of the
HeunC method with other existing methods, we provide numerical comparisons of
energy fluxes ${}_{s}\left<{dE_{22}/dt}\right>_{\infty,{\rm H}}$ of three
perturbation fields obtained from four methods for different floating-point
numbers ${\bf N}$ at $r_{\rm ISCO}$ in Table 3. Meanwhile, the error
comparison of total radiation fluxes is shown in Table 4. Because the relative
errors in Table 3 and Table 4 are not less than $10^{\bf N}$, the error
results of the energy flux are the same as those of the angular momentum flux
(their relationship can be seen in Eq. 27). Analyzing the ROC presented in
Table 3, it can be concluded that the HeunC method is a 10th-order method for
the calculation of electromagnetic and gravitational luminosities, and the
MST-RTE and MST-RWE methods are 5th-order methods. Despite a decrease in ROC
as ${\bf N}$ increases for the scalar perturbation field, the HeunC method
maintains superior accuracy in calculating the energy fluxes compared to the
MST-RTE method. This behavior is evident from the observations presented in
Fig. 6. As $\bf N$ increases, the HeunC method exhibits a notably larger
decrease in relative errors in contrast to both the MST-RTE and MST-RWE
methods. These results show that the convergence speed of HeunC is faster than
that of the MST-RTE and MST-RWE methods.

To explore the precision of energy fluxes across various modes, Table 4
presents the relative errors of total energy flux for three perturbation
fields. To obtain the total energy flux, it is necessary to sum $\ell$ in the
range $2\leq\ell\leq 25$ ($2\leq\ell\leq 11$) for the infinity (horizon) part.
It can be seen from Table 4 that the HeunC method calculates electromagnetic
and gravitational energy fluxes, and obtains more accurate and efficient
results than other methods. Even when calculating the scalar energy fluxes,
the accuracy is attenuated, but it is still superior to other methods. Table 5
exhibits the relative errors of energy fluxes within a region
$r_{0}\in[5M,200M]$, aiming to assess the precision of four distinct methods
as $r_{0}$ is close to the horizon or infinity. Notably, in both near the
horizon and at infinity, the HeunC method showcases exceptional accuracy in
calculating the energy fluxes, surpassing the other three methods by a
significant margin. Moreover, ROC of 23PN method presented in Table 3 and
relative errors in Table 5, reveal that the relative error of the 23PN method
decreases as $r_{0}$ increases, while remaining unaffected by the floating-
point numbers. This also exposes a limitation of the 23 PN method: it exhibits
inaccuracies in the vicinity of the horizon, but when
$r_{0}\rightarrow\infty$, the results obtained from this method become
increasingly accurate.

All numerical results of this test show that the computational accuracy and
efficiency of the HeunC method surpasses that of the other three methods to a
significant extent. This is advantageous for waveform construction.

##  VIII Conclusion

The purpose of this paper is to present an exact method to calculate
gravitational, electromagnetic and scalar radiation fluxes for any type-D BHs.
We first reformulate the radial Teukolsky equations into a general form, which
is classified into multiple types based on $\Delta_{n}$. Then, we focus our
attention on solving the $\Delta_{2}$-type Teukolsky equation which includes
ten kinds of Teukolsky equations shown in Table 1. The general solution (IV.1)
of homogenous form (29) of the $\Delta_{2}$-type Teukolsky equation is
constituted by the linear combination of two linearly independent solutions
which are expressed by confluent Heun functions. This solution, benefiting
from the analytical asymptotic expression (V.4) of the confluent Heun function
at infinity which is obtained in this paper for the first time, can be applied
to the physical models corresponding to various boundary conditions.
Meanwhile, ingoing wave and outgoing wave solutions $R_{\ell
m\omega}^{{\rm{in,up}}}$ can directly be obtained from the exact general
solution according to boundary conditions. In a significant contribution, , we
provide a complete and strict derivation process for expressing $R_{\ell
m\omega}^{{\rm{in,up}}}$ in terms of confluent Heun functions. Consequently,
the exact solution of inhomogeneous GTEs is found utilizing Green’s function.
It’s worth noting that these exact solutions are not subject to any
constraints (such as slow-motion, low-frequency, and weak-field). Next,
Green’s function method is employed to evaluate $Z^{\infty,{\rm H}}_{\ell
m\omega}$, which is used in Eqs. 25 and 26 to compute the energy flux of the
gravitational, electromagnetic and scalar waves. Finally, the HeunC Method is
used for numerical simulations of the radiation fluxes of the Schwarzschild
black hole, and the results are entirely satisfactory in comparison with the
analytical (23th post-Newtonian order expansion) and numerical methods [62,
63]. The findings are summarized in detail as follows:

1\. For calculating the gravitational, electromagnetic, and scalar fluxes for
Schwarzschild black hole, our method can accurately calculate any mode of the
entire spatial region $r_{0}\in[r_{+},\infty)$, and its accuracy is much
higher than the results obtained by other methods. Additionally, the PN
expansion suffers from slow convergence and limitations in the low-frequency
approximation, making its results near the event horizon less precise. Even in
the far distance, our approach produces superior results to those obtained
using the PN expansion, MST-RTE, and MST-RWE methods. And the computational
time of our method is less than half of other numerical methods. Thus, our
method is superior to both the MST and PN methods in terms of its exceptional
precision and efficiency.

2\. Fujita and Tagoshi have demonstrated that the MST method converges very
fast [62, 63]. However, they did not provide a mathematical evaluation of the
convergence specifically for floating-point numbers. The increase in the
floating-point numbers is accompanied by improved accuracy but also increased
computation time. To this end, we introduce ROC for floating-point numbers. By
comparing the ROC of the four methods, it has been determined that HeunC is a
10th-order method, while MST is 5th-order methods. But, the 23PN method is not
convergent about the floating-point numbers. These comparisons mean that for
each additional floating-point number $\bf N$, the relative error of the HeunC
method will be multiplied by $10^{-1}$, whereas the relative error of the MST
and NI methods will be multiplied by $5\times 10^{-1}$.

3\. Due to the complexity of the source term of the neutron star [135], the
radiation fluxes calculated in this paper only considers the gravitational,
electromagnetic, and scalar perturbations. However, if the source term of the
neutron star is known, our method can also find the radiation fluxes of the
neutron star. In theory, our method can be applied to compute the radiation
fluxes of any perturbation fields, as long as their source terms are known.

All in all, our method exhibits high precision, high efficiency, and excellent
applicability. It can be applied to calculate the radiation fluxes of any
perturbation fields quickly and accurately. Accurate and efficient energy flux
calculation is the basis for obtaining GW waveforms. Therefore, in future
work, we will use the HeunC method to derive the radiation fluxes and GW
waveforms for the test particles moving around a rotating black hole in a
general orbit. Furthermore, our general solution is beneficial for solving the
Teukolsky equation in complicated spacetimes ( Kerr-Newman anti-de Sitter BHs
[91, 92, 93] ), modified gravity theories [136], or effective-one-body
theories based on post-Minkowskian approximation [137, 138].

## Acknowledgement

This work was supported by the Grant of NSFC No. 12035005, and National Key
Research and Development Program of China No. 2020YFC2201400.

##  Appendix A Expression of ${}_{s}\tilde{Z}_{\ell m\omega}^{{\rm
H},\infty}$, ${{}_{s}{\beta_{\ell m\omega}}}$, and ${{}_{s}{\alpha_{\ell
m\omega}}}$

Substituting the point sources ${}_{s}T_{\ell m\omega}$ of gravitational,
electromagnetic, and scalar perturbations in Ref. [6] into Eq. 22
respectively, we can obtain the expression of ${}_{s}\tilde{Z}_{\ell
m\omega}^{{\rm H},\infty}$.

${}_{s}\tilde{Z}_{\ell m\omega}^{{\rm H},\infty}=\frac{1}{{W_{\rm
C}}}\int_{-\infty}^{\infty}dt{e^{i\omega t-im\varphi(t)}}({}_{s}{\cal I}_{\ell
m\omega}^{\infty,{\rm H}})_{r=r(t),\theta=\theta(t)},$ (107)

For the gravitational, electromagnetic, and scalar perturbations, we can
incorporate the general formula (107). Each perturbation field introduces
distinct functions ${}_{s}{\cal I}_{\ell m\omega}^{\infty,{\rm H}}$. In this
section, we provide the specific form of ${}_{s}{\cal I}_{\ell
m\omega}^{\infty,{\rm H}}$ associated with different perturbation fields.

###  A.1 Gravitational field

For the general formula (107), ${}_{s}{\cal I}_{\ell m\omega}^{\infty,{\rm
H}}$ of the gravitational perturbation field of the spin-weight $s=-2$, and
the coefficients in fluxes (25) and (26) can be given by

$\displaystyle{}_{-2}{\cal I}_{\ell m\omega}^{\infty,\rm{H}}$
$\displaystyle=\frac{8\pi}{\Upsilon_{t}}\left[{\left({{A_{nn0}}+{A_{\bar{m}n0}}+{A_{\bar{m}\bar{m}0}}}\right)R_{\ell
m\omega}^{{\rm{in,up}}}}\right.$
$\displaystyle-\left({{A_{\bar{m}n1}}+{A_{\bar{m}\bar{m}1}}}\right){\left({R_{\ell
m\omega}^{{\rm{in,up}}}}\right)^{\prime}}{\left.+{A_{\bar{m}\bar{m}2}}{{\left({R_{\ell
m\omega}^{{\rm{in,up}}}}\right)}^{\prime\prime}}\right]},$ (108)

and

$\displaystyle{}_{-2}\alpha_{\ell m\omega}$
$\displaystyle=\frac{{256{{(2M{r_{\rm{H}}})}^{5}}({\omega^{2}}+4{{\tilde{\epsilon}}^{2}})({\omega^{2}}+16{{\tilde{\epsilon}}^{2}}){\omega^{4}}}}{{|{{\bf{C}}_{-2}}{|^{2}}}},$
(109a) $\displaystyle{}_{-2}\beta_{\ell m\omega}$ $\displaystyle=1.$ (109b)

where $\tilde{\epsilon}={(4{r_{\rm{H}}})^{-1}}$ and
${{|{\bf{C}}_{-2}|{{}^{2}}}}={\lambda^{2}}{(\lambda+2)^{2}}+144{\omega^{2}}{M^{2}}.$
And the explicit form of $A_{nn0}$ and other terms can be found and derived in
Ref. [2, 6].

Similarly, ${}_{s}{\cal I}_{\ell m\omega}^{\infty,{\rm H}}$ of the spin-weight
$s=+2$ and the coefficients in fluxes (25) and (26) case can be given by

$\displaystyle{}_{+2}{\cal I}_{\ell m\omega}^{\infty,\rm{H}}$
$\displaystyle=\frac{8\pi}{\Upsilon_{t}}\left[{\left({{A_{ll0}}+{A_{lm0}}+{A_{mm0}}}\right)R_{\ell
m\omega}^{{\rm{in,up}}}}\right.$
$\displaystyle-\left({{A_{lm1}}+{A_{mm1}}}\right){\left({R_{\ell
m\omega}^{{\rm{in,up}}}}\right)^{\prime}}{\left.+{A_{mm2}}{{\left({R_{\ell
m\omega}^{{\rm{in,up}}}}\right)}^{\prime\prime}}\right]},$ (110)

and

$\displaystyle{}_{+2}\alpha_{\ell m\omega}$
$\displaystyle=\frac{{{\omega^{2}}}}{{128r_{\rm{H}}^{3}\left({{\omega^{2}}+4{\epsilon^{2}}}\right)}},$
(111a) $\displaystyle{}_{+2}\beta_{\ell m\omega}$
$\displaystyle=\frac{{16{\omega^{8}}}}{{|{{\bf{C}}_{2}}{|^{2}}}}.$ (111b)

with $|{{\bf
C}_{2}}|^{2}={(\lambda+4)^{2}}{(\lambda+6)^{2}}+144{M^{2}}{\omega^{2}}$.

###  A.2 Electromagnetic field

For the general formula (107), ${}_{s}{\cal I}_{\ell m\omega}^{\infty,{\rm
H}}$ of the electromagnetic perturbation field of the spin-weight $s=-1$, and
the coefficients in fluxes (25) and (26) can be given by

${\cal I}_{\ell
m\omega}^{\infty,\rm{H},-1}={\frac{\mathfrak{A}_{-}}{\Delta}{R_{\ell
m\omega}^{{\rm{in,up}}}}-\frac{\mathfrak{B}_{-}}{\Delta}{\left({R_{\ell
m\omega}^{{\rm{in,up}}}}\right)^{\prime}}},$ (112)

and

$\displaystyle{}_{-1}{\alpha_{\ell m\omega}}$
$\displaystyle=\frac{{64{r_{\rm{H}}}{M^{3}}{\omega^{4}}\left({4{r_{\rm{H}}}^{2}{\omega^{2}}+1}\right)}}{{{\lambda^{2}}}},$
(113a) $\displaystyle{}_{-1}{\beta_{\ell m\omega}}$ $\displaystyle=2.$ (113b)

where the explicit form of ${\mathfrak{A}_{\pm}}$ and ${\mathfrak{B}_{\pm}}$
can be found in Eq. (48) of Ref. [113]

Similarly, ${}_{s}{\cal I}_{\ell m\omega}^{\infty,{\rm H}}$ of the spin-weight
$s=+1$ and the coefficients in fluxes (25) and (26) case can be given by

${}_{+1}{\cal I}_{\ell m\omega}^{\infty,\rm{H}}={{\mathfrak{A}_{+}}R_{\ell
m\omega}^{{\rm{in}},{\rm{up}}}-\frac{\mathfrak{B}_{+}}{\Delta}{\left(\Delta
R_{\ell m\omega}^{{\rm{in}},{\rm{up}}}\right)^{\prime}}},$ (114)

and

$\displaystyle{}_{+1}{\alpha_{\ell m\omega}}$
$\displaystyle=\frac{{{\omega^{2}}}}{{4r_{\rm H}}},$ (115a)
$\displaystyle{}_{+1}{\beta_{\ell m\omega}}$
$\displaystyle=\frac{{8{\omega^{4}}}}{{{{(\lambda+2)}^{2}}}},$ (115b)

###  A.3 Scalar field

For the general formula (107), ${}_{s}{\cal I}_{\ell m\omega}^{\infty,{\rm
H}}$ of the scalar perturbation field of the spin-weight $s=0$, and the
coefficients in fluxes (25) and (26) can be given by

${}_{0}{{\cal I}_{\ell
m\omega}^{\infty,\rm{H}}}=\frac{-4\pi{r^{2}_{0}}}{\Upsilon_{t}}{}_{0}Y_{\ell
m\omega}R_{\ell m\omega}^{{\rm{in}},{\rm{up}}},$ (116)

and

$\displaystyle{}_{0}{\alpha_{\ell m\omega}}$ $\displaystyle=4M{\omega^{4}},$
(117a) $\displaystyle{}_{0}{\beta_{\ell m\omega}}$ $\displaystyle=1.$ (117b)

where ${}_{s}Y_{\ell m\omega}$ is the spin-weighted spherical harmonics.

## References

  * Mino  _et al._ [1997] Y. Mino, M. Sasaki, M. Shibata, H. Tagoshi, and T. Tanaka, Chapter 1. Black Hole Perturbation, Progress of Theoretical Physics Supplement 128, 1 (1997).
  * Sasaki and Tagoshi [2003] M. Sasaki and H. Tagoshi, Analytic black hole perturbation approach to gravitational radiation, Living Reviews in Relativity 6 (2003).
  * Regge and Wheeler [1957] T. Regge and J. A. Wheeler, Stability of a Schwarzschild singularity, Physical Review 108, 1063 (1957).
  * Zerilli [1970] F. J. Zerilli, Gravitational field of a particle falling in a Schwarzschild geometry analyzed in tensor harmonics, Physical Review D 2, 2141 (1970).
  * Bardeen and Press [1973] J. M. Bardeen and W. H. Press, Radiation fields in the Schwarzschild background, Journal of Mathematical Physics 14, 7 (1973).
  * Teukolsky [1973] S. A. Teukolsky, Perturbations of a rotating black hole. I. fundamental equations for gravitational, electromagnetic, and neutrino-field perturbations, The Astrophysical Journal 185, 635 (1973).
  * Abbott  _et al._ [2016a] B. P. Abbott _et al._ , Improved analysis of GW150914 using a fully spin-precessing waveform model, Physical Review X 6, 041014 (2016a).
  * Abbott  _et al._ [2016b] B. P. Abbott _et al._ , Binary Black Hole Mergers in the First Advanced LIGO Observing Run, Physical Review X 6, 041015 (2016b).
  * Abbott  _et al._ [2019] B. P. Abbott _et al._ , GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs, Physical Review X 9, 031040 (2019).
  * Abbott  _et al._ [2021] B. P. Abbott _et al._ , GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo during the First Half of the Third Observing Run, Physical Review X 11, 021053 (2021).
  * Abbott  _et al._ [2023] B. P. Abbott _et al._ , Population of Merging Compact Binaries Inferred Using Gravitational Waves through GWTC-3, Physical Review X 13, 011048 (2023).
  * Amaro-Seoane  _et al._ [2017] P. Amaro-Seoane _et al._ , Laser Interferometer Space Antenna, (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Luo  _et al._ [2016] J. Luo _et al._ , TianQin: a space-borne gravitational wave detector, Classical and Quantum Gravity 33, 035010 (2016).
  * Mei  _et al._ [2020] J. Mei _et al._ , The TianQin project: Current progress on science and technology, Progress of Theoretical and Experimental Physics 2021, 05A107 (2020).
  * Gong  _et al._ [2021] Y. Gong, J. Luo, and B. Wang, Concepts and status of chinese space gravitational wave detection projects, Nature Astronomy 5, 881 (2021).
  * Ruan  _et al._ [2020] W.-H. Ruan, Z.-K. Guo, R.-G. Cai, and Y.-Z. Zhang, Taiji program: Gravitational-wave sources, International Journal of Modern Physics A 35, 2050075 (2020).
  * Amaro-Seoane [2018] P. Amaro-Seoane, Relativistic dynamics and extreme mass ratio inspirals, Living Reviews in Relativity 21 (2018).
  * Amaro Seoane [2020] P. Amaro Seoane, The gravitational capture of compact objects by massive black holes, in  _Handbook of Gravitational Wave Astronomy_, edited by C. Bambi, S. Katsanevas, and K. D. Kokkotas (Springer Singapore, Singapore, 2020) pp. 1–79.
  * Isoyama  _et al._ [2022] S. Isoyama, R. Fujita, A. J. K. Chua, H. Nakano, A. Pound, and N. Sago, Adiabatic Waveforms from Extreme-Mass-Ratio Inspirals: An Analytical Approach, Physical Review Letters 128, 231101 (2022).
  * Poisson  _et al._ [2011] E. Poisson, A. Pound, and I. Vega, The Motion of Point Particles in Curved Spacetime, Living Reviews in Relativity 14, 7 (2011).
  * Pound and Wardell [2020] A. Pound and B. Wardell, Black hole perturbation theory and gravitational self-force, in  _Handbook of Gravitational Wave Astronomy_, edited by C. Bambi, S. Katsanevas, and K. D. Kokkotas (Springer Singapore, Singapore, 2020) pp. 1–119.
  * Babak  _et al._ [2017] S. Babak, J. Gair, A. Sesana, E. Barausse, C. F. Sopuerta, C. P. L. Berry, E. Berti, P. Amaro-Seoane, A. Petiteau, and A. Klein, Science with the space-based interferometer LISA. V. Extreme mass-ratio inspirals, Physical Review D 95, 103012 (2017).
  * Askar  _et al._ [2019] A. Askar _et al._ , Black holes, gravitational waves and fundamental physics: a roadmap, Classical and Quantum Gravity 36, 143001 (2019).
  * Amaro-Seoane  _et al._ [2007] P. Amaro-Seoane, J. R. Gair, M. Freitag, M. C. Miller, I. Mandel, C. J. Cutler, and S. Babak, Intermediate and extreme mass-ratio inspirals—astrophysics, science applications and detection using lisa, Classical and Quantum Gravity 24, R113 (2007).
  * Fan  _et al._ [2020] H.-M. Fan, Y.-M. Hu, E. Barausse, A. Sesana, J.-d. Zhang, X. Zhang, T.-G. Zi, and J. Mei, Science with the TianQin observatory: Preliminary result on extreme-mass-ratio inspirals, Physical Review D 102, 063016 (2020).
  * Zi  _et al._ [2021] T.-G. Zi, J.-D. Zhang, H.-M. Fan, X.-T. Zhang, Y.-M. Hu, C. Shi, and J. Mei, Science with the TianQin Observatory: Preliminary results on testing the no-hair theorem with extreme mass ratio inspirals, Physical Review D 104, 064008 (2021).
  * Leaver [1986a] E. W. Leaver, Solutions to a generalized spheroidal wave equation: Teukolsky’s equations in general relativity, and the two-center problem in molecular quantum mechanics, Journal of Mathematical Physics 27, 1238 (1986a).
  * Leaver and Chandrasekhar [1985] E. W. Leaver and S. Chandrasekhar, An analytic representation for the quasi-normal modes of Kerr black holes, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 402, 285 (1985).
  * Leaver [1986b] E. W. Leaver, Spectral decomposition of the perturbation response of the Schwarzschild geometry, Physical Review D 34, 384 (1986b).
  * Ronveaux and Arscott [1995] A. Ronveaux and F. M. Arscott,  _Heun’s differential equations_ (Oxford University Press, 1995).
  * Slavyanov and Lay [2000] S. Slavyanov and W. Lay,  _Special Functions: Unified Theory Based on Singularities_ (Oxford University Press, 2000).
  * Motygin [2015] O. V. Motygin, On numerical evaluation of the Heun functions, in  _2015 Days on Diffraction (DD)_ (IEEE, 2015).
  * Motygin [2018] O. V. Motygin, On evaluation of the confluent Heun functions, in  _2018 Days on Diffraction (DD)_ (IEEE, 2018).
  * Fiziev [2006] P. P. Fiziev, Exact solutions of Regge–Wheeler equation and quasi-normal modes of compact objects, Classical and Quantum Gravity 23, 2447 (2006).
  * Fiziev [2007] P. P. Fiziev, Exact solutions of Regge-Wheeler equation, Journal of Physics: Conference Series 66, 012016 (2007).
  * Fiziev [2010] P. P. Fiziev, Classes of exact solutions to the Teukolsky master equation, Classical and Quantum Gravity 27, 135001 (2010).
  * Fiziev [2009] P. P. Fiziev, Classes of exact solutions to Regge-Wheeler and Teukolsky equations (2009), arXiv:0902.1277 [gr-qc] .
  * Fiziev and Staicova [2011] P. Fiziev and D. Staicova, Application of the confluent Heun functions for finding the quasinormal modes of nonrotating black holes, Physical Review D 84, 127502 (2011).
  * Borissov and Fiziev [2010] R. S. Borissov and P. P. Fiziev, Exact solutions of Teukolsky master equation with continuous spectrum (2010), arXiv:0903.3617 [gr-qc] .
  * Fiziev and Staicova [2009a] P. P. Fiziev and D. R. Staicova, A new model of the central engine of grb and the cosmic jets (2009a), arXiv:0902.2408 [astro-ph.HE] .
  * Staicova and Fiziev [2010] D. Staicova and P. Fiziev, The spectrum of electromagnetic jets from kerr black holes and naked singularities in the Teukolsky perturbation theory, Astrophysics and Space Science 332, 385 (2010).
  * Gurtas Dogan and Sucu [2019] S. Gurtas Dogan and Y. Sucu, Quasinormal modes of Dirac field in 2+1 dimensional gravitational wave background, Physics Letters B 797, 134839 (2019).
  * Zhang  _et al._ [2014] H. Zhang, D.-J. Liu, and X.-Z. Li, Black holes and gravitational waves in three-dimensional $f(r)$ gravity, Physical Review D 90, 124051 (2014).
  * Hortaçsu [2021] M. Hortaçsu, The radial Teukolsky equation for Kerr-Newman-de Sitter geometry: revisited, The European Physical Journal Plus 136 (2021).
  * Hatsuda [2020] Y. Hatsuda, Quasinormal modes of Kerr–de Sitter black holes via the Heun function, Classical and Quantum Gravity 38, 025015 (2020).
  * Noda and Motohashi [2022] S. Noda and H. Motohashi, Spectroscopy of Kerr-AdS5 spacetime with the Heun function: Quasinormal modes, greybody factor, and evaporation, Physical Review D 106, 064025 (2022).
  * Cook and Zalutskiy [2014] G. B. Cook and M. Zalutskiy, Gravitational perturbations of the Kerr geometry: High-accuracy study, Physical Review D 90, 124021 (2014).
  * Press and Teukolsky [1973] W. H. Press and S. A. Teukolsky, Perturbations of a rotating black hole. II. dynamical stability of the Kerr metric, The Astrophysical Journal 185, 649 (1973).
  * Chrzanowski [1975] P. L. Chrzanowski, Vector potential and metric perturbations of a rotating black hole, Physical Review D 11, 2042 (1975).
  * Tashiro and Ezawa [1981] Y. Tashiro and H. Ezawa, Gravitational radiation of a particle falling towards a black hole. I: The case of non-rotating black hole, Progress of Theoretical Physics 66, 1612 (1981).
  * Chandrasekhar [1998] S. Chandrasekhar,  _The mathematical theory of black holes_ , Vol. 69 (Oxford university press, 1998).
  * Nakamura  _et al._ [1987] T. Nakamura, K. Oohara, and Y. Kojima, General relativistic collapse to black holes and gravitational waves from black holes, Progress of Theoretical Physics Supplement 90, 1 (1987).
  * Sasaki [1994] M. Sasaki, Post-Newtonian expansion of the ingoing-wave Regge-Wheeler function, Progress of Theoretical Physics 92, 17 (1994).
  * Tagoshi and Sasaki [1994] H. Tagoshi and M. Sasaki, Post-newtonian expansion of gravitational waves from a particle in circular orbit around a schwarzschild black hole, Progress of Theoretical Physics 92, 745 (1994).
  * Shibata  _et al._ [1995] M. Shibata, M. Sasaki, H. Tagoshi, and T. Tanaka, Gravitational waves from a particle orbiting around a rotating black hole: Post-Newtonian expansion, Physical Review D 51, 1646 (1995).
  * Tanaka  _et al._ [1996] T. Tanaka, H. Tagoshi, and M. Sasaki, Gravitational waves by a particle in circular orbits around a Schwarzschild black hole: 5.5 Post-Newtonian formula, Progress of Theoretical Physics 96, 1087 (1996).
  * Tagoshi  _et al._ [1996] H. Tagoshi, M. Shibata, T. Tanaka, and M. Sasaki, Postnewtonian expansion of gravitational waves from a particle in circular orbits around a rotating black hole: Up to O$(v^{8})$ beyond the quadrupole formula, Physical Review D 54, 1439 (1996).
  * Poisson and Sasaki [1995] E. Poisson and M. Sasaki, Gravitational radiation from a particle in circular orbit around a black hole. V. Black-hole absorption and tail corrections, Physical Review D 51, 5753 (1995).
  * Mano  _et al._ [1996a] S. Mano, H. Suzuki, and E. Takasugi, Analytic Solutions of the Regge-Wheeler Equation and the Post-Minkowskian Expansion, Progress of Theoretical Physics 96, 549 (1996a).
  * Mano  _et al._ [1996b] S. Mano, H. Suzuki, and E. Takasugi, Analytic solutions of the Teukolsky equation and their low frequency expansions, Progress of Theoretical Physics 95, 1079 (1996b).
  * Casals and Ottewill [2015] M. Casals and A. Ottewill, High-order tail in Schwarzschild spacetime, Physical Review D 92, 124055 (2015).
  * Fujita and Tagoshi [2004] R. Fujita and H. Tagoshi, New numerical methods to evaluate homogeneous solutions of the teukolsky equation, Progress of Theoretical Physics 112, 415 (2004).
  * Fujita and Tagoshi [2005] R. Fujita and H. Tagoshi, New numerical methods to evaluate homogeneous solutions of the teukolsky equation. II: – solutions of the continued fraction equation –, Progress of Theoretical Physics 113, 1165 (2005).
  * Fujita [2012] R. Fujita, Gravitational waves from a particle in circular orbits around a Schwarzschild black hole to the 22nd Post-Newtonian order, Progress of Theoretical Physics 128, 971 (2012).
  * Fujita [2015a] R. Fujita, Gravitational waves from a particle in circular orbits around a rotating black hole to the 11th post-newtonian order, Progress of Theoretical and Experimental Physics 2015, 33E01 (2015a).
  * Fujita and Iyer [2010] R. Fujita and B. R. Iyer, Spherical harmonic modes of 5.5 post-newtonian gravitational wave polarizations and associated factorized resummed waveforms for a particle in circular orbit around a schwarzschild black hole, Physical Review D 82, 044051 (2010).
  * Sago and Fujita [2015] N. Sago and R. Fujita, Calculation of radiation reaction effect on orbital parameters in Kerr spacetime, Progress of Theoretical and Experimental Physics 2015, 073E03 (2015).
  * Hikida  _et al._ [2004] W. Hikida, S. Jhingan, H. Nakano, N. Sago, M. Sasaki, and T. Tanaka, A new analytical method for self-force regularization. I: – charged scalar particles in schwarzschild spacetime –, Progress of Theoretical Physics 111, 821 (2004).
  * Hikida  _et al._ [2005] W. Hikida, S. Jhingan, H. Nakano, N. Sago, M. Sasaki, and T. Tanaka, A new analytical method for self-force regularization. II: – testing the efficiency for circular orbits –, Progress of Theoretical Physics 113, 283 (2005).
  * Casals  _et al._ [2013] M. Casals, S. Dolan, A. C. Ottewill, and B. Wardell, Self-force and Green function in Schwarzschild spacetime via quasinormal modes and branch cut, Physical Review D 88, 044022 (2013).
  * Cutler  _et al._ [1994] C. Cutler, D. Kennefick, and E. Poisson, Gravitational radiation reaction for bound motion around a schwarzschild black hole, Physical Review D 50, 3816 (1994).
  * Finn and Thorne [2000] L. S. Finn and K. S. Thorne, Gravitational waves from a compact star in a circular, inspiral orbit, in the equatorial plane of a massive, spinning black hole, as observed by lisa, Physical Review D 62, 124021 (2000).
  * Glampedakis and Kennefick [2002] K. Glampedakis and D. Kennefick, Zoom and whirl: Eccentric equatorial orbits around spinning black holes and their evolution under gravitational radiation reaction, Physical Review D 66, 044002 (2002).
  * Drasco and Hughes [2006] S. Drasco and S. A. Hughes, Gravitational wave snapshots of generic extreme mass ratio inspirals, Physical Review D 73, 024027 (2006).
  * Hughes  _et al._ [2021] S. A. Hughes, N. Warburton, G. Khanna, A. J. K. Chua, and M. L. Katz, Adiabatic waveforms for extreme mass-ratio inspirals via multivoice decomposition in time and frequency, Physical Review D 103, 104014 (2021).
  * Han [2010] W.-B. Han, Gravitational radiation from a spinning compact object around a supermassive Kerr black hole in circular orbit, Physical Review D 82, 084013 (2010).
  * Harms  _et al._ [2016a] E. Harms, G. Lukes-Gerakopoulos, S. Bernuzzi, and A. Nagar, Asymptotic gravitational wave fluxes from a spinning particle in circular equatorial orbits around a rotating black hole, Physical Review D 93, 044015 (2016a).
  * Harms  _et al._ [2016b] E. Harms, G. Lukes-Gerakopoulos, S. Bernuzzi, and A. Nagar, Spinning test body orbiting around a Schwarzschild black hole: Circular dynamics and gravitational-wave fluxes, Physical Review D 94, 104010 (2016b).
  * Lukes-Gerakopoulos  _et al._ [2017] G. Lukes-Gerakopoulos, E. Harms, S. Bernuzzi, and A. Nagar, Spinning test body orbiting around a Kerr black hole: Circular dynamics and gravitational-wave fluxes, Physical Review D 96, 064051 (2017).
  * Akcay  _et al._ [2020] S. Akcay, S. R. Dolan, C. Kavanagh, J. Moxon, N. Warburton, and B. Wardell, Dissipation in extreme mass-ratio binaries with a spinning secondary, Physical Review D 102, 064013 (2020).
  * Piovano  _et al._ [2020] G. A. Piovano, A. Maselli, and P. Pani, Extreme mass ratio inspirals with spinning secondary: A detailed study of equatorial circular motion, Physical Review D 102, 024041 (2020).
  * Piovano  _et al._ [2021] G. A. Piovano, R. Brito, A. Maselli, and P. Pani, Assessing the detectability of the secondary spin in extreme mass-ratio inspirals with fully relativistic numerical waveforms, Physical Review D 104, 124019 (2021).
  * Skoupý and Lukes-Gerakopoulos [2021a] V. Skoupý and G. Lukes-Gerakopoulos, Gravitational wave templates from Extreme Mass Ratio Inspirals, arXiv e-prints , arXiv:2101.04533 (2021a), arXiv:2101.04533 [gr-qc] .
  * Rahman and Bhattacharyya [2023] M. Rahman and A. Bhattacharyya, Prospects for determining the nature of the secondaries of extreme mass-ratio inspirals using the spin-induced quadrupole deformation, Physical Review D 107, 024006 (2023).
  * Mathews  _et al._ [2022] J. Mathews, A. Pound, and B. Wardell, Self-force calculations with a spinning secondary, Physical Review D 105, 084031 (2022).
  * Skoupý and Lukes-Gerakopoulos [2021b] V. Skoupý and G. Lukes-Gerakopoulos, Spinning test body orbiting around a Kerr black hole: Eccentric equatorial orbits and their asymptotic gravitational-wave fluxes, Physical Review D 103, 104045 (2021b).
  * Skoupý and Lukes-Gerakopoulos [2022] V. Skoupý and G. Lukes-Gerakopoulos, Adiabatic equatorial inspirals of a spinning body into a Kerr black hole, Physical Review D 105, 084033 (2022).
  * Drummond and Hughes [2022a] L. V. Drummond and S. A. Hughes, Precisely computing bound orbits of spinning bodies around black holes. I. General framework and results for nearly equatorial orbits, Physical Review D 105, 124040 (2022a).
  * Drummond and Hughes [2022b] L. V. Drummond and S. A. Hughes, Precisely computing bound orbits of spinning bodies around black holes. II. Generic orbits, Physical Review D 105, 124041 (2022b).