# Gravitational radiation from a particle plunging into a Schwarzschild black
hole:
frequency-domain and semirelativistic analyses

Hector O. Silva  Max Planck Institute for Gravitational Physics (Albert
Einstein Institute), D-14476 Potsdam, Germany  Giovanni Tambalo  Max Planck
Institute for Gravitational Physics (Albert Einstein Institute), D-14476
Potsdam, Germany  Kostas Glampedakis  Departamento de Física, Universidad de
Murcia, Murcia, E-30100, Spain  Theoretical Astrophysics, University of
Tübingen, Auf der Morgenstelle 10, D-72076 Tübingen, Germany  Kent Yagi
Department of Physics, University of Virginia, Charlottesville, Virginia
22904, USA

###### Abstract

We revisit the classic problem of gravitational wave emission by a test
particle plunging into a Schwarzschild black hole both in the frequency-domain
Regge-Wheeler-Zerilli formalism and in the semirelativistic approximation. We
use, and generalize, a transformation due to Nakamura, Sasaki, and Shibata to
improve the falloff of the source term of the Zerilli function. The faster
decay improves the numerical convergence of quantities of interest, such as
the energy radiated at spatial infinity through gravitational waves. As a test
of the method, we study the gravitational radiation produced by test particles
that plunge into the black hole with impact parameters close to the threshold
for scattering. We recover and expand upon previous results that were obtained
using the Sasaki-Nakamura equation. In particular, we study the relative
contributions to the total energy radiated due to waves of axial and polar
parity, and uncover an universal behavior in the waveforms at late times. We
complement our study with a semirelativistic analysis of the problem, and we
compare the two approaches. The generalized Nakamura-Sasaki-Shibata
transformation presented here is a simple and practical alternative for the
analysis of gravitational-wave emission by unbound orbits in the Schwarzschild
spacetime using the frequency-domain Regge-Wheeler-Zerilli formalism.

##  I Introduction

The study of gravitational radiation produced by test particles in black-hole
spacetimes has a long history dating back to the early 1970s, and played a
central role in the early development of the understanding of potential
gravitational-wave sources [1, 2]. In the framework of black-hole perturbation
theory, developed by Regge, Wheeler, and Zerilli [3, 4], the pioneering works
on this problem were done by Zerilli [5], Davies et al. [6, 7, 8], Chung [9],
and Ruffini [10]. These works assumed particles in unbound trajectories that
start from spatial infinity and plunge into a Schwarzschild black hole.

Critical to these calculations is the asymptotic behavior of the source term
that encapsulates how the test particle excites the gravitational
perturbations. As an extreme example, the source term of the Teukolsky
equation, that describes the gravitational perturbations of a Kerr black hole
[11], diverges at spatial infinity for unbound geodesics. In principle, this
jeopardizes the calculation of physical quantities of interest, such as
gravitational waveforms or the energy carried away to infinity by the waves.
To circumvent this problem, one can either develop a regularization scheme
[12, 13, 14, 15], or rewrite the Teukolsky equation to tame the source’s
asymptotic behavior. Pursuing the latter approach, Sasaki and Nakamura [16,
17, 18] found their eponymous equation, widely used in the study of unbound
geodesics both in Schwarzschild and Kerr spacetimes; see, e.g., Refs. [19, 20,
21] and [22, 23, 24], respectively, for early work. A somewhat similar
situation also happens when the perturbations of a Schwarzschild black hole
are described in terms of the Cunningham-Price-Moncrief [25] and Zerilli-
Moncrief [26] functions. Less dramatically, the source term of the Zerilli
equation has a slow falloff at spatial infinity. See Hopper [27] for a
detailed discussion.

In a serendipitous event, in the course of a related investigation, we learned
of a work by Shibata and Nakamura that presents a simple transformation of the
Zerilli function to improve the falloff of the source term of the Zerilli
equation, and that was applied to the radial infall case only
[28].111Reference [28] cites an unpublished work by Nakamura and Sasaski when
introducing this transformation, that we will refer to as the Nakamura-Sasaki-
Shibata transformation. Here we revisit this method and extend its application
range to problems involving particles plunging with nonzero angular momentum.
This allows us to reexamine some aspects of the gravitational radiation
produced by particles that plunge with angular momenta near the threshold for
scattering. This situation is relevant in the context of understanding
ultrarelativistic binary black hole collisions, and that was until now only
studied in detail using the Sasaki-Nakamura equation; see Berti et al. [29]
and references therein. We complement our study of this problem with a
calculation of the energy radiated within the semirelativistic approximation
of Ruffini and Sasaki [30].

This paper is organized as follows. In Sec. II, we review the motion of test
particles plunging into a Schwarzschild black hole. In Sec. III, we provide a
summary of Regge-Wheeler-Zerilli formalism and identify the main issue we want
to resolve. In Sec. IV, we review and generalize the method of Ref. [28]. In
Sec. V, we describe our numerical methods and present our numerical results.
In Sec. VI, we compare our results against an analysis in the semirelativistic
approximation. We summarize our findings in Sec. VII. We use the mostly plus
metric signature and use geometrical units with $c=G=1$, unless stated
otherwise.

##  II Geodesic motion

We consider a particle of mass $\mu$ in geodesic motion the spacetime of a
Schwarzschild black hole of mass $M$, with $\mu/M\ll 1$. We use Schwarzschild-
Droste coordinates $x^{\mu}=\\{t,r,\theta,\phi\\}$ in which the spacetime’s
line element is

${\rm d}s^{2}=-f(r)\,{\rm d}t^{2}+f^{-1}(r)\,{\rm d}r^{2}+r^{2}({\rm
d}\theta^{2}+\sin^{2}\theta\,{\rm d}\phi^{2})\,,$ (1)

where $f=1-2M/r$ and $r=2M$ is the location of the event horizon.

Figure 1:  The effective potential $U(r,{\cal L})$ for various values of
${\cal L}$, which ranges from $3.25M$ (lowermost solid gray line) to $4.25M$
(uppermost solid gray line). For a particle falling from rest from infinity
(${\cal E}=1$) into the hole, there is a critical value of ${\cal L}_{\rm
crit}=4M$ (solid black line) for which the particle becomes trapped in a
marginally stable circular orbit at $r=4M$. This value of ${\cal L}$ separates
plunging (${\cal L}<{\cal L}_{\rm crit}$) from scattering (${\cal L}>{\cal
L}_{\rm crit}$) geodesics.

We assume that the particle starts from rest at infinity with (conserved)
energy ${\cal E}=E/\mu=1$ and angular momentum ${\cal L}=L/\mu$ per unit mass.
If we chose, without loss of generality, that the particle’s motion happens in
the equatorial plane $\theta=\pi/2$, then we can parametrize the particle’s
worldline in terms of the proper time $\tau$ as
$z^{\mu}(\tau)=\\{t_{p},r_{p},\pi/2,\phi_{p}\\}$, and, from the geodesic
equation and the timelike constraint $g_{ab}u^{a}u^{b}=-1$, with $u^{a}={\rm
d}z^{a}/{\rm d}t$, we obtain:

$\dot{t}_{p}={\cal E}/f_{p}\,,\quad\dot{\phi}_{p}={\cal
L}/r_{p}^{2}\,,\quad\dot{r}_{p}^{2}={\cal E}^{2}-U(r_{p},{\cal L})\,,$ (2)

where $\dot{}={\rm d}/{\rm d}\tau$ and

$U=f\,(1+{\cal L}^{2}/r^{2})\,,$ (3)

is the effective radial potential.

The particle’s trajectories are classified according to the number of real
roots of ${\cal E}^{2}-U$. For the special case ${\cal E}=1$ the analysis is
simple, and we find the roots to be,

$r_{\pm}=\frac{{\cal L}}{4M}\,[{\cal L}\pm({\cal L}^{2}-16M^{2})^{1/2}]\,.$
(4)

Hence, a plunging orbit from infinity requires that ${\cal L}<4M$, for
otherwise the turning points are real and positive (i.e., the particle is
scattered.) We will define this special value of the angular momentum as
${\cal L}_{\rm crit}=4M$.

In Fig. 1 we show the effective potential (3) for a range of values ${\cal
L}/M\in\\{3.25,\,4.25\\}$ (light curves). The thicker line corresponds to
$U(r,{\cal L}_{\rm crit})$, which peaks at $r=4M$ with value of 1. Hence, a
particle falling from rest and with angular momentum ${\cal L}={\cal L}_{\rm
crit}$ will be captured in a marginally stable circular orbit. A particle with
${\cal L}$ larger (smaller) than ${\cal L}_{\rm crit}$ will scatter (plunge).

Figure 2:  Timelike geodesics with angular momentum per unit mass ${\cal
L}/M\in\\{3.25,\,3.9996\\}$ plunging into a Schwarzschild black hole (black
disk) starting from rest at spatial infinity, ${\cal E}=1$. The dot-dashed
line represents to the innermost stable circular orbit ($r=6M$) and the dashed
line corresponds to the location of the marginally stable circular orbit
($r=4M$). In the limit ${\cal L}/{\cal L}_{\rm crit}\to 1$, the particle
executes an increasing fractional number of orbits $\phi/(2\pi)$, as seen in
the red curve that corresponds to ${\cal L}=3.9996M$, or $99.99\%$ of ${\cal
L}_{\rm crit}$.

.

It is convenient to rewrite Eq. (2) as first-order in $r$ equations,

$\displaystyle{\rm d}t_{p}/{\rm d}r_{p}$ $\displaystyle=-({\cal
E}/f_{p})\,({\cal E}^{2}-U_{p})^{-1/2}\,,$ (5a) $\displaystyle{\rm
d}\phi_{p}/{\rm d}r_{p}$ $\displaystyle=-({\cal L}/r_{p}^{2})\,({\cal
E}^{2}-U_{p})^{-1/2}\,,$ (5b)

where we have taken $\dot{r}_{p}<0$. We integrate Eqs. (5) with initial
conditions $t_{p}(r_{\rm max})=0$ and $\phi_{p}(r_{\rm max})=0$ at some
arbitrarily large $r_{p}=r_{\rm max}$ down to the horizon, $r_{p}=2M$. In Fig.
2 we show a sequence of trajectories starting from ${\cal L}/M=3.25$ and up to
${\cal L}/M=3.9996$ (i.e., with 99.99% of ${\cal L}_{\rm crit}$). We translate
from Schwarzschild-Droste to Cartesian coordinates using

$x_{p}=r_{p}\cos\phi_{p}\,,\quad y_{p}=r_{p}\sin\phi_{p}\,,\quad z_{p}=0\,.$
(6)

As the ratio ${\cal L}/{\cal L}_{\rm crit}$ approaches unity from below, the
particle executes an increasing fractional number of orbits $\phi_{p}/(2\pi)$,
given approximately by $-(\sqrt{2}\pi)^{-1}\,\log(1-{\cal L}/{\cal L}_{\rm
crit})$ [31].

##  III Black hole perturbations in the Regge-Wheeler-Zerilli gauge

We are interested in calculating the gravitational waves produced by a
particle plunging in a Schwarzschild black hole. The standard treatment of
this problem, in the metric-perturbation formalism, is due to Regge and
Wheeler [3] and Zerilli [4, 5]. The problem reduces to solving two equations
in the time domain:

$\left[-\frac{\partial^{2}}{\partial t^{2}}+\frac{\partial^{2}}{\partial
x^{2}}-V^{(\pm)}_{\ell}(r)\right]X^{(\pm)}_{\ell m}(t,r)=S^{(\pm)}_{\ell
m}(t,r)\,,$ (7)

or, by going to the Fourier domain through

$X^{(\pm)}_{\ell m}(t,r)=\frac{1}{2\pi}\int_{-\infty}^{\infty}{\rm
d}\omega\,e^{-i\omega t}\,X^{(\pm)}_{\ell m\omega}(r)\,,$ (8)

we have alternatively

$\left[\frac{{\rm d}^{2}}{{\rm
d}x^{2}}+\omega^{2}-V^{(\pm)}_{\ell}(r)\right]X^{(\pm)}_{\ell
m\omega}(r)=S^{(\pm)}_{\ell m\omega}(r)\,.$ (9)

In these equations, $x$ is the tortoise coordinate

$x=r+2M\log[r/(2M)-1]\,,$ (10)

