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

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.


I. INTRODUCTION
Consider a small (compact) body of mass µM (with 0 < µ 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 longstanding 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 ∼ 1-100M orbits a massive black hole of mass ∼ 10 5 -10 7 M , 1 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 ( [1][2][3][4]). 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 1 M denotes the solar mass.
unbound (scattering) flyby can excite quasinormal modes (QNMs) of the background black hole. Kojima and Nakamura [5] 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 quasinormal mode". Their figure 3(b) shows an example of the QNM oscillations in the radiated gravitational waves at J + . Burko and Khanna [6] 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 [7] 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 [8] (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 ("wig-arXiv:1906.06791v3 [gr-qc] 9 Nov 2019 gles") 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 [9,10] presented fits of dampedsinusoid 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 [11] calculated the scalarfield self-force and the radiated field at 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 J + , fitting these to a superposition of = m = 1, = m = 2, = m = 3, and = 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 [12] 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. [9,10] 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) O(µ) field perturbations near the Kerr black hole, O(µ) radiation at J + , and O(µ 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, [13,Section 6.10]; [14][15][16][17][18][19] and references therein) are not reliably accurate. Since the timescale for radiation reaction to shrink the orbit is very long (∼ µ −1 M ) while the required resolution near the small body is very high (∼ µM ), a direct "numerical relativity" integration of the Einstein equations (see, for example, [20][21][22][23][24] and references therein) would be prohibitively expensive (and probably insufficiently accurate) for this problem. 2 Instead, we use black hole perturbation theory, treating the small body as an O(µ) 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 [8,33] extended to compute additional field diagnostics. This code uses an effective-source regularization followed by an e imφ Fourier-mode decomposition and a separate 2+1-dimensional timedomain 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 J + .
• Our gravitational-field results were obtained using van de Meent's code [34][35][36][37]. This code obtains the local metric perturbation in the frequency domain by reconstructing the metric perturbation from the Weyl scalar Ψ 4 , followed by -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 Ψ 4 at 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 J + . Except for a few anomalous results for near-extremal Kerr spacetimes (dimensionless Kerr spins 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 A summarizes our notation and our parameterization of bound geodesic orbits in Kerr spacetime. Section II briefly summarizes our calculations of scalarfield (section II A) and gravitational (section II B) perturbations of Kerr spacetime. Section III describes our localand 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 scalarfield and gravitational perturbations of Kerr spacetime, and QNM-model fits to this data. Finally, section VI discusses our results and presents our conclusions.

A. Notation and parameterization of Kerr geodesics
We generally follow the sign and notation conventions of Wald [38], 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. ∇ 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 . := ∇ a ∇ a is the 4-dimensional (scalar) wave operator [39,40].
We use Boyer-Lindquist coordinates (t, r, θ, φ) on Kerr spacetime, defined by the line element where M is the black hole's mass,ã = J/M 2 is the black hole's dimensionless spin (limited to |ã| < 1), Σ := r 2 + M 2ã2 cos 2 θ, and ∆ := r 2 − 2M r + M 2ã2 . We also define a := Mã (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 1 + √ 1 −ã 2 and the inner horizon is the coordinate (See footnote 3 for a discussion of TW's coordinate compactification near the event horizon and J + .) Following Ref. [41], we define the tortoise coordinate r * by and fix the additive constant by choosing u := t − r * is thus an outgoing null coordinate. The particle's worldline is x i = x i 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φ/dt > 0 direction, withã > 0 for prograde orbits andã < 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 (1 + e) and r max = pM (1 − e) , so that the particle orbit is given by for a suitable orbital-phase function χ r . We refer to the combination of a spacetime and a (bound geodesic) particle orbit as a "configuration", and parameterize it with the triplet (ã, p, e). We define T r to be the coordinatetime 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
A. Scalar-field perturbations of Kerr spacetime In this section we summarize TW's scalar-field calculations [8,33]. These authors consider a real scalar field Φ physical satisfying the wave equation in Kerr spacetime, sourced by a point "particle" of scalar charge q which is taken to move on a (pre-specified) equatorial geodesic orbit. Φ physical satisfies outgoing-radiation (retarded) boundary conditions at J + . This system provides a toy model of the full gravitational perturbation problem without the complexity of gauge choice. Because Φ physical diverges on the particle worldline, (2.1) must be regularized. TW use an effective-source regularization of the type introduced by Barack and Golbourn [42] and Vega and Detweiler [43] (see [44] for a review), using the puncture function and effective source described by Wardell et al. [33]. In a neighborhood of the particle worldline, TW define a (real) regularized scalar field Φ regularized = Φ physical − Φ puncture , where Φ puncture is a suitably-chosen approximation to the Detweiler-Whiting singular field [45]. The (4-vector) selfforce acting on the particle is then given by (2.2) TW make an azimuthal Fourier decompositions of all the spacetime scalar fields into complex e imφ modes, where the extra factor of 1/r is introduced for computational convenience and whereφ := φ + f (r) is an "untwisted" azimuthal coordinate, with For each m-mode, TW introduce a finite worldtube surrounding the particle worldline in (t, r, θ) space. For particle orbits with significant eccentricity (e 0.2) the worldtube (now viewed as a region in (r, θ) in each t = 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, θ), with size 10M in r * by approximately 0.25 radians in θ, symmetric about the equatorial plane, and kept centered on the particle to within 0.25M in r * as the particle moves.
(2.5) using a time-domain 2+1-dimensional finite-difference numerical evolution with mesh refinement. Because the (Kerr) background spacetime is axisymmetric, the Fourier modes ϕ m evolve independently -there is no mixing of the modes. Because the physical scalar field Φ is real, only the m ≥ 0 modes need to be explicitly computed; the m < 0 modes may be obtained by symmetry.
Corresponding to the Fourier decomposition (2.3), the self-force (2.2) can be written as a similar sum of e imφ modes, where each F (m) a may be computed from the corresponding (ϕ m ) regularized field near the particle. TW compute a finite set of modes (typically 0 ≤ m ≤ 20) and estimate the m > 20 contributions to the sum (2.6) via a large-m asymptotic series.

B. Gravitational perturbations of Kerr spacetime
In this section we summarize the metric reconstruction approach used by van de Meent [34][35][36][37] to calculate gravitational perturbation of Kerr spacetime generated by particles on bound geodesic orbits. This approach starts from the the spin-(-2) Teukolsky variable, 3 More 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 ≤ r ≤ rmax containing the particle orbit, but which are compactified near the event horizon and J + . T is a Bonditype retarded time coordinate at 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 J + .
where Ψ 4 is one of the Weyl scalars. As shown by Teukolsky [40], 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, where the R lmω are solutions of the radial Teukolsky equation, the S lmω are spin-weight spheroidal harmonics, and 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 [54,55].
Remarkably, Φ −2 contains almost all information about the corresponding perturbation of the metric [56], and in vacuum it is possible to reconstruct the metric perturbation in a radiation gauge [57][58][59][60]. As detailed in Refs. [34][35][36][37], 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 where the R ± lmω (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 [61]), 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 -behavior of the singular field [62,63], we need to re-expand the spheroidal l-modes to spherical -modes, With this re-expansion, the local gravitational self force is given by where, as shown in [64], we can use the Lorenz-gauge B a parameter given in [62,63]. The metric reconstruction formalism can only recover parts of the metric that contribute to Ψ 4 . This means that the reconstructed metric carries an ambiguity, which can be shown [56] 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 [65,66]. 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. [35], the method of extended homogeneous solutions leads to large cancellations between different (low) frequency modes for high-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 selfforce 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 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-l QNMs, then these will contribute (mostly) to the low-l modes of the gravitational metric perturbation. For this purpose, we define the individual lm modes of the Teukolsky variable 12) and the gravitational self-force (2.13) These are much easier to compute, and for low 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 pres-ence of wiggles in the physical scalar field or metric perturbation implies the presence of wiggles in some or all of the e imφ scalar-field modes or (lm) metric-perturbation modes (respectively), and vice versa. 4 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 I 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, θ, φ) = 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 (ϕ m ) regularized . However, this is only defined within the worldtube, so for orbits with significant eccentricity (where the worldtube moves in (r, θ) during the particle orbit) any given watchpoint may lie outside the worldtube (and thus leave (ϕ m ) 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 Φ (lm) −2 on the symmetry axis of the background Kerr spacetime (θ = 0) and the equatorial plane (θ = π/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φ (scalar-field) and (lm) (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 (φ) 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 J + . These have the advantage of being physically observable and of allowing the e imφ (scalar-field) and (lm) (gravitational) mode decompositions to be defined in a weak-field region (for the gravitational case, this also avoids any gauge dependence).
At 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 ∂ tt (ϕ m ) physical evaluated in the equatorial plane (scalar field) 5 and Ψ 4 evaluated either on the z axis or in the equatorial plane (gravitational).

IV. QUASINORMAL-MODE MODELS AND FITS
A. Scalar-field perturbations

Perturbations at a fixed spatial position
Consider first the case of wiggles in an individual e imφ Fourier mode of the regularized scalar field, observed at a fixed "watchpoint" spatial position in the strong-field region. We consider the model where B is a spline function representing the slowlyvarying "background" variation of the scalar field, k indexes the damped-sinusoids included in the model, A (k) , λ (k) , P (k) , and η (k) are respectively the amplitude, damping rate, period, and phase offset of each dampedsinusoid, and the subscript 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) := max k P (k) is the period of the longest-period damped-sinusoid in the model.

Perturbations at the particle position
Consider next the case of wiggles in the radiationreaction 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 φ components for our equatorial orbits), and the field perturbation is being sampled at a time-varying position. Analogously to (4.1), we thus consider the model 5 Recall (cf. footnote 3) that in our scalar-field computations, t is a Bondi-type retarded time coordinate at J + .
where we now parameterize the particle's motion using the retarded time coordinate u, 6 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, are now respectively the 4-vector amplitude, damping rate, period, and 4-vector phase offset of each damped-sinusoid, and the subscript 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φ particle (u) term models the variation in oscillation phase due to particle's (i.e., the sampling point's) motion in φ. The 1/r 3 particle (u) and 1/r 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.) 3. Perturbations at J + Finally, consider the case of wiggles in the radiated (physical) field at J + . Because of the hyperboloidal time slices, we observe these at finite coordinate time t (the J + time has an arbitrary offset relative to the strongfield 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 (4.1) and (4.2), we thus consider the model where the meanings of all terms (and the non-degeneracy condition on the background spline) are the same as in (4.1).

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 P (k) , τ (k) (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 fit ≥ 4P (max) , we randomly choose N 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) . 7 Then we repeat the wiggle-model fit for each subinterval. 8 The ensemble of the N trial sets of Monte-Carlotrial 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 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 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.

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 II) differs from that computed at the next-lower resolution listed in table IV 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 nonwiggle variations, particularly for small-m (longerperiod) wiggles where the spline control points are forced to be quite far apart.
• For wiggles in F a , our wiggle model (4.2) 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. 9 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.

B. Gravitational perturbations
In the gravitational case we search for the QNM frequencies in the individual (spheroidal) (lm)-modes. By looking at individual (lm)-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 -modes.

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 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 (lm)-modes contributing to the local gravitational self-force are modelled by (4.5) 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 J + is very similar to the model for the watchpoints,

Fitting procedure
The only non-linear parameters in the above models are the QNM frequencies ω k and decay constants α k . Consequently, for fixed ω k and α 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α 1 t). Typically, the fits include around 20 terms in the background model and up to 8 QNMs.
The ω k and α 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 ω k and α k are set by the numerically known corresponding QNMs offset by a random 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. This figure shows the phase space of the Kerr spinã = 0.99 scalar-field configurations presented here, plotted in terms of the periapsis radius rmin and the eccentricity e. The shaded region at the left shows orbits with periapsis inside the horizon. The light ring, innermost bound circular orbit (IBCO), innermost stable circular orbit (ISCO), and the locus of marginally stable orbits are also shown. The meanings of the plot symbols are described in detail in Table III. Nameã p e Tr rminφperiapsis r watchpoint Resolution m=1 m=2 m=3 m=4 m=5 m=6  , within this by increasing particle periapsis radius rmin, and within this by decreasing particle orbital eccentricity e. (ã, p, e) describe the configuration. Tr is the coordinate-time period of the orbit's radial motion, rmin is the particle's periapsis r coordinate, andφperiapsis is the particle's angular 3-velocity in φ (i.e., dφ/dt) at the point of periapsis. r watchpoint is the 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 IV. The final sets of columns describe the presence or absence of wiggles in our field diagnostics for m = 1 through m = 6; for each m there are 3 columns labeled "w", "F ", and "J ", describing (respectively) wiggles in the regularized field (ϕm) regularized at a fixed "watchpoint" (located on the equator at r = r watchpoint ), wiggles in the self-force Fa, and wiggles in the second time derivative of the physical field, ∂tt (ϕm) physical , at J + . The meanings of the symbols in these columns are described in detail in table III; briefly, ⊕ or + means that we observed wiggles and fit them to the appropriate model described in section IV A, ∼ means that we observed wiggle-like oscillations but did not attempt a quantitative fit, ♦ means that we observed oscillations which might be wiggles but were not clearly separated from the background variation, and × means that we did not observe wiggles (n.b. this does not prove the absence of wiggles, only that we did not observe them).
Symbol Meaning ⊕ we observed oscillations in the diagnostics, successfully fit the appropriate wiggle model described in section IV A over a t or u range of ≥ 4 wiggle periods, and performed the Monte-Carlo error analysis described in section IV A 5.
+ we observed oscillations in the diagnostics and successfully fit the appropriate wiggle model described in section IV A, but the model was fitted over too short a t or u range (< 4 wiggle periods) for the Monte-Carlo error analysis described in section IV A 5 to be used ∼ 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 × 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 than the range over which the model is fitted. The y coordinates at the spline control points outside the modelfitting 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 [68,69]. 10 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).

