Excitation of Kerr quasinormal modes in extreme--mass-ratio inspirals

TL;DR

This paper investigates the excitation of Kerr quasinormal modes during extreme-mass-ratio inspirals, revealing that these wiggles are a common, observable phenomenon that depends on black hole spin and orbital parameters, with implications for testing general relativity.

Contribution

It provides a comprehensive survey of QNM excitation phenomena in EMRIs, highlighting their observability and dependence on system parameters, and discusses their potential for testing gravity.

Findings

Wiggles are a generic feature across a wide parameter space.

QNM excitations are most prominent with high black hole spin and eccentric prograde orbits.

Wiggles could enable direct observation of Kerr QNMs in gravitational wave signals.

Abstract

If a small compact object orbits a black hole, it is known that it can excite the black hole's quasinormal modes (QNMs), leading to high-frequency oscillations ("wiggles") in the radiated field at $\mathcal{J}^+$, and in the radiation-reaction self-force acting on the object after its orbit passes through periapsis. Here we survey the phenomenology of these wiggles across a range of black hole spins and equatorial orbits. In both the scalar-field and gravitational cases we find that wiggles are a generic feature across a wide range of parameter space, and that they are observable in field perturbations at fixed spatial positions, in the self-force, and in radiated fields at $\mathcal{J}^+$. For a given charge or mass of the small body, the QNM excitations have the highest amplitudes for systems with a highly spinning central black hole, a prograde orbit with high eccentricity, and an orbital periapsis close to the light ring. The QNM amplitudes depend smoothly on the orbital parameters, with only very small amplitude changes when the orbit's (discrete) frequency spectrum is tuned to match QNM frequencies. The association of wiggles with QNM excitations suggest that they represent a situation where the \emph{nonlocal} nature of the self-force is particularly apparent, with the wiggles arising as result of QNM excitation by the compact object near periapsis, and then encountered later in the orbit. Astrophysically, the effects of wiggles at $\mathcal{J}^+$ might allow direct observation of Kerr QNMs in extreme-mass-ratio inspiral (EMRI) binary black hole systems, potentially enabling new tests of general relativity.

Figures

Tables

Table 1. Table 1: This table shows the local- and radiated-field diagnostics in which we study wiggles.

Field perturbation at …	Scalar field	Gravitational
fixed spatial position	$.\left({\phi }_m\right){.}_{\text{regularized}}$	${\Phi }_{-2}^{\left(𝔩m\right)}$
(strong-field)
particle position	$F_a^{\left(m\right)}$	${}_{-2}{}^{}F_a^{\left(𝔩m\right)}$
${𝒥}^{+}$	${\partial }_{tt}(.\left({\phi }_m\right){.}_{\text{physical}})$	${\Psi }_4$

Table 2. Table 2: This table describes the Kerr scalar-field configurations surveyed here. The table rows are ordered by decreasing Kerr spin  a ~ ~ 𝑎 \tilde{a} (horizontal lines separate groups of rows for the same spin), within this by increasing particle periapsis radius r min subscript 𝑟 r_{\min} , and within this by decreasing particle orbital eccentricity  e 𝑒 e . ( a ~ , p , e ) ~ 𝑎 𝑝 𝑒 (\tilde{a},p,e) describe the configuration. T r subscript 𝑇 𝑟 T_{r} is the coordinate-time period of the orbit’s radial motion, r min subscript 𝑟 r_{\min} is the particle’s periapsis r 𝑟 r  coordinate, and ϕ ˙ periapsis subscript ˙ italic-ϕ periapsis \dot{\phi}_{\text{periapsis}} is the particle’s angular 3-velocity in ϕ italic-ϕ \phi (i.e., d ​ ϕ / d ​ t 𝑑 italic-ϕ 𝑑 𝑡 d\phi/dt ) at the point of periapsis. r watchpoint subscript 𝑟 watchpoint r_{\text{watchpoint}} is the r 𝑟 r  coordinate of the watchpoint used for observing wiggles, or blank if no watchpoint data was available for this configuration. “Resolution” is the highest numerical resolution used for simulations of this configuration, and refers to the grid structures described in table  4 . The final sets of columns describe the presence or absence of wiggles in our field diagnostics for m = 1 𝑚 1 m=1 through m = 6 𝑚 6 m=6 ; for each m 𝑚 m there are 3 columns labeled “w”, “ F 𝐹 F ”, and “ 𝒥 𝒥 \mathcal{J} ”, describing (respectively) wiggles in the regularized field . ( φ m ) . regularized \bigl{.}(\varphi_{m})\bigr{.}_{\text{regularized}} at a fixed “watchpoint” (located on the equator at r = r watchpoint 𝑟 subscript 𝑟 watchpoint r=r_{\text{watchpoint}} ), wiggles in the self-force F a subscript 𝐹 𝑎 F_{a} , and wiggles in the second time derivative of the physical field, ∂ t ​ t ( . ( φ m ) . physical ) \partial_{tt}\left(\bigl{.}(\varphi_{m})\bigr{.}_{\text{physical}}\right) , at 𝒥 + superscript 𝒥 \mathcal{J}^{+} . The meanings of the symbols in these columns are described in detail in table  3 ; briefly, ⊕ direct-sum \bm{\oplus} or + \bm{+} means that we observed wiggles and fit them to the appropriate model described in section  IV.1 , ∼ similar-to \bm{\sim}  means that we observed wiggle-like oscillations but did not attempt a quantitative fit, ♢ bold-♢ \bm{\diamondsuit}  means that we observed oscillations which might be wiggles but were not clearly separated from the background variation, and × \bm{\times}  means that we did not observe wiggles (n.b. this does not prove the absence of wiggles, only that we did not observe them).

Name

$\tilde{a}$

$p$

$e$

$T_r$

$r_{\min }$

${\dot{\varphi }}_{\text{periapsis}}$

$r_{\text{watchpoint}}$

Resolution

m=1

m=2

m=3

m=4

m=5

m=6

($M$)

($M$)

($M^{-1}$)

($M$)

w

$F$

$𝒥$

w

$F$

$𝒥$

w

$F$

$𝒥$

w

$F$

$𝒥$

w

$F$

$𝒥$

w

$F$

$𝒥$

w9x5-161

$\mathrm{0.999\hspace{0.17em}99}$

$\mathrm{2.918\hspace{0.17em}315}$

$\mathrm{0.807\hspace{0.17em}519}$

$230.442$

$1.615$

0.344 414

$11.810$

dro12-96

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

w9x4-368

$\mathrm{0.999\hspace{0.17em}9}$

$7.0$

$0.9$

$1513.112$

$3.684$

0.145 444

$46.304$

dro8-64

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\times$

$\times$

$\times$

$\times$

$\times$

w999-278

$0.999$

$5.0$

$0.8$

$400.508$

$2.778$

0.199 135

$21.448$

dro10-80

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

w999-368

$0.999$

$7.0$

$0.9$

$1513.179$

$3.684$

0.145 443

dro8-64

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\sim$

$\times$

$\times$

$\times$

$\times$

ze98a

$0.99$

$\mathrm{2.398\hspace{0.17em}1}$

$0.98$

$3414.259$

$1.211$

0.430 498

$61.223$

dro6-48

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

ze98

$0.99$

$2.4$

$0.98$

$3304.620$

$1.212$

0.430 300

dro10-80

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\sim$

w99-125a

$0.99$

$\mathrm{2.437\hspace{0.17em}5}$

$0.95$

$957.757$

$1.25$

0.421 924

$36.314$

dro8-64

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\sim$

$\sim$

w99-125b

$0.99$

$2.375$

$0.9$

$432.084$

$1.25$

0.421 555

$20.543$

dro8-64

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\sim$

w99-125c

$0.99$

$2.25$

$0.8$

$247.845$

$1.25$

0.420 811

$8.595$

dro6-48

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

w99-125d

$0.99$

$2.0$

$0.6$

$228.354$

$1.25$

0.419 300

$3.003$

dro8-64

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

w99-139

$0.99$

$2.5$

$0.8$

$211.605$

$1.389$

0.389 231

dro8-64

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

w99-139b

$0.99$

$\mathrm{2.222\hspace{0.17em}222}$

$0.6$

$133.939$

$1.389$

0.386 760

$3.926$

dro8-64

$+$

$+$

$+$

$+$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\sim$

$\sim$

w99-139d

$0.99$

$\mathrm{1.944\hspace{0.17em}444}$

$0.4$

$126.071$

$1.389$

0.384 220

$1.984$

dro6-48

$\times$

$♢$

$\times$

$\sim$

$\times$

$\sim$

$\times$

$\times$

$\times$

w99-139c

$0.99$

$\mathrm{1.805\hspace{0.17em}556}$

$0.3$

$135.092$

$1.389$

0.382 960

$1.743$

dro6-48

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

w99-167m

$0.99$

$3.25$

$0.95$

$1329.680$