that maps the region $2M<r<\infty$ to $-\infty<x<\infty$. The superscript
$(\pm)$ denotes variables associated to metric perturbations of polar ($+$) or
axial ($-$) parity, and $V^{(\pm)}_{\ell}$ is an effective potential.
Perturbations of each parity are described by a single master function, known
as the Zerilli $X^{(+)}$ and Regge-Wheeler $X^{(-)}$ functions, respectively.
The effective potentials associated to each of these perturbations bear the
same respective names and are given by

$\displaystyle V_{\ell}^{\,(+)}$
$\displaystyle=\frac{f}{r^{2}\Lambda^{2}}\left[2\lambda^{2}\left(\Lambda+1\right)+\frac{18M^{2}}{r^{2}}\left(\lambda+\frac{M}{r}\right)\right],$
(11a) $\displaystyle V_{\ell}^{\,(-)}$
$\displaystyle=\frac{f}{r^{2}}\left[\ell(\ell+1)-\frac{6M}{r}\right],$ (11b)

where we defined

$\lambda=(\ell+2)(\ell-1)/2\quad\textrm{and}\quad\Lambda=\lambda+3M/r\,.$ (12)

The function $S^{(\pm)}_{\ell m\omega}$ is the source term, responsible for
the excitation of the gravitational perturbations. A detailed derivation of
this source term in the Regge-Wheeler-Zerilli formalism can be found, e.g., in
Refs. [32, 33], whose results we quote in Appendix A.

Because the potential $V_{\ell}^{(\pm)}$ vanishes both at the horizon and at
infinity, and provided that the source vanishes sufficiently fast at both
boundaries, the solutions to Eq. (7) can be written as plane waves for
$x\to\pm\infty$. We will impose that $X^{\,(\pm)}_{\ell m\omega}$ is purely
ingoing at the event horizon and purely outgoing at spatial infinity, that is,

$\displaystyle X^{\,(\pm)}_{\ell
m\omega}\simeq\begin{cases}\,C^{\,(\pm),\,{\rm in}}_{\ell
m\omega}\,e^{-i\omega x}\quad&x\to-\infty\\\ \,C^{\,(\pm),\,{\rm out}}_{\ell
m\omega}\,e^{+i\omega x}\quad&x\to+\infty\end{cases}\,.$ (13)

To solve the inhomogeneous differential equation (7), we use the method of
variation of parameters. The method consists of finding two linearly
independent solutions, say $X^{\,(\pm),\,{\rm in}}_{\ell m\omega}$ and
$X^{\,(\pm),\,{\rm up}}_{\ell m\omega}$, of the homogeneous equations

$\left[\frac{{\rm d}^{2}}{{\rm
d}x^{2}}+\omega^{2}-V^{(\pm)}_{\ell}(r)\right]X^{\,(\pm),\,{\rm up/in}}_{\ell
m\omega}=0\,.$ (14)

The two solutions differ by the boundary conditions we impose when we solve
Eq. (14), namely,

$X^{\,(\pm),\,{\rm in}}_{\ell m\omega}\simeq\begin{cases}\,e^{-i\omega
x}\,&x\to-\infty\\\ \,A^{\,(\pm),\,{\rm in}}_{\ell m\omega}\,e^{-i\omega
x}+A^{\,(\pm),\,{\rm out}}_{\ell m\omega}\,e^{+i\omega
x}\,&x\to+\infty\end{cases}\\\ $ (15)

and

$X^{\,(\pm),\,{\rm up}}_{\ell m\omega}\simeq\begin{cases}\,B^{\,(\pm),\,{\rm
in}}_{\ell m\omega}\,e^{-i\omega x}+B^{\,(\pm),\,{\rm out}}_{\ell
m\omega}\,e^{+i\omega x}\,&x\to-\infty\\\ \,e^{i\omega
x}\,&x\to+\infty\end{cases}\,.$ (16)

We use these solutions to define the Wronskian,

$\displaystyle W^{\,(\pm)}_{\ell m\omega}=X^{\,(\pm),\,{\rm in}}_{\ell
m\omega}\,\frac{{\rm d}{X^{\,(\pm),\,{\rm up}}_{\ell m\omega}}}{{\rm
d}x}-X^{\,(\pm),\,{\rm up}}_{\ell m\omega}\,\frac{{\rm d}{X^{\,(\pm),\,{\rm
in}}_{\ell m\omega}}}{{\rm d}x}\,,$ (17)

which is constant in $x$. We then construct the Green’s function as

$\displaystyle G^{\,(\pm)}_{\ell m\omega}(x,x^{\prime})$
$\displaystyle=\frac{1}{W^{\,(\pm)}_{\ell m\omega}}\left[X^{\,(\pm),\,{\rm
up}}_{\ell m\omega}(x)\,X^{\,(\pm),\,{\rm in}}_{\ell
m\omega}(x^{\prime})\,\Theta(x-x^{\prime})\right.$
$\displaystyle\quad\left.+\,X^{\,(\pm),\,{\rm in}}_{\ell
m\omega}(x)\,X^{\,(\pm),\,{\rm up}}_{\ell
m\omega}(x^{\prime})\,\Theta(x^{\prime}-x)\right],$

where $\Theta(\cdot)$ is the Heaviside step function. Finally, with the
Green’s function (LABEL:eq:greens_fun), we can write the solution of Eq. (7)
as

$\displaystyle X^{\,(\pm)}_{\ell m\omega}$
$\displaystyle=\frac{1}{W^{\,(\pm)}_{\ell m\omega}}\left[X^{\,(\pm),\,{\rm
up}}_{\ell m\omega}\int_{-\infty}^{x}{\rm d}x^{\prime}\,X^{\,(\pm),\,{\rm
in}}_{\ell m\omega}\,S^{(\pm)}_{\ell m\omega}\right.$
$\displaystyle\quad+\left.X^{\,(\pm),\,{\rm in}}_{\ell
m\omega}\int_{x}^{+\infty}{\rm d}x^{\prime}\,X^{\,(\pm),\,{\rm up}}_{\ell
m\omega}\,S^{(\pm)}_{\ell m\omega}\right]\,.$ (19)

Equation (19) has the desired property of being purely ingoing at the horizon
and purely outgoing at infinity:

$\displaystyle\lim_{{x\to-\infty}}$ $\displaystyle X^{\,(\pm)}_{\ell
m\omega}\simeq\frac{e^{-i\omega x}}{W^{\,(\pm)}_{\ell
m\omega}}\int_{-\infty}^{\infty}{\rm d}x^{\prime}\,X^{\,(\pm),\,{\rm
up}}_{\ell m\omega}\,S^{(\pm)}_{\ell m\omega}\,,$ (20a)
$\displaystyle\lim_{{x\to+\infty}}$ $\displaystyle X^{\,(\pm)}_{\ell
m\omega}\simeq\frac{e^{+i\omega x}}{W^{\,(\pm)}_{\ell
m\omega}}\int_{-\infty}^{\infty}{\rm d}x^{\prime}\,X^{\,(\pm),\,{\rm
in}}_{\ell m\omega}\,S^{(\pm)}_{\ell m\omega}\,,$ (20b)

where we used the boundary conditions on $X^{\,(\pm),\,{\rm up/in}}_{\ell
m\omega}$.

We can compare Eq. (20b) with the boundary conditions Eq. (13) imposed on
$X^{\,(\pm)}_{\ell m\omega}$, and make the identification,

$\displaystyle C^{\,(\pm),\,{\rm out}}_{\ell
m\omega}=\frac{1}{W^{\,(\pm)}_{\ell m\omega}}\int_{-\infty}^{+\infty}{\rm
d}x^{\prime}\,X^{\,(\pm),\,{\rm in}}_{\ell m\omega}\,S^{(\pm)}_{\ell
m\omega}\,.$ (21)

This is the main quantity we must numerically calculate to obtain, e.g., the
time-domain gravitational waveform or the energy radiated to infinity. More
specifically, we will be interested in the spectral energy distribution in
each $(\ell,m)$ mode

$\displaystyle\frac{{\rm d}E_{\ell m}}{{\rm
d}\omega}=\frac{\omega^{2}}{64\pi^{2}}\frac{(\ell+2)!}{(\ell-2)!}\left[|C^{\,(+),\,{\rm
out}}_{\ell m\omega}|^{2}+\frac{4}{\omega^{2}}|C^{\,(-),\,{\rm out}}_{\ell
m\omega}|^{2}\right],$

and in the time-domain Regge-Wheeler and Zerilli mode functions

$X_{\ell m}^{\,(\pm)}(t,x)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}{\rm
d}\omega\,C^{\,(\pm),\,{\rm out}}_{\ell m\omega}\,e^{-i\omega(t-x)}\,.$ (23)

In practice, we note that as $x\to\infty$, the Wronskian becomes

$W^{\,(\pm)}_{\ell m\omega}=2i\omega A^{\,(\pm),\,{\rm in}}_{\ell m\omega}\,,$
(24)

and that we can also rewrite Eq. (21) as

$\displaystyle C^{\,(\pm),\,{\rm out}}_{\ell
m\omega}=\frac{1}{W^{\,(\pm)}_{\ell m\omega}}\int_{2M}^{+\infty}\frac{{\rm
d}r^{\prime}}{f}\,X^{\,(\pm),\,{\rm in}}_{\ell m\omega}\,S^{(\pm)}_{\ell
m\omega}\,.$ (25)

Because $f\simeq 1$ and $X^{\,(\pm),\,{\rm in}}_{\ell m\omega}\simeq\exp(\pm
i\omega x)$ as $x\to\infty$, we see that the convergence of this integral
depends critically on the asymptotic properties of $S^{(\pm)}_{\ell m\omega}$.

##  IV Asymptotic behavior of the frequency domain sources

We now review the properties of $S_{\ell m\omega}^{\,(+)}$ in two cases of
interest. We will start by reviewing the case in which the particle falls
radially, ${\cal L}=0$, into the black hole and how Ref. [28] improved the
asymptotic behavior of $S_{\ell m\omega}^{\,(+)}$. We will then show one way
in which the Nakamura-Sasaki-Shibata transformation can be generalized to
general plunging trajectories, ${\cal L}\neq 0$, and discuss the asymptotic
properties of this new source term.

###  IV.1 The case of radial infall

Figure 3: The $\ell=m=2$ source term for the Zerilli equation for a radially
plunging particle starting from rest at spatial infinity, for $M\omega=0.3$.
Left panel: the absolute values of the original Zerilli source term,
$S^{\,(+)}$, and its Nakamura-Sasaki-Shibata-transformed version,
$\tilde{S}^{\,(+)}$. Right panels: the real (“Re”) and imaginary (“Im”) parts
of $S^{\,(+)}$ (top) and $\tilde{S}^{\,(+)}$ (bottom). Both sources vanish at
the event horizon $r=2M$, and the faster asymptotic decay of
$\tilde{S}^{\,(+)}$ improves the convergence of Eq. (21). The source terms are
largest near $r\simeq 3M$, which corresponds approximately to the location of
the peak of the Zerilli potential $r^{(+)}_{\ell=2,\,{\rm peak}}\simeq 3.1M$
and of the light ring $r=3M$.

When ${\cal L}=0$, $S^{\,(-)}_{\ell m\omega}$ vanishes and $S^{\,(+)}_{\ell
m\omega}$ simplifies to [7, 6]

$\displaystyle S^{\,(+)}_{\ell m\omega}$ $\displaystyle=-\,8\pi\mu{\cal
A}_{\ell
m}\frac{f}{r\Lambda}\left[\sqrt{\frac{r}{2M}}-\frac{2i}{\omega}\frac{\lambda}{r\Lambda}\right]e^{i\omega
t_{p}(r)}\,,$ (26)

where ${\cal A}_{\ell m}=Y^{\ast}_{\ell m}(\pi/2,\phi)\exp(im\phi)$, $Y_{\ell
m}$ are the spherical harmonics, the asterisk indicates complex conjugation,
and $t_{p}$ is given by

$\displaystyle\frac{t_{p}}{2M}=-\frac{2}{3}\left(\frac{r}{2M}\right)^{\tfrac{3}{2}}-2\left(\frac{r}{2M}\right)^{\tfrac{1}{2}}+\log\left[\frac{\sqrt{r/(2M)}+1}{\sqrt{r/(2M)}-1}\right].$

We find that the near-horizon and spatial-infinity behaviors of
$S^{\,(+)}_{\ell m\omega}$ are

