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

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]. 1 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.
We assume that the particle starts from rest at infinity with (conserved) energy E = E/µ = 1 and angular momentum L = L/µ per unit mass.If we chose, without loss of generality, that the particle's motion happens in the equatorial plane θ = π/2, then we can parametrize the particle's worldline in terms of the proper time τ as z µ (τ ) = {t p , r p , π/2, ϕ p }, and, from the geodesic equation and the timelike constraint g ab u a u b = −1, with u a = dz a /dt, we obtain: where ˙= d/dτ and is the effective radial potential.The particle's trajectories are classified according to the number of real roots of E 2 − U .For the special case E = 1 the analysis is simple, and we find the roots to be, Hence, a plunging orbit from infinity requires that L < 4M , for otherwise the turning points are real and positive The effective potential U (r, L) for various values of 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 (E = 1) into the hole, there is a critical value of Lcrit = 4M (solid black line) for which the particle becomes trapped in a marginally stable circular orbit at r = 4M .This value of L separates plunging (L < Lcrit) from scattering (L > Lcrit) geodesics.
(i.e., the particle is scattered.)We will define this special value of the angular momentum as L crit = 4M .
In Fig. 1 we show the effective potential (3) for a range of values L/M ∈ {3.25, 4.25} (light curves).The thicker line corresponds to U (r, L crit ), which peaks at r = 4M with value of 1.Hence, a particle falling from rest and with angular momentum L = L crit will be captured in a marginally stable circular orbit.A particle with L larger (smaller) than L crit will scatter (plunge).
It is convenient to rewrite Eq. ( 2) as first-order in r equations, where we have taken ṙp < 0. We integrate Eqs. ( 5) with initial conditions t p (r max ) = 0 and ϕ p (r max ) = 0 at some arbitrarily large r p = r max down to the horizon, r p = 2M .In Fig. 2 we show a sequence of trajectories starting from L/M = 3.25 and up to L/M = 3.9996 (i.e., with 99.99% of L crit ).We translate from Schwarzschild-Droste to Cartesian coordinates using As the ratio L/L crit approaches unity from below, the particle executes an increasing fractional number of or- x/M In the limit L/Lcrit → 1, the particle executes an increasing fractional number of orbits ϕ/(2π), as seen in the red curve that corresponds to L = 3.9996M , or 99.99% of Lcrit. .bits ϕ p /(2π), given approximately by −( √ 2π) −1 log(1 − L/L 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: ℓm (t, r) , (7) or, by going to the Fourier domain through we have alternatively In these equations, x is the tortoise coordinate that maps the region 2M < r < ∞ to −∞ < x < ∞.The superscript (±) denotes variables associated to metric perturbations of polar (+) or axial (−) parity, and V (±) ℓ 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 where we defined The function S (±) ℓmω 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 (±) ℓ 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 → ±∞.We will impose that X (±) ℓmω is purely ingoing at the event horizon and purely outgoing at spatial infinity, that is, 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 (±), in ℓmω and X (±), up ℓmω , of the homogeneous equations The two solutions differ by the boundary conditions we impose when we solve Eq. ( 14), namely, x → +∞ .
We use these solutions to define the Wronskian, which is constant in x.We then construct the Green's function as where Θ(•) is the Heaviside step function.Finally, with the Green's function (16), we can write the solution of Eq. (7) as Equation (17) has the desired property of being purely ingoing at the horizon and purely outgoing at infinity: where we used the boundary conditions on X (±), up/in ℓmω .
We can compare Eq. (18b) with the boundary conditions Eq. ( 13) imposed on X (±) ℓmω , and make the identification, 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 (ℓ, m) mode and in the time-domain Regge-Wheeler and Zerilli mode functions In practice, we note that as x → ∞, the Wronskian becomes and that we can also rewrite Eq. (19) as Because f ≃ 1 and X (±), in ℓmω ≃ exp(±iωx) as x → ∞, we see that the convergence of this integral depends critically on the asymptotic properties of S (±) ℓmω .

IV. ASYMPTOTIC BEHAVIOR OF THE FREQUENCY DOMAIN SOURCES
We now review the properties of ℓmω in two cases of interest.We will start by reviewing the case in which the particle falls radially, L = 0, into the black hole and how Ref. [28] improved the asymptotic behavior of S (+) ℓmω .We will then show one way in which the Nakamura-Sasaki-Shibata transformation can be generalized to general plunging trajectories, L ̸ = 0, and discuss the asymptotic properties of this new source term.

A. The case of radial infall
When L = 0, S (−) ℓmω vanishes and S (+) ℓmω simplifies to [6,7] where A ℓm = Y * ℓm (π/2, ϕ) exp(imϕ), Y ℓm are the spherical harmonics, the asterisk indicates complex conjugation, and t p is given by We find that the near-horizon and spatial-infinity behaviors of S (+) Hence the integral Eq. ( 19) converges slowly at spatial infinity.To improve the convergence, Ref. [28] proposed the substitution2 where The function Q ℓmω vanishes at the event horizon x → −∞ and decays as x −3/2 for x → ∞.Thus, X ℓmω ≃ X (+) ℓmω in the latter limit.
We then insert Eq. ( 27) in Eq. ( 7) to find where S (+) ℓmω is a new source term given by The asymptotic behaviors of the new source term are making the integral in Eq. ( 19) converge faster.One can verify that it is the second x derivative of Q ℓmω in Eq. ( 30) that yields a term that decays as r −1/2 as x → ∞, and that cancels the leading-order asymptotic term in the large-x expansion of the original source term.We note that both X  ) , and its Nakamura-Sasaki-Shibata-transformed version, S (+) .Right panels: the real ("Re") and imaginary ("Im") parts of S (+) (top) and S (+) (bottom).Both sources vanish at the event horizon r = 2M , and the faster asymptotic decay of S (+) improves the convergence of Eq. ( 19).The source terms are largest near r ≃ 3M , which corresponds approximately to the location of the peak of the Zerilli potential r ℓ=2, peak ≃ 3.1M and of the light ring r = 3M .
of the Zerilli equation for a particle falling radially into the black hole, and ℓ = m = 2 and M ω = 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. ( 30) accelerates the convergence of Eq. (19).The value M ω = 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 ≃ 3M , which corresponds approximately to the location of peak of the Zerilli potential r (+) ℓ=2, peak ≃ 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. ( 27) and (28) to the case in which the particle plunges with a nonzero angular momentum?In the next section we give positive answers to both questions.