$1.667$

0.336 304

$47.740$

dro8-64

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

w99-167

$0.99$

$3.0$

$0.8$

$230.442$

$1.667$

0.333 682

$11.810$

dro10-80

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

w99-167k

$0.99$

$\mathrm{2.333\hspace{0.17em}333}$

$0.4$

$100.014$

$1.667$

0.326 143

$2.332$

dro6-48

$♢$

$♢$

$\times$

$\sim$

$\times$

$\sim$

$\times$

$♢$

$\times$

$\times$

$\times$

w99-167d

$0.99$

$\mathrm{2.166\hspace{0.17em}667}$

$0.3$

$97.092$

$1.667$

0.324 141

$2.332$

dro6-48

$\times$

$♢$

$\times$

$♢$

$\times$

$♢$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

w99-167j

$0.99$

$2$

$0.2$

$98.818$

$1.667$

0.322 115

$1.984$

dro6-48

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

w99-200d

$0.99$

$3.6$

$0.8$

$273.551$

$2.0$

0.281 145

$14.784$

dro8-64

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

w99-200

$0.99$

$3.0$

$0.5$

$111.575$

$2.000$

0.274 441

dro6-48

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\times$

$\times$

w99-200b

$0.99$

$2.8$

$0.4$

$99.000$

$2.000$

0.272 062

$3.003$

dro6-48

$\times$

$\sim$

$\times$

$\sim$

$\times$

$\sim$

$\times$

$♢$

$\times$

$\times$

$\times$

$\times$

w99-200c

$0.99$

$2.6$

$0.3$

$91.817$

$2.000$

0.269 607

$2.636$

dro6-48

$\times$

$♢$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

w99-222b

$0.99$

$\mathrm{3.333\hspace{0.17em}333}$

$0.5$

$117.908$

$2.222$

0.245 832

$4.193$

dro6-48

$\sim$

$\sim$

$\sim$

$\sim$

$♢$

$\sim$

$\times$

$\times$

$\times$

w99-222

$0.99$

$\mathrm{3.111\hspace{0.17em}111}$

$0.4$

$102.677$

$2.222$

0.243 319

$3.926$

dro6-48

$\sim$

$\sim$

$\times$

$\sim$

$\times$

$\sim$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

w99-222c

$0.99$

$\mathrm{2.888\hspace{0.17em}889}$

$0.3$

$93.395$

$2.222$

0.240 719

$3.003$

dro6-48

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

e95

$0.99$

$5.0$

$0.95$

$2436.050$

$2.564$

0.220 775

dro8-64

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

w99-278

$0.99$

$5.0$

$0.8$

$401.302$

$2.778$

0.199 076

$20.543$

dro8-64

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\sim$

$\times$

$\sim$

w99-278b

$0.99$

$\mathrm{4.166\hspace{0.17em}667}$

$0.5$

$139.595$

$2.778$

0.191 854

$5.712$

dro8-64

$\oplus$

$\oplus$

$\oplus$

$♢$

$\oplus$

$\oplus$

$\times$

$\oplus$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

w99-278c

$0.99$

$\mathrm{3.888\hspace{0.17em}889}$

$0.4$

$117.842$

$2.778$

0.189 278

$4.193$

dro6-48

$♢$

$\sim$

$\times$

$\sim$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

w99-278d

$0.99$

$\mathrm{3.611\hspace{0.17em}111}$

$0.3$

$103.835$

$2.778$

0.186 607

$3.672$

dro6-48

$\times$

$♢$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

s99

$0.99$

$\mathrm{5.850\hspace{0.17em}762}$

$\mathrm{0.861\hspace{0.17em}866}$

$771.968$

$3.142$

0.174 409

dro6-48

$\sim$

$\sim$

$\sim$

$\sim$

w99-357

$0.99$

$5.0$

$0.4$

$146.751$

$3.571$

0.139 748

dro6-48

$♢$

$\times$

$\times$

$\times$

$\times$

$\times$

w99-360c

$0.99$

$6.48$

$0.8$

$560.918$

$3.6$

0.147 411

dro6-48

$\sim$

$\sim$

$\sim$

$\times$

w99-360b

$0.99$

$6.12$

$0.7$

$330.647$

$3.6$

0.145 260

dro6-48

$\sim$

$\sim$

$\sim$

$\times$

$\times$

$\times$

w99-360a

$0.99$

$5.76$

$0.6$

$232.413$

$3.6$

0.143 040

$11.810$

dro8-64

$\oplus$

$\oplus$

$\oplus$

$\times$

$\oplus$

$\oplus$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

w99-360j

$0.99$

$5.4$

$0.5$

$179.842$

$3.6$

0.140 745

$7.835$

dro6-48

$♢$

$\sim$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

e9

$0.99$

$7.0$

$0.9$

$1513.855$

$3.684$

0.145 429

$25.105$

dro8-64

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\times$

$\oplus$

$\oplus$

$\times$

$\oplus$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

w99-444

$0.99$

$8.0$

$0.8$

$745.170$

$4.444$

0.113 763

dro6-48

$\sim$

$\times$

$\times$

$\times$

$\times$

$\times$

w95-368

$0.95$

$7.0$

$0.9$

$1516.962$

$3.684$

0.145 349

dro6-48

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\times$

$\times$

n96

$0.9$

$5.5$

$0.6$

$227.038$

$\mathrm{3.437\hspace{0.17em}5}$

0.151 199

dro8-64

$\sim$

$\sim$

$\times$

$\times$

$\times$

$\times$

w9-368

$0.9$

$7.0$

$0.9$

$1521.097$

$3.684$

0.145 202

dro6-48

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\oplus$

$\times$

$\times$

$\times$

$\times$

n95

$0.9$

$10.0$

$0.5$

$378.408$

$6.667$

0.062 994

dro8-64

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

w8-368

$0.8$

$7.0$

$0.9$

$1530.314$

$3.684$

0.144 751

$18.300$

dro6-48

$\times$

$\oplus$

$+$

$\times$

$\oplus$

$\oplus$

$\times$

$\oplus$

$\oplus$

$\times$

$+$

$\oplus$

$\times$

$\times$

$\times$

$\times$

e8b

$0.8$

$8.0$

$0.8$

$756.641$

$4.444$

0.113 578

dro6-48

$\sim$

$\times$

$\times$

$\times$

$\times$

$\times$

e8

$0.6$

$8.0$

$0.8$

$771.968$

$4.444$

0.113 000

dro8-64

$\sim$

$\times$

$\times$

$\times$

$\times$

$\times$

w4-368

$0.4$

$7.0$

$0.9$

$1588.133$

$3.684$

0.140 816

dro6-48

$+$

$+$

$\oplus$

$\oplus$

$\oplus$

$+$

$\oplus$

$\times$

$\sim$

$\times$

$\times$

$\times$

w2-368

$0.2$

$7.0$

$0.9$

$1712.163$

$3.684$

0.137 622

dro6-48

$+$

$\times$

$+$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

ze4

$0.2$

$6.15$

$0.4$

$354.628$

$4.393$

0.106 691

dro6-48

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

zze9

$0.0$

$\mathrm{7.800\hspace{0.17em}001}$

$0.9$

$2224.815$

$4.105$

0.120 223

dro8-64

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

ze9

$0.0$

$\mathrm{7.800\hspace{0.17em}1}$

$0.9$

$2112.079$

$4.105$

0.120 222

dro6-48

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

ns5

$0.0$

$7.2$

$0.5$

$405.662$

$4.8$

0.095 855

dro8-64

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

s0

$0.0$

$\mathrm{10.695\hspace{0.17em}207}$

$\mathrm{0.708\hspace{0.17em}941}$

$771.968$

$6.258$

0.070 830

dro4-32

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

n-55

$-0.5$

$10.0$

$0.5$

$505.428$

$6.667$

0.062 012

dro6-48

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

s-6

$-0.6$

$\mathrm{13.083\hspace{0.17em}066}$

$\mathrm{0.609\hspace{0.17em}412}$

$771.968$

$8.129$

0.048 498

dro4-32

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

wm8-631

$-0.8$

$10.1$

$0.6$

$747.545$

$6.313$

0.066 625

dro4-32

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

wm99-605

$-0.99$

$11.5$

$0.9$

$3401.251$

$6.053$

0.072 312

dro4-32

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

s-99

$-0.99$

$\mathrm{14.542\hspace{0.17em}929}$

$\mathrm{0.534\hspace{0.17em}714}$

$771.968$

$9.476$

0.038 531

dro4-32

$\times$

$\times$

$\times$

$\times$

$\times$

$\times$

Table 3. Table 3: This table explains the meanings of the “wiggle symbols” used in table  2 and figure  1 .