$\displaystyle S^{\,(+)}_{\ell m\omega}\simeq\begin{cases}0&\quad
x\to-\infty\\\ x^{-1/2}&\quad x\to+\infty\\\ \end{cases}\,.$ (28)

Hence the integral Eq. (21) converges slowly at spatial infinity. To improve
the convergence, Ref. [28] proposed the substitution222Our notation differs
from that used in Ref. [28]. See also Ref. [13], for a similar substitution in
the context of the Teukolsky equation.

$X^{\,(+)}_{\ell m\omega}=\tilde{X}^{\,(+)}_{\ell m\omega}+Q_{\ell
m\omega}\,,$ (29)

where

$Q_{\ell m\omega}=-\frac{8\pi\mu{\cal A}_{\ell
m}}{\omega^{2}}\frac{f}{r\Lambda}\sqrt{\frac{2M}{r}}e^{i\omega t_{p}}\,.$ (30)

The function $Q_{\ell m\omega}$ vanishes at the event horizon $x\to-\infty$
and decays as $x^{-3/2}$ for $x\to\infty$. Thus, $X^{\,(+)}_{\ell
m\omega}\simeq\tilde{X}^{\,(+)}_{\ell m\omega}$ in the latter limit.

We then insert Eq. (29) in Eq. (7) to find

$\left[\frac{{\rm d}^{2}}{{\rm
d}x^{2}}+\omega^{2}-V^{\,(+)}_{\ell}\right]\tilde{X}^{\,(+)}_{\ell
m\omega}=\tilde{S}^{\,(+)}_{\ell m\omega}\,,$ (31)

where $\tilde{S}^{\,(+)}_{\ell m\omega}$ is a new source term given by

$\tilde{S}^{\,(+)}_{\ell m\omega}=S^{\,(+)}_{\ell m\omega}-\left[\frac{{\rm
d}^{2}}{{\rm d}x^{2}}+\omega^{2}-V^{\,(+)}_{\ell}\right]Q_{\ell m\omega}\,.$
(32)

The asymptotic behaviors of the new source term are

$\displaystyle\tilde{S}^{\,(+)}_{\ell m\omega}\simeq\begin{cases}0&\quad
x\to-\infty\\\ x^{-3/2}&\quad x\to+\infty\\\ \end{cases}\,,$ (33)

making the integral in Eq. (21) converge faster. One can verify that it is the
second $x$ derivative of $Q_{\ell m\omega}$ in Eq. (32) that yields a term
that decays as $r^{-1/2}$ as $x\to\infty$, and that cancels the leading-order
asymptotic term in the large-$x$ expansion of the original source term. We
note that both $X^{\,(+)}_{\ell m\omega}$ and $\tilde{X}^{\,(+)}_{\ell
m\omega}$ satisfy the same homogeneous equation. Hence, we can calculate
$C^{\,(\pm),\,{\rm out}}_{\ell m\omega}$ with the mere replacement
$S^{\,(+)}_{\ell m\omega}\to\tilde{S}^{\,(+)}_{\ell m\omega}$ in Eq. (25).

In Fig. 3, we show both the original (26) and the Nakamura-Sasaki-Shibata-
transformed (32) source terms of the Zerilli equation for a particle falling
radially into the black hole, and $\ell=m=2$ and $M\omega=0.3$. In the left
panel we show their absolute values. In the right panels we show the real
(solid line) and imaginary (dashed line) parts of the original (top panel) and
new (bottom panel) source terms. The faster decay of Eq. (32) accelerates the
convergence of Eq. (21). The value $M\omega=0.3$ is approximately the real
part of the fundamental quasinormal mode frequency of a Schwarzschild black
hole [34] and dominates the emission of energy in the form of gravitational
waves [12]. The source terms are largest near $r\simeq 3M$, which corresponds
approximately to the location of peak of the Zerilli potential
$r^{(+)}_{\ell=2,\,{\rm peak}}\simeq 3.1M$, and of the light ring $r=3M$.

The success of the Nakamura-Sasaki-Shibata transformation opens two questions.
First, can we, even for a radially infalling particle, make the source term
decay more rapidly? Second, can we generalize Eqs. (29) and (30) to the case
in which the particle plunges with a nonzero angular momentum? In the next
section we give positive answers to both questions.

###  IV.2 The case of general plunging trajectories

We seek a minimal extension of Eqs. (29) and (30) for the case of a plunging
particle with general angular momentum ${\cal L}\neq 0$. We focus only on the
source term in the Zerilli equation, because the source term in the Regge-
Wheeler equation already has a fast falloff:

$S^{\,(-)}_{\ell m\omega}\simeq x^{-2}\,,\quad x\to\infty\,.$ (34)

We would like to make the source term in the Zerilli equation to decay at
least as rapidly as the source term for the Regge-Wheeler equation.

Figure 4: The $(\ell,m)=(2,2)$ source term for the Zerilli equation for a
plunging particle with ${\cal L}=0.9999\,{\cal L}_{\rm crit}$ starting from
rest at spatial infinity, for $M\omega=0.3$. Left panel: the absolute values
of the original Zerilli source term, $S^{\,(+)}$, and its Nakamura-Sasaki-
Shibata-transformed version, $\tilde{S}^{\,(+),\,(N)}$, for $N=0$, 1, and 2 in
the expansion (35). Right panels: the real (“Re”) and imaginary (“Im”) parts
of $S^{\,(+)}$ (top) and $\tilde{S}^{\,(+),\,(N)}$ (three remaining panels)
for increasing values of $N$. The source term is largest near $r\simeq 4M$,
that corresponds to the location of the marginally stable circular orbit at
which the particle particle executes an increasing number of revolutions in
the limit ${\cal L}/{\cal L}_{\rm crit}\to 1$.

Our strategy is to retain Eq. (29), but generalize $Q_{\ell m\omega}$ to the
form

$\displaystyle Q_{\ell m\omega}$ $\displaystyle=-\frac{8\pi\mu{\cal A}_{\ell
m}}{\omega^{2}}\frac{f}{\Lambda
r}\left[\,\sum_{n=0}^{N}a_{n}\,\left({2M}/{r}\right)^{{(n+1)}/{2}}\,\right]$
$\displaystyle\quad\times e^{i(\omega t_{p}-m\phi_{p})}\,.$ (35)

Proceeding as in the previous section, we arrive at the same Eq. (32), where
$S^{\,(+)}_{\ell m\omega}$ is now given by a more complicated formula, whose
form can be found in Eq. (59), Appendix A. We then expand Eq. (32) in a power
series in $r$ for $r\to\infty$, and collect powers in $r$. Next, we fix the
coefficients $a_{n}$ in such a way to cancel each successive power of $r$ in
this power series. This yields

$a_{0}=1\,,\quad a_{1}=\frac{im}{\ell(\ell+1)}\frac{{\cal L}}{M}\,,\quad
a_{2}=1-\frac{{\cal L}^{2}}{8M^{2}}\,,$ (36)

by working up to $N=2$. Higher-order terms can be obtained by following the
same procedure just outlined. We notice there is no dependence on ${\cal L}$
at order $N=0$ and that we recover Eq. (30) in this case. The next correction
to the radial infall case occurs at $n=2$.

In Fig. 4 we show the original source term of the Zerilli equation
$S^{\,(+)}_{\ell m\omega}$ and its transformed version
$\tilde{S}^{\,(+),\,(N)}_{\ell m\omega}$ (for $N=$ 0, 1, and 2) for
$\ell=m=2$, $M\omega=0.3$ and a near-critical geodesic with ${\cal
L}=0.9999\,{\cal L}_{\rm crit}$. We chose this value of ${\cal L}$ as a
“stress test” for the method, due to the large amplitude of the source term
around the location of the marginally stable circular orbit, $r=4M$. This peak
dominates over the peak at light ring $r=3M$ present in Fig. 3. The left panel
shows the absolute values of the various source terms, whereas in the right
panels we show their real (solid lines) and imaginary parts (dashed lines).

We see that the original Nakamura-Sasaki-Shibata transformation ($N=0$)
results in a slower decaying source term, $r^{-1}$, when applied to case in
which ${\cal L}\neq 0$, when compared to the radial infall case, $r^{-3/2}$.
This can be understood by examining the asymptotic behavior of
$S^{\,(+)}_{\ell m\omega}$:

$\displaystyle S^{\,(+)}_{\ell m\omega}$ $\displaystyle=-\frac{4\pi\mu{\cal
A}_{\ell m}}{\lambda
M}\left[\sqrt{\frac{2M}{r}}+\frac{im}{\ell(\ell+1)}\frac{{\cal
L}}{M}\frac{2M}{r}\right.$ $\displaystyle\left.+\left(\frac{{\cal
L}^{2}}{8M^{2}}+\frac{3+2\lambda}{2\lambda}\right)\left(\frac{2M}{r}\right)^{3/2}\right]+{\cal
O}(r^{-2})\,.$ (37)

This expansion shows that the first term containing the angular momentum
${\cal L}$ appears at order $r^{-1}$, and it cannot be canceled with Eq. (30).

##  V Numerical methods and results

In this section we describe the numerical methods we use to compute Eq. (25),
and show a few applications of solving this equation using the generalized
Nakamura-Sasaki-Shibata transformation (35).

###  V.1 Numerical methods

To evaluate our main quantity of interest, namely the amplitude
$C^{(\pm),\,\rm out}_{\ell m\omega}$ we proceed as follows.

  1. 1.

Choose a value of $\ell$, $|m|\leq\ell$, and $M\omega$.

  2. 2.

Choose a value of the angular momentum ${\cal L}$, and solve the geodesic
equations (5), to obtain $t_{p}(r)$ and $\phi_{p}(r)$. We use initial
conditions $t_{p}(r_{\rm min})=0$ and $\phi_{p}(r_{\rm min})=0$, and integrate
from $r_{\rm min}=2\,(1+10^{-4})M$ up to $r_{\rm max}=500/\omega$.

  3. 3.

If $\ell+m$ is even, integrate the homogeneous Zerilli equation. If $\ell+m$
is odd, integrate the homogeneous Regge-Wheeler equation; see Eqs. (14) and
(11).

  4. 4.

Integrate the equation from the previous step with “in” boundary conditions,
from $r_{\rm min}=2\,(1+10^{-4})M$ to $r=r_{\rm max}$. From $X^{(\pm)}_{\ell
m\omega}$ and ${\rm d}X^{(\pm)}_{\ell m\omega}/{\rm d}r$ at $r_{\rm max}$
calculate $A^{(\pm),\,{\rm in},\,{\rm out}}_{\ell m\omega}$. See Appendix C.

  5. 5.

Calculate the Wronskian (24) and the source $S^{(\pm)}_{\ell m\omega}$ to
evaluate the integral in Eq. (25).

  6. 6.

Repeat steps 1 through 5, scanning the range $M\omega\in[5\times
10^{-3},\,1.5]$ in steps of size $\Delta(M\omega)=5\times 10^{-3}$, for
$\ell=2$ to 6, covering all $|m|\leq\ell$ in between.

Our code is written in Mathematica. In step 3, we validated our integration of
the homogeneous equations (14) by comparison against the integrators available
in the Black Hole Perturbation Toolkit [35]. Recursion relations for the near-
horizon and far-field expansions of the Regge-Wheeler and Zerilli functions
are summarized in Appendix B. In step 5, we calculate the Wronskian as
explained in Appendix C, and it is useful to rewrite Eq. (25) as a
differential equation [36, 37]:

$\frac{{\rm d}C^{\,(\pm),\,{\rm out}}_{\ell m\omega}}{{\rm
d}r}=\frac{1}{W^{\,(\pm)}_{\ell m\omega}}\,f^{-1}\,X^{\,(\pm),\,{\rm
in}}_{\ell m\omega}\,{S}^{\,(\pm)}_{\ell m\omega}\,.$ (38)