B. The case of general plunging trajectories
We seek a minimal extension of Eqs. ( 27) and (28) for the case of a plunging particle with general angular momentum L ̸ = 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: 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.
Our strategy is to retain Eq. ( 27), but generalize Q ℓmω to the form Proceeding as in the previous section, we arrive at the same Eq.( 30), where S (+) ℓmω is now given by a more complicated formula, whose form can be found in Eq. (A1), Appendix A. We then expand Eq. ( 30) in a power series in r for r → ∞, 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 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 L at order N = 0 and that we recover Eq. ( 28) 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  ) , and its Nakamura-Sasaki-Shibata-transformed version, S (+), (N ) , for N = 0, 1, and 2 in the expansion (33).Right panels: the real ("Re") and imaginary ("Im") parts of S (+) (top) and S (+), (N ) (three remaining panels) for increasing values of N .The source term is largest near r ≃ 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 L/Lcrit → 1.
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 L ̸ = 0, when compared to the radial infall case, r −3/2 .This can be understood by examining the asymptotic behavior of S (+) ℓmω : This expansion shows that the first term containing the angular momentum L appears at order r −1 , and it cannot be canceled with Eq. ( 28).

V. NUMERICAL METHODS AND RESULTS
In this section we describe the numerical methods we use to compute Eq. ( 23), and show a few applications of solving this equation using the generalized Nakamura-Sasaki-Shibata transformation (33).