Symbol	Meaning
$\oplus$	we observed oscillations in the diagnostics, successfully fit the appropriate wiggle model described in section IV.1 over a $t$ or $u$ range of $\ge 4$ wiggle periods, and performed the Monte-Carlo error analysis described in section IV.1.5.
$+$	we observed oscillations in the diagnostics and successfully fit the appropriate wiggle model described in section IV.1, but the model was fitted over too short a $t$ or $u$ range ( wiggle periods) for the Monte-Carlo error analysis described in section IV.1.5 to be used
$\sim$	we observed oscillations in the diagnostics which visually appeared to be wiggles, with physically reasonable period and decay rate, but we did not attempt to quantitatively fit a wiggle model
$♢$	we observed oscillations in the diagnostics which might have been wiggles, but these oscillations were not clearly separated from the background variation in the diagnostics
$\times$	we did not observe wiggle-like oscillations in the diagnostics; this could mean either that no wiggles are present, or alternatively that wiggles were present but they were at too low an amplitude and/or insufficiently separated from the background variation to be observed
(blank)	diagnostics were not computed, not computed sufficiently accurately for studying wiggles, or were computed but not assessed

Table 4. Table 4: This table gives the resolution in  r ∗ subscript 𝑟 r_{*} and  θ 𝜃 \theta of the coarsest and finest grids for each of our standard grid structures. The finest grid covers a neighborhood of the worldtube in each slice, while successively coarser grids cover correspondingly larger regions using standard Berger-Oliger 2 :  1 : 2 1 2\,{:}\,1  mesh refinement.

	coarsest grid		finest grid
	$R_{\ast }$	$\theta$	$R_{\ast }$	$\theta$
	($M$)	(radians)	($M$)	(radians)
dro12-96	$1/12$	$\pi /216$	$1/96$	$\pi /1728$
dro10-80	$1/10$	$\pi /180$	$1/80$	$\pi /1440$
dro8-64	$1/8$	$\pi /144$	$1/64$	$\pi /1152$
dro6-48	$1/6$	$\pi /108$	$1/48$	$\pi /864$
dro4-32	$1/4$	$\pi /72$	$1/32$	$\pi /576$

Table 5. Table 5: Numerical fits of the wiggles in the gravitational (spheroidal) ( 𝔩 ​ m ) 𝔩 𝑚 (\mathfrak{l}m) -mode field diagnostics for the ( a ~ , p , e ) = ( 0.99 , 3 , 0.8 ) ~ 𝑎 𝑝 𝑒 0.99 3 0.8 (\tilde{a},p,e)=(0.99,3,0.8) configuration. In each row the top values are the fitted frequency/decay rate; the bottom values are the known values for (least damped) gravitational QNMs. Fitting errors are indicated in parentheses, e.g. , 0.870891 ​ ( 7 ) 0.870891 7 0.870891(7) means 0.870891 ± 7 × 10 − 6 plus-or-minus 0.870891 7 superscript 10 6 0.870891\pm 7\times 10^{-6} .

$𝔩$	$m$	${\omega }_{𝔩mk}$	${\alpha }_{𝔩mk}$
2	0	0.421(3) 0.423	0.069(4) 0.072
2	1	0.57278(6) 0.57274	0.04624(5) 0.04620
3	1	0.7891(1) 0.7891	0.0618(3) 0.0618
2	2	0.870891(7) 0.870892	0.029386(4) 0.029390
3	2	0.999142(9) 0.999140	0.037128(6) 0.037129
4	2	1.1771(1) 1.1770	0.0510(4) 0.0511
3	3	1.32308(3) 1.32308	0.02940(1) 0.02940
4	3	1.43371(3) 1.43372	0.03323(3) 0.03323
5	3	1.58353(3) 1.58351	0.04306(3) 0.04308
4	4	1.7731(2) 1.7731	0.02946(6) 0.02947
5	4	1.87254(4) 1.87254	0.03160(5) 0.03160
5	5	2.2216(2) 2.2217	0.0295(1) 0.0295

Equations

$$d{s}^2=$$

$$+\frac{\Sigma }{\Delta }d{r}^2+\Sigma d{\theta }^2$$

$$+\left(r^2+M^2{\tilde{a}}^2+2M^3{\tilde{a}}^2\frac{r{\sin }^2\theta }{\Sigma }\right){\sin }^2\theta d{\varphi }^2\text{ ,}$$

$$\frac{d{r}_{\ast }}{dr}=\frac{r^2+M^2{\tilde{a}}^2}{\Delta }.$$

$$r_{\ast }=$$

$$-2M\frac{r_{-}}{r_{+}-r_{-}}\ln \left(\frac{r-r_{-}}{2M}\right)\text{ .}$$

$$r_{\text{particle}}(t)=\frac{pM}{1+e\cos {\chi }_r(t)}\text{ ,}$$

$$\Box\Phi_{\text{physical}}=-4\pi q\int\frac{\delta^{4}\bigl{(}x^{a}-x^{a}_{\text{particle}}(t)\bigr{)}}{\sqrt{-g}}\,d\tau\text{~{},}$$

$$F_{a}=\bigl{.}(\nabla_{a}\Phi_{\text{regularized}})\bigr{|}_{x^{i}=x^{i}_{\text{particle}}(t)}\text{~{}.}$$

$$\Phi (t,r,\theta ,\varphi )=\sum _{m=-\infty }^{\infty }\frac{1}{r}{\phi }_m(t,r,\theta )e^{im\tilde{\varphi }}\text{ ,}$$

$$f(r)=\frac{\tilde{a}}{2\sqrt{1-{\tilde{a}}^2}}\ln \left|\frac{r-r_{+}}{r-r_{-}}\right|\text{ .}$$

$$\bigl{.}(\varphi_{m})\bigr{.}_{\text{numerical}}=\begin{cases}\bigl{.}(\varphi_{m})\bigr{.}_{\text{regularized}}&\text{inside the worldtube}\\ \bigl{.}(\varphi_{m})\bigr{.}_{\text{physical}}&\text{outside the worldtube}\end{cases}\text{~{}.}$$

$$F_a=q\sum _{m=0}^{\infty }{F}_a^{(m)}\text{ ,}$$

$${\Phi }_{-2}:=(r-ia\cos \theta {)}^4{\Psi }_4,$$

$${\Phi }_{-2}=\sum _{\mathfrak{l}m\omega }{R}_{\mathfrak{l}m\omega }(r)S_{\mathfrak{l}m\omega }(\theta )e^{im\varphi -i\omega t}\text{ ,}$$

$$F_a=\sum _{\begin{array}{c}\mathfrak{l}m\omega \\ nk\pm \end{array}}{\mathcal{C}}_{am\omega nk}^{\pm }{R}_{\mathfrak{l}m\omega }^{(n)\pm }(r_0)S_{\mathfrak{l}m\omega }^{(k)}({\theta }_0)e^{im{\varphi }_0-i\omega t_0}\text{ ,}$$

$$F_a^{(\ell )}=\sum _{\begin{array}{c}\mathfrak{l}m\omega \\ n\pm \end{array}}{\tilde{\mathcal{C}}}_{a\mathfrak{l}m\omega n}^{\ell \pm }{R}_{\mathfrak{l}m\omega }^{(n)\pm }(r_0)Y_{\ell m}({\theta }_0)e^{im{\varphi }_0-i\omega t_0}\text{ .}$$

$$F_a=\sum _{\ell }{F}_a^{(\ell )}-B_a\text{ ,}$$

$${\Phi }_{-2}^{(\mathfrak{l}m)}=\sum _{\omega }{R}_{\mathfrak{l}m\omega }(r)S_{\mathfrak{l}m\omega }(\theta )e^{im\varphi -i\omega t},$$

$$F_a^{(\mathfrak{l}m)}:=\sum _{\begin{array}{c}\omega \\ nk\pm \end{array}}{\mathcal{C}}_{am\omega nk}^{\pm }{R}_{\mathfrak{l}m\omega }^{(n)\pm }(r_0)S_{\mathfrak{l}m\omega }^{(k)}({\theta }_0)e^{im{\varphi }_0-i\omega t_0}.$$

$$\bigl{.}(\varphi_{m})\bigr{.}_{\text{regularized}}=B(t)+\sum_{k}A^{(k)}\exp\left(-\lambda^{(k)}(t-t_{\text{ref}})\right)\,\sin\left(2\pi\frac{t-t_{\text{ref}}}{P^{(k)}}+\eta^{(k)}\right)\,\,\text{,}$$

$$F_{a}(u)=\frac{B_{a}(u)}{r^{3}_{\text{particle}}(u)}+\sum_{k}\frac{A^{(k)}_{a}}{r_{\text{particle}}(u)}\exp\left(-\lambda^{(k)}(u-u_{\text{ref}})\right)\sin\left(2\pi\frac{u-u_{\text{ref}}}{P^{(k)}}-m\bigl{(}\phi_{\text{particle}}(u)-\phi_{\text{ref}}\bigr{)}+\eta^{(k)}_{a}\right)\,\,\text{,}$$