We integrate Eq. (38), with initial condition $C^{\,(\pm),\,\rm out}_{\ell
m\omega}(r_{\rm min})=0$ up to $r_{\rm max}$, using the Regge-Wheeler source
(64) or the Nakamura-Sasaki-Shibata-transformed Zerilli source (32), depending
on whether $\ell+m$ is odd or even, respectively. In the context of Eq. (38),
the new Zerilli source term reduces the value of $r_{\rm max}$ at which ${{\rm
d}C^{\,(\pm),\,{\rm out}}_{\ell m\omega}}/{{\rm d}r}$ becomes zero to the
desired accuracy. We decided here only, with respect to the rest this work, to
integrate the geodesic equations from near the horizon outwards. In
particular, when integrating from spatial infinity, $t_{p}$ acquires a large
value around $r\approx 2M$. This causes the factor $\exp(i\omega t_{p})$,
which appears in the source, to become highly oscillatory, and thus sensitive
to our choice of $r_{\rm min}$ and $\omega$. Consequently, the phase (but not
the amplitude) of the $C^{\,(\pm),\,{\rm out}}_{\ell m\omega}$ failed to
converge in our numerical calculations. We solved this issue by fixing the
initial conditions at the horizon instead of spatial infinity. At last, we
performed steps 1 to 6 for the range of angular momentum values ${\cal L}=$
$\\{0$, 0.25, 0.5, 0.75, 0.9, 0.99, 0.999, $0.9999\\}{\cal L}_{\rm crit}$.

###  V.2 Energy spectrum

Figure 5: The energy spectra for a particle plunging with angular momentum
${\cal L}=0.9\,{\cal L}_{\rm crit}$ (left panels) and ${\cal L}=0.9999\,{\cal
L}_{\rm crit}$ (right panels) into a Schwarzschild black hole. Top panels: the
energy spectra from the multipoles $\ell=m$ from 2 to 6, (colored lines), and
their sum (dashed line). Bottom panels: the energy spectra for the quadrupole
perturbation $\ell=2$ and all $|m|\leq 2$ in between. Note the different
ranges in the ordinates across the panels. Our results agree with Ref. [29],
which solved the Sasaki-Nakamura equation instead.

We first consider the energy spectrum for the particle in the near-critical
limit ${\cal L}\approx{\cal L}_{\rm crit}$. To our knowledge, this situation
was studied only in the frequency domain using the Sasaki-Nakamura equation
[29] or using the Zerilli-Moncrief [26] and Cunningham-Price-Moncrief [25]
master functions in Ref. [38]. To validate our calculations, we focus on two
cases, ${\cal L}/{\cal L}_{\rm crit}=0.9$ and $0.9999$, which were examined in
detail by Berti et al. [29].

In Fig. 5, we show the energy spectrum (LABEL:eq:dEdw) for ${\cal
L}=0.9\,{\cal L}_{\rm crit}$ (left panels) and ${\cal L}=0.9999\,{\cal L}_{\rm
crit}$ (right panels). The top panels show the spectral energy for multipoles
$\ell=m$ ranging from 2 to 6 (different line colors) and their sum (dashed
line). The bottom panels show the quadrupole mode $\ell=2$ and all azimuthal
contributions $|m|\leq 2$ (different line styles.) Our results are in
excellent agreement with those of Ref. [29]; cf. Figs. 10 and 11 therein.

For ${\cal L}=0.9\,{\cal L}_{\rm crit}$, the energy radiated in each multipole
$\ell=m$ has a single maximum. The location $M\omega$ of these peaks coincide
approximately with the real part of the fundamental ($n=0$) quasinormal mode
frequency of a Schwarzschild black hole $M\omega_{\ell m0}$, as first observed
by Detweiler and Szedenits [12]. For reference, these values are ${\rm
Re}[M\omega_{\ell m0}]\simeq\\{0.373,\,0.599,\,0.809,\,1.012,\,1.212\\}$, for
$\ell\in\\{2,\,6\\}$ [39, 40, 41]. (Due to the spherical symmetry of the
Schwarzschild solution, all quasinormal modes with $|m|\leq\ell$, at fixed
$\ell\geq 2$, are degenerate.)

In the near-critical limit ${\cal L}=0.9999~{}{\cal L}_{\rm crit}$, the
spectral energy distribution has a peak at $M\omega<{\rm Re}[M\omega_{\ell
m0}]$. This peak corresponds to $m$ times the particle’s orbital frequency
$M\Omega_{\rm orb}=1/8$ at the marginally stable circular orbit at $r=4M$ (cf.
Fig. 2), that is

$M\omega=m\,M\,\Omega_{\rm orb}=m/8\,.$ (39)

Therefore, the energy emitted at these frequencies is dominated by the
particle’s geodesic motion. This is not surprising given that we saw that the
maximum amplitude of the source term of the Zerilli equation is located at
$r\simeq 4M$, as we discussed in Fig. 4. For moderate values of $\ell$, the
contributions to the energy driven by the particle’s orbital motion and the
quasinormal mode excitation overlap, while as $\ell\gg 1$ their contributions
separate sufficiently to result in two peaks in the energy spectra; see the
$\ell=m=5$ and 6 in Fig. 5. For $\ell=m=5$ the peaks are located at
$M\omega\simeq 0.637$ and $0.919$, while for $\ell=m=6$ they are found at
$M\omega\simeq 0.763$ and $1.125$. We can estimate this separation as follows.
First, we use the geometrical-optics (eikonal) limit [42, 43, 44], to
approximate the real part of the quasinormal mode frequency as

${\rm Re}[M\omega_{\rm eik}]\simeq\ell\,\Omega_{\rm LR}=\ell/(3\sqrt{3}),$
(40)

where $M\Omega_{\rm LR}=1/(3\sqrt{3})$ is the orbital frequency of a null
geodesic at the light ring. Equation (40) can also be obtained from the
$\ell\gg 1$ limit of a Wentzel–Kramers–Brillouin approximation to the
calculation of black-hole quasinormal modes [45, 46]. We can then estimate the
separation between the “quasinormal-mode” and “geodesic” peaks by taking the
difference of Eqs. (39) and (40):

$(\Delta M\omega)_{\rm
peak}\simeq\frac{1}{3\sqrt{3}}\left(\ell-\frac{3\sqrt{3}}{8}m\right),$ (41)

which is valid for $\ell\gg 1$, ${\cal L}/{\cal L}_{\rm crit}\simeq 1$, and
${\cal E}=1$. For $\ell=m=6$, we find $(\Delta M\omega)_{\rm peak}\simeq
0.404$, in fair agreement with the numerical result $\simeq 0.362$, obtained
from the difference between the two peak locations.

Although not studied here, it is interesting to analyze the ultrarelativistic
limit, in which the particle plunges with an energy ${\cal E}\to\infty$. In
this limit, the “geodesic” peak moves rightwards, eventually overlapping with
the “quasinormal-mode” peak; see Ref. [29], Fig. 10. This occurs because as
${\cal E}\to\infty$, the particle’s marginally stable circular orbit coincides
with the light ring, hence $\Omega_{\rm orb}\simeq\Omega_{\rm LR}$ [31]. In
this limit, we then have

$(\Delta M\omega)_{\rm peak}=(\ell-m)\,(3\sqrt{3}),$ (42)

which is valid for $\ell\gg 1$, ${\cal L}/{\cal L}_{\rm crit}\simeq 1$, and
${\cal E}\to\infty$. The peak separation vanishes when $\ell=m$, reproducing
the results of Ref. [29].

###  V.3 Total energy

From the energy spectrum, we can compute the total energy $\Delta E$ emitted
in form of gravitational waves by integrating the spectrum density and summing
the contributions from all multipoles:

$\Delta E=\sum_{\ell=2}^{\infty}\,\Delta
E_{\ell}=\sum_{\ell=2}^{\infty}\sum_{m=-\ell}^{\ell}\,\int_{0}^{+\infty}{\rm
d}\omega\,\frac{{\rm d}E_{\ell m}}{{\rm d}\omega}\,.$ (43)

To understand the relative contribution due to perturbations of each parity,
we also define $\Delta E^{\,(\pm)}$, where the subscript $(+)$ means we add
only the contributions to the energy coming from multipoles for which $\ell+m$
is even and $(-)$ when $\ell+m$ is odd.

For particles that plunge from infinity initially from rest, Oohara and
Nakamura [20], based on an earlier observation by Davis et al. [6] for the
case ${\cal L}=0$, proposed the empirical relation

$(M/\mu^{2})\,\Delta E_{\ell}=a\,\exp(-b\,\ell\,),$ (44)

between energy and multipole $\ell$. We fit Eq. (44) to the outcome of
computing $\Delta E_{\ell}$, truncating the $\ell$ sum at $\ell_{\rm max}=6$.
Table 1 shows the fitting coefficients and the total energy emitted, including
the individual polar and axial contributions. The results for the total energy
agree with those of Ref. [29], Table II, with less than 1% difference, when
comparison is possible. As observed by Oohara and Nakamura, the value of $a$
increases while that of $b$ decreases as ${\cal L}$ approaches ${\cal L}_{\rm
crit}$. In other words, the fraction of the total energy radiated contained in
each $\ell$ pole tends to “even out” as ${\cal L}$ increases.

Figure 6 illustrates these observations. For example, we see in the leftmost
panel that $\Delta E_{2}$ varies from being about four orders of magnitude
larger than $\Delta E_{6}$, to being some more than ten times larger as ${\cal
L}$ approaches criticality. The individual contributions of the polar and
axial perturbations to the energy are shown in the middle and rightmost
panels, respectively. We see that the energy emitted through the “axial
radiative channel” is always subdominant relative to the polar one, even for
near-critical values of ${\cal L}$; the extreme case is, of course, when
${\cal L}=0$ where $\Delta E^{\,(-)}$ vanishes.

${\cal L}/{\cal L}_{\rm crit}$ | $a$ | $b$ | $a^{\,(+)}$ | $b^{\,(+)}$ | $a^{\,(-)}$ | $a^{\,(-)}$ | $\Delta E$ | $\Delta E^{\,(+)}$ | $\Delta E^{\,(-)}$
---|---|---|---|---|---|---|---|---|---
|  |  |  |  |  |  | $\cdot 10^{-2}$ | $\cdot 10^{-2}$ | $\cdot 10^{-4}$
0.0000 | 0.45 | 1.99 | 0.45 | 1.99 | – | – | 1.04 | 1.04 | –
0.2500 | 0.29 | 1.56 | 0.27 | 1.56 | 0.01 | 1.59 | 1.71 | 1.63 | 7.59
0.5000 | 0.31 | 1.18 | 0.27 | 1.17 | 0.04 | 1.25 | 4.48 | 4.01 | 47.4
0.7500 | 0.47 | 0.91 | 0.40 | 0.91 | 0.07 | 0.98 | 13.0 | 11.4 | 169
0.9000 | 0.76 | 0.79 | 0.66 | 0.79 | 0.09 | 0.83 | 28.7 | 25.4 | 321
0.9900 | 1.75 | 0.77 | 1.63 | 0.77 | 0.12 | 0.74 | 71.5 | 66.2 | 525
0.9990 | 2.88 | 0.80 | 2.75 | 0.80 | 0.13 | 0.72 | 109 | 103 | 609
0.9999 | 4.02 | 0.82 | 3.89 | 0.82 | 0.15 | 0.71 | 144 | 137 | 669
Table 1: Coefficients $a$ and $b$ in the Oohara-Nakamura relation (44) and the
total energy $(\mu^{2}/M)\,\Delta E$ for different values of angular momentum
${\cal L}$. The superscripts refer to polar $(+)$ and axial-parity $(-)$
perturbations. Figure 6: Energy radiated to infinity in the plunge process.
From left to right, the panels show the total energy radiated, only through
the parity even (“Zerilli”) perturbations, and only through the parity odd
(“Regge-Wheeler”) perturbations. The markers correspond to Eq. (43), whereas
the straight lines corresponds to the fit (44) suggested by Oohara and
Nakamura [20]. The lines starting from the bottom correspond to ${\cal L}=0$
and increase to ${\cal L}=0.9999\,{\cal L}_{\rm crit}$ as one moves upwards.
The energy emitted in the dominant quadrupole mode in the axial-parity
radiative channel is always smaller than that emitted through the polar one,
at fixed ${\cal L}$. When ${\cal L}=0$, there is no energy radiated in the
axial-parity channel.