A. Numerical methods
To evaluate our main quantity of interest, namely the amplitude C (±), out ℓmω we proceed as follows.
2. Choose a value of the angular momentum L, and solve the geodesic equations (5), to obtain t p (r) and ϕ p (r).We use initial conditions t p (r min ) = 0 and ϕ p (r min ) = 0, and integrate from r min = 2 (1 + 10 −4 )M up to r max = 500/ω.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. ( 23) as a differential equation [36,37]: We integrate Eq. ( 36), with initial condition C (±), out ℓmω (r min ) = 0 up to r max , using the Regge-Wheeler source (A6) or the Nakamura-Sasaki-Shibata-transformed Zerilli source (30), depending on whether ℓ + m is odd or even, respectively.In the context of Eq. ( 36), the new Zerilli source term reduces the value of r max at which dC (±), out ℓmω /dr 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 ≈ 2M .This causes the factor exp(iωt p ), which appears in the source, to become highly oscillatory, and thus sensitive to our choice of r min and ω.Consequently, the phase (but not the amplitude) of the C (±), out ℓmω 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 L = {0, 0.25, 0.5, 0.75, 0.9, 0.99, 0.999, 0.9999}L crit .

B. Energy spectrum
We first consider the energy spectrum for the particle in the near-critical limit L ≈ L 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, L/L crit = 0.9 and 0.9999, which were examined in detail by Berti et al. [29].
In Fig. 5, we show the energy spectrum (20) for L = 0.9 L crit (left panels) and L = 0.9999 L crit (right panels).The top panels show the spectral energy for multipoles ℓ = m ranging from 2 to 6 (different line colors) and their sum (dashed line).The bottom panels show the quadrupole mode ℓ = 2 and all azimuthal contributions |m| ≤ 2 (different line styles.)Our results are in excellent agreement with those of Ref. [29]; cf.Figs. 10 and 11 therein.
For L = 0.9 L crit , the energy radiated in each multipole ℓ = m has a single maximum.The location M ω of these peaks coincide approximately with the real part of the fundamental (n = 0) quasinormal mode frequency of a Schwarzschild black hole M ω ℓm0 , as first observed by Detweiler and Szedenits [12].For reference, these values are Re[M ω ℓm0 ] ≃ {0.373, 0.599, 0.809, 1.012, 1.212}, for ℓ ∈ {2, 6} [39][40][41].(Due to the spherical symmetry of the Schwarzschild solution, all quasinormal modes with |m| ≤ ℓ, at fixed ℓ ≥ 2, are degenerate.) In the near-critical limit L = 0.9999 L crit , the spectral energy distribution has a peak at M ω < Re[M ω ℓm0 ].This peak corresponds to m times the particle's orbital frequency M Ω orb = 1/8 at the marginally stable circular orbit at r = 4M (cf.Fig. 2), that is 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 ≃ 4M , as we discussed in Fig. 4. For moderate values of ℓ, the contributions to the energy driven by the particle's orbital motion and the quasinormal mode excitation overlap, while as ℓ ≫ 1 their contributions separate sufficiently to result in two peaks in the energy spectra; see the ℓ = m = 5 and 6 in Fig. 5.For ℓ = m = 5 the peaks are located at M ω ≃ 0.637 and 0.919, while for ℓ = m = 6 they are found at M ω ≃ 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 where is the orbital frequency of a null geodesic at the light ring.Equation (38) can also be obtained from the ℓ ≫ 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 "quasinormalmode" and "geodesic" peaks by taking the difference of Eqs. ( 37) and ( 38): which is valid for ℓ ≫ 1, L/L crit ≃ 1, and E = 1.For ℓ = m = 6, we find (∆M ω) peak ≃ 0.404, in fair agreement with the numerical result ≃ 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 E → ∞.In this limit, the "geodesic" peak moves rightwards, eventually overlapping with the "quasinormal-mode" peak; see Ref. [29], Fig. 10.This occurs because as E → ∞, the particle's marginally stable circular orbit coincides with the light ring, hence Ω orb ≃ Ω LR [31].In this limit, we then have which is valid for ℓ ≫ 1, L/L crit ≃ 1, and E → ∞.The peak separation vanishes when ℓ = m, reproducing the results of Ref. [29].
The energy spectra for a particle plunging with angular momentum L = 0.9 Lcrit (left panels) and L = 0.9999 Lcrit (right panels) into a Schwarzschild black hole.Top panels: the energy spectra from the multipoles ℓ = m from 2 to 6, (colored lines), and their sum (dashed line).Bottom panels: the energy spectra for the quadrupole perturbation ℓ = 2 and all |m| ≤ 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.