B. Gravitational field
We now turn our attention to gravitational perturbations. We first consider the same (ã, 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 J + . When looking at the local self-force the wiggles are most pronounced in the F φ 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 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 = 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 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 l = m = 2 mode.
One of the advantages of using a frequency domain approach is that we can easily isolate individual lm-modes (as defined in section II B). Figure 10 shows different aspects of the l = m = 2 mode of the gravitational perturbation generated by a particle on our standard  . The lower subfigure shows the physical scalar field and its second time derivative, on the equator at J + ; the fields at other viewing angles are very similar in shape but with smaller amplitudes. Due to the orbital precession, the J + fields are not "periodic with period Tr". Note that the zero point of the J + time coordinate does not correspond to periapsis. . The lower subfigure shows the physical scalar field and its second time derivative, on the equator at J + ; the fields at other viewing angles are very similar in shape but with smaller amplitudes. Due to the orbital precession, the J + fields are not "periodic with period Tr". Note that the zero point of the J + time coordinate does not correspond to periapsis.  mode of the field observed at 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 Φ (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 J + , and field at watchpoints) following the methodology set out in section IV B. Table V summarizes the results for some low order lm-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.    In this section we study how the strength of the QNM excitations in the gravitational field depends on the parameters of the orbit. For this investigation we leverage the ease with which the (lm) modes of the gravitational field can be computed using our frequency domain code. As a measure of the strength of the excitations we take the amplitudes A  The three non-zero components of the gravitational self-force at the particle. The F φ component has very prominent wiggles. The wiggles in the other two components are present but only visible after zooming in. (middle) The value of ψ4 observed at J + when the system is viewed "face on" (θ observer = 0). We observe very clean exponentially decaying wiggles right up to the next burst produced when the particle approaches the central black hole again. (bottom) The value of ψ4 observed at J + when the system is viewed "edge on" (θ observer = π/2). In this case the wiggle pattern is much more complex as multiple (l, m)-modes contribute.
the QNMs exhibits a strong dependence on the alignment of the discrete orbital frequency spectrum of the orbit with the QNM mode. If strong localized "resonances" between the orbit frequency spectrum and QNM excitations were to exist, these could have significant impact on waveform modelling strategies, as they would hamper an attempt to apply reduced order modelling to build efficient waveforms. On the other hand, such a phenomenon might lead to interesting and rich dynamics. To quantify the alignment between a QNM and the particle orbit we define i.e., δ (k) is the distance between the QNM frequency ω (k) and the nearest line in the orbit's frequency spectrum, normalized such that δ (k) = 0 corresponds to maximal alignment and δ (k) = 1/2 to maximal misalignment.
In figure 11 we examine the dependence of the QNM amplitude (for the least-damped (k = 1) l = m = 2 QNM) on the spin a of the background Kerr spacetime while keeping eccentricity and ratio of the orbital frequencies (matching the eccentricity and frequency ratio for the (ã, p, e) = (0.99, 3, 0.8) configuration). We see that the dependence of the QNM excitation amplitude on the Kerr spin a is very smooth, with no noticeable dependence on the spectrum misalignment parameter δ. No- tice that the QNM excitations persist for negative spins (i.e., retrograde orbits), although they become exceedingly weak. Figure 12 explores the dependence of the amplitude |A (1) | of the least damped (k = 1) l = m = 2 QNM on the particle's orbital eccentricity. For this exploration we keep the spin of the background Kerr spacetime fixed at a = 0.95M , and we fix the periapsis distance at r min = 1.85M . The relationship between |A (1) | and e appears almost linear by eye, with slight deviations both at high and low eccentricity. If we subtract off the dominant trend in the form of a quintic fit in e, we see what appears to be a systematic trend where orbits with δ (1) = 0 have a slightly larger amplitude |A (1) | than orbits with δ (1) = 1/2. This difference becomes stronger for low eccentricity orbits. This latter effect is consistent with the frequency spectrum of the orbit becoming sparser at lower eccentricities. However, we stress that this effect is very small, with the variation of the QNM amplitude |A (1) | due to changing δ (1) being only about 1 part in 10 3 .
Finally, figure 13 explores the relation between the amplitude |A (1) | of the least damped (k = 1) l = m = 2 QNM and the particle's inverse periapsis distance M/r min , keeping the Kerr spin (a = 0.95M ) and particle eccentricity (e = 0.8) fixed. As is to be expected, the amplitude |A (1) | drops off sharply as we increase the particle periapsis radius.
We emphasize that the overall shape of the plots in one should not read too much into the shapes themselves. However, there are three main lessons that we learn from this investigation that do not depend on the choice of u ref : • The amplitudes of the wiggles depend smoothly on the Kerr spin ( figure 11) and orbital parameters (figures 12 and 13). In particular, no fine-tuning is needed for wiggles to appear.
• The wiggles are strongest for high spin and prograde particle orbits with high eccentricity and low periapsis distance. However, there is no indication that they will completely disappear in any region of the parameter space (although they may become very difficult to separate from the rest of the field due to low amplitudes, high damping rates, and/or longer periods).
• The effect of aligning the orbital frequencies with the QNM frequencies is very small, and decreases still further when the orbital spectrum becomes denser for more eccentric orbits.

VI. DISCUSSION AND CONCLUSIONS
In this paper we study an interesting class of features first observed in the scalar self-force for point particles in orbit in Kerr spacetime [8]. That study identified the feature, introduced the term "wiggles", and argued that it was in some (unspecified) manner "caused by the particle's close passage by the large black hole", but did not attempt to attribute it to any particular physical origin. More recently Refs. [9][10][11] have shown further examples of wiggles, demonstrated that wiggles' complex frequencies agree with known Kerr quasinormal-mode (QNM) frequencies, and concluded that wiggles are in fact "just" a sampling at the measurement point(s) of Kerr QNMs excited by the particle.
Here we survey the phenomenology of wiggles for both the scalar-field and gravitational cases, across a range of Kerr spins and particle orbits. In both the scalar-field and gravitational cases we find that wiggles are essentially a generic phenomenon, i.e., they occur over a wide range of configuration space without any "fine-tuning" of parameters. Wiggles are observable in field perturbations at fixed spatial positions, in the radiation-reaction "self-force", and in the radiated fields at J + .
In both the scalar-field and gravitational cases we find that at all observed locations in spacetime, wiggles can be quantitatively fit by models of QNMs sampled at the observation points. In particular, in both the scalar-field and gravitational cases our fitted wiggle frequencies agree well [in both real (oscillatory) and imaginary (damping) parts] with Kerr QNM frequencies calculated by Berti, Cardoso, and Starinets [68,69] (figures 6-8 and table V). The appearance of pronounced wiggles appears to rely on three key aspects of the configuration of the system: • A highly spinning central (Kerr) black hole (the closer toã = 1, the more pronounced the effect).
• A highly eccentric prograde orbit for the particle (the closer to e = 1, the more pronounced the effect).
• A close periapsis passage by the particle (the closer to the light ring, the more pronounced the effect). 11 This is not surprising: highly spinning black holes have much longer-lived QNMs than those with low spin. Increasing eccentricity of the particle orbit does three things: it increases the strength of the perturbation at the periapsis, it widens the frequency spectrum of the perturbation (increasing the overlap with the QNMs), and it provides a natural "quiet" period when the particle approaches its apoapsis, during which the QNMs can more easily be observed. Finally, bringing the particle periapsis closer to the light ring allows the perturbation to deposit more energy in the QNMs. (QNMs 11 We refer here to the prograde light ring; we observe only very small QNM excitation when the particle periapsis is close to the retrograde light ring. in Kerr spacetime are readily excited by orbits near the light ring [68,70,71].) Interestingly, we find that the amplitude with which wiggles are excited does not depend sensitively on the particle's precise orbital motion near periapsis. Notably, we find that the wiggle amplitude varies smoothly and monotonically with the particle periapsis radius and orbital eccentricity (and with the Kerr spin).
We have not attempted to carefully delineate the exact boundaries of the region in configuration space where wiggles occur (even assuming that there are, in fact, configurations with no QNM excitation, which is not obvious). It is likely that different modelling/fitting schemes could observe and fit low-amplitude and/or rapidlydamped wiggles even in some cases where we fail to observe them (e.g., the × cases in table II). For example, figure 11 strongly suggests that although the wiggle amplitude is very small in some cases, wiggles are present for all Kerr spins along this sequence of orbits, including retrograde as well as prograde orbits.
Our scalar-field wiggle modelling/fitting scheme is (deliberately) quite conservative in requiring visual observation of a wiggle in a time-series plot of the original diagnostic. This requirement reduces the risk of false positives (where we would misidentify a fitting or background-spline artifact as a wiggle), at the cost of reducing our sensitivity to low-amplitude and/or rapidlydamped wiggles.
An interesting example of these factors at play is the (ã, p, e) = (0.99, 8, 0.8) scalar-field configuration, for which Nasipak, Osburn, and Evans [11] observed and fitted = m = 1, = m = 2, = m = 3, and = m = 4 wiggle (QNM) modes. Their figures 8 and 9 show the = m = 4 wiggle as having an amplitude approximately 10 9 times smaller than the = m = 1 wiggle; this is only detectable by virtue of the high accuracy and low numerical noise level of their frequency-domain code. In contrast, for this configuration we observed wiggles for m = 1 but not for m ≥ 2; this is likely because even the m = 2 wiggles are already too low in amplitude to be visually observable in the original time series.
Existing astrophysical models of extreme mass ratio binaries [72,73] and observations of highly spinning black holes [74] suggest that it is quite reasonable to expect some fraction of EMRIs to fall within the region of parameter space where wiggles are excited with significant amplitude. (Both the magnitude of this fraction and the absolute numbers of such systems are still very uncertain.) Given that the QNM excitations appear not just in the local self-force, but also in the gravitational waveform, a natural question is whether they could be experimentally observed by LISA or other detectors. While this is certainly possible in principle, there are two considerations which make it less likely in practice. Most importantly, the effect is quite weak in all but the most extreme cases. In most of the gravitational-field cases investigated here, it was necessary to zoom in on plots in order to see the wiggles visually, reflecting the fact that their magnitude represents at most a few percent of the total signal. A second consideration in terms of detectability is that the dynamical evolution of EMRIs may tend to avoid the wiggles region of parameter space (e.g., if most EMRIs evolve to low orbital eccentricities while still at relatively large periastron radii). This would imply that the event rate for detectable EMRI wiggles would be quite low.
Despite these concerns, it would be worthwhile to conduct a more thorough study to quantitatively address the question of detectability of QNM wiggles by LISA or future next-generation gravitational-wave detectors. It may even be the case that advanced data analysis techniques could be used to boost the detectability. For example, although an individual wiggle is weak, it will repeat for each orbit throughout the entire inspiral. As noted by Ref. [11], wiggles will appear with almost the same frequency throughout the inspiral (the QNM frequency only depends on the mass and spin parameters of the larger black hole, and these change very little during the inspiral). Moreover, this frequency is much higher than than the main orbital frequency, potentially making it easier to separate these signal components in data analysis.
The analysis done here has been somewhat post-hoc, in that we first identified a feature in the signal and then fit this feature to a damped sinusoid representing a QNM ringdown. Our intuitive interpretation of this QNM ring-down is that it is a result of strong QNM excitation near periapsis, which is then encountered over an extended period later in the orbit. The self-force in curved spacetimes arises from nonlocal interactions of the object with its self-field, which was generated in the object's past and scattered off the spacetime curvature. The association of wiggles with QNM excitations suggest that they represent a situation where this nonlocal nature of the self-force is particularly apparent. To more explicitly develop this interpretation, it may be informative to attempt the analysis in the other direction, first by starting with a model for a QNM excitation from a burst of radiation generated near a periapsis passage, and then comparing such a model to the observed signal. This approach would allow one to pinpoint where in the orbit the QNM excitation occurs, would give a deeper understanding of the effect, and may even provide a link to geometric features such as caustics and the propagation of waves on black hole spacetimes. A Green function approach [75,76] would be a natural choice for such a study, but is quite distinct from the methodology used in this paper so we leave it for future work.
In this work we have focused mostly on systems with somewhat realistic Kerr spins J/M 2 0.999. Initial investigations of the near-extremal regime suggest a rich phenomenology, involving many different QNMs at fixed l and m. One puzzling result is the l = m = 2 mode generated at J + by a particle orbiting a black hole with spin J/M 2 = 0.999 99 and the same orbital frequencies as the (ã, p, e) = (0.99, 3, 0.8) orbit, shown in figure 14. One of the puzzling aspects of this waveform is that the decaying wiggles in the highlighted area have a frequency between 0.93M −1 and 0.98M −1 (depending on where it is measured), while the nearest QNMs all have frequency close to 0.995M −1 . Whether this is the result of some complicated collective behaviour of the QNMs or some new physical effect is currently unclear, and should be investigated in future works. 12