$$\biggl{.}\partial_{tt}\left(\bigl{.}(\varphi_{m})\bigr{.}_{\text{regularized}}\right)\biggr{|}_{\mathcal{J}^{+}}=B(t)+\sum_{k}A^{(k)}\exp\left(-\lambda^{(k)}(t-t_{\text{ref}})\right)\,\sin\left(2\pi\frac{t-t_{\text{ref}}}{P^{(k)}}+\eta^{(k)}\right)\,\,\text{,}$$

$${\Phi }_{-2}^{(\mathfrak{l}m)}=\sum _n{B}^{(n)}{t}^n+\sum _k\left\{A_s^{(k)}\sin \left[{\omega }_k(t-t_{\text{ref}})\right]+A_c^{(k)}\cos \left[{\omega }_k(t-t_{\text{ref}})\right]\right\}{e}^{-{\alpha }_k(t-t_{\text{ref}})}.$$

$$F_a^{(\mathfrak{l}m)}=\sum _n{B}_a^{(n)}{u}^n+\sum _k\left\{A_{s,a}^{(k)}\sin \left[{\omega }_k(u-u_{\text{ref}})-m{\varphi }_{\text{particle}}\right]+A_{c,a}^{(k)}\cos \left[{\omega }_k(u-u_{\text{ref}})-m{\varphi }_{\text{particle}}\right]\right\}\frac{e^{-{\alpha }_k(u-u_{\text{ref}})}}{r_{\text{particle}}}.$$

$$\underset{r\to \infty }{\lim }r{\Psi }_4^{(\mathfrak{l}m)}=\sum _n{B}^{(n)}{u}^n+\sum _k\left\{A_s^{(k)}\sin \left[{\omega }_k(u-u_{\text{ref}})\right]+A_c^{(k)}\cos \left[{\omega }_k(u-u_{\text{ref}})\right]\right\}{e}^{-{\alpha }_k(u-u_{\text{ref}})}\text{.}$$

$$|A^{(k)}|:=\sqrt{{\left(A_c^{(k)}\right)}^2+{\left(A_s^{(k)}\right)}^2}.$$

$${\delta }^{(k)}:=\underset{n\in \text{mathbbm}Z}{\min }\frac{|{\omega }^{(k)}-m{\Omega }_{\varphi }-n{\Omega }_r|}{{\Omega }_r},$$

Full text

I Introduction

Consider a small (compact) body of mass $\mu M$ (with $0<\mu\ll 1$) moving freely near a Schwarzschild or Kerr black hole of mass $M$. This system emits gravitational radiation, and there is a corresponding radiation-reaction influence on the small body’s motion. Calculating the resulting perturbed spacetime (including the small body’s motion and the emitted gravitational radiation) is a long-standing research question in general relativity.

There is also an astrophysical motivation for this calculation: If a neutron star or stellar-mass black hole of mass ${\sim}\,1$–$100M_{\odot}$ orbits a massive black hole of mass ${\sim}\,10^{5}$–$10^{7}M_{\odot}$,111$M_{\odot}$ denotes the solar mass. the resulting “extreme–mass-ratio inspiral” (EMRI) system is expected to be a strong astrophysical gravitational-wave (GW) source detectable by the planned Laser Interferometer Space Antenna (LISA) space-based gravitational-wave detector. LISA is expected to observe many such systems, some of them at quite high signal-to-noise ratios (Gair et al. (2004); Barack and Cutler (2004); Amaro-Seoane et al. (2007); Gair (2009)). The data analysis for, and indeed the detection of, such systems will generally require matched-filtering the detector data stream against appropriate precomputed GW templates. The problem of computing such templates provides an astrophysical motivation for EMRI modelling.

In the test-particle limit it has long been known that an unbound (scattering) flyby can excite quasinormal modes (QNMs) of the background black hole. Kojima and Nakamura Kojima and Nakamura (1984) studied this process, finding that “A scattered particle excites the quasi-normal mode under the condition that twice the angular velocity at the periapsis is greater than the real part of the frequency of the quasi-normal mode”. Their figure 3(b) shows an example of the QNM oscillations in the radiated gravitational waves at $\mathcal{J}^{+}$.

Burko and Khanna Burko and Khanna (2007) found small oscillations in the total radiated energy flux from a test particle making a parabolic (unbound) flyby of a Kerr black hole. They attributed these oscillations to the particle encountering scattered gravitational waves emitted during the particle’s inbound motion.

O’Sullivan and Hughes O’Sullivan and Hughes (2016) observed “small-amplitude high-frequency oscillations” in their calculations of the horizon shear of a Kerr black hole orbited by a test particle. Because they did not find corresponding oscillations in the horizon’s tidal distortion field, and their measured oscillation frequencies did not match known Kerr QNM frequencies, they concluded that the horizon-shear oscillations they observed “cannot be related to the [Kerr black] hole’s quasi-normal modes”.

Thornburg and Wardell Thornburg and Wardell (2017) (hereinafter TW) calculated the scalar-field self-force for eccentric equatorial particle orbits in Kerr spacetime. For some systems where the Kerr black hole was highly spinning and the particle orbit was prograde and highly eccentric, TW found that the self-force exhibits large oscillations (“wiggles”) on the outgoing leg of the orbit shortly after periapsis passage. TW suggested that wiggles “are in some way caused by the particle’s close passage by the large black hole”. Thornburg Thornburg (2016, 2017) presented fits of damped-sinusoid models to these wiggles for a range of Kerr spins and particle orbits, found close agreement of the model complex-frequencies with those of known Kerr QNMs, and argued that this agreement shows that wiggles are, in fact, caused by Kerr QNMs excited by the particle’s close periapsis passage.

Nasipak, Osburn, and Evans Nasipak et al. (2019) calculated the scalar-field self-force and the radiated field at $\mathcal{J}^{+}$ for eccentric inclined particle orbits in Kerr spacetime. For one particular (equatorial) orbit configuration they confirmed TW’s finding of wiggles in the self-force and also found wiggles in the radiated scalar field at $\mathcal{J}^{+}$, fitting these to a superposition of $\ell\,{=}\,m\,{=}\,1$, $\ell\,{=}\,m\,{=}\,2$, $\ell\,{=}\,m\,{=}\,3$, and $\ell\,{=}\,m\,{=}\,4$ Kerr quasinormal modes (QNMs). They concluded that wiggles are caused by Kerr QNMs excited by the particle’s close periapsis passage.

Rifat, Khanna, and Burko Rifat et al. (2019) recently studied wiggles for EMRIs where the central Kerr BH is nearly extremal (dimensionless spin up to $0.999\,99$), finding that in such systems many Kerr QNMs can be simultaneously excited.

Here we extend Refs. Thornburg (2016, 2017) and survey wiggles’ phenomenology over a wide range in parameter space for eccentric equatorial orbits in Kerr spacetime, for both the scalar-field model and the full gravitational field. We focus on leading-order radiation and radiation-reaction effects computed via 1st-order perturbations of Kerr spacetime, i.e., (for the gravitational case) $\mathcal{O}(\mu)$ field perturbations near the Kerr black hole, $\mathcal{O}(\mu)$ radiation at $\mathcal{J}^{+}$, and $\mathcal{O}(\mu^{2})$ radiation-reaction “self-forces” acting on the small body. We fit models of Kerr QNMs to all these diagnostics.

Our focus is the case where the perturbing body’s orbit is highly relativistic, so post-Newtonian methods (see, for example, (Damour, 1987, Section 6.10); Poisson and Will (2014); Will (2018); Blanchet (2014); Futamase and Itoh (2007); Blanchet (2011); Schäfer (2011) and references therein) are not reliably accurate. Since the timescale for radiation reaction to shrink the orbit is very long ($\sim\mu^{-1}M$) while the required resolution near the small body is very high ($\sim\mu M$), a direct “numerical relativity” integration of the Einstein equations (see, for example, Pretorius (2007); Hannam et al. (2009); Hannam (2009); Hannam and Hawke (2011); Campanelli et al. (2010) and references therein) would be prohibitively expensive (and probably insufficiently accurate) for this problem.222A number of researchers have attempted direct numerical-relativity binary black hole simulations for systems with “intermediate” mass ratios up to $100\,{:}\,1$ ($\mu=0.01$), (see, for example, Bishop et al. (2003, 2005); Sopuerta et al. (2006); Sopuerta and Laguna (2006); Lousto et al. (2010); Lousto and Zlochower (2011); Husa et al. (2016)). However, it has not (yet) been possible to extend these results to the extreme-mass-ratio case nor to accurately evolve any systems with mass ratios more extreme than $18\,{:}\,1$ for a radiation-reaction time scale Husa et al. (2016) (Hinder Hinder (2019) reports ongoing efforts to extend this to $32\,{:}\,1$). Instead, we use black hole perturbation theory, treating the small body as an $\mathcal{O}(\mu)$ perturbation on the background spacetime.