C. Total energy
From the energy spectrum, we can compute the total energy ∆E emitted in form of gravitational waves by integrating the spectrum density and summing the contributions from all multipoles: To understand the relative contribution due to perturbations of each parity, we also define ∆E (±) , where the subscript (+) means we add only the contributions to the energy coming from multipoles for which ℓ + m is even and (−) when ℓ + 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 L = 0, proposed the empirical relation between energy and multipole ℓ.We fit Eq. ( 42) to the outcome of computing ∆E ℓ , truncating the ℓ sum at ℓ max = 6.Table I 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 L approaches L crit .In other words, the fraction of the total energy radiated contained in each ℓ pole tends to "even out" as L increases.Figure 6 illustrates these observations.For example, we see in the leftmost panel that ∆E 2 varies from being about four orders of magnitude larger than ∆E 6 , to being some more than ten times larger as 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 L; the extreme case is, of course, when L = 0 where ∆E (−) vanishes.
Interestingly, the axial and polar contributions to the total energy ∆E are nonmonotonic with respect to the angular momentum.In Fig. 7, we show the ratio ∆E (±) /∆E as a function of the logarithm of 1 − L/L crit .As we in-  Regge-Wheeler FIG. 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. ( 41), whereas the straight lines corresponds to the fit (42) suggested by Oohara and Nakamura [20].The lines starting from the bottom correspond to L = 0 and increase to L = 0.9999 Lcrit 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 L. When L = 0, there is no energy radiated in the axial-parity channel.
crease the particle's angular momentum (i.e., moving right to left along the abscissa), we see that ∆E (−) /∆E has a local maximum at L ≈ 0.75 L crit , yet with only approximately 14% of the total energy budget.We interpret this result as follows.When L = 0, by symmetry, all energy must be radiated through the polar channel.As we increase L, axial perturbations become increasingly excited, and a nonzero (albeit small) percentage of the energy is emitted through them.As L approaches L 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.

