# Black Hole Greybody Factors from Korteweg-de Vries Integrals: Computation

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

###### Abstract

It has recently been shown that the dynamics of perturbed non-rotating black
holes (BHs) admits an infinite number of symmetries that are generated by the
flow of the Korteweg-de Vries (KdV) equation. These symmetries lead to an
infinite number of conserved quantities that can be obtained as integrals of
differential polynomials in the potential appearing in the gauge-invariant
master equations describing the BH perturbations, the KdV integrals. These
conserved quantities are the same for all the possible potentials, which means
that they are invariant under Darboux transformations, and they fully
determine the BHs transmission amplitudes, or greybody factors, via a moment
problem. In this paper we introduce a new semi-analytical method to obtain the
greybody factors associated with BH scattering processes by solving the moment
problem using only the KdV integrals. The method is based on the use of Padé
approximants and we check it first by comparing with results from the case of
a Pöschl-Teller potential, for which we have analytical expressions for the
greybody factors. Then, we apply it to the case of a Schwarzschild BH and
compare with results from computations based on the Wentzel–Kramers–Brillouin
(WKB) approximation. It turns out that the new method provides accurate
results for the BH greybody factors for all frequencies. The method is also
computationally very efficient.

##  I Introduction

Black Holes (BHs) are one of the most radical predictions of the General
Theory of Relativity [1, 2, 3]. They represent the most extreme gravitational
systems and are defined by the existence of an event horizon that acts as a
membrane causally separating the spacetime in the sense that no physical
signals can get out from the interior of the horizon, which constitutes a
past-future asymmetry [4]. Classical BHs have been widely studied and there
are many textbooks describing their main properties and physical consequences
[5, 6, 7, 8]. On the other hand, there is accumulating evidence of the
existence of astrophysical systems that are compatible with the general
relativistic model of a BH. Actually, in some cases the BH model of General
Relativity turns out to be the most conservative description in the sense that
alternative models require the introduction of exotic forms of matter for
which there is no observational evidence. Observations of BH candidates are
being made using different types of astronomical techniques, from radio
observations to the more recent observations by ground-based gravitational-
wave detectors [9, 10, 11], including X-rays (see, e.g. [12]). Binary BHs are
probably the main source of gravitational waves for current detectors, but
also for future ones like third-generation detectors [13, 14] and space-based
detectors like LISA [15, 16, 17, 18].

BHs are becoming central objects for a number of research areas within
astrophysics and cosmology, and also in fundamental physics, where progress
towards theories that encompass quantum effects and relativistic gravitation
use the study of physical phenomena around BHs that goes beyond the current
established knowledge. This is partly a consequence of the fact that BHs can
have any mass, from the smallest to the largest scales. The only thing we need
for their existence is a viable physical mechanism for their formation.

Scattering processes around BHs, together with quasinormal mode oscillations
of BHs, are one of main physical processes that can be described using
relativistic perturbation theory of BHs (BHPT). Particles and/or fields
scattered by BHs are of great relevance for astrophysics and fundamental
physics as they can provide us with key information about the properties of
BHs. For this reason, there is a number of studies about BH scattering (see,
e.g. [19, 20, 21, 22, 23, 24, 25, 26, 27]), both within classical general
relativity and also in semiclassical gravity, where we can gain significant
insight into quantum gravitational phenomena.

In the case of non-rotating BHs, the perturbations can be expanded in scalar,
vector, and tensor spherical harmonics in such a way that the perturbative
equations for each harmonic mode, and for each of the two possible parities,
decouple from each other. Moreover, for each mode, the perturbative equations
decouple so that the perturbations are governed by a set of master wave-like
equations with a potential that tells us the response of the BH to excitations
[28, 29, 30, 31, 32, 33, 34] (see also Ref. [8]). In the case of scattering
processes, the transmission and reflection probabilities can be obtained from
the analysis of the master equations. Some low- and high-frequency limits for
the greybody factors were obtained in [35, 36, 37, 20] and in [21, 22, 23]
respectively. Other methods to obtain greybody factors are:
Wentzel–Kramers–Brillouin (WKB) approximations [38, 39, 40] (see also [41, 42,
43]); bounding of Bogoliubov coefficients [44, 45]; monodromy techniques [46,
26, 47], etc.

Recently [48], we adopted a different point of view to the computation of BH
greybody factors, which profits from two previous studies, one on the
structure of the space of master functions and master equations [49], and the
other one on the symmetries that one can introduce in that space [50]. The
picture that emerged from these studies is that the space of possible master
functions and equations is much bigger than what it was previously known,
where we can distinguish two different branches of equations. One branch
corresponds to master equations with the well-known potentials: the Regge-
Wheeler [28] (odd-parity perturbations) and Zerilli [32] (even-parity
perturbations) potentials. The second branch contains an infinite number of
master equations with new potentials. In [50] it was shown that the space of
master functions admits different classes of symmetries. On the one hand, all
the master equations and functions are connected via Darboux transformations,
which shows that all the master functions have the same spectral properties,
in particular the same set of quasinormal modes and the same reflection and
transmission coefficients. On the other hand, one can also introduce
symmetries using the techniques developed for inverse scattering and
integrable systems [51, 52, 53, 54, 55, 56]. We showed [50] that the master
equations admit an infinite number of symmetries generated by the Korteweg-de
Vries (KdV) equation [57] and the associated infinite hierarchy of non-linear
partial differential equations (PDEs). These symmetries translate into
conservation laws from which we can obtain an infinite set of conserved
quantities [58, 59, 60], the KdV integrals. Moreover, these symmetries can be
seen as isospectral deformations of the master equations that preserve both
the scattering transmission coefficient and the quasi-normal modes. Finally,
the KdV integrals are the same for all the possible master equations, or in
other words, they are invariant under Darboux transformations.

Building on these ingredients, in [48] it was shown that the BH greybody
factors are completely determined by the KdV integrals associated with the BH
potential barrier. More specifically, the KdV integrals determine, up to a
multiplicative factor, the moments of a (probability) distribution function
associated with the transmission probability coefficient. This relationship
between moments and distribution is usually known as a moment problem, the
problem of finding the distribution from the moments. It was also argued in
[48] that the BH moment problem, where the KdV integrals are determined by the
BH potential, is determinate in the sense that a solution exists and it is
unique. Some other theoretical aspects of this problem in the context of BH
perturbations were also discussed.

Executive Summary. In this paper, we consider the BH moment problem formulated
in [48] and develop computational techniques based on Padé approximants to
solve the problem. We first apply these techniques to a case that can be
solved analytically, the case of a Pöschl-Teller potential [61]. Although this
potential does not belong to the class of potentials describing BH
perturbations, it shares some properties with the BH potential barrier that
make it a good touchstone to test methods to describe BH perturbations. We
compare our results with exact results of the Pöschl-Teller potential and show
that they coincide to a high degree of approximation. We then apply the
technique to the case of the Regge-Wheeler potential and show some error
estimation indicators. In this sense, it is important to mention that all the
possible BH potentials are physically equivalent in regard to the description
of scattering processes (the S matrix), and hence any potential would provide
the same results. The error indicators show that our technique to solve the BH
moment problem gives accurate results and is very competitive as compared with
other techniques, also from the point of view of computational cost.

Structure of the Paper. In Section II we briefly review BHPT for the case of a
Schwarzschild BH, introducing the perturbative master equations. In Section
III we review the important role of KdV symmetries of the master equations
describing the BH perturbations and how the associated KdV integrals uniquely
determine the BH greybody factors via a moment problem. In Section IV, we
introduce a new method of computing the greybody factors by inverting the BH
moment problem using Padé approximants. In Sec. V we discuss the performance
of our Padé-based method by comparing results in the case of a Pöschl-Teller
potential [61], for which we have exact expressions for the greybody factors.
In Sec. VI we apply our Padé-based approximation to perturbations of
Schwarzschild BHs, where in all the computations we use the Regge-Wheeler
potential, but the results are valid for any potential describing
perturbations of a Schwarzschild BH. In Section VII we summarize the main
results of the paper and discuss further applications of our method and
possible avenues for future developments. The paper contains four appendices:
In Appendix A we introduce some basic elements of Padé approximants; in
Appendix B we quote results on the computation of BH greybody factors using
the WKB approximation; in Appendix C we show some results on the Pöschl-Teller
potential; and in Appendix D we list the first KdV integrals for the Pöschl-
Teller potential.

Throughout this paper, otherwise stated, we use geometric units in which
$G=c=1\,$.

##  II Master Equations for BH Perturbations

The scattering of test fields, including the gravitational field, and
particles by a BH, and also quasinormal oscillations of BHs, are physical
processes that can be described by BHPT. In the case of a non-rotating BH, the
static background is the Schwarzschild metric

$d{s}{}^{2}={g}_{\mu\nu}dx^{\mu}dx^{\nu}=-f(r)\,dt^{2}+\frac{dr^{2}}{f(r)}+r^{2}d\Omega^{2}\,,$
(1)

in which

$f(r)=1-\frac{r_{s}}{r}\,,$ (2)

where $r_{s}$ is the location of the event horizon, the Schwarzschild radius:
$r_{s}=2GM/c^{2}=2M\,$. Perturbations around a Schwarzschild BH can be
decomposed into scalar, vector, and tensor spherical harmonics thanks to the
spherical symmetry of the background. The perturbative Einstein equations
decouple for each harmonic $(\ell,m)$ and for each parity 111A given harmonic
component, ${\cal O}^{\ell m}$, is said to be of the even-parity type if,
under a parity transformation $(\theta,\phi)$ $\rightarrow$
$(\pi-\theta,\phi+\pi)$, transforms as ${\cal O}^{\ell
m}\rightarrow(-1)^{\ell}{\cal O}^{\ell m}$, while it is said to be of the odd-
parity type when it transforms as ${\cal O}^{\ell
m}\rightarrow(-1)^{\ell+1}{\cal O}^{\ell m}$.. The result is a set of ten
linear but coupled PDEs for the ten metric perturbations of each harmonic and
parity. Gauge choices can help to simplify the problem, but the most important
realization is that we can find master functions, $\Psi^{\rm even/odd}_{\ell
m}$, i.e. linear combinations of the metric perturbations and their first-
order derivatives, that decouple the system of PDEs for the perturbations. In
other words, these master functions satisfy master equations where no other
combinations of perturbations appear. It turns out that these master equations
are wave-type equations (see, e.g. [62, 63, 49]):

$\left(-\frac{\partial^{2}}{\partial t^{2}}+\frac{\partial^{2}}{\partial
x^{2}}-V^{\rm even/odd}_{\ell}\right)\Psi^{\rm even/odd}_{\ell m}=0\,,$ (3)

where: $x$ is the tortoise coordinate, defined by $dx/dr=1/f\,$; $V^{\rm
even/odd}_{\ell}(r)$ is a potential constructed from the Schwarzschild
background and depends on the harmonic number $\ell$ and the parity; the
master functions $\Psi^{\rm even/odd}_{\ell m}(t,r)$ depend on both harmonic
numbers, the parity and the coordinates of the time-radial sector of the
metric, which constitutes a true Lorentzian metric 222In what follows, for the
sake of simplicity, we drop the harmonic numbers $(\ell,m)$ from the master
functions..

In [49], all the possible master equations with master functions linear in the
metric perturbations and their first-order derivatives were determined. In
this space of possible master functions and equations one can distinguish two
branches: The standard branch, which is characterized by having the known
potentials, namely the Regge-Wheeler potential [28] for odd-parity
perturbations

$V_{\ell}^{\rm
RW}(r)=\left(1-\frac{r_{s}}{r}\right)\left(\frac{\ell(\ell+1)}{r^{2}}-\frac{3r_{s}}{r^{3}}\right)\,,$
(4)

and the Zerilli potential [32] for even-parity perturbations