Interestingly, the axial and polar contributions to the total energy $\Delta
E$ are nonmonotonic with respect to the angular momentum. In Fig. 7, we show
the ratio $\Delta E^{\,(\pm)}/\Delta E$ as a function of the logarithm of
$1-{\cal L}/{\cal L}_{\rm crit}$. As we increase the particle’s angular
momentum (i.e., moving right to left along the abscissa), we see that $\Delta
E^{\,(-)}/\Delta E$ has a local maximum at ${\cal L}\approx 0.75\,{\cal
L}_{\rm crit}$, yet with only approximately 14% of the total energy budget. We
interpret this result as follows. When ${\cal L}=0$, by symmetry, all energy
must be radiated through the polar channel. As we increase ${\cal L}$, axial
perturbations become increasingly excited, and a nonzero (albeit small)
percentage of the energy is emitted through them. As ${\cal L}$ approaches
${\cal L}_{\rm crit}$, the particle orbits the black hole an increasing number
of times around the _circular_ orbit $r=4M$. Hence, again by symmetry, we
expect the total energy fraction emitted through the polar channel to
increase, and it happens to the extent of thereby decreasing the fraction
emitted through the axial channel.

Figure 7: Fraction of the total energy radiated by the polar and axial
perturbations in the plunge process. We show both the polar [$(+)$, left
ordinate] and axial [$(-)$, right ordinate] total energies as functions of the
particle’s angular momentum. The fraction of the total energy radiated via
each radiative channel is nonmonotonic in ${\cal L}$, with the axial channel
having a peak at ${\cal L}/{\cal L}_{\rm crit}\approx 0.75$. The polar
contribution is at least $\approx 85\%$ of the total energy budget for all
values of ${\cal L}$ considered.

###  V.4 Time-domain waveforms

We also computed the time-domain Regge-Wheeler and Zerilli mode functions
using Eq. (23). In Fig. 8, we show the dominant modes for the Regge-Wheeler
and Zerilli functions, $(\ell,m)=(2,1)$ and $(2,2)$, respectively, as a
function of the retarded time $(t-x)/M$. The top panel corresponds to a
particle that plunges with ${\cal L}=0.9~{}{\cal L}_{\rm crit}$, while the
bottom panel to ${\cal L}=0.9999~{}{\cal L}_{\rm crit}$. In the former case,
we see that the waveform has the characteristic “precursor,” “sharp burst,”
and “ringing tail” phases, as first observed for radial infalling particles by
Davis et al. [8] and for wave scattering by Vishveshwara [47]. As ${\cal
L}\to{\cal L}_{\rm crit}$, the large values of ${\cal L}$ cause the Zerilli
(but not the Regge-Wheeler) function to have an intermediate
quasimonochromatic phase, related to the particle’s sweep around $r=4M$. This
behavior is qualitatively the same as that seen in Refs. [48, 49] for test
particles plunging from the innermost stable circular orbit $r=6M$. We see
that the amplitude of the Regge-Wheeler function is always smaller than
Zerilli’s, even in the nearest-critical angular-momentum values studied by us.

In Fig. 9, we take a closer look on the waveforms in the limit ${\cal
L}\to{\cal L}_{\rm crit}$, for the same multipole moments. In this limit, the
waveforms become increasingly similar around their peak amplitudes,
$t-x\simeq-40\,M$, although differences are clearly visible at the “precursor”
phase. Similar results hold for the other multipoles we examined. This result
suggests that a quasiuniversal description of the plunge from the marginally
stable circular orbit $r=4M$ exists. Such a treatment for particle’s plunging
from the innermost stable circular orbit exists [50].

Figure 8: Time-domain Zerilli (solid lines) and Regge-Wheeler (dashed lines)
mode functions excited by a test particle plunging with ${\cal L}=0.9\,{\cal
L}_{\rm crit}$ (top panel) and ${\cal L}=0.9999\,{\cal L}_{\rm crit}$ (bottom
panel) into a Schwarzschild black hole. We focus on the dominant quadrupole
multipole, associated with the perturbations of each parity: $m=2$ and $m=1$
for the Zerilli and Regge-Wheeler modes, respectively. The former is always
larger in amplitude than the latter.  Figure 9: Time-domain Zerilli (top
panel) and Regge-Wheeler (bottom panel) mode functions excited by a test
particle plunging with ${\cal L}=0.99\,{\cal L}_{\rm crit}$ (solid lines),
${\cal L}=0.999\,{\cal L}_{\rm crit}$ (dashed lines) and ${\cal
L}=0.9999\,{\cal L}_{\rm crit}$ (dot-dashed lines) into a Schwarzschild black
hole. As in Fig. 8, we consider the dominant quadrupolar Zerilli and Regge-
Wheeler modes. In the limit ${\cal L}\to{\cal L}_{\rm crit}$, the solutions
become quasiuniversal around and after $t-x\gtrsim-40\,M$.

##  VI The semirelativistic approximation

A complementary way of studying the gravitational wave emission by an
infalling particle is via the so-called “semirelativistic” approximation [30]
(often used in the “kludge” approximation [51, 52, 53]). In this approach the
particle is assumed to move along a fully relativistic geodesic trajectory of
the black-hole spacetime while the gravitational wave emission itself is
treated approximately, using the weak-gravity quadrupole formula.

Despite its inherent inconsistency, the semirelativistic approximation is
known to perform surprisingly well, when compared against more rigorous
results obtained from black-hole perturbation theory. The price one has to pay
for the conceptual and technical simplicity of this approach is that its
accuracy deteriorates as soon as the particle’s trajectory enters the near
horizon, strong field regime where radiation backscattering by the spacetime
curvature becomes an important factor. Unfortunately, this condition is
clearly met by a plunging particle so we expect the semirelativistic method to
provide accurate results only for the early time portion of the trajectory
(i.e., the low frequency part of the gravitational wave spectrum.)
Nevertheless, the less accurate $\omega M\gtrsim 1$ part of the spectrum is of
some interest in its own right as it allows us to understand (and separate)
the effects due to the particle’s motion and due to backscattering.

The quadrupole-order gravitational-wave formalism underpinning the
semirelativistic approximation can be found in many general relativity
textbooks; here we follow and expand the analysis of Ref. [54], Sec. 4.3.1,
for radially infalling particles.

We start by recalling that the appropriately averaged gravitational-wave
luminosity is given by

$L=\tfrac{1}{5}\,\langle\,\dddot{M}_{ij}\dddot{M}^{\,ij}-\tfrac{1}{3}\dddot{M}^{2}\,\rangle,$
(45)

where the (mass) quadrupole moment $M_{ij}$ is defined as

$M^{ij}(t)=\int{\rm d}^{3}x\,\rho(t,\bm{x})\,x^{i}x^{j}\,,$ (46)

for a mass density $\rho(t,\bm{x})$, and with trace $M=M^{i}{}_{i}$, which is
distinguishable from the black hole’s mass $M$ from the context. The total
energy emitted in gravitational waves is given by the integral

$\Delta E=\int_{-\infty}^{t_{\rm max}}{\rm d}t\,L(t)\,,$ (47)

where the luminosity is to be evaluated without any averaging (also, we can
set $t_{\max}=\infty$ at the particle’s horizon crossing time). The same
quantity can be evaluated in the frequency domain

$\Delta E=\int_{0}^{+\infty}{\rm d}\omega\,\frac{{\rm d}E}{{\rm d}\omega}\,,$
(48)

where ${\rm d}E/{\rm d}\omega$ is the spectral energy distribution.

For a point particle moving along a trajectory $\bm{x}_{p}(t)$ we have
$\rho(t,\bm{x})=\mu\,\delta^{(3)}[\bm{x}-\bm{x}_{p}(t)]$, and we find

$M^{ij}(t)=\mu\,x_{p}^{i}(t)\,x_{p}^{j}(t)\,.$ (49)

As in Sec. II, the geodesics under consideration can be taken to lie in the
equatorial $x$-$y$ plane and the Cartesian coordinates can be related to the
Schwarzschild-Droste coordinates (1) through Eqs. (6). In this setup, the only
nonvanishing components of the quadrupole moment (46) are $M_{11}$, $M_{22}$,
$M_{12}$, and the trace is $M=M_{11}+M_{12}$. Then, a short calculation leads
to

$\displaystyle\Delta E$
$\displaystyle=\frac{2}{15}\int_{-\infty}^{+\infty}{\rm
d}t\,[\dddot{M}_{11}^{2}+\dddot{M}_{22}^{2}-\dddot{M}_{11}\dddot{M}_{22}+3\dddot{M}_{12}^{2}]\,,$
$\displaystyle=\frac{2}{15\pi}\int_{0}^{+\infty}{\rm
d}\omega\,\omega^{6}\,[\,|\tilde{M}_{11}|^{2}+|\tilde{M}_{22}|^{2}+3|\tilde{M}_{12}|^{2}$
$\displaystyle\quad-\,{\rm Re}[\tilde{M}_{11}\tilde{M}^{*}_{22}]\,]\,,$ (50)

where the Fourier transform of $M_{ij}$ is defined as333For clarity, we adopt
a slightly different notation for frequency-domain quantities in this section,
mirroring Ref. [54].

$\tilde{M}_{ij}(\omega)=\int_{-\infty}^{+\infty}{\rm d}t\,e^{i\omega
t}\,M_{ij}(t)\,.$ (51)

The fact that $M_{ij}$ is real implies the useful property

$\tilde{M}_{ij}^{*}(\omega)=\tilde{M}_{ij}(-\omega)\,.$ (52)

As discussed in Ref. [54], the integral (51) is divergent at $t\to-\infty$,
since $x_{p}\to+\infty$. Therefore, some regularization procedure is required.
This is achieved by working with the Fourier transform
$\tilde{\ddot{M}}_{ij}$, which is well behaved at spatial infinity. In fact,
this procedure is equivalent to the regularization of $\tilde{M}_{ij}(\omega)$
via integrations by parts, where the produced divergent boundary terms are
discarded (see, e.g., Detweiler and Szedenits [12] for a similar
regularization of the solution of the Teukolsky equation sourced by a plunging
particle). The outcome of this exercise is the regularized quadrupole moment

$\tilde{M}_{ij}^{\rm reg}=-\tilde{\ddot{M}}_{ij}/\omega^{2}\,,$ (53)

and from Eq. (50) we read off the regularized formula:

$\displaystyle\frac{{\rm d}E}{{\rm d}\omega}$
$\displaystyle=\frac{2\,\omega^{2}}{15\pi}[\,|\tilde{\ddot{M}}_{11}|^{2}+|\tilde{\ddot{M}}_{22}|^{2}+3|\tilde{\ddot{M}}_{12}|^{2}-{\rm
Re}[\tilde{\ddot{M}}_{11}\tilde{\ddot{M}}^{*}_{22}]\,]\,.$

For the actual numerical evaluation of $\tilde{\ddot{M}}_{ij}$ (and the
subsequent one of ${\rm d}E/{\rm d}\omega$), it is advantageous to convert the
time integral into a radial integral using the geodesic equations (2), leading
to

$\tilde{\ddot{M}}_{ij}=\int^{r_{\rm max}}_{2M}\frac{{\rm
d}r}{f}\left[\frac{2M}{r}+f\,\frac{{\cal
L}^{2}}{r^{2}}\right]^{-\tfrac{1}{2}}\ddot{M}_{ij}\,e^{i\omega t_{p}(r)}\,.$
(55)

For the problem at hand, the individual $\ddot{M}_{ij}$ components required
for the calculation of Eq. (53) are

$\displaystyle\ddot{M}_{11}$
$\displaystyle=2\mu\left(x_{p}\ddot{x}_{p}+\dot{x}^{2}_{p}\right),$ (56a)
$\displaystyle\ddot{M}_{22}$
$\displaystyle=2\mu\left(y_{p}\ddot{y}_{p}+\dot{y}^{2}_{p}\right),$ (56b)
$\displaystyle\ddot{M}_{12}$
$\displaystyle=\mu\left(2\dot{x}_{p}\dot{y}_{p}+x_{p}\ddot{y}_{p}+y_{p}\ddot{x}_{p}\right).$
(56c)

With the help of Eq. (2) both the accelerations $\ddot{x}_{p}$, $\ddot{y}_{p}$
and the velocities $\dot{x}_{p}$, $\dot{y}_{p}$ can be rewritten solely as
functions of $r$, $\phi_{p}(r)$ and the orbital constant ${\cal L}$. At the
same time $t_{p}(r)$ explicitly appears in the exponential term. Both
$\phi_{p}(r)$ and $t_{p}(r)$ are obtained numerically as in Sec. II. We chose
trajectories with $\phi_{p}(0)=0$ at some initial radius $r_{\rm max}=r(0)$.
The same radius serves as the “infinity” upper limit in the
$\tilde{\ddot{M}}_{ij}$ integral. The initial time $t=0$ can be chosen
arbitrarily because the addition of a constant in $t_{p}(r)$ does not affect
the numerical value of Eq. (55).