We present results obtained from two separate numerical codes which were programmed independently, using completely different theoretical formalisms and numerical methods:

• •

Our scalar-field results were obtained using TW’s code Wardell et al. (2012); Thornburg and Wardell (2017) extended to compute additional field diagnostics. This code uses an effective-source regularization followed by an $e^{im\phi}$ Fourier-mode decomposition and a separate $2{+}1$-dimensional time-domain numerical evolution of each Fourier mode. The main outputs of this code are the regularized scalar field at selected (fixed) spatial positions, the regularized scalar field and self-force at the particle, and the physical radiated scalar field at $\mathcal{J}^{+}$.

• •

Our gravitational-field results were obtained using van de Meent’s code van de Meent and Shah (2015); van de Meent (2016, 2017, 2018). This code obtains the local metric perturbation in the frequency domain by reconstructing the metric perturbation from the Weyl scalar $\Psi_{4}$, followed by $\ell$-mode regularization to obtain the regular part. The main outputs of this code are the regularized outgoing–radiation-gauge metric perturbation and self-force at the particle, and the physical radiated $\Psi_{4}$ at $\mathcal{J}^{+}$ and selected fixed positions in the spacetime.

For both codes we take the particle orbit to be a bound equatorial geodesic; we do not consider changes in the orbit induced by the self-force.

To briefly summarize our main results, we observe wiggles across a wide range of Kerr spins and particle orbits. Wiggles are present in all of our field diagnostics in the strong-field region and at $\mathcal{J}^{+}$. Except for a few anomalous results for near-extremal Kerr spacetimes (dimensionless Kerr spins $\gtrsim 0.9999$), our results are all consistent with the interpretation of wiggles as Kerr QNMs. Wiggles are stronger and more readily observable for high Kerr spins, highly eccentric prograde particle orbits, and particle periapsis radii close to the light ring.

The remainder of this paper is organized as follows:

Section I.1 summarizes our notation and our parameterization of bound geodesic orbits in Kerr spacetime. Section II briefly summarizes our calculations of scalar-field (section II.1) and gravitational (section II.2) perturbations of Kerr spacetime. Section III describes our local- and radiated-field diagnostics. Section IV describes our QNM models and how we fit these to time series of our field diagnostics. Section V presents our data for scalar-field and gravitational perturbations of Kerr spacetime, and QNM-model fits to this data. Finally, section VI discusses our results and presents our conclusions.

I.1 Notation and parameterization of Kerr geodesics

We generally follow the sign and notation conventions of Wald Wald (1984), with $G=c=1$ units and a $(-,+,+,+)$ metric signature. We use the Penrose abstract-index notation, with indices $abcd$ running over spacetime coordinates, and $ijk$ running over the spatial coordinates. $\nabla_{a}$ is the (spacetime) covariant derivative operator and $g$ is the determinant of the spacetime metric. $X:=Y$ means that $X$ is defined to be $Y$. $\Box:=\nabla_{a}\nabla^{a}$ is the 4-dimensional (scalar) wave operator Brill et al. (1972); Teukolsky (1973).

We use Boyer-Lindquist coordinates $(t,r,\theta,\phi)$ on Kerr spacetime, defined by the line element

[TABLE]

where $M$ is the black hole’s mass, $\tilde{a}=J/M^{2}$ is the black hole’s dimensionless spin (limited to $|\tilde{a}|<1$), $\Sigma:=r^{2}+M^{2}\tilde{a}^{2}\cos^{2}\theta$, and $\Delta:=r^{2}-2Mr+M^{2}\tilde{a}^{2}$. We also define $a:=M\tilde{a}$ (this is unrelated to the use of $a$ as an abstract tensor index). In these coordinates the event horizon is the coordinate sphere $r=r_{+}=M\left(1+\sqrt{1-\tilde{a}^{2}}\right)$ and the inner horizon is the coordinate sphere $r=r_{-}=M\left(1-\sqrt{1-\tilde{a}^{2}}\right)$. (See footnote 3 for a discussion of TW’s coordinate compactification near the event horizon and $\mathcal{J}^{+}$.)

Following Ref. Sundararajan et al. (2007), we define the tortoise coordinate $r_{*}$ by

[TABLE]

and fix the additive constant by choosing

[TABLE]

$u:=t-r_{*}$ is thus an outgoing null coordinate.

The particle’s worldline is $x^{i}=x^{i}_{\text{particle}}(t)$; we consider this to be known in advance, i.e., we do not consider changes to the particle’s worldline induced by the self-force. For present purposes we consider only particle orbits in the Kerr spacetime’s equatorial plane; this restriction is for computational convenience and is not fundamental. We take the particle to orbit in the $d\phi/dt>0$ direction, with $\tilde{a}>0$ for prograde orbits and $\tilde{a}<0$ for retrograde orbits.

We define $r_{\min}$ and $r_{\max}$ to be the particle’s periapsis and apoapsis $r$ coordinates, respectively. We parameterize bound geodesic equatorial particle orbits by the usual (dimensionless) semi-latus rectum $p$ and eccentricity $e$ (defined by r_{\min}=pM\big{/}(1+e) and r_{\max}=pM\big{/}(1-e)), so that the particle orbit is given by

[TABLE]

for a suitable orbital-phase function $\chi_{r}$. We refer to the combination of a spacetime and a (bound geodesic) particle orbit as a “configuration”, and parameterize it with the triplet $(\tilde{a},p,e)$. We define $T_{r}$ to be the coordinate-time period of the particle’s radial motion; we usually refer to $T_{r}$ as the particle’s “orbital period”.

II Calculations of scalar-field and gravitational perturbations

of Kerr spacetime

II.1 Scalar-field perturbations of Kerr spacetime

In this section we summarize TW’s scalar-field calculations Wardell et al. (2012); Thornburg and Wardell (2017). These authors consider a real scalar field $\Phi_{\text{physical}}$ satisfying the wave equation in Kerr spacetime,

[TABLE]

sourced by a point “particle” of scalar charge $q$ which is taken to move on a (pre-specified) equatorial geodesic orbit. $\Phi_{\text{physical}}$ satisfies outgoing-radiation (retarded) boundary conditions at $\mathcal{J}^{+}$. This system provides a toy model of the full gravitational perturbation problem without the complexity of gauge choice.

Because $\Phi_{\text{physical}}$ diverges on the particle worldline, (5) must be regularized. TW use an effective-source regularization of the type introduced by Barack and Golbourn Barack and Golbourn (2007) and Vega and Detweiler Vega and Detweiler (2008) (see Vega et al. (2011) for a review), using the puncture function and effective source described by Wardell et al. Wardell et al. (2012). In a neighborhood of the particle worldline, TW define a (real) regularized scalar field $\Phi_{\text{regularized}}=\Phi_{\text{physical}}-\Phi_{\text{puncture}}$, where $\Phi_{\text{puncture}}$ is a suitably-chosen approximation to the Detweiler-Whiting singular field Detweiler and Whiting (2003). The (4-vector) self-force acting on the particle is then given by

[TABLE]

TW make an azimuthal Fourier decompositions of all the spacetime scalar fields into complex $e^{im\tilde{\phi}}$ modes,

[TABLE]

where the extra factor of $1/r$ is introduced for computational convenience and where $\tilde{\phi}:=\phi+f(r)$ is an “untwisted” azimuthal coordinate, with

[TABLE]

For each $m$-mode, TW introduce a finite worldtube surrounding the particle worldline in $(t,r,\theta)$ space. For particle orbits with significant eccentricity ($e\gtrsim 0.2$) the worldtube (now viewed as a region in $(r,\theta)$ in each $t=\text{constant}$ slice) moves inward and outward in $r$ during each orbit so as to always contain the particle. All the results reported here were obtained using a worldtube which is rectangular in $(r,\theta)$, with size $10M$ in $r_{*}$ by approximately $0.25$ radians in $\theta$, symmetric about the equatorial plane, and kept centered on the particle to within $0.25M$ in $r_{*}$ as the particle moves.

TW numerically solve for the piecewise function

[TABLE]

using a time-domain $2{+}1$-dimensional finite-difference numerical evolution with mesh refinement. Because the (Kerr) background spacetime is axisymmetric, the Fourier modes $\varphi_{m}$ evolve independently – there is no mixing of the modes. Because the physical scalar field $\Phi$ is real, only the $m\geq 0$ modes need to be explicitly computed; the $m<0$ modes may be obtained by symmetry.

Corresponding to the Fourier decomposition (7), the self-force (6) can be written as a similar sum of $e^{im\tilde{\phi}}$ modes,

[TABLE]

where each $F_{a}^{(m)}$ may be computed from the corresponding \bigl{.}(\varphi_{m})\bigr{.}_{\text{regularized}} field near the particle. TW compute a finite set of modes (typically $0\leq m\leq 20$) and estimate the $m>20$ contributions to the sum (10) via a large-$m$ asymptotic series.