D. Time-domain waveforms
We also computed the time-domain Regge-Wheeler and Zerilli mode functions using Eq. ( 21).In Fig. 8, we show the dominant modes for the Regge-Wheeler and Zerilli functions, (ℓ, 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 L = 0.9 L crit , while the bottom panel to L = 0.9999 L 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 L → L crit , the large values of 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 L → L crit , for the same multipole moments.In this limit, the waveforms become increasingly similar around their peak amplitudes, t − x ≃ −40 M , although differences are clearly visible at the "precursor" phase.Similar results hold for the other multipoles we examined.
FIG. 8. Time-domain Zerilli (solid lines) and Regge-Wheeler (dashed lines) mode functions excited by a test particle plunging with L = 0.9 Lcrit (top panel) and L = 0.9999 Lcrit (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.
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].

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 ωM ≳ 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 where the (mass) quadrupole moment M ij is defined as for a mass density ρ(t, 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 where the luminosity is to be evaluated without any averaging (also, we can set t max = ∞ at the particle's horizon crossing time).The same quantity can be evaluated in the frequency domain where dE/dω is the spectral energy distribution.For a point particle moving along a trajectory x p (t) we have ρ(t, x) = µ δ (3) [x − x p (t)], and we find 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 (44) are M 11 , M 22 , M 12 , and the trace is Then, a short calculation leads to where the Fourier transform of M ij is defined as The fact that M ij is real implies the useful property As discussed in Ref. [54], the integral ( 49) is divergent at t → −∞, since x p → +∞.Therefore, some regularization procedure is required.This is achieved by working with the Fourier transform Mij , which is well behaved at spatial infinity.In fact, this procedure is equivalent to the regularization of Mij (ω) 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 and from Eq. ( 48) we read off the regularized formula: 3 For clarity, we adopt a slightly different notation for frequencydomain quantities in this section, mirroring Ref. [54].
For the actual numerical evaluation of Mij (and the subsequent one of dE/dω), it is advantageous to convert the time integral into a radial integral using the geodesic equations ( 2), leading to r) .(53) For the problem at hand, the individual Mij components required for the calculation of Eq. ( 51) are With the help of Eq. ( 2) both the accelerations ẍp , ÿp and the velocities ẋp , ẏp can be rewritten solely as functions of r, ϕ p (r) and the orbital constant L. At the same time t p (r) explicitly appears in the exponential term.Both ϕ p (r) and t p (r) are obtained numerically as in Sec.II.We chose trajectories with ϕ p (0) = 0 at some initial radius r max = r(0).The same radius serves as the "infinity" upper limit in the Mij 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. ( 53).
For any value L < L crit , the integral in the M11 component is the slowest converging one at large r (the absolute value of the integrand decays as ∼ 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 B, namely, we subtract the slowly decaying part prior to the numerical integration and then add back the (analytically obtained) asymptotic contribution from that part.
The semirelativistic energy spectrum, calculated from Eqs. ( 52) and (53), is shown in Fig. 10 and has a characteristic single-hump profile.When comparing against the full black-hole perturbation theory result for ℓ = 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, ωM ≪ 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 L 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, ωM ≫ 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 L and the , where the radiation is due to the particle's quasi-Newtonian motion.For L = 0 and L = 0.75 Lcrit, the kludge calculation underestimates the location of the spectrum's peak, which is related to the hole's fundamental quasinormal mode frequency [12].As L → Lcrit, 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 ω ≲ 0.25.The semirelativistic approximation predicts a slowly decaying tail for the spectrum.
modal numbers ℓ and m.We can approximate the highfrequency "tail" of dE ℓm /dω 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. (23).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 ℓmω ≃ 2iω.These manipulations lead to a scaling The same procedure applied to the semirelativistic spectrum leads to a slower decaying tail The difference can be traced back to the additional highly oscillatory function X (±), in ℓmω ≃ exp(iωx) in the full blackhole perturbation theory expression, which suppresses the integral at large ω.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 L from spatial infinity into a Schwarzschild black hole.In particular, we focused on the limit in which L approximates from below the critical value L 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 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. ( 28) and (30)].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 plungingparticle problem in the semirelativistic approximation.We found that the two energy spectra agree best with in the limit L → L 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].
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.
(r ≫ 2M or x → ∞) 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 We now consider the asymptotic expansion of X (±) ℓmω .We assume they can be factorized as, where J (±) ℓmω ("Jost function") approaches a constant as x → ∞ in order to recover a purely outgoing behavior.We then substitute Eq. (B5) in the homogeneous equation (14).This results in a differential equation for J where a (±) n = 0 for negative values of n [34,60].The derivations we just made also apply to modes that behave as ≃ exp(−iωx) at spatial infinity.In practice, we can replace σ → −σ in Eqs.(B8) and (B9).from near the horizon [making use of the recursion relations (B2) and (B4)] out to some large r max .This gives us two constants X (±), in ℓmω (r max ) and dX where we used the factorization (B5) in the second line.
We evaluate the Jost functions J (±) ℓmω using the recursion relations (B8) and (B9), with a (±) 0 = 1 and adding terms in the series until a sufficient level of accuracy is reached.Equation (C1) and its derivative with respect to r provide two equations that depend on X (±), in ℓmω (and its derivative) at r max ; the outcomes of the numerical integration of X (±), in ℓmω .We use these two equations to solve for the two amplitudes A (±), in ℓmω and A (±), out ℓmω .The former is then used to calculate the Wronskian (22).

FIG. 2 . 1 .
FIG.2.Timelike geodesics with angular momentum per unit mass L/M ∈ {3.25, 3.9996} plunging into a Schwarzschild black hole (black disk) starting from rest at spatial infinity, 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 L/Lcrit → 1, the particle executes an increasing fractional number of orbits ϕ/(2π), as seen in the red curve that corresponds to L = 3.9996M , or 99.99% of Lcrit..

FIG. 3 .
FIG.3.The ℓ = m = 2 source term for the Zerilli equation for a radially plunging particle starting from rest at spatial infinity, for M ω = 0.3.Left panel: the absolute values of the original Zerilli source term, S (+) , and its Nakamura-Sasaki-Shibata-transformed version, S (+) .Right panels: the real ("Re") and imaginary ("Im") parts of S (+) (top) and S (+) (bottom).Both sources vanish at the event horizon r = 2M , and the faster asymptotic decay of S (+) improves the convergence of Eq.(19).The source terms are largest near r ≃ 3M , which corresponds approximately to the location of the peak of the Zerilli potential r

FIG. 4 .
FIG.4.The (ℓ, m) = (2, 2) source term for the Zerilli equation for a plunging particle with L = 0.9999 Lcrit starting from rest at spatial infinity, for M ω = 0.3.Left panel: the absolute values of the original Zerilli source term, S (+) , and its Nakamura-Sasaki-Shibata-transformed version, S (+), (N ) , for N = 0, 1, and 2 in the expansion(33).Right panels: the real ("Re") and imaginary ("Im") parts of S (+) (top) and S (+), (N ) (three remaining panels) for increasing values of N .The source term is largest near r ≃ 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 L/Lcrit → 1.
to evaluate the integral in Eq.(23).

FIG. 7 .
FIG.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 L, with the axial channel having a peak at L/Lcrit ≈ 0.75.The polar contribution is at least ≈ 85% of the total energy budget for all values of L considered.

FIG. 9 .
FIG.9.Time-domain Zerilli (top panel) and Regge-Wheeler (bottom panel) mode functions excited by a test particle plunging with L = 0.99 Lcrit (solid lines), L = 0.999 Lcrit (dashed lines) and L = 0.9999 Lcrit (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 L → Lcrit, the solutions become quasiuniversal around and after t − x ≳ −40 M .

FIG. 10 .
FIG.10.The energy spectra for a particle plunging with angular momentum L = 0, L = 0.75 Lcrit, and L = 0.9999 Lcrit into a Schwarzschild black hole in the semirelativistic approximation ("SR," dashed lines) and black-hole perturbation theory, for ℓ = m = 2 ("BHPT," solid lines).Both calculations agree at low frequencies M ω ≪ 1, where the radiation is due to the particle's quasi-Newtonian motion.For L = 0 and L = 0.75 Lcrit, the kludge calculation underestimates the location of the spectrum's peak, which is related to the hole's fundamental quasinormal mode frequency[12].As L → Lcrit, 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 ω ≲ 0.25.The semirelativistic approximation predicts a slowly decaying tail for the spectrum.

(
±), in ℓmω (r max )/dr.For large r, the mode function is approximately given by X

TABLE I .
Coefficients a and b in the Oohara-Nakamura relation (42) and the total energy (µ 2 /M ) ∆E for different values of angular momentum L. The superscripts refer to polar (+) and axial-parity (−) perturbations.