For any value ${\cal L}<{\cal L}_{\rm crit}$, the integral in the
$\tilde{\ddot{M}}_{11}$ component is the slowest converging one at large $r$
(the absolute value of the integrand decays as $\sim r^{-1/2}$). In order to
improve the convergence of this component and reduce the associated numerical
error we employ a procedure similar to that in Sec. IV.2, namely, we subtract
the slowly decaying part prior to the numerical integration and then add back
the (analytically obtained) asymptotic contribution from that part.

Figure 10: The energy spectra for a particle plunging with angular momentum
${\cal L}=0$, ${\cal L}=0.75\,{\cal L}_{\rm crit}$, and ${\cal
L}=0.9999\,{\cal L}_{\rm crit}$ into a Schwarzschild black hole in the
semirelativistic approximation (“SR,” dashed lines) and black-hole
perturbation theory, for $\ell=m=2$ (“BHPT,” solid lines). Both calculations
agree at low frequencies $M\omega\ll 1$, where the radiation is due to the
particle’s quasi-Newtonian motion. For ${\cal L}=0$ and ${\cal L}=0.75\,{\cal
L}_{\rm crit}$, the kludge calculation underestimates the location of the
spectrum’s peak, which is related to the hole’s fundamental quasinormal mode
frequency [12]. As ${\cal L}\to{\cal L}_{\rm crit}$, the peak of the spectrum
is dominated by the particle’s orbit around $r=4M$, and the two calculations
agree qualitatively with each other up to $M\omega\lesssim 0.25$. The
semirelativistic approximation predicts a slowly decaying tail for the
spectrum.

The semirelativistic energy spectrum, calculated from Eqs.
(LABEL:eq:dEdw_kludge) and (55), is shown in Fig. 10 and has a characteristic
single-hump profile. When comparing against the full black-hole perturbation
theory result for $\ell=m=2$ shown in Fig. 5, we notice a moderate
disagreement in the location of the emission peak. This is not surprising
because the exact location of this peak depends on the properties of the BH
potential near $3M$; as pointed out earlier, this information is missing
altogether in the semirelativistic approximation. As expected, the two spectra
are in good agreement in the low frequency, $\omega M\ll 1$, end of the
spectrum which is the part that corresponds to quasi-Newtonian motion, that
is, when the particle is still moving far away from the black-hole horizon.
Interestingly, the agreement improves to some degree as we approach the
critical value ${\cal L}_{\rm crit}$ for scattering. In this case, the
particle spends a large amount of time orbiting around $r=4\,M$ and the bulk
of gravitational-wave emission is generated there, as we have discussed in
Sec. V.

The opposite high-frequency end of the spectrum, $\omega M\gg 1$, has also
some interest of its own. In the fully relativistic calculation, this part of
the spectrum appears to be independent of the angular momentum ${\cal L}$ and
the modal numbers $\ell$ and $m$. We can approximate the high-frequency “tail”
of ${\rm d}E_{\ell m}/{\rm d}\omega$ in the special case of radial infall,
taking into account that this part of the spectrum corresponds to the near-
horizon region of integration in Eq. (25). Under these circumstances we can
recast the integrand into a single exponential, and expand the exponent up to
linear order in $r/(2M)-1$. The integral can then be computed analytically,
and it is found to be dominated by the lower limit of integration. Moreover,
in the same high-frequency limit the Wronskian can be approximated as $W_{\ell
m\omega}\simeq 2i\omega$. These manipulations lead to a scaling

$\frac{{\rm d}E}{{\rm d}\omega}\propto\frac{e^{-8\pi
M\omega}}{M\omega}\quad\textrm{(perturbation theory)}\,.$ (57)

The same procedure applied to the semirelativistic spectrum leads to a slower
decaying tail

$\frac{{\rm d}E}{{\rm d}\omega}\propto\frac{e^{-4\pi
M\omega}}{M\omega}\quad\textrm{(semirelativistic)}\,.$ (58)

The difference can be traced back to the additional highly oscillatory
function $X^{(\pm),\,{\rm in}}_{\ell m\omega}\simeq\exp(i\omega x)$ in the
full black-hole perturbation theory expression, which suppresses the integral
at large $\omega$. This analysis thus explains the difference in the high-
frequency tails shown in Fig. 10.

##  VII Conclusions

We reviewed, and generalized, a transformation by Nakamura, Sasaki, and
Shibata [28] that makes the source term of the Zerilli function to have a
faster falloff at spatial infinity, thereby improving the numerical
convergence of the convolution integral that arises in the calculation of
gravitational radiation by particles in unbound geodesic motion in the
Schwarzschild spacetime.

As an application, we studied the gravitational radiation produced by test
particles that plunge from rest and with angular momentum ${\cal L}$ from
spatial infinity into a Schwarzschild black hole. In particular, we focused on
the limit in which ${\cal L}$ approximates from below the critical value
${\cal L}_{\rm crit}$ for scattering. To our knowledge, this is the first time
this calculation was done using the original Regge-Wheeler and Zerilli master
functions. Our results are in agreement with the work by Berti et al. [29]
that used the Sasaki-Nakamura equation. We studied in detail the relative
contributions to the energy radiated in gravitational waves due to
perturbations of polar and axial parity. We found that the former always
dominates. We also observed an quasiuniversal late-time behavior of the
waveforms in limit in which ${\cal L}$ approaches the critical value for
scattering, $4M$.

The main merit of the Nakamura-Sasaki-Shibata transformation is that it only
requires minimal modifications to the source term of the Zerilli function [cf.
Eqs. (30) and (32)]. The new source term can be easily computed analytically
with any symbolic algebra software. In contrast, the Sasaki-Nakamura formalism
requires the numerical integration of an auxiliary second-order differential
equation for the calculation of the source term [16, 17, 18]. However, an
advantage of the Sasaki-Nakamura formalism is that it applies to the Kerr
spacetime, while the method presented here is not. Our work is alternative to
that of Hopper [38] that addresses the sources of the Zerilli-Moncrief and
Cunningham-Price-Moncrief master functions.

We complemented our calculations in black-hole perturbation theory, with an
analysis of the same plunging-particle problem in the semirelativistic
approximation. We found that the two energy spectra agree best with in the
limit ${\cal L}\to{\cal L}_{\rm crit}$, when the energy is dominated by the
particle’s motion around the marginally stable circular orbit. We also studied
the high-frequency limit of the spectrum, expanding upon the discussion in
Ref. [54].

The method presented here can be used in other problems involving unbound
motion in the Schwarzschild spacetime, including ultrarelativistic plunges
[10, 55], scattering [27], or in the spacetime of relativistic stars [56, 57,
58, 59] using the Regge-Wheeler-Zerilli formalism.

## Acknowledgments

We thank Emanuele Berti, Benjamin Leather, Caio F. B. Macedo, Raj Patil,
Masaru Shibata, Jan Steinhoff, Nicholas C. Stone and Helvi Witek for
discussions. H.O.S acknowledges funding from the Deutsche
Forschungsgemeinschaft (DFG) – Project No. 386119226. K.G. acknowledges
support from research Grant No. PID2020-1149GB-I00 of the Spanish Ministerio
de Ciencia e Innovación. K.Y. acknowledges support from NSF Grants No.
PHY-2207349 and No. PHY-2309066, a Sloan Foundation Research Fellowship, and
the Owens Family Foundation. This work makes use of the Black Hole
Perturbation Toolkit [35]. Some of our calculations were performed in the
Hypatia cluster at the Max Planck Institute for Gravitational Physics.

##  Appendix A The sources of the Regge-Wheeler and Zerilli equations

In this Appendix we present the sources for the Regge-Wheeler and Zerilli
equations (7), quoting from the formulas presented in Refs. [32, 33]; the
original derivation is due Zerilli [5].

The source $S^{\,(+)}_{\ell m\omega}$, which excites the perturbations of
polar parity, is given by

$\displaystyle S^{\,(+)}_{\ell m\omega}$ $\displaystyle=-\,if\frac{{\rm
d}{}}{{\rm d}r}\left[\frac{f^{2}}{\Lambda}\left(\frac{ir}{f}\,\tilde{C}_{1\ell
m\omega}+\tilde{C}_{2\ell m\omega}\right)\right]$
$\displaystyle\quad+i\frac{f^{2}}{r\Lambda^{2}}\left[i\frac{\lambda
r^{2}-3\lambda Mr-3M^{2}}{rf}\,\tilde{C}_{1\ell m\omega}\right.$
$\displaystyle\quad\left.+\frac{\lambda(\lambda+1)r^{2}+3\lambda
Mr+6M^{2}}{r^{2}}\,\tilde{C}_{2\ell m\omega}\right],$ (59)

where444We note that in the equation for $\tilde{B}_{\ell m\omega}$, the term
proportional to $A^{(1)}_{\ell m\omega}$ in Eq. (A42) of Ref. [32] has a typo,
which is corrected in Eq. (4.12) of Ref. [33].

$\displaystyle\tilde{B}_{\ell m\omega}$ $\displaystyle=\frac{8\pi
r^{2}f}{\Lambda}\left[A_{\ell
m\omega}+\frac{1}{\sqrt{\ell(\ell+1)/2}}\,B_{\ell m\omega}\right]$
$\displaystyle\quad-4\pi\frac{\sqrt{2}}{\Lambda}\frac{M}{\omega}A^{(1)}_{\ell
m\omega}\,,$ (60) $\displaystyle\tilde{C}_{1\ell m\omega}$
$\displaystyle=\frac{8\pi}{\sqrt{2}\omega}A^{(1)}_{\ell
m\omega}+\frac{1}{r}\tilde{B}_{\ell m\omega}$ $\displaystyle\quad-{16\pi
r}\left[\frac{1}{2}\frac{(\ell+2)!}{(\ell-2)!}\right]^{-\tfrac{1}{2}}\,F_{\ell
m\omega}\,,$ (61) $\displaystyle\tilde{C}_{2\ell m\omega}$
$\displaystyle=\frac{8\pi
i}{\omega\sqrt{\ell(\ell+1)/2}}\frac{r}{f}\,B^{(0)}_{\ell
m\omega}-\frac{i}{f}\,\tilde{B}_{\ell m\omega}$
$\displaystyle\quad+\frac{16\pi
ir^{2}}{f}\left[\frac{1}{2}\frac{(\ell+2)!}{(\ell-2)!}\right]^{-\tfrac{1}{2}}F_{\ell
m\omega}\,.$ (62)

Here, $f=1-2M/r$, and $A_{\ell m\omega}$, $A^{(1)}_{\ell m\omega}$,
$B^{(0)}_{\ell m\omega}$, $B_{\ell m\omega}$, $\tilde{B}_{\ell m\omega}$ and
$F_{\ell m\omega}$ are the Fourier-domain projections of the particle’s
energy-momentum tensor onto the tensor harmonic basis; cf. Ref. [32], Table I,
for their expressions in time-domain and generic orbit. These functions encode
information about the particle’s geodesic motion, and when particularized for
plunging geodesics they read

$\displaystyle A_{\ell m\omega}$
$\displaystyle=\mu\,\frac{V}{r^{2}f^{2}}\,Y^{\ast}_{\ell m}\,e^{i\omega
t_{p}}\,,$ (63a) $\displaystyle A^{(1)}_{\ell m\omega}$
$\displaystyle=-i\sqrt{2}\,\mu\,\frac{{\cal E}}{r^{2}f}\,Y^{\ast}_{\ell
m}\,e^{i\omega t_{p}}\,,$ (63b) $\displaystyle B^{(0)}_{\ell m\omega}$
$\displaystyle=i\mu\,\frac{{\cal E}{\cal
L}}{Vr^{3}}\frac{1}{\sqrt{\ell(\ell+1)/2}}\partial_{\phi}Y^{\ast}_{\ell
m}\,e^{i\omega t_{p}}\,,$ (63c) $\displaystyle B_{\ell m\omega}$
$\displaystyle=-\mu\,\frac{{\cal
L}}{r^{3}f}\frac{1}{\sqrt{\ell(\ell+1)/2}}\,\partial_{\phi}Y^{\ast}_{\ell
m}\,e^{i\omega t_{p}}\,,$ (63d) $\displaystyle F_{\ell m\omega}$
$\displaystyle=\mu\,\frac{{\cal
L}^{2}}{Vr^{4}}\left[\frac{1}{2}\frac{(\ell+2)!}{(\ell-2)!}\right]^{-\tfrac{1}{2}}\partial_{\phi\phi}Y^{\ast}_{\ell
m}\,e^{i\omega t_{p}},$ (63e)