TW use a Zenginoğlu-type hyperboloidal compactification Zenginoğlu (2008a, b, 2011); Zenginoğlu and Khanna (2011); Zenginoğlu and Kidder (2010); Zenginoğlu and Tiglio (2009); Bernuzzi et al. (2011, 2012) so they also have direct access to far-field quantities at $\mathcal{J}^{+}$ (where the coordinate $t$ becomes a Bondi-type retarded time coordinate).333More precisely, TW define compactified coordinates $(T,R_{*})$ which are identical to (respectively) the Boyer-Lindquist $t$ and the tortoise coordinate $r_{*}$ throughout a neighborhood of the region $r_{\min}\leq r\leq r_{\max}$ containing the particle orbit, but which are compactified near the event horizon and $\mathcal{J}^{+}$. $T$ is a Bondi-type retarded time coordinate at $\mathcal{J}^{+}$. For present purposes the details of the compactification are not important, so for convenience of exposition we refer to $T$ as $t$ hereinafter when describing our diagnostics at $\mathcal{J}^{+}$.

II.2 Gravitational perturbations of Kerr spacetime

In this section we summarize the metric reconstruction approach used by van de Meent van de Meent and Shah (2015); van de Meent (2016, 2017, 2018) to calculate gravitational perturbation of Kerr spacetime generated by particles on bound geodesic orbits. This approach starts from the the spin-(-2) Teukolsky variable,

[TABLE]

where $\Psi_{4}$ is one of the Weyl scalars. As shown by Teukolsky Teukolsky (1973), linear perturbations to this variable satisfy an equation of motion that decouples from all other degrees of freedom. Moreover, unlike the linearized Einstein equation, solutions to the Teukolsky equation can be decomposed into Fourier-harmonic modes,

[TABLE]

where the $R_{\mathfrak{l}m\omega}$ are solutions of the radial Teukolsky equation, the $S_{\mathfrak{l}m\omega}$ are spin-weight spheroidal harmonics, and $\mathfrak{l}$ is the spheroidal mode number. In van de Meent’s code the radial Teukolsky equation can be solved to arbitrarily high precision using a numerical implementation of the semi-analytical methods of Mano, Suzuki, and Takasugi Mano and Takasugi (1997); Mano et al. (1996).

Remarkably, $\Phi_{-2}$ contains almost all information about the corresponding perturbation of the metric Wald (1973), and in vacuum it is possible to reconstruct the metric perturbation in a radiation gauge Cohen and Kegeles (1974); Kegeles and Cohen (1979); Chrzanowski (1975); Wald (1978). As detailed in Refs. van de Meent and Shah (2015); van de Meent (2016, 2017, 2018), this procedure can be used to calculate the backreaction of the metric perturbation on the particle, the gravitational self-force, which then takes the form

[TABLE]

where the $R_{\mathfrak{l}m\omega}^{\pm}(r_{0})$ are vacuum solutions of the radial Teukolsky equation analytically extended to the particle position $r_{0}$ from either infinity ($+$) or the black hole horizon ($-$) (method of extended homogeneous solutions Barack et al. (2008)), and the $(n)$ and $(k)$ superscripts on a function denote derivatives with respect to the argument. The indices $n$ and $k$ run from 0 to 3.

Although each individual term in the sum above is finite, the sum as a whole does not converge. This is a simple consequence of the fact that it was built from the retarded field perturbation rather than the Detweiler-Whiting regular field. To obtain the regular field we still need to subtract the Detweiler-Whiting singular field. In principle this can be done mode-by-mode. To match previous analytical calculations of the large $\ell$-behavior of the singular field Barack and Ori (2003); Barack (2009), we need to re-expand the spheroidal $\mathfrak{l}$-modes to spherical $\ell$-modes,

[TABLE]

With this re-expansion, the local gravitational self force is given by

[TABLE]

where, as shown in Pound et al. (2014), we can use the Lorenz-gauge $B_{a}$ parameter given in Barack and Ori (2003); Barack (2009).

The metric reconstruction formalism can only recover parts of the metric that contribute to $\Psi_{4}$. This means that the reconstructed metric carries an ambiguity, which can be shown Wald (1973) to consist of perturbations of the background within the class of Kerr metrics and pure gauge contributions. These ambiguities can be uniquely fixed based on physical considerations Merlin et al. (2016); van De Meent (2017). The corrections needed to fix these ambiguities are known and straightforward to add. They contribute only to the low frequency envelope of the self-force. Hence, to facilitate identification and extraction of the wiggles in the gravitational self-force, we omit them in this work.

Frequency domain calculations of the type used here are ideally suited for calculating metric perturbations with a sparse discrete frequency spectrum, such as those sourced by a particle on a low eccentricity geodesic. That spectrum becomes denser at higher eccentricities, necessitating the calculation of more frequency modes and making the calculation more time-consuming. Moreover, as discussed in detail in Ref. van de Meent (2016), the method of extended homogeneous solutions leads to large cancellations between different (low) frequency modes for high-$\ell$ modes, causing a large loss of precision. In this work we tackle this problem by harnessing the full power of the arbitrary precision implementation of our code and simply throw more precision at the computation than we lose in the cancellation.

For this work we calculated the full gravitational self-force for orbits with eccentricities up to $e=0.8$, which involves dealing with cancellations of around 30 orders of magnitude. These calculations are fairly resource intensive, taking $\mathcal{O}(10^{4})$ CPU hours (or a few days on 400 CPUs) to compute.

However, for many aspects of our investigation here, we do not need the full local regular metric perturbation. If we want to look for the dominant low-$\mathfrak{l}$ QNMs, then these will contribute (mostly) to the low-$\mathfrak{l}$ modes of the gravitational metric perturbation. For this purpose, we define the individual $\mathfrak{l}m$ modes of the Teukolsky variable

[TABLE]

and the gravitational self-force

[TABLE]

These are much easier to compute, and for low $\mathfrak{l}$ do not suffer from the large loss of precision due to the method of extended homogeneous solutions, thus allowing for very high accuracy calculations without excessive computational cost.

III Field diagnostics

We consider several different types of local- and radiated-field diagnostics, and attempt to fit the wiggles in these diagnostics to QNM models. Clearly the presence of wiggles in the physical scalar field or metric perturbation implies the presence of wiggles in some or all of the $e^{im\tilde{\phi}}$ scalar-field modes or $(\mathfrak{l}m)$ metric-perturbation modes (respectively), and vice versa.444While it would be theoretically possible for multiple modes to have wiggles of the same frequency whose amplitudes sum to approximately zero (leading to an absence of wiggles in the physical fields), in practice we have never observed this. Because many fewer QNMs are present at significant amplitude (usually only one, or in a few cases two), it is much simpler to analyze wiggles in the individual modes.

Table 1 summarizes our local- and radiated-field diagnostics for studying wiggles. We consider (time series of) diagnostics at three locations in spacetime:

• •

Diagnostics of the local field at selected fixed spatial “watchpoint” coordinate positions $(r,\theta,\phi)=\text{constant}$. These diagnostics directly sample any QNMs that may be present, but the diagnostics may be contaminated by the direct field when the particle is close to the watchpoint position.

For the scalar-field case, we avoid any such possible contamination by considering the regularized field mode \bigl{.}(\varphi_{m})\bigr{.}_{\text{regularized}}. However, this is only defined within the worldtube, so for orbits with significant eccentricity (where the worldtube moves in $(r,\theta)$ during the particle orbit) any given watchpoint may lie outside the worldtube (and thus leave \bigl{.}(\varphi_{m})\bigr{.}_{\text{regularized}} undefined) during some parts of the orbit. To minimize this effect, for many of the analyses reported here we use watchpoint positions which are near the orbit’s apoapsis, where the particle (and hence the worldtube) motion is relatively slow and hence a suitable watchpoint can remain within the worldtube for a relatively long time in each orbit. All our scalar-field watchpoints are in the equatorial plane.

For the gravitational case, the regularized field is not readily available, so instead we have the code output the retarded $\Phi_{-2}^{(\mathfrak{l}m)}$ on the symmetry axis of the background Kerr spacetime ($\theta=0$) and the equatorial plane $(\theta=\pi/2)$ at coordinate radii corresponding to the particle’s periapsis and apoapsis distances.

• •

Diagnostics of the local field at the particle. Here we consider the $e^{im\tilde{\phi}}$ (scalar-field) and $(\mathfrak{l}m)$ (gravitational) modes of the self-force itself. The main complication here is that these diagnostics sample the field perturbation at a time-dependent position (the particle position), so our fitting model for the QNM effects must include corrections for the spatial variation of the QNM eigenfunctions as the particle (sampling point) moves. For the azimuthal ($\phi$) particle motion this is straightforward (described in section IV) but for the radial ($r$) motion we include this correction only approximately.

• •