$\displaystyle V^{\rm Z}_{\ell}(r)$ $\displaystyle=$
$\displaystyle\frac{f}{\lambda^{2}}\left[\frac{(\ell-1)^{2}(\ell+2)^{2}}{r^{2}}\left(\ell(\ell+1)+\frac{3r_{s}}{r}\right)\right.$
(5) $\displaystyle+$
$\displaystyle\left.\frac{9r^{2}_{s}}{r^{4}}\left((\ell-1)(\ell+2)+\frac{r_{s}}{r}\right)\right]\,,$

where $\lambda$ is a function of $r$ given by

$\displaystyle\lambda(r)=(\ell-1)(\ell+2)+\frac{3r_{s}}{r}\,.$ (6)

Regarding master functions, it was found [49] that in the odd-parity case the
most general master function is a general linear combination of the Regge-
Wheeler [28] and the Cunningham-Price-Moncrief [29, 30, 31] master functions.
For the even-parity sector, it was found that the most general master function
is a general linear combination of the Zerilli-Moncrief master function [32,
34], and another master function that appears to be new [49].

The second branch of master functions and equations, named the Darboux branch,
is described in [49]. It was shown that it contains an infinite number of
different master equations of the form of Eq. (3). As a consequence, there is
an infinite family of possible potentials that have to satisfy a non-linear
ordinary differential equation. A remarkable aspect of this branch is that the
master functions depend explicitly on the potential. Indeed, for each
potential, the master functions are written in terms of the metric
perturbations and an integral containing the potential (see [49] for the
explicit expressions).

The structure of this infinite landscape of master equations and functions was
investigated in [50]. It was shown that all the pairs $(V,\Psi)$ are connected
by Darboux transformations (DTs). From a more physical point of view, DTs are
isospectral transformations, which means that they preserve the spectrum of
the frequency domain operator associated with the master equation, as well as
the reflection and transmission coefficients [64, 50]. Indeed, let us consider
single frequency solutions, i.e. let us take the following form of the master
function: $\Psi(t,r)=e^{ikt}\psi(x;k)$. Then, the master equation (3) becomes
a time-independent Schrödinger equation of the form

$\psi_{,xx}-V\psi=-k^{2}\psi\,.$ (7)

It was shown in [50] (see also [48]) that DTs map a time-independent
Schrödinger equation to a physically equivalent one, with different potential
barrier and master function, but with the same spectrum.

##  III BH Greybody factors and the Korteweg-de Vries Integrals

We have mentioned that the different possible master equations describing the
first-order perturbations of non-rotating BHs are connected by DTs, which turn
out to be isospectral transformations. In addition, there is another type of
transformations of the master equations that are isospectral. They consist in
deformations [50, 48] of the time-independent Schrödinger equation that follow
the Korteweg-de Vries (KdV) equation and the associated infinite hierarchy of
KdV equations. The KdV deformations consist in introducing, in the time-
independent Schrödinger equation (7), a dependence on a parameter $\tau$, that
is

$\psi(x)\rightarrow\psi(\tau,x)\,,\;\;V(x)\rightarrow
V(\tau,x)\,,\;\;k\rightarrow k(\tau)\,,$ (8)

in such a way that the potential $V(\tau.x)$ follows the KdV equation

$V_{,\tau}-6VV_{,x}+V_{,xxx}=0\,,$ (9)

which is a non-linear PDE. It is precisely the fact that the deformation of
the potential follows the KdV equation that makes the spectrum conserved [50]
(see also [48]) under the KdV flow, that is: $(k^{2})_{,\tau}=0\,$. This
result holds for the continuous spectrum, bound-states and resonances (or
quasinormal modes), although in the case of the Schwarzschild BHs, where the
potential is positive everywhere, there is no discrete spectrum.

The KdV deformations constitute symmetries of our problem. These symmetries
lead to conservation laws which, in turn, lead to conserved quantities [58,
59, 60], the so-called KdV integrals. In [48], we have showed different paths
to these conserved quantities. A particularly interesting approach is the
Hamiltonian formulation of the KdV equation, which was first studied by
Gardner [65]. Soon afterwards, Zakharov and Faddeev [59] showed that the
results by Gardner and collaborators [52, 65] can be seen from the point of
view of action-angle variables. The action-angle variables associated to the
Hamiltonian formulation of the KdV equation appear naturally when we look at
the scattering problem associated with the time-independent Schrödinger
equation. In this way, Zakharov and Faddeev [59] were able to show that the
KdV equation constitutes a completely integrable Hamiltonian system. Following
this line of thought, we can introduce an infinite hierarchy of KdV evolution
equations, associated with the infinite chain of KdV conservation laws, using
the Hamiltonian formulation of the KdV equation (see also [66]). The conserved
quantities, i.e. the KdV integrals, can be obtained in terms of integrals of
differential polynomials of the potential $V$ (that is, polynomials in $V$ and
its derivatives). The explicit expressions of the infinite series of KdV
integrals is given by (see, e.g. [59]):

$\mathcal{K}_{n}=\int^{\infty}_{-\infty}dx\,\kappa_{n}(x)\,.$ (10)

where the densities $\kappa_{n}(x)$ are obtained from the following recurrence
relation

$\displaystyle\kappa_{1}(x)$ $\displaystyle=$ $\displaystyle V(x)\,,$ (11)
$\displaystyle\kappa_{n}(x)$ $\displaystyle=$
$\displaystyle-\frac{d}{dx}\kappa_{n-1}(x)-\sum_{k=1}^{n-1}\kappa_{n-k-1}(x)\kappa_{k}(x)\,,$
(12) $\displaystyle(n=2,\ldots)\,.$

In summary, the potential of a time-independent Schrödinger equation (7) has
associated an infinite sequence of conserved quantities, the KdV integrals,
generated by the KdV flows. We have also seen that BH perturbations admit an
infinite number of descriptions in terms of master functions and equations
[50]. The obvious question is how the KdV integrals change for the different
master equations (for fixed $\ell$). The answer to this question is very
simple: The KdV integrals are the same for all the possible potentials that
appear in the master equations describing BH perturbations. In other words,
the KdV integrals are invariant under DTs [50] (see also [48] for more
details).

This interplay between BHPT, inverse scattering theory, and integrable
systems, in particular the key role of DTs and the KdV equation, is what has
motivated the research program initiated in [50], and which has led to the
result [48] that the KdV integrals determine completely the BH greybody
factors. This is what has been called the BH moment problem, that is, the
problem of finding the BH greybody factors from the KdV integrals.

BH Greybody factors are defined in the context of BH scattering processes, and
they are nothing but the modulus square of the transmission coefficient
through the BH potential barrier. To see this in more detail, let us consider
the scattering of a wave coming from $x\rightarrow-\infty$ (the BH horizon) by
the BH potential barrier. After interacting with the potential barrier, part
of the wave is transmitted, goes to $x\rightarrow\infty$, and part is
reflected, going back to $x\rightarrow-\infty$. In mathemtical terms, this can
expressed in the following way:

$\displaystyle\psi(x,k)=\left\\{\begin{array}[]{lc}a(k)e^{ikx}+b(k)e^{-ikx}&\mbox{for~{}}x\to-\infty\,,\\\\[14.22636pt]
e^{ikx}&\mbox{for~{}}x\to\infty\,,\end{array}\right.$ (15)

The coefficients $a(k)$ and $b(k)$ are usually known as the Bogoliubov
coefficients. They completely determine the BH scattering matrix as well as
the reflection and transmission coefficients as follows

$\displaystyle t(k)=\frac{1}{a(k)}\,,\quad r(k)=\frac{b(k)}{a(k)}\,.$ (16)

The transmission and reflection probabilities, or greybody factors, are given
by the modulus square of the corresponding coefficients, i.e.

$T(k)=|t(k)|^{2}\,,\quad R(k)=|r(k)|^{2}\,.$ (17)

For real $k$ they satisfy the unitarity condition

$\displaystyle T(k)+R(k)=1\,.$ (18)

In Ref. [48], it was shown that the KdV integrals and the BH greybody factors
are related by an infinite set of integral relations, the so-called trace
identities [59], which play an important role in the spectral theory of the
time-independent Schrödinger equations (see, e.g. [67] and references
therein). These identities can be written in the following form

$\displaystyle\mu_{2j}=\int_{-\infty}^{\infty}dk\,k^{2j}\,p(k)\,.$ (19)

Here, the $\mu_{j}$ are positive constants that coincide with the standard
definition of the moments of the (distribution) function $p(k)$. The
remarkable fact is that that they are proportional to the KdV integrals as
follows

$\displaystyle\mu_{2j}=(-1)^{j}\,\frac{\mathcal{K}_{2j+1}}{2^{2j+1}}\,,$ (20)

and the distribution function $p(k)$ only depends on the transmission
amplitude $T(k)$

$p(k)=-\frac{\ln T(k)}{2\pi}\,,$ (21)

where the minus sign guarantees that we are expressing everything in terms of
positive quantities. Therefore, the problem of finding the greybody factors of
a given potential barrier (without bound states) is equivalent to a moment
problem [68, 69, 70].

The moment problem is a long standing problem in mathematics and it appears in
a large variety of physical situations. Roughly speaking, the moment problem
consists in inverting equation (19) to find the probability distribution
$p(k)$ (only) from the knowledge of the moments $\mu_{i}$. In our case, the
moments are proportional to the KdV integrals associated with the BH
potential, and the distribution is a logarithmic function of the greybody
factors [see Eqs. (20) and (21)].

To be more precise, the moment problem we are facing is known in the
mathematical literature as a Hamburger moment problem (see [48] for more
details), since the interval in which the distribution function $p(k)$ is
defined corresponds to the whole real line, and the the moments are computed
as integrals over it [see Eq. (20)]. Furthermore, Eq. (19) defines a special
case of the Hamburger moment problem, namely the symmetric Hamburger moment
problem [70], as the probability distribution $p(k)$ appears to be a symmetric
function of $k$, that is $p(-k)=p(k)$, so that all the odd moments vanish:
$\mu_{2j+1}=0\,$. This is an interesting property of the integral equations
(19) whose nature can be related to the fact that only the odd integrals (10)
correspond to the true first integrals of the KdV equation (see [48] for more
details), that is, to the first integrals that arise from the Hamiltonian
formulation of the KdV equation.

The symmetric Hamburger moment problem can be uniquely mapped into an
associated Stieltjes moment problem (there is actually a bijective
correspondence [70] between the two moment problems). The Stieljes moment
problem has the particularity that the distribution function $p(k)$ is defined
in the positive real line. Indeed, we can rewrite Eq. (19) in the following
form

$\displaystyle\mu_{2j}=2\int_{0}^{\infty}\\!\\!\\!dk\,k^{2j}p(k)\,.$ (22)

We can now introduce the following change of variable:

$\xi=\sigma^{2}k^{2}\,,$ (23)

where $\sigma$ is constant with dimensions of length that makes $\xi$ to be a
dimensionless quantity. Then, we can write

$\displaystyle\mu_{2j}=\frac{1}{\sigma^{2j+1}}\int_{0}^{\infty}\\!\\!\\!d\xi\,\xi^{j}\tilde{p}(\xi)\equiv\tilde{\mu}_{j}\,,$
(24)

where

$\tilde{p}(\xi)=\frac{p(\xi)}{\sqrt{\xi}}=\frac{p(k)}{\sigma\,k}\,.$ (25)

The associated Stieltjes distribution has non-vanishing odd and even moments
$\tilde{\mu}_{j}$, both defined by the even moments $\mu_{2j}$ of the
Hamburger moment problem. We can also introduce dimensionless moments from Eq.
(24) in the following way

$\displaystyle
m_{j}\equiv\int_{0}^{\infty}d\xi\,\xi^{j}\,\tilde{p}(\xi)=\sigma^{2j+1}\tilde{\mu}_{j}\,,$
(26)

which provides a formulation of the moment problem in terms of dimensionless
quantities only. Finally, it is worth mentioning that it is possible to show
[70] that when the symmetric Hamburger moment problem (20) has a unique
solution, so does the associated Stieltjes problem of Eq. (24).

##  IV Solving the Moment Problem using Padé Approximants

The connection between the scattering trace formulas and the moment problem
opens a new avenue to develop novel methods to calculate the transmission (and
reflection) coefficients of any potential barrier that does not have bound
states (as in the case of positive potential barriers as it happens for non-
rotating BHs). In Ref. [48], the moment problem approach has been used to show
that the KdV integrals, which are proportional to the moments, completely and
uniquely determine the greybody factors of the BH potential barrier. In other
words, the BH moment problem associated to the calculation of the BH greybody
factors has a unique solution.

Moreover, in Ref. [48] it was described how, when analytical expressions are
available, the moment problem can be solved with the use of the Stieltjes-
Perron inversion formula [71] (see also [68, 69, 70]). This was illustrated by
considering the particular case of a Pöschl-Teller potential [61].

However, since it is quite unrealistic to think that we can have analytic
expressions for a general potential, as we do for the Pöschl-Teller case, we
present here a general approach to construct approximations to the BH greybody
factors. In this sense, it is important to remark that in the literature there
is no a unique general solution method for the moment problem, but various
approximate and/or numerical approaches that have been explored during the
last century like, for instance: The method of moments and its generalizations
[72, 73]; kernel density functions approximation [74]; maximum entropy methods
[75, 76] (see also [77, 78] for some reviews); etc.

Since we are interested in exploiting the connection between KdV integrals and
greybody factors in a semi-analytical way, we choose to follow the approach of
Ref. [79], where the authors show how to find a probability distribution
function starting from the asymptotic form of the corresponding moment
generating function (MGF) in combination with the use of Padé approximants
(see also Ref. [80]). This choice is further motivated by the particular
features of the Padé approximants to a Stieltjes series. The reason is that
for Stieltjes series the convergence of diagonal and subdiagonal Padé
approximants is guaranteed [80, 81]. In what follows we use this method to
obtain semi-analytical expressions for the BH greybody factors in terms only
of the KdV integrals associated to the BH potential barrier.

The procedure we follow consists in determining the Padé approximants
corresponding to the asymptotic series of the MGF and then to take the inverse
Laplace transform. We first illustrate how this can be done in practice and
then we apply it to the well-known Pöschl-Teller potential barrier, in Sec. V,
where we compare our results with exact expressions. Once we show that the
method works properly we apply it, in Sec. VI, to the BH potential barrier,
the potential that appears in the master equation for perturbations of the
Schwarzschild metric.

Let us start with the formulation of the Stieltjes moment problem in terms of
dimensionless quantities [see Eq. (26)]. We have computed analytically a few
KdV integrals both for the Pöschl-Teller potential [61] and also for the
Regge-Wheeler potential (and hence for any potential in any of the master
equations describing Schwarzschild perturbations). The ones for the Pöschl-
Teller potential are given in Appendix D while the ones for the Regge-Wheeler
potential are given in Appendix E of Ref. [48]. The constant $\sigma$, in the
case of the Pöschl-Teller potential can be taken to be $\sigma=1/\alpha$,
while for the Regge-Wheeler potential we can choose it simply as
$\sigma=r_{s}\,$.

The MGF for the moments $m_{j}$ corresponding to the distribution $\tilde{p}$
is usually defined as the Laplace transform of the (probability) distribution
[82]

$M(t)=\int_{0}^{\infty}d\xi\,e^{-t\xi}\,\tilde{p}(\xi)\equiv\mathcal{L}(\tilde{p})\,,$
(27)

where $\mathcal{L}$ denotes the Laplace transform operator. By expanding the
exponential inside the integral in a Taylor series in $t$, we get the
asymptotic behavior of $M(t)$

$M(t)=\sum_{n=0}^{\infty}\frac{m_{n}}{n!}(-t)^{n}\,.$ (28)

Notice that the $n$-th moment can be found as the $n$-th derivative of $M(t)$
at $t=0$, i.e.

$m_{n}=(-1)^{n}M^{(n)}(0)\,,$ (29)

which is the reason why this function is called the MGF. The series in Eq.
(28) is a formal expression that does not necessarily need to converge for our
procedure to work [80, 81].

The step by step procedure to find the probability distribution starting from
its moments that we are going to use in this work follows the method proposed
in Ref. [79] (see also [80]). This method consists in two basic steps: First,
to construct the Padé approximants for the asymptotic expansion of the MGF
[Eq. (28)]. And second, to take the Laplace inverse of the Padé approximation
to find $\tilde{p}(\xi)$.

As is well known, the Padé approximants are approximations by means of
rational functions to a power series, in such a way that the order of accuracy
depends on the degree of the polynomials in the numerator and denominator of
the rational function (see Appendix A for the basic ingredients of Padé
approximants that are relevant for this work). Then, the Padé approximants,
for the series of $M(t)$ in Eq. (28), of order $K+L$, where $K$ and $L$ are
positive integers, are written as follows

$\left[K/L\right](t)=\frac{P(t)}{Q(t)}=\sum_{n=0}^{K+L}\frac{m_{n}}{n!}(-t)^{n}\,.$
(30)

where $P(t)$ and $Q(t)$ are polynomials of order $K$ and $L$ respectively,
that is,

$P(t)=\sum_{n=0}^{K}P_{n}t^{n}\,,\quad Q(t)=\sum_{n=0}^{L}Q_{n}t^{n}\,.$ (31)

As expected [see Eqs. (57) and (58) in Appendix A], all the coefficients of
the polynomials $P_{n}$ and $Q_{n}$ are completely determined by the KdV
integrals. In this work we only consider the diagonal ($K=L$) and sub-diagonal
($K=L-1$) Padé approximants because they play a particular role in the
convergence of the Padé sequence as, for example, they provide upper and lower
bounds to any Stieltjes series and if the problem is determinate they converge
to the unique solution [81]. In order to put ourselves in a situation that
would allow us to perform easily the Laplace inversion of Eq. (27), we need to
first find the poles of the Padé approximants, $t=-t_{i}\,$. In this way, we
can decompose the Padé approximants in partial fractions, that is, we can
write them in the form

$\left[K/L\right](t)=\sum_{i=1}^{L}\frac{\lambda_{i}}{t+t_{i}}\,.$ (32)

Here $K=L,L-1$ and the $\lambda_{i}$’s are the residues of the poles, i.e.

$\lambda_{i}=\frac{P(-t_{i})}{Q^{\prime}(-t_{i})}\,,$ (33)

where here the prime denotes differentiation with respect to $t\,$. The poles
of the Padé approximants are expected to be located on the negative half-plane
in $t\in\mathbb{C}\,$, therefore we expect $\Re(t_{i})>0\,$ [80]. Then, as a
general rule, whenever we encounter a Padé approximant for which
$\Re(t_{i})<0$ for some $i$ (a pole in the positive half-plane), we should
discard it as an unacceptable approximation since it clearly has the wrong
analytic structure [80]. Moreover, since the distribution function is real, we
must have either real poles or complex ones coming in pairs together with the
complex conjugate.

The last step is to carry out the inverse Laplace transform of the Padé
approximant. This is a simple task since we know that (see, e.g. [83])

$\mathcal{L}^{-1}\left(\frac{1}{t+t_{i}}\right)=e^{-t_{i}\,\xi}\,.$ (34)

Therefore, we have the following approximated expression for the distribution
function associated with the BH greybody factors

$\tilde{p}(\xi)\simeq\sum_{i=1}^{L}\lambda_{i}\,e^{-t_{i}\,\xi}\,,$ (35)

from where it is straightforward to obtain an expression for the transmission
probability $T(k)$. Indeed, taking into account Eqs. (25) and (35) we get

$p(k)\simeq\sigma\,k\;\sum_{i=1}^{L}\lambda_{i}\,e^{-t_{i}\,\sigma^{2}k^{2}}\,.$
(36)

Finally, using Eq. (21) we obtain the following approximate form for $T(k)$

$\displaystyle T(k)$ $\displaystyle\simeq$
$\displaystyle\exp{\left(-2\,\pi\,\sigma\,k\;\sum_{i=1}^{L}\lambda_{i}e^{-t_{i}\,\sigma^{2}k^{2}}\right)}$
(37) $\displaystyle\equiv$ $\displaystyle T_{[K/L]}(k)\,.$

From these expressions it becomes clear that if we allow for poles with
positive real part, the solution would be unphysical as it would either
describe an infinitely growing greybody factor or decaying one, depending on
the sign of the residue. The expression in Eq. (37) constitutes a semi-
analytic approximation for the transmission probability associated with a
potential barrier, a BH barrier in our context, exclusively in terms of their
KdV integrals. Indeed, all the coefficients in the approximation are evaluated
solely from the KdV integrals (10). To show this explicitly we evaluate the
first two Padé approximants. The first is the Padé approximant $[0/1]$:

$[0/1](t)=\frac{P_{0}}{1+Q_{1}\,t}\,,$ (38)

It depends only on the first two moments as the coefficients $P_{0}$ and
$Q_{1}$ are given by

$P_{0}=m_{0}\,,\quad Q_{1}=\frac{m_{1}}{m_{0}}\,.$ (39)

There is a single real pole, which is given by

$t_{1}=\frac{m_{0}}{m_{1}}\,,\quad\lambda_{1}=\frac{m_{0}^{2}}{m_{1}}\,,$ (40)

where $\lambda_{1}$ is the associated residue. Therefore, the greybody factor
(37) has the following simple expression

$T_{[0/1]}(k)=\exp{\left(-2\pi\frac{m^{2}_{0}}{m_{1}}\sigma\,k\,e^{-\frac{m_{0}}{m_{1}}\sigma^{2}k^{2}}\right)}\,.$
(41)

The second Padé approximant is $[1/1]$:

$[1/1](t)=\frac{P_{0}+P_{1}\,t}{1+Q_{1}\,t}\,,$ (42)

which contains only the first three moments

$P_{0}=m_{0}\,,\quad P_{1}=\frac{m_{0}m_{2}-2\,m_{1}^{2}}{2m_{1}}\,,\quad
Q_{1}=\frac{m_{2}}{2\,m_{1}}\,.$ (43)

Again, there is a single pole which, together with the corresponding residue,
are given by

$t_{1}=\frac{2m_{1}}{m_{2}}\,,\quad\lambda_{1}=\frac{4m^{3}_{1}}{m^{2}_{2}}\,,$
(44)

and the greybody factor (37) adopts the following expression

$\displaystyle T_{[1/1]}(k)$ $\displaystyle=$
$\displaystyle\exp\left(-2\pi\frac{4m^{3}_{1}}{m^{2}_{2}}\sigma\,k\,e^{-\frac{2m_{1}}{m_{2}}\sigma^{2}k^{2}}\right)\,.$
(45)

These expressions do not provide the most accurate results but may be helpful
in cases in which a simple analytical expression is needed.

We can summarize the method we just have shown for constructing approximations
for the BH greybody factors in a schematic way as an algorithmic list of
steps:

  1. 1.

Evaluate the first $n$ KdV integrals for the chosen potential.

  2. 2.

Obtain the moments from the KdV integrals and construct the MGF by using the
expansion (28) at order $n$.

  3. 3.

Construct Padé approximants of order $[K/L]$, with $K+L\leq n$.

  4. 4.

Evaluate the poles and residues [see Eq. (33)] of the Padé approximants and
discard those with poles that have positive real part.

  5. 5.

Apply the Laplace inversion formula (34) to finally obtain the approximations
[see Eq. (37)] for the greybody factors.

##  V Greybody Factors for the Pöschl-Teller Potential

In order to test our method for the computation of greybody factors, it is
convenient to compare our results with a case in which we have analytical
expressions for them. In the case of BH perturbations we do not have exactly
solvable potentials, and for this reason we consider the well-known case of
the Pöschl-Teller potential [61]:

$\displaystyle V_{\rm PT}(x)=\frac{U_{0}}{\cosh^{2}{(\alpha
x)}}=\frac{\alpha^{2}(\beta^{2}+\frac{1}{4})}{\cosh^{2}(\alpha x)}\,,$ (46)

where $\beta$ is a dimensionless constant given by

$\beta=\sqrt{\frac{U_{0}}{\alpha^{2}}-\frac{1}{4}}>0\,.$ (47)

This is one of the few cases of a potential, defined on the whole real line,
in which the time-independent Schrödinger equation (7) is exactly solvable
[84] (see Appendix C for a summary of the main relevant results). The greybody
factors for the Pöschl-Teller potential are given by Eq. (70), i.e.

$\displaystyle T_{\rm PT}=\frac{\sinh^{2}\left(\frac{\pi
k}{\alpha}\right)}{\cosh^{2}\left(\frac{\pi
k}{\alpha}\right)+\sinh^{2}\left(\pi\beta\right)}\,.$ (48)

Exact solvability, together with the similarity in shape333We can adjust the
Pöschl-Teller potential parameters to have a potential barrier that looks
quite similar to the Schwarzschild BH potential barrier. Nevertheless, both
potentials differ in the decay to zero near spatial infinity, which has
important consequences for the time-dependent response to perturbations: In
the BH case the long-term behavior is dominated by power law tails [85] in
contrast with the Pöschl-Teller case, where there are no tails. with the BH
potential barrier, makes it a perfect playground to get deeper insight on the
BH scattering (see e.g. Ref. [86, 87]). In this section, we exploit this
comparison to test our approximation based on Padé approximants.

Let us now apply the method we described in Sec. IV (the steps involved are
outlined at the end of the section) to find approximations for the greybody
factors of the Pöschl-Teller potential barrier. We start with the computation
of the KdV integrals. The first ten non-zero integrals for the the Pöschl-
Teller potential are given in Appendix D. From these integrals, we find the
moments, construct the Padé approximants together with their poles and
residues, and finally apply the inverse Laplace transform. If we consider the
first six non-zero KdV integrals, the $T_{[2/3]}(k)$ approximation to the
transmission coefficient is given by (we set the Pöschl-Teller parameter
$\beta=5$):

$\displaystyle T_{[2/3]}$ $\displaystyle\simeq$
$\displaystyle\exp\left\\{-\frac{2\pi
k}{\alpha}\left[14.895\,e^{-0.256\frac{k^{2}}{\alpha^{2}}}\sin\left(0.093\,\frac{k^{2}}{\alpha^{2}}\right)\right.\right.$
(49) $\displaystyle-$ $\displaystyle
5.586\,e^{-0.256\,\frac{k^{2}}{\alpha^{2}}}\cos{\left(0.093\,\frac{k^{2}}{\alpha^{2}}\right)}$
$\displaystyle+$
$\displaystyle\left.\left.19.167\,{e^{-0.743\,\frac{k^{2}}{\alpha^{2}}}}\right]\right\\}\,,$

where the trigonometric functions appear due to the presence of pairs of
complex conjugate poles. A similar structure is shared by the other
$T_{[K/L]}$ functions. Different comparisons of the approximations obtained
with our method with the exact solution for $T(k)$ [see Eq. (48) and Appendix
C] are given in Fig. 1 for different values of the parameter $\beta$. As we
can see, the approximation (37) proves to be highly accurate, even when we use
low-order Padé approximants. Indeed, we can see in Fig. 1 that the
approximations based on the first Padé approximants are almost
indistinguishable from the exact solution.

Figure 1: Greybody factors for the Pöschl-Teller potential barrier [61] with
respect to the wave number $k$ obtained from the Padé approximation of Eq.
(37) and compared with the exact solution in Eq. (48) for $\beta=1.8$ (top
panel) and $\beta=5$ (bottom panel). The notation $T[N/M]$ means that the
greybody factor is computed from the $[N/M]$ Padé approximant to the MGF
asymptotic series.

The main drawback of the semi-analytic approximation represented by Eq. (37)
is the not so good behavior at small frequencies (small $k$). However, as one
can see in Fig. 2, the region in which the behavior deviates from the expected
one is pushed towards the origin as we include more KdV integrals in the
approximation. Or in other words, the region were the approximation gets worse
is shrinked towards the origin as we use high-order Padé approximants. This
suggests that the approximation given by Eq. (37) converges also at low
frequencies in the limit in which we include infinite moments (an infinite
number of KdV integrals).

Figure 2: Plots of the behavior of the greybody factors for the Pöschl-Teller
potential in the low-frequency (low $k$) regions obtained by zooming near
$k=0$ in the plots of Fig. 1 for both $\beta=1.8$ (top panel) and $\beta=5$
(bottom panel).

In order to obtain a better quantitative understanding of the behavior of our
approximation technique, we need to find appropriate error estimates. The fact
that we are using the Pöschl-Teller potential to test our method allows us to
make direct comparisons with the exact greybody factor in Eq. (48) and have a
better insight of the convergence properties of the procedure. In this sense,
let us consider the following $k$-dependent error function

$\delta T_{\rm PT}(k)=\left|T_{\rm PT}(k)-T_{[K/L]}(k)\right|\,,$ (50)

that is, the absolute difference between the exact solution and our
approximation. Plots of this error estimate are shown in Figs. 3 and 4 for
different values of the parameter $\beta$. As we can see from these plots,
both diagonal and subdiagonal Padé approximants to the MGF seem to produce a
converging approximation to the exact value. Apart from the region near the
origin $k=0$, the error estimate in the intermediate region appears to be
slightly worse than in the rest of the $k$ line. Nevertheless, we still obtain
a very good degree of precision in that intermediate region, and this improves
with the number of KdV integrals that we include in the approximation.

Regarding the error coming from the low-frequency behavior (low $k$), the
results shown in Figs. 3 and 4 clarify that the deviations become significant
only in a small area close to the origin $k=0$. The size of this region
decreases as we increase the order of the approximation (or the number of KdV
integrals used), so that we can say that we obtain an accurate description
also at small frequencies (small $k$).

The comparison between Figs. 3 and 4 further shows that the higher the
potential barrier is the better is the approximation of Eq. (37) to the exact
value $T_{\rm PT}$, both at low and high frequencies, while it provides
slightly worse results in the intermediate region (see also Figs. 1 and 2). On
the other hand, it is worth mentioning that due to the shape of $T_{\rm
PT}(k)$ and $T_{[K/L]}(k)$, the relative error would not provide a fully
meaningful estimate. Indeed, since $T_{\rm PT}$ goes to zero when $k$ does,
the contributions from small frequencies to the relative error grow near $k=0$
and hence, we would obtain a misleading measure of the error with respect to
the exact solution.

Figure 3: Plots of the error estimate $\delta T_{\rm PT}(k)$ [see Eq. (50)],
in logarithmic scale of base 10, for the diagonal (top panel) and sub-diagonal
(bottom panel) Padé-based approximation of Eq. (37) for the greybody factors
with $\beta=1.8\,$.

Figure 4: Plots of the error estimate $\delta T_{\rm PT}(k)$ [see Eq. (50)],
in logarithmic scale of base 10, for the diagonal (top panel) and sub-diagonal
(bottom panel) Padé-based approximation of Eq. (37) for the greybody factors
with $\beta=5\,$.

The integration of $\delta T_{\rm PT}(k)$ over $k$ provides a global error
indicator that can also inform us about the improvement of the approximation.
Then, we define

$\Delta T_{\rm PT}=\int_{k_{0}}^{k_{\infty}}dk\,\delta T_{\rm PT}(k)\,,$ (51)

where $k_{0}$ is a cut-off at low frequencies and $k_{\infty}$ at high
frequencies. The low-frequency cut-off $k_{0}$ is introduced to separate the
very low-frequency contributions, where our approximations may misbehave, but
without changing the global qualitative behavior. The high-frequency cut-off
$k_{\infty}$ is introduced instead for computational convenience, since the
greybody factors become almost constant for large $k$ and hence, integrating
all the way to $k\rightarrow\infty$ may produce meaningless results from the
numerical point of view. In Fig. 5 we show results of the computation of
$\Delta T_{\rm PT}$, where we can see that, irrespective of the details of the
approximate greybody factors [see Eq. (37)] at each range of frequencies, the
convergence of both diagonal and sub-diagonal Padé approximants is guaranteed.

Figure 5: Plots of the error $\Delta T_{\rm PT}$ [Eq. (51)] as computed from
diagonal (blue line) and sub-diagonal (black line) Padé approximations with
cut-offs $(k_{0}\,,\,k_{\infty})=(0.25\,,\,20)$ and with Pöschl-Teller
parameter $\beta=1.8$ (top panel) and $\beta=5$ (bottom panel).

The diagonal and subdiagonal Padé approximants provide two different kinds of
approximations for the greybody factors in terms of convergence of our Padé-
based approximation to the real function (48). Each of them follow its own
pattern of convergence. This can clearly be seen when studying Stieltjes
sequences [80, 81], as already mentioned in Sec. IV. However, for our
purposes, it can be useful to study the improvements in the approximation when
considering the addition of just one more KdV integral. This turns out to be
crucial when the exact solution is not provided, as in the relevant case of a
BH potential barrier (see Sec. VI). To that end, we can introduce an error
estimator for successive approximations as follows

$\delta T_{[K/L]}(k)=\left|T_{[(K+d)/K]}(k)-T_{[K/L]}(k)\right|\,,$ (52)

where either $L=K$ and $d=-1$ or $L=K+1$ and $d=0$. This is a measure of how
much the Padé-based approximation of Eq. (37) improves every time that we
include in the analysis one more moment/KdV integral. The results are shown in
Fig. 6, where we can see how the convergence of the approximation with $n+1$
moments improves up to more than one order of magnitude with respect to the
one calculated from $n$ moments.

Figure 6: Plots of the error estimate in Eq. (53) for successive
approximations to the greybody factor of the Pöschl-Teller potential, in
logarithmic scale of base 10, for $\beta=1.8$ (top panel) and for $\beta=5$
(bottom panel).

We can introduce a global error estimate associated with the error estimate in
Eq. (52) for succesive Padé-base approximations [see Eq. (37)] by integrating
this error function over $k$:

$\Delta T_{[K/L]}=\int_{k_{0}}^{k_{\infty}}dk\,\delta T_{[K/L]}(k)\,,$ (53)

where we have introduced again the cut-offs $(k_{0},k_{\infty})$ for the same
reasons described above. This global error estimate is shown in Fig. 7. By
looking at Figs. 6 and 7 we see that $\delta T_{[K/L]}(k)$ and $\Delta
T_{[K/L]}$ have the same qualitative behavior as $\delta T_{\rm PT}(k)$ and
$\Delta T_{\rm PT}$, thus the same conclusions can be drawn here. There error
estimates are going to be particularly useful in Sec. VI, where we deal with
the Regge-Wheeler potential.

Figure 7: Plots for the global error estimator $\Delta T_{[K/L]}$ [see Eq.
(53)] for successive Padé-based approximations of the Pöschl-Teller gredybody
factors with increasing number of moments/KdV integrals, with cut-offs
$(k_{0}\,,\,k_{\infty})=(0.25\,,\,20)\,$, for $\beta=1.8$ (top panel) and
$\beta=5$ (bottom panel).

It is also worth mentioning that from these results it appears that the cases
in which we have Padé approximants that present positive poles, these poles do
not completely spoil the convergence. Instead, they just appear to slow down
the convergence rate of the approximation.

In addition to the comparisons with the exact expressions that we have for the
Pöschl-Teller potential and between successive approximations, there is a
different type of test of our approximation method that we can carry out,
namely the comparison with WKB approximations for the computations of greybody
factors [38, 39, 40, 43]. The expressions of the WKB approximations that we
use in this work are shown in Appendix B. The results of the comparisons with
third-order WKB approximations (the fourth order happens to vanish for the
Pöschl-Teller potential) are shown in Fig. 8. It is clear from these plots
that the accuracy of the WKB approximation depends on the parameter $\beta$,
which controls the strength of the Pöschl-Teller potential. Indeed, it is well
known that the WKB approximation works better for higher values of the angular
momentum in the case of Regge-Wheeler potential [38, 39, 40, 43]. Since the
higher the angular momentum is the higher the maximum of the potential barrier
becomes, we can apply the same argument to the case of the Pöschl-Teller
potential in terms of the parameter $\beta$. In this sense, our Padé-based
approximation seems to be less sensitive to the potential strength.

On the other hand, our approach differs from the WKB approximation in the way
that we build more accurate estimations of the greybody factors. In the case
of the WKB method, it is an accumulative process in which we compute higher
orders of an approximation series. The size of the expression for each order
of approximation in the WKB method grows considerably as we can see in the
Appendix B. In our approach, we can start from an arbitrary number of KdV
integrals (moments) and compute the desired order of approximation from the
associated Padé approximant. The computational cost of the KdV integrals is
quite low as we can see from their explicit expressions in Appendix D (see
also Appendix E in Ref. [48]). Moreover, the computational cost of our Padé-
based approximation is quite affordable since the evaluation of high-order
Padé approximants only require a few seconds using a computer algebra system
like Maxima [88] on a personal laptop computer.

Finally, it is interesting to note that both approximations, our Padé-based
approximation and the WKB approximation, fail in a similar way in describing
the intermediate frequency region. In fact, in Fig. 1 we see that, for high
barriers (in the case of the figure $\beta=5$), our Padé -based approximation,
from low to high frequencies, grows up to a maximum slightly greater than one,
in the intermediate region, before settling down to a constant value equal to
one. The same happens in the case of the WKB approximation in the low barrier
case as one can see by looking closely at Fig. 8 (in the case of the figure
$\beta=1.8$), and it becomes more evident for lower barriers (see also Figs. 9
and 13 for the Regge-Wheeler case in the next section.

Figure 8: Plots of the comparison of computations of the greybody factors for
the Pöschl-Teller potential using the Padé-baed approximation of Eq. (37), the
WKB approximation method of Eq. (59), and the exact values in Eq. (70) for
$\beta=1.8$ (top panel) and $\beta=5$ (bottom panel).

##  VI Greybody Factors for Schwarzschild Black Holes

The greybody factors associated with gravitational perturbations around a
Schwarzschild BH background can be evaluated in an analogous way as we have
done in the previous section for the case of a Pöschl-Teller potential [61].
We just need to follow the step-by-step procedure outlined at the end of Sec.
IV. There is however a significant difference between the two potentials,
beyond exact solvability, and it has to do with the decay towards infinity of
the potentials. Indeed, while both potentials decay exponentially as we go
towards $x\rightarrow-\infty$ (the BH horizon in the case of the BH potential
barrier), the approach towards $x\rightarrow\infty$ (spatial infinity in the
case of the BH potential barrier) is very different. In the case of the
Pöschl-Teller potential, it also decays exponentially, while in the case of
the BH potential barrier, take the Regge-Wheeler potential as an example, it
decays quadratically, that is $V_{\rm RW}\sim L/r^{2}\sim L/x^{2}$, where
$L=\ell(\ell+1)$. As we can see, the so-called angular momentum barrier
dominates at $x\rightarrow\infty$ in the BH caes. The different decay towards
$x\rightarrow\infty$ has important consequences for the properties of the
solutions of the time-independent Schrödinger equation (7), and even more in
the time-domain. Indeed, it turns out that in the case of a non-rotating BH,
the structure of the Green function present a branch cut that is responsible
for the power-law tails that we observe in the time-evolution of BH
perturbations (see [85] for more details), and which are not present in the
case of the Pöschl-Teller potential in Eq. (46). Nevertheless, for the
computation of the KdV integrals that we need for the greybody factors, both
potentials decay sufficiently fast so that the corresponding integrals are
well defined (finite).

Returning to our procedure to compute the BH greybody factors, the first step
is to find the KdV integrals. To that end, due to the Darboux invariance of
the KdV integrals (see Refs. [50, 48] for details on this symmetry of the
space of master functions and equations), we can use any of the potentials
that appear in the master equations that describe the perturbations of a
Schwarzschild BH, for instance, the Regge-Wheeler potential (4). The first ten
non-vanishing KdV integrals for this potential can be found in Appendix E of
Ref. [48]. The second step is to build the asymptotic series for the MGF [Eq.
(27)]. For this we need the moments, which are proportional to the KdV
integrals according to Eq. (19). The third step is to construct the Padé
approximant to the MGF series to the desired degree of precision. The order of
the Padé approximation, $(K,L)$ [see Eq. (37)] is directly connected to the
number of KdV integrals (moments), $n=K+L$, that we use in the computation.
Then, the fourth step consists in finding the poles and residues of the Padé
approximant. This is done numerically and in our case we use the computer
algebra system Maxima [88]. This leads to the semi-analytical approximation of
Eq. (32), in terms of partial fractions, for the Padé approximant. The last
step is to perform the Laplace inversion, as in Eq. (36), to finally obtain
the result represented by Eq. (37). In Fig. 9 we show the results for the
transmission coefficient using different Padé approximants and for two values
of the harmonic number $\ell$: $\ell=2$ (top panel) and $\ell=10$ (bottom
panel).

Figure 9: Plots of the greybody factors for the Regge-Wheeler potential
barrier calculated using the Padé approximation of Eq. (37) for $\ell=2$ (top
panel) and $\ell=10$ (bottom panel). The notation $T[N/M]$ refers to the fact
that the greybody factor is calculated from the $[N/M]$ Padé approximant to
the MGF asymptotic series of Eq. (28).

The behavior of the Padé approximation for the greybody factors shown in Fig.
9 appears to be similar to what we found for the Pöschl-Teller potential in
the previous section. Actually, we observe that the higher the harmonic number
$\ell$ (angular momentum), and hence the higher the strength of the potential
is, the higher is the number of moments that we need to include in order to
obtain a good accuracy for intermediate frequencies (intermediate values of
$k$).

On the other hand, in Fig. 10 we show the behavior of the Regge-Wheeler
greybody factor in the low-frequency (low $k$) region (by zooming into the
plots of Fig. 9), and in Fig. 11 we plot the error $\delta T_{[K/L]}$ for
successive Padé-based approximations, introduced in Eq. (52). By looking at
these plots, we can notice again that the higher the harmonic number $\ell$
is, the better the approximation near the origin becomes. We can also see that
the general behavior of the approximation of Eq. 37 is that it improves as we
increase the number of moments involved in the computation, as expected.
Therefore, following the Pöschl-Teller case as a guidance, our approximation
seems to be very good as we can see by looking at the order of magnitude of
$\delta T_{[K/L]}$ in Fig. 11

Figure 10: The plot shows the behavior of the Regge-Wheeler potential greybody
factor in the low-frequency (low $k$) regions obtained by zooming near $k=0$
in plots of Fig. 9 for $\ell=2$ (top panel) and $\ell=10$ (bottom panel).

Figure 11: Plot of the Regge-Wheeler greybody factor error $\delta
T_{[K/L]}(k)$ [Eq. (52)] for various orders of the Padé approximant used, or
equivalently, for different numbers of the KdV integrals involved in the
computation. The top panel shows plots for $\ell=2$ and the bottom panel for
$\ell=10\,$.

In Fig. 12 we show two plots (for two values of the harmonic number
$\ell=2,10$) for the global error estimate (53) for the Regge-Wheeler greybody
factor. Using this error estimate we have a unique measure of the improvement
of the approximation with the increase of KdV integrals/moments. Fig. 12 shows
the improvement in the approximation as we increase the order of the Padé
approximants used in our computations.

Figure 12: Plot of the Regge-Wheeler greybody factor global error $\Delta
T_{[K/L]}$ [Eq. (53)] for various orders of the Padé approximant used, or
equivalently, for different numbers of the KdV integrals involved in the
computation. The top panel shows plots for $\ell=2$ and cutoffs
$(k_{0}\,,\,k_{\infty})=(0.05\,,\,20)$, while the bottom panel shows plots for
$\ell=10$ and cutoffs $(k_{0}\,,\,k_{\infty})=(0.2\,,\,20)\,$.

In Fig. 13 we plot our Padé-based approximation against the fourth-order WKB
approximation (see Appendix B for details). The comparison shows that the two
approximations become almost indistinguishable for high $\ell$ but differ more
significantly for lower values of $\ell$. Fig. 13 also shows that by
considering more moments (KdV integrals) in the construction of the Padé-based
approximation (37) we can reach very good accuracy (see also Fig. 9).

To summarize, we have found that the Padé-based approximation we have
introduced in Eq. (37) provides accurate estimates for the greybody factors
associated with gravitational scattering processes around a Schwarzschild BH.
The accuracy of the greybody factors computed in this way is comparable to the
one obtained from the WKB approximation scheme, provided we use sufficient KdV
integrals. In this sense, it is important to mention that while increasing the
WKB order of approximation implies a rapid growth of the number of terms that
appear at each order [42] (see also Appendix B), the increase in the Padé
order (or in other words, in the number of KdV integrals/moments used in the
computation) appears to be quite simple and requires a relative quite low
increase in computational cost. Therefore, our Padé-based approximation can be
pushed to very high-order without neither modifying the procedure nor
increasing significantly the computational cost associated with it.

Figure 13: Comparison between the Regge-Wheeler greybody factors as computed
from the Padé-based approximation with those computed using the fourth-order
WKB approximation scheme. The top planel shows the $\ell=2$ case while the
bottom panel shows the $\ell=10$ case.

##  VII Conclusions and Discussion

In this paper we have introduced a new method to compute the greybody factors
of a potential barrier, in particular of the gravitational Schwarzschild
potential barrier. The method is based on the main result of Ref. [48], which
tells us that the greybody factors are uniquely determined by the KdV
integrals associated with the BH potential barrier via a moment problem, what
we called the BH moment problem. This result originates from the study of the
space of all the possible master functions and equations describing BH
perturbations that was carried out in Refs. [50, 49]. In those works it was
found that all the master equations are connected by Darboux Transformations,
which preserve the continuum spectrum (the greybody factors) and the
quasinormal spectrum (the resonances), showing their physical equivalence. An
additional set of symmetries was found, consisting in deformations that follow
the hierarchy of KdV equations and generate an infinite sequence of
conservation laws with their associated conserved quantities, the KdV
integrals (they are integrals of differential polynomials of the potential
that characterize the master equation). It turns out that the KdV integrals
are the same for all the possible potentials (for all the possible master
equations), that is, they are invariant under Darboux Transformations. In this
sense, the KdV integrals characterize the physics behind the master equations,
as we have seen with the case of greybody factors and their connection with
the KdV integrals via a moment problem.

There is not a unique way of solving/inverting the moment problem to obtain
the greybody factors once the KdV integrals are given. In this work, we have
made the choice of using Padé approximants mainly because of their intimate
connection to the moment problem, which goes back to the first studies on
continued fractions [71]. Actually, an important results in this line is that
the Padé approximants of a Stieltjes series converge to the corresponding
exact Stieltjes function when the problem is determinate [81, 80], i.e. it
admits a unique solution (see also [48]). Moreover, the Padé approximants have
been also proposed as a tool to carry out the Laplace inversion [89, 90], so
that the moment problem can be solved by inverting the MGF [79], which is in
itself a Laplace transform, to find the distribution function associated with
the moments. As a final remark on the use of Padé approximants, it is
important to mention that they have been widely used and investigated in the
literature and, as a result of this, a number of useful results are known. For
instance, the use of Shanks transformations, which improves their convergence
rate [81], and/or the multipoint Padé approximants (see, e.g. [80]) to improve
the approximation near zero frequency.

The results of this paper have to be considered as a first step towards
solving the BH moment problem to find the BH greybody factors in this new
approach to the problem, introduced in Ref. [48]. Despite the good precision
we have achieved in the computation of the transmission coefficient, there is
still significant room for improvement in the way we use Padé approximants, or
even for introducing more sophisticated techniques, which may trigger the
exploration of better methods and techniques to solve the BH moment problem.
For instance, one possibility would be to look for alternative methods to
invert the Laplace transform [90, 80]. Leaving aside the approach we use in
this paper of introducing the MGF as a main object in the computation, a
number of numerical or semi-analytical methods have been studied to invert the
moment problem in other contexts (see, for instance [72, 73, 74, 75, 76, 77,
78]).

Finally, in this paper, we have applied our techniques to the case of
scattering of gravitational perturbations off the potential barrier of a non-
rotating BH. By looking at the details of the whole procedure, it seems
reasonable to think that with little modifications, the BH moment problem can
be adapted to other physical situations of interest, in particular to the case
of BH perturbations of different character (spin): scalar, vector and neutrino
perturbations. Other relevant scenarios where it would be worth applying these
techniques include: Perturbations of spinning BHs in General Relativity;
perturbations of compact objects, in particular of those exotic cases that can
mimic BHs; perturbations of BHs in alternative theories of gravity;
perturbations of BHs in higher dimensions; etc.

###### Acknowledgements.

ML and CFS are supported by contract
PID2019-106515GB-I00/AEI/10.13039/501100011033 (Spanish Ministry of Science
and Innovation) and 2017-SGR-1469 (AGAUR, Generalitat de Catalunya). ML is
also supported by Juan de la Cierva contract FJC2021-047289-I funded by
program MCIN/AEI/10.13039/501100011033 (Spanish Ministry of Science and
Innovation) and by NextGenerationEU/PRTR (European Union). This work was also
partially supported by the program Unidad de Excelencia María de Maeztu
CEX2020-001058-M (Spanish Ministry of Science and Innovation). We have used
the Computer Algebra System Maxima [88] to carry out most of the computations
necessary for this paper.

##  Appendix A Basics of Padé Approximants

Padé approximants (see Ref. [80]) are a powerful way of reproducing a power
series, in a such a way that the Padé approximant exhibits in general better
convergence properties. Let us consider the following power series:

$f(z)=\sum_{i=0}^{\infty}c_{i}z^{i}\,.$ (54)

The Padé approximants are a sequence of rational functions of the form

$\left[K/L\right](z)=\frac{\sum_{i=0}^{K}P_{i}z^{i}}{\sum_{j=0}^{L}Q_{j}z^{j}}\,,$
(55)

such that each term of the sequence is equal to the power series expansion up
to the order $K+L+1$, where the $(P_{i},Q_{j})$ are just constant
coefficients. The coefficient $Q_{0}$ in particular can be fixed to one
without loosing generality. The rest of the coefficients are found by matching
the expansion of $\left[K/L\right](z)$ up to order $K+L+1$ with the first
$K+L+1$ coefficients of Eq. (54), i..e.

$\left(\sum_{i=0}^{K+L}c_{i}z^{i}\right)\left(\sum_{j=0}^{L}Q_{j}z^{j}\right)-\sum_{i=0}^{K}P_{i}z^{i}=O\left(z^{K+L+1}\right)\,.$
(56)

Then, the coefficients in the denominator of Eq. (55) are related to the
coefficients of the power series in the following way

$\sum_{j=1}^{L}C_{ij}Q_{j}=c_{K+i}\,,\quad 1\leq j\leq L\,,$ (57)

where $C_{ij}=c_{K+i-j}$ are the element of a $L\times L$ matrix. The
coefficients in the numerator are instead found as follows

$P_{n}=\sum_{j=0}^{n}c_{n-j}Q_{j}\,,\quad 0\leq j\leq K\,,$ (58)

where $Q_{j}=0$ if $j>L$.

##  Appendix B WKB Approximations for BH Greybody Factors

The application of the WKB approximation to the study of perturbations of BHs
[38] provides useful and accurate expressions for the estimation of
quasinormal mode frequencies and greybody factors. In the particular case of
greybody factors, the result provided by the WKB approximation can be written
as follows [39]:

$T(k)=\frac{1}{\left|1+e^{2\pi i\left(\nu(k)+1/2\right)}\right|}\,,$ (59)

where

$\nu(k)+1/2=i\frac{\left(k^{2}-V_{0}\right)}{\sqrt{-2V_{0}^{{}^{\prime\prime}}}}-\sum_{i=2}\Lambda_{i}\,,$
(60)

and

$V_{0}=V(x_{0})\,,\quad
V^{(n)}_{0}=\left.\frac{d^{n}V}{dx^{n}}\right|_{x=x_{0}}\,,$ (61)

where $x_{0}$ is the location of the maximum of the potential barrier $V(x)$,
that is, $x_{0}$ is defined by the following relation

$V^{\prime}(x_{0})=0\,.$ (62)

The quantities $\Lambda_{i}$ for $i=1,\ldots,13$ can be found in Refs. [38,
39, 40, 41, 42]. In particular, the quantities $\Lambda_{2}\,$,
$\Lambda_{3}\,$, and $\Lambda_{4}$, the ones used in this work, are given by

$\displaystyle\Lambda_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{-2V_{0}^{{}^{\prime\prime}}}}\left[\frac{1}{8}\left(\frac{V_{0}^{(4)}}{V_{0}^{{}^{\prime\prime}}}\right)\left(\frac{1}{4}+\gamma^{2}\right)-\frac{1}{288}\left(\frac{V_{0}^{(3)}}{V_{0}^{{}^{\prime\prime}}}\right)^{2}\left(7+60\gamma^{2}\right)\right]\,,$
(63) $\displaystyle\Lambda_{3}$ $\displaystyle=$
$\displaystyle\frac{\gamma}{2V_{0}^{{}^{\prime\prime}}}\bigg{[}\frac{5}{6912}\left(\frac{V_{0}^{(3)}}{V_{0}^{{}^{\prime\prime}}}\right)^{4}\left(77+188\gamma^{2}\right)-\frac{1}{384}\left(\frac{V_{0}^{(3)\,2}V_{0}^{(4)}}{V_{0}^{{}^{\prime\prime}\,3}}\right)\left(51+100\gamma^{2}\right)+\frac{1}{2304}\left(\frac{V_{0}^{(4)}}{V_{0}^{{}^{\prime\prime}}}\right)^{2}\left(67+68\gamma^{2}\right)$
(64) $\displaystyle+$
$\displaystyle\frac{1}{288}\left(\frac{V_{0}^{(3)}V_{0}^{(5)}}{V_{0}^{{}^{\prime\prime}\,2}}\right)\left(19+28\gamma^{2}\right)-\frac{1}{288}\left(\frac{V_{0}^{(6)}}{V_{0}^{{}^{\prime\prime}}}\right)\left(5+4\gamma^{2}\right)\bigg{]}\,,$
$\displaystyle\Lambda_{4}$ $\displaystyle=$
$\displaystyle\frac{1}{597196800\sqrt{2}V_{0}^{{}^{\prime\prime}\,7}\sqrt{V_{0}^{{}^{\prime\prime}}}}\left\\{2536975V_{0}^{(3)\,6}-9886275V_{0}^{{}^{\prime\prime}}V_{0}^{(3)\,4}V_{0}^{(4)}+5319720V_{0}^{{}^{\prime\prime}\,2}V_{0}^{(3)\,3}V_{0}^{(5)}-225V_{0}^{{}^{\prime\prime}\,2}V_{0}^{(3)\,2}\right.\times$
(65) $\displaystyle\times$
$\displaystyle\left(-40261V_{0}^{(4)\,2}+9688V_{0}^{{}^{\prime\prime}}V_{0}^{(6)}\right)+3240V_{0}^{{}^{\prime\prime}\,3}V_{0}^{(3)}\left(-1889V_{0}^{(4)}V_{0}^{(5)}+220V_{0}^{{}^{\prime\prime}}V_{0}^{(7)}\right)-729V_{0}^{{}^{\prime\prime}\,3}\left[1425V_{0}^{(4)\,3}\right.$
$\displaystyle-$
$\displaystyle\left.\left.1400V_{0}^{{}^{\prime\prime}}V_{0}^{(4)}V_{0}^{(6)}+8V_{0}^{{}^{\prime\prime}}\left(-123V_{0}^{(5)\,2}+25V_{0}^{{}^{\prime\prime}}V_{0}^{(8)}\right)\right]\right\\}$
$\displaystyle+$
$\displaystyle\frac{\gamma^{2}}{4976640\sqrt{2}V_{0}^{{}^{\prime\prime}\,7}\sqrt{V_{0}^{{}^{\prime\prime}}}}\left\\{348425V_{0}^{(3)\,6}-1199925V_{0}^{{}^{\prime\prime}}V_{0}^{(3)\,4}V_{0}^{(4)}+57276V_{0}^{{}^{\prime\prime}\,2}V_{0}^{(3)\,3}V_{0}^{(5)}-45V_{0}^{{}^{\prime\prime}\,2}V_{0}^{(3)\,2}\times\right.$
$\displaystyle\times$
$\displaystyle\left(-20671V_{0}^{(4)\,2}+4552V_{0}^{{}^{\prime\prime}}V_{0}^{(6)}\right)+1080V_{0}^{{}^{\prime\prime}\,3}V_{0}^{(3)}\left(-489V_{0}^{(4)}V_{0}^{(5)}+52V_{0}^{{}^{\prime\prime}}V_{0}^{(7)}\right)-27V_{0}^{{}^{\prime\prime}\,3}\left[2845V_{0}^{(4)\,3}\right.$
$\displaystyle-$
$\displaystyle\left.\left.2360V_{0}^{{}^{\prime\prime}}V_{0}^{(4)}V_{0}^{(6)}+56V_{0}^{{}^{\prime\prime}}\left(-31V_{0}^{(5)\,2}+5V_{0}^{(2)}V_{0}^{(8)}\right)\right]\right\\}$
$\displaystyle+$
$\displaystyle\frac{\gamma^{4}}{2488320\sqrt{2}V_{0}^{{}^{\prime\prime}\,7}\sqrt{V_{0}^{{}^{\prime\prime}}}}\left\\{192925V_{0}^{(3)\,6}-581625V_{0}^{{}^{\prime\prime}}V_{0}^{(3)\,4}V_{0}^{(4)}+234360V_{0}^{{}^{\prime\prime}\,2}V_{0}^{(3)\,3}V_{0}^{(5)}-45V_{0}^{{}^{\prime\prime}\,2}V_{0}^{(3)\,2}\times\right.$
$\displaystyle\times$
$\displaystyle\left(-8315V_{0}^{(4)\,2}+1448V_{0}^{{}^{\prime\prime}}V_{0}^{(6)}\right)+1080V_{0}^{{}^{\prime\prime}\,3}V_{0}^{(3)}\left(-161V_{0}^{(4)}V_{0}^{(5)}+12V_{0}^{{}^{\prime\prime}}V_{0}^{(7)}\right)-27V_{0}^{{}^{\prime\prime}\,3}\left[625V_{0}^{(4)\,3}\right.$
$\displaystyle-$
$\displaystyle\left.\left.440V_{0}^{{}^{\prime\prime}}V_{0}^{(4)}V_{0}^{(6)}+8V_{0}^{{}^{\prime\prime}}\left(-63V_{0}^{(5)\,2}+5V_{0}^{(2)}V_{0}^{(8)}\right)\right]\right\\}\,.$

where $\gamma=n+1/2$ and $n=0,1,2...$ is the overtone number.

##  Appendix C Some Formulas for the Pöschl-Teller Potential

Let us consider the time-independent Schrödinger equation (7) for the case of
a Pöschl-Teller potential barrier [see Eq. (46). The general solution reads
(Refs. [84, 91, 48]):

$\displaystyle\psi(x,k)$ $\displaystyle=$ $\displaystyle
A\,2^{ik/\alpha}\left[1-\tanh^{2}(\alpha
x)\right]^{-ik/2\alpha}{}_{2}F_{1}\left(\lambda-\frac{ik}{\alpha},1-\lambda-\frac{ik}{\alpha},1-\frac{ik}{\alpha};\frac{1-\tanh(\alpha
x)}{2}\right)$ (66) $\displaystyle+$ $\displaystyle B\,\left[1-\tanh(\alpha
x)\right]^{ik/2\alpha}\left[1+\tanh(\alpha
x)\right]^{-ik/2\alpha}{}_{2}F_{1}\left(\lambda,1-\lambda,1+\frac{ik}{\alpha};\frac{1-\tanh(\alpha
x)}{2}\right)\,,$

where

$\lambda=i\beta+\frac{1}{2}\,.$ (67)

By analyzing the asymptotic behavior of the general solution, both at
$x\rightarrow\infty$ and at $x\rightarrow-\infty$, we can find analytical
expressions for the reflection and transmission coefficients. First, one can
obtain the following expressions for the Bogoliubov coefficients

$\displaystyle a(k)$ $\displaystyle=$
$\displaystyle\frac{\Gamma\left(1-ik\right)\Gamma\left(-ik\right)}{\Gamma\left(\frac{1}{2}-i(k-\beta)\right)\Gamma\left(\frac{1}{2}-i(k+\beta)\right)}\,,$
(68) $\displaystyle b(k)$ $\displaystyle=$
$\displaystyle\frac{\Gamma\left(1-ik\right)\Gamma\left(ik\right)}{\Gamma\left(\frac{1}{2}+i\beta)\Gamma(\frac{1}{2}-i\beta\right)}\,,$
(69)

where for the sake of simplicity we set $\alpha=1$. Then, we can find the
transmission and reflection coefficients from Eq. (16). The transmission
probability can be evaluated by taking the inverse modulus square of Eq. (68)
for real $k$ and reads

$\displaystyle T(k)=\frac{\sinh^{2}\left(\pi k\right)}{\cosh^{2}\left(\pi
k\right)+\sinh^{2}\left(\pi\beta\right)}\,,$ (70)

where we used the properties of the gamma functions [92].

##  Appendix D KdV Integrals for the Pöschl-Teller Potential

Here we list the first non-vanishing KdV integrals [see Eq. (10)] for the
Pöschl-Teller potential [see Eq. (46)]:

$\displaystyle I_{1}$ $\displaystyle=$
$\displaystyle\frac{\alpha}{2}\left(4\beta^{2}+1\right)\,,$ (71)
$\displaystyle I_{3}$ $\displaystyle=$
$\displaystyle-\frac{\alpha^{3}}{12}\left(4\beta^{2}+1\right)^{2}\,,$ (72)
$\displaystyle I_{5}$ $\displaystyle=$
$\displaystyle\alpha^{5}\left(4\beta^{2}+1\right)^{2}\left(\frac{2}{15}\beta^{2}-\frac{1}{6}\right)\,,$
(73) $\displaystyle I_{7}$ $\displaystyle=$
$\displaystyle-\alpha^{7}\left(4\beta^{2}+1\right)^{2}\left(\frac{2}{7}\beta^{4}-\frac{41}{105}\beta^{2}-\frac{577}{840}\right)\,,$
(74) $\displaystyle I_{9}$ $\displaystyle=$
$\displaystyle\alpha^{9}\left(4\beta^{2}+1\right)^{2}\left(\frac{32}{45}\beta^{6}-\frac{232}{315}\beta^{4}-\frac{1406}{315}\beta^{2}-\frac{2683}{630}\right)\,,$
(75) $\displaystyle I_{11}$ $\displaystyle=$
$\displaystyle-\alpha^{11}\left(4\beta^{2}+1\right)^{2}\left(\frac{64}{33}\beta^{8}-\frac{64}{99}\beta^{6}-\frac{12968}{693}\beta^{4}-\frac{168668}{3465}\beta^{2}-\frac{61767}{1540}\right)\,,$
(76) $\displaystyle I_{13}$ $\displaystyle=$
$\displaystyle\alpha^{13}\left(4\beta^{2}+1\right)^{2}\left(\frac{512}{91}\beta^{10}+\frac{12160}{3003}\beta^{8}-\frac{956096}{15015}\beta^{6}-\frac{1123952}{3465}\beta^{4}-\frac{2879438}{4095}\beta^{2}-\frac{49126459}{90090}\right)\,,$
(77) $\displaystyle I_{15}$ $\displaystyle=$
$\displaystyle-\alpha^{15}\left(4\beta^{2}+1\right)^{2}\left(\frac{256}{15}\beta^{12}+\frac{1408}{39}\beta^{10}-\frac{1057744}{6435}\beta^{8}-\frac{826544}{495}\beta^{6}-\frac{303716683}{45045}\beta^{4}-\frac{406249013}{30030}\beta^{2}\right.$
(78) $\displaystyle\left.-\frac{813135229}{80080}\right)\,,$ $\displaystyle
I_{17}$ $\displaystyle=$
$\displaystyle\alpha^{17}\left(4\beta^{2}+1\right)^{2}\left(\frac{8192}{153}\beta^{14}+\frac{157696}{765}\beta^{12}-\frac{436736}{3315}\beta^{10}-\frac{755645056}{109395}\beta^{8}-\frac{116271200}{2431}\beta^{6}-\frac{134553980488}{765765}\beta^{4}\right.$
(79)
$\displaystyle\left.-\frac{51719984938}{153153}\beta^{2}-\frac{126841684787}{510510}\right)\,,$
$\displaystyle I_{19}$ $\displaystyle=$
$\displaystyle-\alpha^{19}\left(4\beta^{2}+1\right)^{2}\left(\frac{16384}{95}\beta^{16}+\frac{4947968}{4845}\beta^{14}+\frac{10694656}{4845}\beta^{12}-\frac{2947115008}{146965}\beta^{10}-\frac{119163316352}{440895}\beta^{8}\right.$
(80)
$\displaystyle\left.-\frac{720103571584}{440895}\beta^{6}-\frac{115858335632}{20349}\beta^{4}+\frac{14090850111224}{1322685}\beta^{2}-\frac{8185650267425}{1058148}\right)\,.$

## References

  * Einstein [1915a] A. Einstein, Zur Allgemeinen Relativitätstheorie, Sitzungsber. Preuss. Akad. Wiss. Berlin (Math. Phys.) 1915, 778 (1915a), [Addendum: Sitzungsber.Preuss.Akad.Wiss.Berlin (Math.Phys.) 1915, 799–801 (1915)].
  * Einstein [1915b] A. Einstein, The Field Equations of Gravitation, Sitzungsber. Preuss. Akad. Wiss. Berlin (Math. Phys.) 1915, 844 (1915b).
  * Einstein [1916] A. Einstein, The Foundation of the General Theory of Relativity, Annalen Phys. 49, 769 (1916).
  * Finkelstein [1958] D. Finkelstein, Past-Future Asymmetry of the Gravitational Field of a Point Particle, Phys. Rev. 110, 965 (1958).
  * Hawking and Ellis [2011] S. W. Hawking and G. F. R. Ellis,  _The Large Scale Structure of Space-Time_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2011).
  * Wald [1984] R. Wald,  _General Relativity_ (The University of Chicago Press, Chicago, 1984).
  * Novikov and Frolov [1989] I. D. Novikov and V. P. Frolov,  _Physics of Black Holes_ (Kluwer Academic, Dordrecht, Netherlands, 1989).
  * Chandrasekhar [1992] S. Chandrasekhar,  _The Mathematical Theory of Black Holes_ (Oxford University Press, New York, 1992).
  * Abbott  _et al._ [2019] B. P. Abbott _et al._ (LIGO Scientific, Virgo), GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs, Phys. Rev. X 9, 031040 (2019), arXiv:1811.12907 [astro-ph.HE] .
  * Abbott  _et al._ [2021a] R. Abbott _et al._ (LIGO Scientific, VIRGO), GWTC-2.1: Deep Extended Catalog of Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run, (2021a), arXiv:2108.01045 [gr-qc] .
  * Abbott  _et al._ [2021b] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run, (2021b), arXiv:2111.03606 [gr-qc] .
  * Corral-Santana  _et al._ [2016] J. M. Corral-Santana, J. Casares, T. Muñoz-Darias, F. E. Bauer, I. G. Martínez-Pais, and D. M. Russell, BlackCAT: A catalogue of stellar-mass black holes in X-ray transients, Astron. Astrophys. 587, A61 (2016), arXiv:1510.08869 [astro-ph.HE] .
  * Sathyaprakash  _et al._ [2012] B. Sathyaprakash _et al._ , Scientific Objectives of Einstein Telescope, Class. Quant. Grav. 29, 124013 (2012), [Erratum: Class.Quant.Grav. 30, 079501 (2013)], arXiv:1206.0331 [gr-qc] .
  * Evans  _et al._ [2021] M. Evans _et al._ , A Horizon Study for Cosmic Explorer: Science, Observatories, and Community, (2021), arXiv:2109.09882 [astro-ph.IM] .
  * Amaro-Seoane  _et al._ [2017] P. Amaro-Seoane _et al._ (LISA), Laser Interferometer Space Antenna, (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Amaro-Seoane  _et al._ [2022] P. Amaro-Seoane _et al._ , Astrophysics with the Laser Interferometer Space Antenna, (2022), arXiv:2203.06016 [gr-qc] .
  * Auclair  _et al._ [2022] P. Auclair _et al._ (LISA Cosmology Working Group), Cosmology with the Laser Interferometer Space Antenna, (2022), arXiv:2204.05434 [astro-ph.CO] .
  * Arun  _et al._ [2022] K. G. Arun _et al._ (LISA), New horizons for fundamental physics with LISA, Living Rev. Rel. 25, 4 (2022), arXiv:2205.01597 [gr-qc] .
  * Vishveshwara [1970] C. Vishveshwara, Scattering of Gravitational Radiation by a Schwarzschild Black-hole, Nature 227, 936 (1970).
  * Starobinskii and Churilov [1974] A. A. Starobinskii and S. M. Churilov, Amplification of electromagnetic and gravitational waves scattered by a rotating ”black hole”, Sov. Phys. JETP 65, 1 (1974).
  * Sanchez [1976] N. G. Sanchez, Scattering of scalar waves from a Schwarzschild black hole, J. Math. Phys. 17, 688 (1976).
  * Sanchez [1977] N. G. Sanchez, The Wave Scattering Theory and the Absorption Problem for a Black Hole, Phys. Rev. D 16, 937 (1977).
  * Sanchez [1978] N. G. Sanchez, Elastic Scattering of Waves by a Black Hole, Phys. Rev. D 18, 1798 (1978).
  * Futterman  _et al._ [2012] J. A. H. Futterman, F. A. Handler, and R. A. Matzner,  _Scattering from Black Holes_ (Cambridge University Press, 2012).
  * Andersson and Jensen [2000] N. Andersson and B. P. Jensen, Scattering by black holes. Chapter 0.1, , 1607 (2000), arXiv:gr-qc/0011025 .
  * Castro  _et al._ [2013] A. Castro, J. M. Lapan, A. Maloney, and M. J. Rodriguez, Black Hole Scattering from Monodromy, Class. Quant. Grav. 30, 165005 (2013), arXiv:1304.3781 [hep-th] .
  * Folacci and Ould El Hadj [2019] A. Folacci and M. Ould El Hadj, Regge pole description of scattering of gravitational waves by a Schwarzschild black hole, Phys. Rev. D 100, 064009 (2019), arXiv:1906.01441 [gr-qc] .
  * Regge and Wheeler [1957] T. Regge and J. A. Wheeler, Stability of a Schwarzschild singularity, Phys. Rev. 108, 1063 (1957).
  * Cunningham  _et al._ [1978] C. T. Cunningham, R. H. Price, and V. Moncrief, Radiation from collapsing relativistic stars. I - Linearized odd-parity radiation, Astrophys. J. 224, 643 (1978).
  * Cunningham  _et al._ [1979] C. T. Cunningham, R. H. Price, and V. Moncrief, Radiation from collapsing relativistic stars. II. Linearized even parity radiation, Astrophys. J. 230, 870 (1979).
  * Cunningham  _et al._ [1980] C. T. Cunningham, R. H. Price, and V. Moncrief, Radiation from collapsing relativistic stars. III - Second order perturbations of collapse with rotation, Astrophys. J. 236, 674 (1980).
  * Zerilli [1970a] F. J. Zerilli, Gravitational Field of a Particle Falling in a Schwarzschild Geometry Analyzed in Tensor Harmonics, Phys. Rev. D 2, 2141 (1970a).
  * Zerilli [1970b] F. J. Zerilli, Effective potential for even parity Regge-Wheeler gravitational perturbation equations, Phys. Rev. Lett. 24, 737 (1970b).
  * Moncrief [1974] V. Moncrief, Gravitational perturbations of spherically symmetric systems. I. The exterior problem, Ann. Phys. (N.Y.) 88, 323 (1974).
  * Page [1976] D. N. Page, Particle Emission Rates from a Black Hole: Massless Particles from an Uncharged, Nonrotating Hole, Phys. Rev. D 13, 198 (1976).
  * Unruh [1976] W. G. Unruh, Absorption Cross-Section of Small Black Holes, Phys. Rev. D 14, 3251 (1976).
  * Starobinskii [1973] A. A. Starobinskii, Amplification of waves reflected from a rotating ”black hole”., Sov. Phys. JETP 37, 28 (1973).
  * Schutz and Will [1985] B. F. Schutz and C. M. Will, Black Hole Normal Modes: A Semianalytic Approach, Astrophys. J. Lett. 291, L33 (1985).
  * Iyer and Will [1987] S. Iyer and C. M. Will, Black Hole Normal Modes: A WKB Approach. 1. Foundations and Application of a Higher Order WKB Analysis of Potential Barrier Scattering, Phys. Rev. D 35, 3621 (1987).
  * Iyer [1987] S. Iyer, Black Hole Normal Modes: A WKB Approach. 2. Schwarzschild Black Holes, Phys. Rev. D 35, 3632 (1987).
  * Konoplya [2003] R. A. Konoplya, Quasinormal behavior of the d-dimensional Schwarzschild black hole and higher order WKB approach, Phys. Rev. D 68, 024018 (2003), arXiv:gr-qc/0303052 .
  * Matyjasek and Opala [2017] J. Matyjasek and M. Opala, Quasinormal modes of black holes. The improved semianalytic approach, Phys. Rev. D 96, 024011 (2017), arXiv:1704.00361 [gr-qc] .
  * Konoplya  _et al._ [2019] R. A. Konoplya, A. Zhidenko, and A. F. Zinhailo, Higher order WKB formula for quasinormal modes and grey-body factors: recipes for quick and accurate calculations, Class. Quant. Grav. 36, 155002 (2019), arXiv:1904.10333 [gr-qc] .
  * Visser [1999] M. Visser, Some general bounds for 1-D scattering, Phys. Rev. A 59, 427 (1999), arXiv:quant-ph/9901030 .
  * Boonserm and Visser [2008] P. Boonserm and M. Visser, Bounding the greybody factors for Schwarzschild black holes, Phys. Rev. D 78, 101502 (2008), arXiv:0806.2209 [gr-qc] .
  * Neitzke [2003] A. Neitzke, Greybody factors at large imaginary frequencies, (2003), arXiv:hep-th/0304080 .
  * Harmark  _et al._ [2010] T. Harmark, J. Natario, and R. Schiappa, Greybody Factors for d-Dimensional Black Holes, Adv. Theor. Math. Phys. 14, 727 (2010), arXiv:0708.0017 [hep-th] .
  * Lenzi and Sopuerta [2022] M. Lenzi and C. F. Sopuerta, Black Hole Greybody Factors from Korteweg-de Vries Integrals: Theory, (2022), arXiv:2212.03732 [gr-qc] .
  * Lenzi and Sopuerta [2021a] M. Lenzi and C. F. Sopuerta, Master functions and equations for perturbations of vacuum spherically symmetric spacetimes, Phys. Rev. D 104, 084053 (2021a), arXiv:2108.08668 [gr-qc] .
  * Lenzi and Sopuerta [2021b] M. Lenzi and C. F. Sopuerta, Darboux covariance: A hidden symmetry of perturbed Schwarzschild black holes, Phys. Rev. D 104, 124068 (2021b), arXiv:2109.00503 [gr-qc] .
  * Miura [1968] R. M. Miura, Korteweg-de Vries Equation and Generalizations. I. A Remarkable Explicit Nonlinear Transformation, J. Math. Phys. 9, 1202 (1968).
  * Gardner  _et al._ [1967] C. S. Gardner, J. M. Greene, M. D. Kruskal, and R. M. Miura, Method for solving the Korteweg-de Vries equation, Phys. Rev. Lett. 19, 1095 (1967).
  * Faddeev [1976] L. D. Faddeev, Inverse problem of quantum scattering theory. II, J. Sov. Math. 5, 334 (1976).
  * Deift and Trubowitz [1979] P. Deift and E. Trubowitz, Inverse scattering on the line, Commun. Pure Appl. Math. 32, 121 (1979).
  * Novikov  _et al._ [1984] S. Novikov, S. V. Manakov, L. P. Pitaevsky, and V. E. Zakharov,  _Theory of Solitons. The Inverse Scattering Method_ (Consultants Bureau, New York, 1984).
  * Ablowitz and Segur [1981] M. J. Ablowitz and H. Segur,  _Solitons and the Inverse Scattering Transform_ (Society for Industrial and Applied Mathematics, Philadelphia, 1981).
  * Miura [1976] R. M. Miura, The Korteweg–de Vries Equation: A Survey of Results, SIAM Review 18, 412 (1976), https://doi.org/10.1137/1018076 .
  * Miura  _et al._ [1968] R. M. Miura, C. S. Gardner, and M. D. Kruskal, Korteweg-de Vries Equation and Generalizations. II. Existence of Conservation Laws and Constants of Motion, J. Math. Phys. 9, 1204 (1968).
  * Zakharov and Faddeev [1971] V. E. Zakharov and L. D. Faddeev, Korteweg-de Vries equation: A completely integrable Hamiltonian system, Functional Analysis and Its Applications 5, 280 (1971).
  * Lax [1968] P. D. Lax, Integrals of Nonlinear Equations of Evolution and Solitary Waves, Commun. Pure Appl. Math. 21, 467 (1968).
  * Pöschl, G. and Teller, E. [1933] Pöschl, G. and Teller, E., Bemerkungen zur Quantenmechanik des anharmonischen Oszillators, Z. Phys. 83, 143 (1933).
  * Martel and Poisson [2005] K. Martel and E. Poisson, Gravitational perturbations of the Schwarzschild spacetime: A practical covariant and gauge-invariant formalism, Phys. Rev. D 71, 104003 (2005), gr-qc/0502028 .
  * Nagar and Rezzolla [2005] A. Nagar and L. Rezzolla, Gauge-invariant non-spherical metric perturbations of Schwarzschild black-hole spacetimes, Class. Quant. Grav. 22, R167 (2005), [Erratum: Class. Quant. Grav. 23, 4297 (2006)], arXiv:gr-qc/0502064 .
  * Glampedakis  _et al._ [2017] K. Glampedakis, A. D. Johnson, and D. Kennefick, Darboux transformation in black hole perturbation theory, Phys. Rev. D 96, 024036 (2017), arXiv:1702.06459 [gr-qc] .
  * Gardner [1971] C. S. Gardner, Korteweg-de Vries equation and generalizations. IV. The Korteweg-de Vries equation as a Hamiltonian system, J. Math. Phys. 12, 1548 (1971).
  * Faddeev and Takhtajan [1987] L. D. Faddeev and L. A. Takhtajan,  _Hamiltonian Methods in the Theory of Solitons_ (Springer-Verlag, Berlin, 1987).
  * Craig [1989] W. Craig, The trace formula for Schrödinger operators on the line, Communications in Mathematical Physics 126, 379 (1989).
  * Shohat and Tamarkin [1943] J. A. Shohat and J. D. Tamarkin,  _The problem of moments_ (American Mathematical Society, Providence, Rhode Island, 1943).
  * Akhiezer [1965] N. Akhiezer,  _The Classical Moment Problem and Some Related Questions in Analysis_, University mathematical monographs (Oliver & Boyd, 1965).
  * Schmüdgen [2017] K. Schmüdgen,  _The Moment Problem_, Graduate Texts in Mathematics (Springer International Publishing, 2017).
  * Stieltjes [1894] T. J. Stieltjes, Recherches sur les fractions continues, Annales de la Faculté des sciences de Toulouse : Mathématiques 8, J1 (1894).
  * Hulburt and Katz [1964] H. Hulburt and S. Katz, Some problems in particle technology: A statistical mechanical formulation, Chemical Engineering Science 19, 555 (1964).
  * Marchisio  _et al._ [2003] D. L. Marchisio, J. T. Pikturna, R. O. Fox, R. D. Vigil, and A. A. Barresi, Quadrature method of moments for population-balance equations, AIChE Journal 49, 1266 (2003).
  * Gavriliadis and Athanassoulis [2012] P. Gavriliadis and G. Athanassoulis, The truncated stieltjes moment problem solved by using kernel density functions, Journal of Computational and Applied Mathematics 236, 4193 (2012).
  * Mead and Papanicolaou [1984] L. R. Mead and N. Papanicolaou, Maximum entropy in the problem of moments, Journal of Mathematical Physics 25, 2404 (1984), https://doi.org/10.1063/1.526446 .
  * Tagliani [1999] A. Tagliani, Hausdorff moment problem and maximum entropy: A unified approach, Applied Mathematics and Computation 105, 291 (1999).
  * John  _et al._ [2007] V. John, I. Angelov, A. Öncül, and D. Thévenin, Techniques for the reconstruction of a distribution from a finite number of its moments, Chemical Engineering Science 62, 2890 (2007).
  * Lebaz  _et al._ [2016] N. Lebaz, A. Cockx, M. Sperandio, and J. Morchain, Reconstruction of a distribution from a finite number of its moments: A comparative study in the case of depolymerization process, Computers & Chemical Engineering 84, 326 (2016).
  * Amindavar and Ritcey [1994] H. Amindavar and J. A. Ritcey, Pade approximations of probability density functions, IEEE Transactions on Aerospace and Electronic Systems 30, 416 (1994).
  * Baker and Graves-Morris [1996] G. A. Baker and P. Graves-Morris,  _Padé Approximants_, 2nd ed., Encyclopedia of Mathematics and its Applications (Cambridge University Press, 1996).
  * Bender and Orszag [1978] C. Bender and S. Orszag,  _Advanced Mathematical Methods for Scientists and Engineers_ (McGraw-Hill, 1978).
  * Billingsley [1995] P. Billingsley,  _Probability and Measure_ (John Wiley & Sons, New York, 1995).
  * Project  _et al._ [1954] B. M. Project, H. Bateman, A. Erdélyi, and U. S. O. of Naval Research,  _Tables of Integral Transforms: Based, in Part, on Notes Left by Harry Bateman_, California Institute of Technology Bateman manuscript project No. v. 1 (McGraw-Hill, 1954).
  * Landau and Lifshitz [1981] L. Landau and E. Lifshitz,  _Quantum Mechanics: Non-Relativistic Theory_, Course of Theoretical Physics (Elsevier Science, 1981).
  * Leaver [1986] E. W. Leaver, Spectral decomposition of the perturbation response of the Schwarzschild geometry, Phys. Rev. D 34, 384 (1986).
  * Ferrari and Mashhoon [1984a] V. Ferrari and B. Mashhoon, New approach to the quasinormal modes of a black hole, Phys. Rev. D 30, 295 (1984a).
  * Ferrari and Mashhoon [1984b] V. Ferrari and B. Mashhoon, Oscillations of a Black Hole, Phys. Rev. Lett. 52, 1361 (1984b).
  * Maxima [2022] Maxima, Maxima, a computer algebra system. version 5.46.0 (2022).
  * Longman [1972] I. Longman, Numerical Laplace transform inversion of a function arising in viscoelasticity, Journal of Computational Physics 10, 224 (1972).
  * Cohen [2007] A. Cohen,  _Numerical Methods for Laplace Transform Inversion_, Numerical Methods and Algorithms (Springer US, 2007).
  * Çevik  _et al._ [2016] D. Çevik, M. Gadella, Ş. Kuru, and J. Negro, Resonances and antibound states for the Pöschl–Teller potential: Ladder operators and SUSY partners, Phys. Lett. A 380, 1600 (2016), arXiv:1601.05134 [math-ph] .
  * Abramowitz and Stegun [1972] M. Abramowitz and I. A. Stegun,  _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_ (Dover, New York, 1972).