where, for brevity, we wrote $Y^{\ast}_{\ell m}=Y^{\ast}_{\ell
m}(\pi/2,\phi_{p})$ and defined $V=\sqrt{{\cal E}^{2}-U}$, where $U$ is given
by Eq. (3). Here, $\mu$ is the particle’s mass and $t_{p}$ and $\phi_{p}$ are
functions of $r$, obtained by integrating the geodesic equation (5).

The source $S^{\,(-)}_{\ell m\omega}$, that excites the perturbations of axial
parity, is given by

$\displaystyle S^{\,(-)}_{\ell m\omega}$ $\displaystyle=\frac{8\pi
if}{r}\left[\frac{1}{2}\frac{(\ell+2)!}{(\ell-2)!}\right]^{-\tfrac{1}{2}}\left[-r^{2}\frac{{\rm
d}}{{\rm d}r}(fD_{\ell m\omega})\right.$
$\displaystyle\quad\left.+\,\sqrt{2\lambda}\,rfQ_{\ell
m\omega}\vphantom{\frac{1}{1}}\right],$ (64)

where, analogously to the polar case, $D_{\ell m\omega}$ and $Q_{\ell
m\omega}$ are the Fourier-domain projections of the particle’s energy-momentum
tensor onto the tensor harmonic basis; cf. Ref. [32], Table I, for their time
domain, general orbit forms. When particularized for plunging trajectories,
$D_{\ell m\omega}$ and $Q_{\ell m\omega}$ become

$\displaystyle D_{\ell m\omega}$ $\displaystyle=i\mu\,\frac{{\cal
L}^{2}}{Vr^{4}}\left[\frac{1}{2}\frac{(\ell+2)!}{(\ell-2)!}\right]^{-\tfrac{1}{2}}X^{\ast}_{\ell
m}\,e^{i\omega t_{p}}\,,$ (65a) $\displaystyle Q_{\ell m\omega}$
$\displaystyle=-i\mu\,\frac{{\cal
L}}{fr^{3}}\frac{1}{\sqrt{\ell(\ell+1)/2}}\,\partial_{\theta}Y^{\ast}_{\ell
m}\,e^{i\omega t_{p}}\,,$ (65b)

where we introduced the shorthand notation,

$\displaystyle X_{\ell
m}=2\partial_{\phi}(\partial_{\theta}-\cot\theta)Y_{\ell m}\,.$ (66a)

In the case of radial infall, ${\cal L}=0$, both $D_{\ell m\omega}$ and
$Q_{\ell m\omega}$ are zero, and consequently $S^{\,(-)}_{\ell m\omega}$,
vanishes.

##  Appendix B The near-horizon and far-field expansions of the Regge-Wheeler
and Zerilli functions

When integrating numerically the homogeneous Regge-Wheeler and Zerilli
equations, we use series representations of the solutions to these equations
in the near-horizon [$r/(2M)-1\ll 1$ or $x\to-\infty$] and asymptotic ($r\gg
2M$ or $x\to\infty$) regions of the Schwarzschild spacetime. The coefficients
in these series satisfy certain recursion relations that we summarize here.

We first consider the near-horizon limit. We assume that the Regge-Wheeler and
Zerilli mode functions $X^{\,(\pm)}_{\ell m\omega}$ can be factorized as,

$X^{\,(\pm)}_{\ell m\omega}=e^{-i\omega
x}\sum_{n=0}^{\infty}c_{n}^{\,(\pm)}\,z(r)^{n}\,\,,\quad z=r/(2M)-1\,.$ (67)

We substitute Eq. (71) in the homogeneous equation (14) and derive a recursion
relation between the coefficients $c_{n}^{\,(\pm)}$. For the Zerilli equation
this relation is

$\displaystyle\mathfrak{a}\,c^{\,(+)}_{n}$
$\displaystyle=\mathfrak{b}\,c_{n-1}^{\,(+)}+\mathfrak{c}\,c_{n-2}^{\,(+)}+\mathfrak{d}\,c_{n-3}^{\,(+)}+\mathfrak{e}\,c_{n-4}^{\,(+)}+\mathfrak{f}\,c_{n-5}^{\,(+)}\,,$

where

$\displaystyle\mathfrak{a}$ $\displaystyle=n(3+2\lambda)^{2}(n-4i\sigma)\,,$
$\displaystyle\mathfrak{b}$
$\displaystyle=(3+2\lambda)\,\\{4\lambda^{2}+i(40\lambda+36)(n-1)\sigma$
$\displaystyle\quad+\lambda[2(9-4n)n-6]-3(n-2)(2n-1)\\}\,,$
$\displaystyle\mathfrak{c}$
$\displaystyle=24\lambda^{3}-\lambda^{2}(24n^{2}-108n+72)-\lambda(36n^{2}-168n+174)$
$\displaystyle\quad-9n^{2}+i(160\lambda^{2}+288\lambda+108)(n-2)\sigma+45n-54\,,$
$\displaystyle\mathfrak{d}$
$\displaystyle=24\lambda^{3}-\lambda^{2}(16n^{2}-108n+144)-\lambda(12n^{2}-84n+144)$
$\displaystyle\quad+i(160\lambda^{2}+192\lambda+36)(n-3)\sigma\,,$
$\displaystyle\mathfrak{e}$
$\displaystyle=i\lambda\left(80\lambda+48\right)(n-4)\sigma-4\lambda^{2}(n-6)(n-3)+8\lambda^{3}\,,$
$\displaystyle\mathfrak{f}$ $\displaystyle=16i\lambda^{2}(n-5)\sigma\,,$ (69)

and $\sigma=M\omega$ [39]. For the Regge-Wheeler equation we obtain instead

$\displaystyle n(n+4i\sigma)\,c_{n}^{\,(-)}$
$\displaystyle=-[2n^{2}-12i(n-1)\sigma-5n+\ell(\ell+1)$
$\displaystyle\quad-6]\,c_{n-1}^{\,(-)}-[(n-\ell-3)(n+\ell-2)$
$\displaystyle\quad-12i(n-2)\sigma]\,c_{n-2}^{\,(-)}+4i(n-3)\sigma\,c_{n-3}^{\,(-)}.$

In Eqs. (LABEL:eq:rec_z_nh) and (LABEL:eq:rec_rw_nh), $c^{\,(\pm)}_{n}=0$ for
negative values of $n$.

We now consider the asymptotic expansion of $X^{\,(\pm)}_{\ell m\omega}$. We
assume they can be factorized as,

$X^{\,(\pm)}_{\ell m\omega}=J^{\,(\pm)}_{\ell m\omega}(r)\,e^{+i\omega x}\,,$
(71)

where $J^{\,(\pm)}_{\ell m\omega}$ (“Jost function”) approaches a constant as
$x\to\infty$ in order to recover a purely outgoing behavior. We then
substitute Eq. (71) in the homogeneous equation (14). This results in a
differential equation for $J^{\,(\pm)}_{\ell m\omega}$,

$\left[f\frac{{\rm d}^{2}}{{\rm
d}r^{2}}+\left(\frac{2M}{r^{2}}+2i\omega\right)\frac{{\rm d}}{{\rm
d}r}-\frac{V^{\,(\pm)}_{\ell}}{f}\right]J^{\,(\pm)}_{\ell m\omega}=0\,,$ (72)

which we solve with Frobenius’ method. We assume a series expansion of the
form,

$J^{\,(\pm)}_{\ell m\omega}=\sum_{n=0}^{\infty}{a^{\,(\pm)}_{n}}/(\omega
r)^{n}\,.$ (73)

We then substitute this expansion in Eq. (72), and derive a recursion relation
between the coefficients $a^{\,(\pm)}_{n}$. When using the Zerilli potential
(11a), we obtain

$\displaystyle 2i\lambda^{2}na^{\,(+)}_{n}$
$\displaystyle=\lambda[\lambda(n-1)n-12i\sigma(n-1)-2\lambda(\lambda+1)]\,a^{\,(+)}_{n-1}$
$\displaystyle\quad+2\sigma[\lambda(3-\lambda)(n-2)(n-1)-(\lambda^{2}+9i\sigma)$
$\displaystyle\quad\times(n-2)-3\lambda^{2}]\,a^{\,(+)}_{n-2}+3\sigma^{2}[(3-4\lambda)(n-3)]$
$\displaystyle\quad\times(n-2)-4\lambda(n-3)-6\lambda]\,a^{\,(+)}_{n-3}$
$\displaystyle\quad-18\sigma^{3}(n-3)^{2}\,a^{\,(+)}_{n-4}\,,$ (74)

while using the Regge-Wheeler potential (11b) we find

$\displaystyle 2ina^{\,(-)}_{n}$
$\displaystyle=-2\sigma[(n+1)(n-3)]\,a^{\,(-)}_{n-1}$
$\displaystyle\quad-[\ell(\ell+1)-n(n-1)]\,a^{\,(-)}_{n-2}\,,$ (75)

where $a^{\,(\pm)}_{n}=0$ for negative values of $n$ [34, 60]. The derivations
we just made also apply to modes that behave as $\simeq\exp(-i\omega x)$ at
spatial infinity. In practice, we can replace $\sigma\to-\sigma$ in Eqs. (B)
and (B).

##  Appendix C Calculation of the Wronskian

To calculate the Wronskian (24), we need first to calculate the mode amplitude
$A_{\ell m\omega}^{\,(\pm)\,{\rm in}}$. We do this calculation following the
same strategy outlined in Ref. [29], Appendix A1, in the context of the
Sasaki-Nakamura formalism.

We first integrate $X^{\,(\pm),\,{\rm in}}_{\ell m\omega}$ from near the
horizon [making use of the recursion relations (LABEL:eq:rec_z_nh) and
(LABEL:eq:rec_rw_nh)] out to some large $r_{\rm max}$. This gives us two
constants $X^{\,(\pm),\,{\rm in}}_{\ell m\omega}(r_{\rm max})$ and ${\rm
d}X^{\,(\pm),\,{\rm in}}_{\ell m\omega}(r_{\rm max})/{\rm d}r$. For large $r$,
the mode function is approximately given by

$\displaystyle X^{\,(\pm),\,{\rm in}}_{\ell m\omega}$ $\displaystyle\simeq
A^{\,(\pm),\,{\rm in}}_{\ell m\omega}\,e^{-i\omega x}+A^{\,(\pm),\,{\rm
out}}_{\ell m\omega}\,e^{+i\omega x}$ $\displaystyle=A^{\,(\pm),\,{\rm
in}}_{\ell m\omega}\,J^{\,(\pm),\,{\rm in}}_{\ell m\omega}\,e^{-i\omega
x}+A^{\,(\pm),\,{\rm out}}_{\ell m\omega}\,J^{\,(\pm),\,{\rm out}}_{\ell
m\omega}\,e^{+i\omega x},$

where we used the factorization (71) in the second line.

We evaluate the Jost functions $J_{\ell m\omega}^{\,(\pm)}$ using the
recursion relations (B) and (B), with $a_{0}^{\,(\pm)}=1$ and adding terms in
the series until a sufficient level of accuracy is reached. Equation
(LABEL:eq:modeIn_jost) and its derivative with respect to $r$ provide two
equations that depend on $X^{\,(\pm),\,{\rm in}}_{\ell m\omega}$ (and its
derivative) at $r_{\rm max}$; the outcomes of the numerical integration of
$X^{\,(\pm),\,{\rm in}}_{\ell m\omega}$. We use these two equations to solve
for the two amplitudes $A_{\ell m\omega}^{\,(\pm),\,{\rm in}}$ and $A_{\ell
m\omega}^{\,(\pm),\,{\rm out}}$. The former is then used to calculate the
Wronskian (24).