Diagnostics of the radiated field at $\mathcal{J}^{+}$. These have the advantage of being physically observable and of allowing the $e^{im\tilde{\phi}}$ (scalar-field) and $(\mathfrak{l}m)$ (gravitational) mode decompositions to be defined in a weak-field region (for the gravitational case, this also avoids any gauge dependence).

At $\mathcal{J}^{+}$ we only have the physical (retarded) fields available, so it is more difficult to separate wiggles from the overall radiation pattern. To help in making this separation, we observe that wiggles are of relatively high (temporal) frequency relative to other major features in the radiated fields, so that taking time derivatives of the radiated fields increases the wiggles’ amplitude relative to that of the other features. We have found that good results are obtained by using as diagnostics the second time derivatives \partial_{tt}\left(\bigl{.}(\varphi_{m})\bigr{.}_{\text{physical}}\right) evaluated in the equatorial plane (scalar field)555Recall (cf. footnote 3) that in our scalar-field computations, $t$ is a Bondi-type retarded time coordinate at $\mathcal{J}^{+}$. and $\Psi_{4}$ evaluated either on the z axis or in the equatorial plane (gravitational).

IV Quasinormal-mode models and fits

IV.1 Scalar-field perturbations

IV.1.1 Perturbations at a fixed spatial position

Consider first the case of wiggles in an individual $e^{im\phi}$ Fourier mode of the regularized scalar field, observed at a fixed “watchpoint” spatial position in the strong-field region. We consider the model

[TABLE]

where $B$ is a spline function representing the slowly-varying “background” variation of the scalar field, $k$ indexes the damped-sinusoids included in the model, $A^{(k)}$, $\lambda^{(k)}$, $P^{(k)}$, and $\eta^{(k)}$ are respectively the amplitude, damping rate, period, and phase offset of each damped-sinusoid, and the subscript ${}_{\text{ref}}$ denotes a “reference” time chosen for convenience. To avoid degeneracy between the spline and the damped-sinusoid we require that the spacing in $t$ between adjacent spline control points always be at least $1.5P^{(\max)}$, where $P^{(\max)}:=\displaystyle\max_{k}P^{(k)}$ is the period of the longest-period damped-sinusoid in the model.

IV.1.2 Perturbations at the particle position

Consider next the case of wiggles in the radiation-reaction self-force (which is implicitly defined at the particle position). This introduces two complications: the self-force is a 4-vector (with nontrivial $t$, $r$, and $\phi$ components for our equatorial orbits), and the field perturbation is being sampled at a time-varying position. Analogously to (18), we thus consider the model

[TABLE]

where we now parameterize the particle’s motion using the retarded time coordinate $u$,666Heuristically, the choice of $u$ rather than $t$ as an independent variable in the model is suggested by the QNM signals propagating outward along light cones after being excited near periapsis. However, it is not clear that this is a good approximation to actual QNM dynamics, so we experimented with models using both $t$ and $u$ as independent variables. We found the latter to give better fits to our numerical data. $B_{a}$ is now a 4-vector spline function representing the background variation of the self-force along the particle worldline, $k$ again indexes the damped-sinusoids included in the model, $A_{a}^{(k)}$, $\lambda^{(k)}$, $P^{(k)}$, and $\eta_{a}^{(k)}$ are now respectively the 4-vector amplitude, damping rate, period, and 4-vector phase offset of each damped-sinusoid, and the subscript ${}_{\text{ref}}$ again denotes a “reference” time chosen for convenience. The non-degeneracy condition on the background spline now applies to the spacing in $u$ between adjacent spline control points.

Notice that the damping rate and oscillation period of each damped-sinusoid are common across all tensor components of the self-force. The $-m\phi_{\text{particle}}(u)$ term models the variation in oscillation phase due to particle’s (i.e., the sampling point’s) motion in $\phi$. The $1/r^{3}_{\text{particle}}(u)$ and $1/r_{\text{particle}}(u)$ factors model the expected far-field variation in self-force and in the oscillation eigenfunction amplitude with position. (Actual QNM eigenfunctions have a much more complicated variation of amplitude with spatial position, but for simplicity we omit this from our model.)

IV.1.3 Perturbations at $\mathcal{J}^{+}$

Finally, consider the case of wiggles in the radiated (physical) field at $\mathcal{J}^{+}$. Because of the hyperboloidal time slices, we observe these at finite coordinate time $t$ (the $\mathcal{J}^{+}$ time has an arbitrary offset relative to the strong-field coordinate time). As noted in section III, here it is somewhat difficult to separate wiggles from the overall radiation pattern, so we consider second time derivatives of the physical scalar-field modes. Analogously to (18) and (19), we thus consider the model

[TABLE]

where the meanings of all terms (and the non-degeneracy condition on the background spline) are the same as in (18).

IV.1.4 Fitting the models

For each of the models (18), (19), and (20), we visually inspect plots of our time-series data to identify a suitable range of the independent variable for fitting and to choose initial guesses for the background and wiggle parameters, then use the nonlinear least-squares fitting subroutine LMDIF1 from the MINPACK library Moré et al. (1980) to fit the model to the data. To make the model closer to linear (which improves the convergence of the nonlinear fitting), we fit the wiggle amplitudes and phases as cosine- and sine-component amplitudes (i.e., $A\sin(X+\eta)$ is actually fitted as $A^{(\cos)}\cos(X)+A^{(\sin)}\sin(X)$). In most cases we used uniform weighting for the fits, but in a few cases we used weights proportional to $r^{3}$ so as to improve the fit at late times (close to apoapsis).

IV.1.5 Uncertainties in the fitted wiggle parameters

The residuals from our wiggle fits are not random, but rather are dominated by low-amplitude oscillations of similar frequency to the wiggles themselves (this can be seen in figures 4 and 5). This means that formal uncertainty estimates for the fitted parameters \bigl{(}P^{(k)},\tau^{(k)}\bigr{)} (derived assuming uncorrelated Gaussian residuals) are not realistic. Because of the oscillatory nature of the residuals, the fitted parameters are slightly dependent on the precise choice of fitting interval; this is, in fact, usually the dominant uncertainty in the fitted parameters.

We use a Monte-Carlo procedure to estimate realistic uncertainties in the fitted parameters: Given a fit of one of the above wiggle models to our data in some interval $I$ (in either $t$ or $u$) of length $L_{\text{fit}}\geq 4P^{(\max)}$, we randomly choose $N_{\text{trial}}=300$ subintervals of $I$ (randomly sampling each lower and upper interval endpoint from a uniform distribution) subject to the constraint that each subinterval must have a minimum length of $L_{\min}\,{=}\,3P^{(\max)}$.777The minimum-length requirement for the subintervals ensures that each subinterval is long enough to allow a reasonable estimate of the wiggle decay rate and period (the fitting errors should scale roughly inversely with $L_{\min}$). The minimum-length requirement for the full fitting interval ensures that different subintervals can sample significantly different regions of the data. Then we repeat the wiggle-model fit for each subinterval.888Each subinterval fit uses a subset of the original fit’s background spline control points which just span that subinterval, plus one point outside the subinterval on each of the interval’s left and right endpoints. The ensemble of the $N_{\text{trial}}$ sets of Monte-Carlo–trial fitted parameters then provides an estimate of the uncertainty in the fitted parameters from the full-interval fit.

After allowing initial transients to decay, our numerical calculations extend over a number of particle orbits. Because the particle orbit precesses strongly, each orbit places the particle in a different position with respect to any fixed watchpoint or $\mathcal{J}^{+}$ observer. For each orbit we repeat the entire fitting procedure (including the full set of Monte-Carlo sub-interval trials). Our final estimate for the uncertainty in the fitted parameters is obtained from the union of all the Monte-Carlo trials over several (typically 3) distinct orbits.

This procedure has two main limitations:

• •

The procedure is not applicable to cases where the overall fitting interval is too short (length $L_{\text{fit}}<4P^{(\max)}$). (Footnote 7 outlines the reasons for this.)

• •

If a wiggle is rapidly damped, then the wiggle amplitude becomes very small at the right (large $t$ or $u$) end of a long fitting interval, so a subinterval of near-minimum length ($3P^{(\max)}$) which is close to the right end of the overall fitting interval will have a poorly-constrained fit, yielding a large scatter in the fitted parameters.

These limitations are most severe when the wiggles have low amplitude and are rapidly damped, as is the case for low Kerr spins.

IV.1.6 Other error sources

There are a number of other error sources not taken into account in our Monte-Carlo error estimates:

• •

Our (TW’s) numerical code only computes the diagnostics to finite accuracy. Comparing diagnostics between calculations done with different numerical resolutions, we have generally excluded any data where the diagnostic computed at our highest resolution (that shown in table 2) differs from that computed at the next-lower resolution listed in table 4 by more than a few percent.

• •

Wiggles are not perfectly separable from the background variation of the diagnostics. That is, the actual frequency spectra of the diagnostics are almost certainly continuous, and cannot be unambiguously separated into low-frequency (background) and high-frequency (wiggle) components.

• •

Our models for the background variation are imperfect. Our constraint that background spline control points must be spaced at least $1.5$ wiggle periods apart keeps the background and wiggles from being degenerate, but at the cost of leaving the background model unable to accurately fit some non-wiggle variations, particularly for small-$m$ (longer-period) wiggles where the spline control points are forced to be quite far apart.

• •

For wiggles in $F_{a}$, our wiggle model (19) doesn’t accurately include the actual spatial variation of the wiggle (QNM) eigenfunctions.

• •

There may be multiple wiggle modes present simultaneously in the diagnostics for a single $m$. Although our wiggle models and fitting software support simultaneously fitting an arbitrary number of wiggles, we have generally not done this, i.e., we have generally only attempted to fit a single-wiggle model for each diagnostic time series.999In a few cases where our best-fitting single-wiggle model’s residuals showed strong systematics, we then proceeded to fit 2-wiggle models. These improved the residuals by at least an order of magnitude. However, all the results presented in section V.1 are based on single-wiggle fits.

We believe that all these other error sources are small, but it is difficult to quantify them.

For each wiggle fit we visually assess the fit residuals to look for obvious systematics. For all results reported here the fit residuals are at least a factor of $10$ smaller than the wiggle amplitude; in most cases they are a factor of $30$ to $100$ smaller. This suggests that our fits are indeed accurately modelling at least the dominant wiggle features of the diagnostics.

IV.2 Gravitational perturbations

In the gravitational case we search for the QNM frequencies in the individual (spheroidal) $(\mathfrak{l}m)$-modes. By looking at individual $(\mathfrak{l}m)$-modes we minimize the number of QNMs that need to be modelled (fitted) simultaneously. Since QNMs appear naturally as spheroidal modes, using spheroidal modes minimizes the amount of “crosstalk” mixed in from neighboring modes. Note that since the spheroidicity of the spheroidal harmonics in the Teukolsky equation depends on the frequency of the modes, the QNMs have complex spheroidicity and will not project perfectly on the corresponding (real) spheroidal modes that appear in the field solutions. Consequently, “crosstalk” between the modes cannot be fully eliminated. Nonetheless, the crosstalk in the spheroidal modes should be significantly smaller than if one were to use the spherical $\ell$-modes.

IV.2.1 Fit models

The fit models used in the gravitational case are very similar to the ones used in the scalar case. For the “watchpoint” diagnostics we use

[TABLE]

In this case the smooth background of the signal is modelled by a simple polynomial in $t$. We maximize the number of linear fit parameters by writing the model as a sum of sines and cosines.

Similarly, the $(\mathfrak{l}m)$-modes contributing to the local gravitational self-force are modelled by

[TABLE]

As in the scalar case, the main shortcoming in this model is the inaccurate modelling of the QNMs’ radial profiles.

Finally, the model for the gravitational waveform at $\mathcal{J}^{+}$ is very similar to the model for the watchpoints,

[TABLE]

IV.2.2 Fitting procedure

The only non-linear parameters in the above models are the QNM frequencies $\omega_{k}$ and decay constants $\alpha_{k}$. Consequently, for fixed $\omega_{k}$ and $\alpha_{k}$ the remaining parameters can be determined efficiently through a linear least squares procedure. This is implemented by using Mathematica’s LinearModelFit routine for each diagnostic on a suitable time window of data. To reduce the impact of un-modelled higher order QNMs, these fits are weighted by $\exp(2\alpha_{1}t)$. Typically, the fits include around 20 terms in the background model and up to 8 QNMs.

The $\omega_{k}$ and $\alpha_{k}$ are then determined by maximizing the sum of the adjusted $R^{2}$ values of all the component fits. This is implemented using Mathematica’s FindMaximum with the PrincipleAxis method. The initial values for $\omega_{k}$ and $\alpha_{k}$ are set by the numerically known corresponding QNMs offset by a random $\mathcal{O}(1\%)$ amount. An indication of the modelling error is obtained by varying the fit window and number of background terms, and determining the spread of the best fits.

V Data and quasinormal-mode fits

V.1 Scalar field

We have surveyed a large number of configurations for Kerr spin $\tilde{a}=0.99$, together with a smaller number of configurations for other Kerr spins; for selected configurations we have fitted (or attempted to fit) wiggle models as described in section IV.1. Tables 2 and 3 describe all the configurations surveyed here, and figure 1 shows the $(\text{periapsis radius},\text{orbital eccentricity})$ phase space of the $\tilde{a}=0.99$ configurations.

Figures 2 and 3 show the wiggles in the scalar-field diagnostics for the $(\tilde{a},p,e)=(0.99,3,0.8)$ configuration, for $m=1$ and $m=4$ respectively. Notice that the wiggles are visible in all the field diagnostics. Notice also the much higher frequency and smaller amplitude of the $m=4$ wiggles.

Figures 4 and 5 show our model fits to these wiggles for $m=1$ and $m=4$ respectively. Notice that in each case the spline control points span a wider range of $t$ or $u$ than the range over which the model is fitted. The $y$ coordinates at the spline control points outside the model-fitting range are still adjusted by the least-squares fitting algorithm, but have only small influences on the model within the fitting range.

Figures 6–8 show the fitted complex frequencies and their Monte-Carlo error estimates, compared to Kerr QNM frequencies calculated by Berti, Cardoso, and Starinets Berti et al. (2009, 2006).101010Data tables downloaded from https://pages.jh.edu/~eberti2/ringdown/ on 19 April 2019.

In each case the fitted frequencies agree with the calculated QNM frequencies, lending further support to the identification of wiggles with QNMs (more precisely, QNMs sampled at the observation points).

V.2 Gravitational field

We now turn our attention to gravitational perturbations. We first consider the same $(\tilde{a},p,e)=(0.99,3,0.8)$ configuration studied in figures 2–7 in the scalar field case. Figure 9 displays both the gravitational self-force at the particle location and the waveform observed at $\mathcal{J}^{+}$. When looking at the local self-force the wiggles are most pronounced in the $F_{\phi}$ component. However, faint traces of wiggles can be found by zooming in on the $F_{t}$ and $F_{r}$ components. We note that the relative amplitudes of the wiggles in the gravitational self-force are much smaller than those in the scalar case for the same orbit (shown in figures 2 and 3).

The waveform observed at $\mathcal{J}^{+}$ depends on the viewing angle. When the system is viewed “face on” (middle panel of figure 9) the waveform is determined by the $m=2$ modes with the $m=\mathfrak{l}=2$ dominating. In this case the wiggles appear as a clear exponentially damped sinusoid. When the system is viewed “edge on” (bottom panel of figure 9), the wiggles have a much more irregular shape, consistent with a much larger collection of $m$’s and $\mathfrak{l}$’s contributing to the wiggles. Also note that while the overall waveform has a much larger amplitude when viewed edge on (due to contributions from higher modes), the observed wiggles are actually stronger when the system is viewed “face on”. This is consistent with the wiggles being dominated by the $\mathfrak{l}=m=2$ mode.

One of the advantages of using a frequency domain approach is that we can easily isolate individual $\mathfrak{l}m$-modes (as defined in section II.2). Figure 10 shows different aspects of the $\mathfrak{l}=m=2$ mode of the gravitational perturbation generated by a particle on our standard $(\tilde{a},p,e)=(0.99,3,0.8)$ configuration. The $\mathfrak{l}=m=2$ mode of the gravitational self-force (top panel) shows the same qualitative features as the full GSF; the $F_{\phi}$ components show the most obvious wiggles with weak wiggles visible in the other components. The $\mathfrak{l}=m=2$ mode of the field observed at $\mathcal{J}^{+}$ shows a clean exponentially decaying sinusoid wiggle just as the full field. In addition, the bottom panel of figure 10 shows the local Teukolsky variable $\Phi^{(22)}_{-2}$ at two watch points located on the background Kerr spacetime’s symmetry axis at radii corresponding to the periapsis and apoapsis of the particle orbit. These show the cleanest wiggles of any of our diagnostics.

To test our hypothesis that the observed wiggles are, in fact, QNM excitations, we perform a global fit of our three field diagnostics (local gravitational self-force, field at $\mathcal{J}^{+}$, and field at watchpoints) following the methodology set out in section IV.2. Table V.2 summarizes the results for some low order $\mathfrak{l}m$-modes. In each, case we recover the principal QNM frequency and damping time of the gravitational field within the estimated numerical precision of the fits. This provides yet more evidence for our hypothesis that the observed wiggles are QNM excitations. Note that while our fits include multiple QNMs, we do not conclusively recover any of the modes beyond the principal mode. (More precisely, we find that the estimated numerical errors of the recovered complex frequencies are comparable to the variation of the initial seed for the optimization.) Not including the higher modes, however, led to observable bias in the recovery of the principal QNMs.