## References

  * Rees  _et al._ [1974] M. Rees, R. Ruffini, and J. A. Wheeler,  _Black Holes, Gravitational Waves and Cosmology: An Introduction to Current Research_ (1974).
  * Smarr [1979] L. L. Smarr, ed.,  _Proceedings, Sources of Gravitational Radiation: Seattle, WA, USA, July 24 - August 4, 1978_ (Cambridge Univ. Press, Cambridge, 1979).
  * Regge and Wheeler [1957] T. Regge and J. A. Wheeler, Stability of a Schwarzschild singularity, Phys. Rev. 108, 1063 (1957).
  * Zerilli [1970a] F. J. Zerilli, Effective potential for even-parity Regge-Wheeler gravitational perturbation equations, Phys. Rev. Lett. 24, 737 (1970a).
  * Zerilli [1970b] F. J. Zerilli, Gravitational field of a particle falling in a Schwarzschild geometry analyzed in tensor harmonics, Phys. Rev. D 2, 2141 (1970b).
  * Davis  _et al._ [1971] M. Davis, R. Ruffini, W. H. Press, and R. H. Price, Gravitational radiation from a particle falling radially into a Schwarzschild black hole, Phys. Rev. Lett. 27, 1466 (1971).
  * Davis and Ruffini [1971] M. Davis and R. Ruffini, Gravitational radiation in the presence of a Schwarzschild black hole. A boundary value search, Lett. Nuovo Cim. 2S2, 1165 (1971).
  * Davis  _et al._ [1972] M. Davis, R. Ruffini, and J. Tiomno, Pulses of gravitational radiation of a particle falling radially into a Schwarzschild black hole, Phys. Rev. D 5, 2932 (1972).
  * Chung [1973] K. P. Chung, A four-dimensional Green’s function approach to the calculation of gravitational radiation from a particle falling into a black hole, Nuovo Cim. B 14, 293 (1973).
  * Ruffini [1973] R. Ruffini, Gravitational radiation from a mass projected into a Schwarzschild black hole, Phys. Rev. D 7, 972 (1973).
  * Teukolsky [1973] S. A. Teukolsky, Perturbations of a rotating black hole. 1. Fundamental equations for gravitational electromagnetic and neutrino field perturbations, Astrophys. J. 185, 635 (1973).
  * Detweiler and Szedenits [1979] S. L. Detweiler and E. Szedenits, Black holes and gravitational waves. II. Trajectories plunging into a nonrotating hole, Astrophys. J. 231, 211 (1979).
  * Tashiro and Ezawa [1981] Y. Tashiro and H. Ezawa, Gravitational Radiation of a Particle Falling Towards a Black Hole. I. The Case of Nonrotating Black Hole, Prog. Theor. Phys. 66, 1612 (1981).
  * Poisson [1997] E. Poisson, Gravitational radiation from infall into a black hole: Regularization of the Teukolsky equation, Phys. Rev. D 55, 639 (1997), arXiv:gr-qc/9606078 .
  * Campanelli and Lousto [1997] M. Campanelli and C. O. Lousto, Regularization of the Teukolsky equation for rotating black holes, Phys. Rev. D 56, 6363 (1997), arXiv:gr-qc/9707017 .
  * Sasaki and Nakamura [1981] M. Sasaki and T. Nakamura, The Regge-Wheeler equation with sources for both even and odd parity perturbations of the Schwarzschild geometry, Phys. Lett. A 87, 85 (1981).
  * Sasaki and Nakamura [1982a] M. Sasaki and T. Nakamura, A Class of New Perturbation Equations for the Kerr Geometry, Phys. Lett. A 89, 68 (1982a).
  * Sasaki and Nakamura [1982b] M. Sasaki and T. Nakamura, Gravitational Radiation From a Kerr Black Hole. 1. Formulation and a Method for Numerical Analysis, Prog. Theor. Phys. 67, 1788 (1982b).
  * Oohara and Nakamura [1983a] K. i. Oohara and T. Nakamura, Energy, momentum and angular momentum of gravitational radiation from a particle plunging into a nonrotating black hole, Phys. Lett. A 94, 349 (1983a).
  * Oohara and Nakamura [1983b] K. I. Oohara and T. Nakamura, Energy, momentum and angular momentum of gravitational waves induced by a particle plunging into a Schwarzschild black hole, Prog. Theor. Phys. 70, 757 (1983b).
  * Oohara and Nakamura [1984] K. I. Oohara and T. Nakamura, Gravitational waves from a particle scattered by a schwarzschild black hole, Prog. Theor. Phys. 71, 91 (1984).
  * Kojima and Nakamura [1983a] Y. Kojima and T. Nakamura, Gravitational radiation from a particle with zero orbital angular momentum plunging into a Kerr black hole, Phys. Lett. A 96, 335 (1983a).
  * Kojima and Nakamura [1983b] Y. Kojima and T. Nakamura, Gravitational radiation by a particle with nonzero orbital angular momentum falling into a Kerr black hole, Phys. Lett. A 99, 37 (1983b).
  * Sasaki and Nakamura [1990] M. Sasaki and T. Nakamura, Gravitational Radiation from Extreme Kerr Black Hole, Gen. Rel. Grav. 22, 1351 (1990).
  * 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).
  * Moncrief [1974] V. Moncrief, Gravitational perturbations of spherically symmetric systems. I. The exterior problem., Annals Phys. 88, 323 (1974).
  * Hopper and Cardoso [2018] S. Hopper and V. Cardoso, Scattering of point particles by black holes: gravitational radiation, Phys. Rev. D 97, 044031 (2018), arXiv:1706.02791 [gr-qc] .
  * Shibata and Nakamura [1992] M. Shibata and T. Nakamura, Metric perturbations induced by a particle falling into a Schwarzschild black hole 1: Formulation, Prog. Theor. Phys. 87, 1139 (1992).
  * Berti  _et al._ [2010] E. Berti, V. Cardoso, T. Hinderer, M. Lemos, F. Pretorius, U. Sperhake, and N. Yunes, Semianalytical estimates of scattering thresholds and gravitational radiation in ultrarelativistic black hole encounters, Phys. Rev. D 81, 104048 (2010), arXiv:1003.0812 [gr-qc] .
  * Ruffini and Sasaki [1981] R. Ruffini and M. Sasaki, On a Semirelativistic Treatment of the Gravitational Radiation From a Mass Thrusted Into a Black Hole, Prog. Theor. Phys. 66, 1627 (1981).
  * Berti  _et al._ [2009] E. Berti, V. Cardoso, L. Gualtieri, F. Pretorius, and U. Sperhake, Comment on ‘Kerr Black Holes as Particle Accelerators to Arbitrarily High Energy’, Phys. Rev. Lett. 103, 239001 (2009), arXiv:0911.2243 [gr-qc] .
  * Sago  _et al._ [2003] N. Sago, H. Nakano, and M. Sasaki, Gauge problem in the gravitational self-force. 1. Harmonic gauge approach in the Schwarzschild background, Phys. Rev. D 67, 104017 (2003), arXiv:gr-qc/0208060 .
  * Nakano  _et al._ [2003] H. Nakano, N. Sago, and M. Sasaki, Gauge problem in the gravitational self-force. 2. First post-Newtonian force under Regge-Wheeler gauge, Phys. Rev. D 68, 124003 (2003), arXiv:gr-qc/0308027 .
  * Chandrasekhar and Detweiler [1975] S. Chandrasekhar and S. L. Detweiler, The quasi-normal modes of the Schwarzschild black hole, Proc. Roy. Soc. Lond. A 344, 441 (1975).
  * [35] Black Hole Perturbation Toolkit, (bhptoolkit.org).
  * Warburton and Wardell [2014] N. Warburton and B. Wardell, Applying the effective-source approach to frequency-domain self-force calculations, Phys. Rev. D 89, 044046 (2014), arXiv:1311.3104 [gr-qc] .
  * Leather and Warburton [2023] B. Leather and N. Warburton, Applying the effective-source approach to frequency-domain self-force calculations for eccentric orbits, Phys. Rev. D 108, 084045 (2023), arXiv:2306.17221 [gr-qc] .
  * Hopper [2018] S. Hopper, Unbound motion on a Schwarzschild background: Practical approaches to frequency domain computations, Phys. Rev. D 97, 064007 (2018), arXiv:1706.05455 [gr-qc] .
  * Chandrasekhar [1985] S. Chandrasekhar,  _The mathematical theory of black holes_ (1985).
  * Leaver [1985] E. W. Leaver, An Analytic representation for the quasi normal modes of Kerr black holes, Proc. Roy. Soc. Lond. A 402, 285 (1985).
  * [41] https://pages.jh.edu/eberti2/ringdown.
  * Press [1971] W. H. Press, Long Wave Trains of Gravitational Waves from a Vibrating Black Hole, Astrophys. J. Lett. 170, L105 (1971).
  * Goebel [1972] C. J. Goebel, Comments on the “vibrations” of a black hole., Astrophys. J. 172, L95 (1972).
  * Cardoso  _et al._ [2009] V. Cardoso, A. S. Miranda, E. Berti, H. Witek, and V. T. Zanchin, Geodesic stability, Lyapunov exponents and quasinormal modes, Phys. Rev. D 79, 064016 (2009), arXiv:0812.1806 [hep-th] .
  * 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).
  * Vishveshwara [1970] C. V. Vishveshwara, Scattering of Gravitational Radiation by a Schwarzschild Black-hole, Nature 227, 936 (1970).
  * Hadar  _et al._ [2011] S. Hadar, B. Kol, E. Berti, and V. Cardoso, Comparing numerical and analytical calculations of post-ISCO ringdown amplitudes, Phys. Rev. D 84, 047501 (2011), arXiv:1105.3861 [gr-qc] .
  * Folacci and Ould El Hadj [2018] A. Folacci and M. Ould El Hadj, Multipolar gravitational waveforms and ringdowns generated during the plunge from the innermost stable circular orbit into a Schwarzschild black hole, Phys. Rev. D 98, 084008 (2018), arXiv:1806.01577 [gr-qc] .
  * Hadar and Kol [2011] S. Hadar and B. Kol, Post-ISCO Ringdown Amplitudes in Extreme Mass Ratio Inspiral, Phys. Rev. D 84, 044019 (2011), arXiv:0911.3899 [gr-qc] .
  * Moreno-Garrido  _et al._ [1995] C. Moreno-Garrido, E. Mediavilla, and J. Buitrago, Gravitational radiation from point masses in elliptical orbits: spectral analysis and orbital parameters, MNRAS 274, 115 (1995).
  * Glampedakis  _et al._ [2002] K. Glampedakis, S. A. Hughes, and D. Kennefick, Approximating the inspiral of test bodies into Kerr black holes, Phys. Rev. D 66, 064005 (2002), arXiv:gr-qc/0205033 .
  * Babak  _et al._ [2007] S. Babak, H. Fang, J. R. Gair, K. Glampedakis, and S. A. Hughes, ‘Kludge’ gravitational waveforms for a test-body orbiting a Kerr black hole, Phys. Rev. D 75, 024005 (2007), [Erratum: Phys.Rev.D 77, 04990 (2008)], arXiv:gr-qc/0607007 .
  * Maggiore [2007] M. Maggiore,  _Gravitational waves: Volume 1: Theory and experiments_ (OUP Oxford, 2007).
  * Ferrari and Ruffini [1981] V. Ferrari and R. Ruffini, On the Structure of Gravitational Wave Bursts: Implosion With Finite Kinetic Energy, Phys. Lett. B 98, 381 (1981).
  * Borrelli [1997] A. Borrelli, Gravitational radiation emitted when a mass falls onto a compact star, Nuovo Cim. B 112, 225 (1997).
  * Ferrari  _et al._ [1999] V. Ferrari, L. Gualtieri, and A. Borrelli, Stellar pulsations excited by a scattered mass, Phys. Rev. D 59, 124020 (1999), arXiv:gr-qc/9901060 .
  * Tominaga  _et al._ [1999] K. Tominaga, M. Saijo, and K.-i. Maeda, Gravitational waves from a test particle scattered by a neutron star: Axial mode case, Phys. Rev. D 60, 024004 (1999), arXiv:gr-qc/9901040 .
  * Ruoff  _et al._ [2001] J. Ruoff, P. Laguna, and J. Pullin, Excitation of neutron star oscillations by an orbiting particle, Phys. Rev. D 63, 064019 (2001), arXiv:gr-qc/0005002 .
  * Hopper and Evans [2010] S. Hopper and C. R. Evans, Gravitational perturbations and metric reconstruction: Method of extended homogeneous solutions applied to eccentric orbits on a Schwarzschild black hole, Phys. Rev. D 82, 084010 (2010), arXiv:1006.4907 [gr-qc] .