Repeated faint quasinormal bursts in extreme-mass-ratio inspiral waveforms: Evidence from frequency-domain scalar self-force calculations on generic Kerr orbits
Zachary Nasipak, Thomas Osburn, Charles R. Evans
TL;DR
This paper develops a frequency-domain code to compute scalar self-force on particles in Kerr spacetime, revealing quasinormal bursts in waveforms that could enhance gravitational wave data analysis for extreme-mass-ratio inspirals.
Contribution
The authors extend spectral source integration techniques to Kerr orbits and confirm the presence of quasinormal bursts in waveforms, a novel observation in frequency-domain calculations.
Findings
QNBs appear after each periastron passage in waveforms.
QNBs are superpositions of least-damped quasinormal modes.
QNBs could be observable in high SNR EMRIs, aiding parameter estimation.
Abstract
We report development of a code to calculate the scalar self-force on a scalar-charged particle moving on generic bound orbits in the Kerr spacetime. The scalar self-force model allows rapid development of computational techniques relevant to generic gravitational extreme-mass-ratio inspirals (EMRIs). Our frequency-domain calculations are made with arbitrary numerical precision code written in \textsc{Mathematica}. We extend spectral source integration techniques to the Kerr spacetime, increasing computational efficiency. We model orbits with nearly arbitrary inclinations and eccentricities up to . This effort extends earlier work by Warburton and Barack where motion was restricted to the equatorial plane or to inclined spherical orbits. Consistent with a recent discovery by Thornburg and Wardell \cite{ThorWard17} in time-domain calculations, we…
| model | ||||
|---|---|---|---|---|
| base | ||||
| large | ||||
| large | ||||
| large |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Repeated faint quasinormal bursts in extreme-mass-ratio-inspiral
waveforms: Evidence from frequency-domain scalar self-force calculations on generic Kerr orbits
Zachary Nasipak
Department of Physics and Astronomy, University of North Carolina, Chapel Hill, North Carolina 27599, USA
Thomas Osburn
Department of Physics and Astronomy, State University of New York at Geneseo, New York 14454, USA
Department of Physics and Astronomy, University of North Carolina, Chapel Hill, North Carolina 27599, USA
Charles R. Evans
Department of Physics and Astronomy, University of North Carolina, Chapel Hill, North Carolina 27599, USA
Abstract
We report development of a code to calculate the scalar self-force on a scalar-charged particle moving on generic bound orbits in the Kerr spacetime. The scalar self-force model allows rapid development of computational techniques relevant to generic gravitational extreme-mass-ratio inspirals (EMRIs). Our frequency-domain calculations are made with arbitrary numerical precision code written in Mathematica. We extend spectral source integration techniques to the Kerr spacetime, increasing computational efficiency. We model orbits with nearly arbitrary inclinations and eccentricities up to . This effort extends earlier work by Warburton and Barack where motion was restricted to the equatorial plane or to inclined spherical orbits. Consistent with a recent discovery by Thornburg and Wardell (2017) in time-domain calculations, we observe self-force oscillations during the radially outbound portion of highly eccentric orbits around a rapidly rotating black hole. As noted previously, these oscillations reflect coupling into the self-force by quasinormal modes excited during pericenter passage. Our results confirm the effect with a frequency-domain code. More importantly, we find that quasinormal bursts (QNBs) appear directly in the waveform following each periastron passage. These faint bursts are shown to be a superposition of the least-damped overtone (i.e., fundamental) of at least four () quasinormal modes. Our results suggest that QNBs should appear in gravitational waveforms, and thus provide a gauge-invariant signal. Potentially observable in high signal-to-noise ratio EMRIs, QNBs would provide high-frequency components to the parameter estimation problem that would complement low-frequency elements of the waveform.
pacs:
04.25.dg, 04.30.-w, 04.25.Nx, 04.30.Db
I Introduction
Recent direct detections of gravitational waves have inaugurated a new branch of multimessenger astronomy. These observations of compact binary mergers by advanced LIGO and advanced Virgo Abbott, B. P. et al. (2016a) (The LIGO Scientific Collaboration and VIRGO Collaboration); Abbott, B. P. et al. (2016b) (The LIGO Scientific Collaboration and VIRGO Collaboration); Abbott, B. P. et al. (2017a) (The LIGO Scientific Collaboration and VIRGO Collaboration); Abbott, B. P. et al. (2017b) (The LIGO Scientific Collaboration and VIRGO Collaboration); Abbott, B. P. et al. (2017c) (The LIGO Scientific Collaboration and VIRGO Collaboration); Abbott, B. P. et al. (2017d) (The LIGO Scientific Collaboration and VIRGO Collaboration) have led to discovery of a new class of heavy stellar mass black holes, confirmed the primary site of the r-process for creation of heavy elements, provided strong-field tests of general relativity Abbott, B. P. et al. (2016c) (The LIGO Scientific Collaboration and VIRGO Collaboration), placed limits on the astrophysical environments of compact binaries Abbott, B. P. et al. (2016d) (The LIGO Scientific Collaboration and VIRGO Collaboration), and made connection with other parts of astronomy Abbott, B. P. et al. (2017e) (The LIGO Scientific Collaboration and VIRGO Collaboration); Abbott, B. P. et al. (2017f) (The LIGO Scientific Collaboration and VIRGO Collaboration). Detection rates are poised to increase following recent sensitivity enhancements in LIGO and Virgo, eventual completion of KAGRA KAG , and development of LIGO-India Unnikrishnan (2013). Ground-based detectors will be complemented by the LISA mission Amaro-Seoane et al. (2013, 2017); NASA (2011); ESA (2012) recently approved by the European Space Agency, which will be sensitive to gravitational waves in a lower frequency band ( Hz).
A prime target for LISA will be extreme-mass-ratio inspirals (EMRIs) consisting of a small compact object of mass (neutron star or black hole) in orbit about a supermassive black hole () Berry et al. (2019). With a small mass ratio , a gradual, adiabatic inspiral occurs, which provides a natural application of black hole perturbation theory (BHPT) and attendant gravitational self-force (GSF) calculations. Once an EMRI crosses into the detector passband, its orbital motion will accumulate a total change in phase of order radians prior to merger, with the implication that the small black hole will skim close to the event horizon hundreds of thousands of times and provide an unprecedented test of general relativity Vigeland and Hughes (2010); Barack and Cutler (2007); Brown et al. (2007); Amaro-Seoane et al. (2014); Babak et al. (2017). LISA will also serve as a cosmological probe, detecting EMRIs out to redshifts of Amaro-Seoane et al. (2007); Barack (2009); Amaro-Seoane et al. (2014); Babak et al. (2017).
Waveform templates produced from self-force calculations will be useful in aiding signal detection of EMRIs and be essential for parameter estimation, supplanting kludge waveforms derived from adiabatic inspiral calculations Babak et al. (2017). Long term self-force inspiral calculations of Schwarzschild EMRIs are well advanced Warburton et al. (2012); Osburn et al. (2016); Warburton et al. (2017); van de Meent and Warburton (2018), tracking the accumulated orbital or gravitational wave phase to accuracies better than due to all first-order-in-the-mass-ratio effects at post-1-adiabatic order Hinderer and Flanagan (2008), lacking only the orbit-averaged dissipative part of the second-order self-force. Progress is also being made on understanding and calculating the second-order GSF Pound (2012a, b); Pound and Miller (2014); Wardell and Warburton (2015); Miller et al. (2016); Pound (2017); Moxon and Flanagan (2018). In the case of Kerr EMRIs, steady developments have been made in GSF calculations for circular and bound equatorial orbits Shah et al. (2012); Isoyama et al. (2014); van de Meent and Shah (2015); van de Meent (2016); Merlin et al. (2016); Fujita et al. (2017); Barack and Giudice (2017). Progress has now been reported van de Meent (2018) in calculating the GSF on generic Kerr orbits. In principle this latest self-force result could serve as the basis for long-term inspiral models of astrophysically relevant EMRIs, but prospects are dimmed at present by high computational costs of these GSF calculations.
In the past, the scalar field self-force analogue Quinn (2000), where a scalar point charge orbiting a black hole sources a scalar wave that acts back on the charge, has frequently been used as a simplified model to provide understanding and to develop tools for use in the gravitational case. The scalar self-force (SSF) has been computed in Schwarzschild spacetime Burko (2000); Wiseman (2000); Burko (2000); Barack and Burko (2000); Detweiler et al. (2003); Diaz-Rivera et al. (2004); Haas (2007); Vega and Detweiler (2008); Vega et al. (2009); Casals et al. (2009); Cañizares and Sopuerta (2009); Canizares et al. (2010); Wardell et al. (2012); Diener et al. (2012); Vega et al. (2013); Wardell et al. (2014) and in Kerr spacetime using frequency domain (FD) Warburton and Barack (2010, 2011); Warburton (2015) and time domain (TD) calculations Thornburg and Wardell (2017a).
The present work generalizes the previous FD SSF calculations of Warburton and Barack to arbitrary eccentric inclined orbits in Kerr spacetime. Part of our procedure involves calculating modes with the Mano-Suzuki-Takasugi (MST) analytic function expansion approach Mano et al. (1996a, b); Sasaki and Tagoshi (2003) using Mathematica. The code serves as a test bed for developing more advanced physical and numerical techniques to aid downstream work in making generic Kerr GSF calculations more practical. For example, in this paper we adapted spectral source integration (SSI) Hopper et al. (2015) to the Kerr generic-orbit source problem, significantly optimizing computational efficiency. Physically, we are able to explore rapidly the SSF in interesting high eccentricity and high black hole spin systems and follow-up work will examine the behavior of resonant-orbit configurations.
A primary physical result in this paper is confirming with our FD calculations the existence of quasinormal mode excitations in the self-force, which was discovered by Thornburg in TD SSF simulations of highly eccentric Kerr orbits. This finding was discussed in a series of talks Thornburg and Wardell (2014, 2016, 2017b) by Thornburg and reported in a paper by Thornburg and Wardell Thornburg and Wardell (2017a). Oscillations are observed in the self-force during the outbound portion of certain highly eccentric orbits following periastron passage near a rapidly rotating black hole. These oscillations were confirmed to fit the least-damped overtone of the quasinormal mode. We see precisely the same behavior in the self-force in our FD calculations (see V.2) of a similar highly eccentric () equatorial orbit about a rapidly rotating () primary.
More interestingly, we decided to take a look at the waveform in this same model to see if the excitation is imprinted in an asymptotically accessible signal. Confirming our expectation, it is indeed possible to discern repeated (albeit faint) quasinormal bursts (QNBs) in the waveform following each periastron passage. Fig. 1 shows the asymptotic waveform over a period of two radial librations at several observer angles. Without further processing, no quasinormal oscillations are directly apparent. However, by high-pass filtering or otherwise enhancing high frequencies in the signal, we can make the low-level QNBs evident. One particular way of enhancing high frequencies is shown in Fig. 2 where the log (base 10) of the absolute value of the second time derivative of the waveform is plotted. (Computing the second derivative is reminiscent of some numerical relativity codes where, to extract gravitational radiation, is first obtained, from which the waveforms are derived by integrating twice or by Fourier processing.) Now the QNBs are revealed, superimposed on the lower frequency waveform components. Use of a high-pass filter has similar effect (see V.2.1). We show in that later section that the QNBs are in fact a superposition of (at least) four least-damped quasinormal modes, with through .
Our scalar self-force results suggest that comparable QNBs may appear in the gravitational waveform, which would provide a gauge-invariant signal of the effect. These bursts are faint and might be fainter still in the gravitational case where will be the first mode excited. On the other hand, we have not yet conducted a thorough parameter survey to find where the excitation is maximized. Furthermore, it is entirely possible that even faint QNBs might be detected and measured using template matching. QNBs in EMRIs provide the exciting possibility of measuring black hole properties by repeatedly “tickling the dragon’s tail,” as opposed to settling for the single final excitation of quasinormal modes seen in LIGO/Virgo mergers. Finally, QNBs might reveal the presence of EMRIs in systems with heavy primaries, where the usual, low-frequency parts of the signal are difficult to detect but the periodic, higher-frequency QNBs lie in LISA’s area of best sensitivity.111During discussion at a recent (May 2019) LISA Waveform Working Group meeting we were made aware that this signal has been observed in a gravitational self-force code; M. van de Meent, private communication.
The layout of this paper is as follows. Section II covers the formalism, with Sec. II.1 discussing the general nature of the scalar self-force, Sec. II.2 reviewing the generic geodesic motion problem and setting our notation, and Sec. II.3 outlining the Fourier-harmonic decomposition of the scalar Teukolsky equation. Section III gives key details about the techniques we developed and adapted to efficiently handle each phase of the generic Kerr SSF problem including spectral solution of the orbital motion (Sec. III.1), optimizations of the MST method for solving the homogeneous wave equation (Sec. III.2), and efficient spectral source integration for solutions to the inhomogeneous wave equation (Sec. III.3). Section IV discusses the regularization procedures and computation of all four components of the scalar self-force. We also discuss there the split between conservative and dissipative parts of the self-force on generic, nonresonant Kerr orbits. In Sec. V we present our results, including the QNB-in-waveform discovery highlighted above, and discuss various implications. For this paper we use units such that , use metric signature , and the sign conventions of Misner, Thorne, and Wheeler Misner et al. (1973).
II Review of the Formalism for the Scalar Self-Force Problem
II.1 The scalar self-force
The SSF model we consider assumes a point particle of mass and scalar-charge in bound motion about a Kerr black hole of mass and spin parameter . Perturbations in the gravitational field and the associated GSF are neglected. Instead the particle motion generates a scalar field , whose local behavior acts back on the scalar charge to produce the SSF. Absent the SSF, the motion of the particle is a geodesic in the Kerr spacetime. The scalar field satisfies the curved-space Klein-Gordon equation (i.e., the spin-0 Teukolsky equation Teukolsky (1973))
[TABLE]
where is the scalar (point) charge density and is the (inverse) Kerr metric. Causal boundary conditions are selected, making the resulting solution the retarded field . The particle’s timelike worldline is and its four-velocity is , where is proper time. Formally, the SSF will make the motion nongeodesic and the SSF will in principle depend upon the entire past inspiral. However, if is sufficiently small and the SSF weak, the inspiral will be adiabatic, mimicking the GSF case with EMRIs. Making this assumption here, we take the past worldline as some (arbitrary) bound geodesic and calculate the SSF along that fixed motion, the result being the (approximate) geodesic self-force. While not a topic of this paper, once the geodesic SSF is obtained in this way, it might be used in an osculating elements calculation to determine the inspiral as is done with the GSF Pound and Poisson (2008); Gair et al. (2011); Warburton et al. (2012); Osburn et al. (2016); Warburton et al. (2017). The multiple periodicity of the background geodesic makes it possible to solve the field equation in the FD, which we do in this paper.
The retarded field diverges at the point charge, necessitating a regularization procedure Quinn (2000) to compute the SSF. Detweiler and Whiting Detweiler and Whiting (2003) gave one particular separation of the retarded field into regular and singular pieces , where satisfies the same inhomogeneous wave equation (1) as but with (different) boundary conditions such that not only satisfies the source-free wave equation but is the part of the field solely responsible for the SSF
[TABLE]
Because the SSF is not orthogonal to the four-velocity Quinn (2000), all four components of must be determined. Substitution of or in (2) in place of produces corresponding forces, and , both of which are divergent on the particle worldline. Thus even though one might write
[TABLE]
the expression is not immediately useful given the divergences. Instead, one practical procedure is mode-sum regularization Barack and Ori (2000, 2003), wherein the retarded, singular, and regular fields (as well as their associated forces) are decomposed into angular harmonics (typically using scalar spherical harmonics for everything including components of vectors). The individual mode amplitudes are finite and if the subtraction in (3) is taken before summing (over ), the finite SSF is recovered
[TABLE]
The singular part, , can be obtained by local analytic expansion in an -dependent series with -independent regularization parameters. The lower-order parameters are known Barack and Ori (2003). The structure of higher-order terms is also understood Detweiler et al. (2003) and analytic expressions have been given for certain restricted motions on Schwarzschild Heffernan et al. (2012, 2018) and Kerr Heffernan et al. (2014) backgrounds. We fit numerically Detweiler et al. (2003) for higher-order parameters in our more general application (Sec. IV).
With an assumed fixed background geodesic, the SSF can be decomposed Warburton and Barack (2011) into dissipative () and conservative () pieces
[TABLE]
though assembling these pieces of the SSF from symmetries of the retarded field is more subtle for generic, nonresonant orbits on Kerr. The dissipative part is responsible for the secular orbital decay producing the inspiral, while serves to perturb the orbital parameters. The dissipative self-force does not require regularization, as it is derived from the difference between retarded and advanced fields, . The regularization procedure is necessary to determine . This decomposition is further discussed in Sec. IV.2.
II.2 Bound geodesic motion in Kerr spacetime
We briefly review the generic geodesic motion problem in Kerr spacetime to set our notation for use later in the paper. In Boyer-Lindquist coordinates a Kerr black hole of mass and spin has the line element
[TABLE]
where
[TABLE]
We define the conserved specific energy and -component of the specific angular momentum
[TABLE]
using the Killing vectors and , and define the (scaled) Carter constant
[TABLE]
associated with the Killing tensor Walker and Penrose (1970).
The geodesic equations are then Carter (1968); Misner et al. (1973); Drasco et al. (2005)
[TABLE]
where the separate -dependent and -dependent functions appearing in the last two equations are
[TABLE]
Instead of parametrizing the orbit by , , and , alternative constants of the motion can directly characterize the size, shape, and orientation of the orbit. The potential in (II.2) is a quartic polynomial in and has four roots, which we denote by the following ordering: . For a bound, stable orbit the two largest roots are finite and give the limits of radial motion. Analogous to Keplerian orbits, these extrema serve to define an eccentricity and a semi-latus rectum
[TABLE]
Having fixed in this way, it is useful to follow Drasco et al. (2005) in defining the dimensionless quantities and from the final two roots for use in later expressions
[TABLE]
An inclination angle is then defined from and Hughes (2001)
[TABLE]
It is straightforward to choose , , and as orbital parameters and then solve for , , and Schmidt (2002); Drasco and Hughes (2004).
A key additional reparameterization is to switch from to Mino time Mino (2003)
[TABLE]
which allows the and motions to separate
[TABLE]
These equations yield solutions that are functions of the new curve parameter , e.g., , with confusion over the slight abuse of notation avoided by explicit reference to the new curve parameter. The subscript continues to mean “on the worldline.” With this in mind, further reparameterizations are made by introducing Darwin-like Darwin (1961) angular coordinates and Schmidt (2002); Drasco and Hughes (2004)
[TABLE]
where . In the last equation, are roots of , with ordering of roots taken to be and associated with the turning points.
Equations (25) and (26) may be combined to find differential equations relating and to , or vice versa with functions and satisfying
[TABLE]
The definitions of and in (26) are made to improve the behavior of the differential equations at what would otherwise be turning points for and . The solutions for and can be expressed as integrals
[TABLE]
where and are integration constants, with providing initial conditions for orbits that do not simultaneously pass through and . The effect of choosing a nonzero value for , for example, is demonstrated in Fig. 4. The integrals in Eqs. (30) and (31) may be reexpressed in terms of elliptic integrals Drasco and Hughes (2004); Fujita and Hikida (2009) and thereby regarded as solved. We adopt an alternate approach in this paper, based on results in Hopper et al. (2015) and the observation that the integrands in Eqs. (30) and (31) are smooth and periodic functions. This allows functions like in Eq. (II.2) to be represented by exponentially convergent Fourier series that can be accurately truncated at some
[TABLE]
Term-by-term integration of (30) then gives
[TABLE]
with a similar expression for (31). The Fourier series coefficients () are ostensibly derived themselves from integrals, but it proves possible in a numerical calculation to replace the Fourier series representation with the discrete Fourier transform (DFT). The coefficients in the DFT are then rapidly and accurately obtained using the fast Fourier transform (FFT) algorithm. Section III.1 details this new application of spectral integration to Kerr orbits; reference Hopper et al. (2015) demonstrates the application to integrating Schwarzschild geodesics.
The periods of motion in and measured in Mino time are
[TABLE]
and the corresponding frequencies with respect to Mino time are
[TABLE]
Eqs. (16) and (15) can be reexpressed in terms of Mino time derivatives and the evolution of and in terms of have the following formal dependence
[TABLE]
with and constants. In these expressions the average rates of accumulation of and in are, respectively
[TABLE]
while and are oscillatory functions with period and and are oscillatory functions with period . These oscillatory functions are described by similar integrals, and we obtain their numerical solution via spectral integration in like fashion to Eq. (33) (see Sec. III.1). The average motion of and , along with the Mino time frequencies, then provide the fundamental (coordinate time) frequencies
[TABLE]
The motion of the particle can then be described by a discrete frequency spectrum
[TABLE]
with , , and being integers.
II.3 Scalar wave equation
The charge density , which acts as the source of the wave equation (1), is that of a point charge following the timelike orbital motion
[TABLE]
where and and are given by Eqs. (19) and (20), respectively. The wave equation is equivalent to the TD spin-0 Teukolsky equation Teukolsky (1973)
[TABLE]
II.3.1 Separation of variables
Equation (II.3) is amenable to solution via separation of variables Brill et al. (1972); Teukolsky (1973)
[TABLE]
Here is the Teukolsky radial function, is the spheroidal Legendre function with and multipole indices and spheroidicity (hence the subscripts). In the above equation and henceforth, the following condensed notion is introduced to represent the sum over modes
[TABLE]
Following Warburton and Barack Warburton and Barack (2011), we use for the spheroidal harmonic index and reserve for the spherical harmonic index used in the mode-sum regularization. The FD decomposition in (44) assumes bound motion, with a resulting discrete frequency spectrum that allows the field to be represented by a multiple Fourier series.
We follow Hughes (2000a); Sasaki and Nakamura (1982); Brill et al. (1972) in connecting the Teukolsky function, , to a new radial function,
[TABLE]
(Warburton and Barack Warburton and Barack (2010, 2011); Warburton (2015) make a different transformation.) Both and are used in what follows (see especially Sec. III). Inserting Eqs. (44) and (46) into Eq. (II.3), we arrive at two ordinary differential equations for and
[TABLE]
where is the angular eigenvalue, or separation constant, and is the tortoise coordinate
[TABLE]
which follows from integrating
[TABLE]
Here are the outer and inner horizon radii (roots of ). Our definition of agrees with e.g., Drasco et al. (2005); Sasaki and Tagoshi (2003) but differs from Warburton and Barack (2010, 2011); Warburton (2015). The radial potential is
[TABLE]
and gives the radial behavior of the source in the FD
[TABLE]
II.3.2 Radial solutions and time domain reconstruction
General solution of Eq. (48) requires two independent homogeneous solutions, and , that hold throughout the region and have respective asymptotic dependence
[TABLE]
Here is the wave number at the horizon, with denoting the angular velocity of the event horizon. The solution is “outgoing,” while the solution is “downgoing.” These two can be combined to construct the causal Green function for the radial equation (48), associated ultimately with the retarded solution in the TD. The solution of the inhomogeneous Eq. (48) is then found to be
[TABLE]
where
[TABLE]
is the (constant) Wronskian.
An attempt to use from (55) with (46) in (44) to make a (time domain) Fourier reconstruction of the field at points within the libration region is fraught with difficulty due to Gibbs oscillations caused by the delta function source. In this region, at points away from the worldline, the convergence in and is slow, while derivatives (needed for the SSF) may not even converge at the particle. The usual path around this problem, at least in a background spacetime with spherical symmetry, is the method of extended homogeneous solutions (EHS) Barack et al. (2008). In that case the four-dimensional wave equation separates into two-dimensional wave equations in and for each spherical harmonic order , . Extended homogeneous solutions are found mode by mode, which are finite at the particle as needed for mode-sum regularization. Unfortunately, in Kerr spacetime the angular decomposition in spheroidal harmonics is inseparably linked to the transformation into the FD. As Warburton and Barack Warburton and Barack (2011) have shown however, it is still possible to define functions on the spherical harmonic basis that can be extended to the particle location and are finite there.
This procedure begins with determining normalization coefficients, , which are found by evaluating at the limits of the radial libration region
[TABLE]
and which are used to define the properly normalized extended homogeneous radial modes in the FD
[TABLE]
These solutions in turn may be used in (44) to define extended solutions in the full time and space domain
[TABLE]
from which the retarded solution to (II.3), at least off the worldline, can be given as
[TABLE]
While the functions (61) converge exponentially in and and their use eliminates the Gibbs behavior near the particle in the libration region, the full reconstruction (62) is not of immediate use in calculating the SSF. The approach taken by Warburton and Barack relies upon using the representation Hughes (2000b) of spheroidal angular harmonics in terms of spherical harmonics
[TABLE]
While the spheroidal harmonics of order couple to an infinite number of spherical harmonics, the coupling coefficients decay in size rapidly as the difference in orders grows Warburton and Barack (2010), the rate dependent upon the spheroidicity . In a numerical calculation, the number of spherical harmonics needed for a given accuracy can be determined. The coupling coefficients are determined by a three-term recurrence relation that results from inserting (63) into (47).
Substituting (63) into (61), the five-fold summation may be reordered to leave and for last. This allows the extended functions to be defined,
[TABLE]
where in a practical numerical calculation the sum over will be finite in number, as will the sums over and given their exponential convergence. The remaining sums allow to be recovered
[TABLE]
The functions are not modes in the fullest sense, since there are no wave equations in and that they satisfy. However, they do derive from linear combinations of extended (homogeneous) radial modes in the FD, they provide a decomposition of , and they are finite at the location of the particle. These properties are all that is essential for employing mode-sum regularization, as shown by Warburton and Barack (2010, 2011); Warburton (2015) and as outlined in Sec. IV. Our generalization here to eccentric inclined orbits introduces no qualitatively new element in the Kerr SSF regularization, only a further dimension in the mode calculations.
The homogeneous solutions are often obtained by numerical integration Warburton and Barack (2010, 2011); Warburton (2015) of (48). In this work, however, we use a Mathematica code employing the MST method (Sec. III.2) to derive the mode functions. The resulting code is very accurate but slow. We are concurrently developing and testing a faster, complementary C-code based on numerical integration of (48).
III Analytic and Numerical Solution Techniques
This section describes some of the analytic and numerical techniques we use to solve the SSF problem in the generic Kerr case. It provides some details on spectral solution of the orbit equations, on efficient use of the MST method to obtain certain mode functions, and especially on spectral integration of source terms in the scalar case. The computational roadmap is as follows:
- (i)
After specifying the parameters , , , and , we solve for the geodesic motion on Kerr using spectral integration techniques (Sec. III.1). From the geodesic, we determine the fundamental frequencies of the orbit, , , and . 2. (ii)
We calculate the radial and polar mode functions for each frequency and multipole. The polar mode functions (spheroidal harmonics) are constructed using Eq. (63). We calculate the homogeneous radial mode functions, , using the MST function expansion formalism, with Sec. III.2 serving primarily to discuss an efficient approach to finding the near-horizon modes. 3. (iii)
Finally, we discuss in Sec. III.3 means to evaluate the normalization constants , which determine the scalar field via the EHS method, using spectral source integration techniques. In the scalar case, it proves possible to decompose the source integration (59) into products of one-dimensional integrals.
III.1 Spectral integration of the geodesic equations
As an alternative to using initial value integration, or to using special functions Drasco and Hughes (2004); Fujita and Hikida (2009), we employ a spectral (Fourier) integration technique to find the Kerr geodesics numerically. Spectral integration of the orbital motion problem in Schwarzschild spacetime was previously carefully laid out in Hopper et al. (2015). In this subsection we generalize that approach to generic bound geodesics in Kerr spacetime.
We first consider the dependence of on Darwin angles and . The integration for is given as an example, but the same approach applies to . As discussed in Sec. II.2, the function can be written as a cosine series because it is smooth, even, and periodic
[TABLE]
Because is , (66) converges exponentially with the number of harmonics, and for a given accuracy may be truncated at some as in (32).
Fourier series coefficients like are derived from integrals, so computing many of these by, for example, adaptive stepsize integration is no improvement over simply integrating (II.2) itself. Instead, an efficient alternative is to use the discrete Fourier transform (DFT). To do so, we use (32) to sample at evenly spaced points . The sampled values are the DFT of Fourier coefficients, . Up to a normalization factor, the DFT coefficients (with no tilde) converge exponentially to the Fourier series coefficients (with tilde) as the sample number increases. Since is even, we discretely sample the arc of half the radial motion and represent the function with the type-I discrete cosine transform (DCT-I)
[TABLE]
[TABLE]
[TABLE]
The DFT (or in this case DCT) may be computed numerically using a fast Fourier transform (FFT) algorithm, efficiently finding all of the Fourier coefficients . The angular sampling of is made over equally spaced points. The required radial and angular sample numbers are independent and subject only to desired numerical accuracy goals.
Returning to the radial motion example, once is adequately represented, then is found by substituting (III.1) into (30) and integrating term-by-term
[TABLE]
an expression that can be evaluated at any . The same is done for . The Mino time periods, and , are related to the leading Fourier coefficients
[TABLE]
Taken all together, these solutions for and end up accurately relating motion in and with . This approach models that found in Sec. II of Hopper et al. (2015).
We proceed next to find the motion in and . With (15) and (16) reexpressed in terms of Mino time, the periodicity of (17)-(20), and the ability to express those functions in terms of , suggests a Mino-time Fourier decomposition of , , , and Drasco and Hughes (2004)
[TABLE]
where in keeping with the left-hand sides being real functions the coefficients will satisfy crossing relations (e.g., ). As before, the series might be truncated (here with some upper and lower bounds on and ). The Fourier coefficients are found from integrals over ; for example
[TABLE]
with similar integrals for , , and .
If we introduced sufficiently fine, evenly spaced divisions of the respective periods in , each of the Fourier coefficient integrals, like (76), could be accurately replaced with a finite sum. Unfortunately, the functions being integrated depend on or (e.g., above), which are known from the previous analysis as functions sampled on evenly spaced grids in or . Rather than resample them to an evenly spaced grid in , we instead convert the integrals and integrate over or . For example
[TABLE]
with similar expressions for and . Despite the transformations, all of these integrands are still periodic functions of (now) or . As was shown in Hopper et al. (2015) (Sec. III.B.3), smooth reparameterizations of this sort still allow exponentially convergent approximations to be made by replacing an integral with a finite sum on an evenly spaced grid in the new coordinate (either or )
[TABLE]
Similar expressions again hold for and . Because (80) and (82) are not evaluated on an evenly spaced, periodic grid in , they do not represent DFT sums (the argument of the exponential is nonlinear in or ). Accordingly, the coefficients cannot be computed with the FFT algorithm, but instead are evaluated directly, which is an process.
Once the Fourier coefficients are known, the average accumulation rates, and , are found from the leading coefficients
[TABLE]
The remaining parts that determine the advance of and in (36) and (37), the periodic functions and , may be expressed as functions of by integrating (72)-(75) term-by-term
[TABLE]
Here and are assumed to be even and the restricted range of the sums reflects use of the crossing relations.
III.2 Analytic mode functions from MST formalism
The MST formalism Mano et al. (1996a, b) ultimately provides radial mode function solutions to (48) subject to the boundary conditions (53) and (54). The formalism more traditionally yields the radial Teukolsky functions (in our case spin weight equal zero), from which follow . A comprehensive review of the MST formalism is given in Sasaki and Tagoshi (2003). Our presentation here primarily focuses on efficient calculation of one set of these solutions. The calculation first starts by determining the separation constant . We make use of the Black Hole Perturbation Toolkit’s BHP Mathematica package SpinWeightedSpheroidalHarmonics to evaluate .
The Teukolsky functions and are the solutions to the radial Teukolsky equation with boundary conditions
[TABLE]
that correspond to the conditions (54) and (53), respectively, on . Here and are asymptotic amplitudes. By introducing the renormalized angular momentum and rescaling the radial coordinate in two convenient ways
[TABLE]
the functions and are expressed as series of hypergeometric functions,
[TABLE]
where here we have adopted (spin weight of the scalar case), which we maintain henceforth in this paper. Other parameters are
[TABLE]
In the expressions above, is the Gauss hypergeometric function and is the irregular confluent hypergeometric function.
The series coefficients are the minimal solution to a three-term recurrence relation that allows the series to converge once is determined. The second set of coefficients are completely determined by via
[TABLE]
making the “up” series convergent also. Here is the Pochhammer symbol. For , we calculate and using Mathematica’s built-in functions Hypergeometric2F1 and HypergeometricU, respectively. For , we construct both types of hypergeometric functions using their respective three-term recursion relations (provided in Sasaki and Tagoshi (2003)).
The eigenvalue is often determined by solving for the root of a complex equation with coefficients that are built from continued fractions Sasaki and Tagoshi (2003). An alternative method, employed in this paper, relates to the eigenvalue of the monodromy matrix defined for the irregular singular point of the Teukolsky equation at Castro et al. (2013a, b). Then is determined numerically by calculating Stokes multipliers Rodriguez ; Castro et al. (2013b).
An accuracy goal in determining the radial functions is met in part by terminating the hypergeometric series at a sufficiently large value of (where is not necessarily the same for both series). The MST technique provides precise, semianalytic solutions, but it can be computationally expensive, especially when programmed in Mathematica. As the frequency increases, the hypergeometric series expansions must range over an increasing number of terms to meet a pre-defined accuracy goal. Computational costs are exacerbated by roundoff errors from near cancellations in the sums. Roundoff errors are circumvented by making internal Mathematica calculations at working precisions significantly higher than desired accuracy in final results.
We found empirically that, for the radial positions considered in this work, the series of confluent hypergeometric functions converges more rapidly than the series of Gauss hypergeometric functions (used in the “in” solution). Further study showed that computational costs can be mitigated on the horizon side in calculating by using an alternative expression given in the MST literature (see Eq. (166) in Sasaki and Tagoshi (2003))
[TABLE]
where is expressed as a series of regular confluent hypergeometric functions ,
[TABLE]
Here is a new set of series coefficients (given below) and is a (phase) factor that involves summing over the prior series coefficients and . The exact form of in our case is given by
[TABLE]
The new series coefficients , , and can be expressed in terms of the prior coefficients and by
[TABLE]
The review by Sasaki and Tagoshi Sasaki and Tagoshi (2003) discusses (95) as a complement of (91) that provides convergent coverage of the entire domain , but does not mention its computational efficiency. Rapid convergence was our focus in comparing these expressions and settling on use of (95). While writing this paper, we sought other MST users’ experiences with the potential practical virtues of (95) and (96). Casals Buss and Casals (2018) and Wardell (private communication) were aware of the benefits of (95) and make use of it in their work, though have not previously discussed this particular issue in detail. Use of both (95) and (91) are described by Throwe Throwe (2010), with his observation that both formulae have their own regions in which they are numerically more suitable. Elsewhere BHP Eq. (91) is exclusively used.
A side benefit in our approach is that the series of regular confluent hypergeometric functions given in (96) converges with similar rapidness as the series of irregular confluent hypergeometric functions given in (92). Thus the same value of can be used to truncate both series.
While use of has benefits, it is not straightforward to construct the underlying functions numerically. The functions satisfy a three-term recurrence relation Sasaki and Tagoshi (2003) but evaluating the functions by stepping through the recurrence formula is numerically unstable in the increasing- direction. There are several ways to circumvent this problem: increase the code’s internal precision, calculate directly using Mathematica’s built-in function Hypergeometric1F1, or translate the three-term recurrence relation into a continued fraction, which does not suffer from the same cancellation errors in the increasing- direction. Alternatively, since the recurrence relation does not suffer the same instability when moving down in , one can begin the summation of (96) at and evaluate the terms as decreases down to . The value of is conveniently determined by evaluating first. A mixture of these strategies is employed to maximize computational efficiency. Ultimately the improved convergence of (95) and (96), compared to (91), offsets the computational cost of summing two series instead of one.
Using these expressions for and , we can construct the unit-normalized functions by comparing (88) and (89) with (46), (53), and (54)
[TABLE]
where . The asymptotic amplitudes can be found by expanding both solutions near the horizon and at large
[TABLE]
with
[TABLE]
III.3 Optimized source integration
We consider next the optimized calculation of the normalization coefficients defined in (59). That reduction begins with a review of the derivation of the FD source function , exploiting the orthogonality of the harmonics in and , and the spheroidal Legendre functions found in (52). Integrating the product of (52) and over azimuth angle and using the delta function in , we find
[TABLE]
We next remove the linear phase factor , which makes the remaining expression
[TABLE]
biperiodic with fundamental frequencies and , since up to an irrelevant constant.
We next reduce (106) to a single sum over by using orthogonality of the factor . To do so, we convert to Mino time Fourier series, with , using results in Drasco and Hughes (2004)
[TABLE]
where the function is
[TABLE]
which can be thought of as a function of and . The final step in deriving is multiplying the above expression by and integrating over
[TABLE]
With the FD source function in hand, we may calculate the normalization constants by substituting (109) into (59)
[TABLE]
The order of integration is exchanged, allowing the integral to be evaluated first
[TABLE]
where , implicitly a function of and , is given by
[TABLE]
The double integral in (111) may be computed directly using adaptive-step-size integration Drasco and Hughes (2006). We refer to this method henceforth as the “2D-integral” approach, which can be shown to deliver numerical results that converge algebraically (i.e., as a power law). Given the number of modes in the Kerr generic-orbit problem, this is a computationally expensive method that compelled us to search for more efficient alternatives in evaluating (111).
A first alternative is to exploit the integrand’s smoothness and bi-periodicity to make a discrete, evenly spaced sampling in two dimensions that is analogous to the approach we took with the orbit equations. Just as in that case, where an equally spaced sum over samples of a smooth periodic integrand converged exponentially, we find “spectral” convergence in the two-dimensional integral as well. Using the discrete sampling locations of (79) and (81), we calculate
[TABLE]
where we have changed the integration variables from and to and , adopted and , and let the arguments of reflect the discrete sampling. The integration approach in (113) is referred to here as the “2D-SSI” method (i.e., the two-dimensional generalization of the SSI technique Hopper et al. (2015)). Figure 3 demonstrates the increased efficiency of the 2D-SSI method compared to the 2D-integral scheme. The 2D-SSI method has been independently adopted by van de Meent van de Meent (2018) in his GSF FD calculations on inclined eccentric orbits in Kerr spacetime. We also understand that the code used in Drasco and Hughes (2006) has been upgraded to use the 2D-SSI method (Hughes, private communication).
The explicit dependence on and found in (112) allows for further optimization. Because can be written in the following form
[TABLE]
[TABLE]
[TABLE]
the double integral in (111) can be calculated from products of one-dimensional integrals
[TABLE]
If we use these equations and just compute the integrals (118)-(121) with a straightforward adaptive integrator, we get an algebraically convergent method that we refer to as the “1D-integral” approach. Despite its algebraic convergence, it is much faster at any required level of accuracy than the 2D-integral approach, by as much as two orders of magnitude at conventional double precision (as seen in Fig. 3). At that accuracy level it is also faster than 2D-SSI, though the faster convergence rate of 2D-SSI would ultimately win at higher accuracies.
Finally, the 1D integrals are just as amenable to the SSI method as the double integral and it is possible to make an exponentially convergent discrete representation for (118)-(121)
[TABLE]
When (122)-(125) are used to evaluate (117), we refer to it as the “1D-SSI” method. Figure 3 shows that the 1D-SSI method is the most efficient and most rapidly convergent technique. Switching to 2D-SSI from 2D adaptive-step integration is nearly two orders of magnitude faster at double precision accuracies. Switching from 2D-SSI to 1D-SSI yields another factor of 30.
The 1D-SSI method is possible because the two-dimensional source integrations decompose as shown in (117) into products of 1D integrals. Unfortunately a similar decomposition does not occur in any obvious way for gravitational perturbations in Kerr spacetime due to leading factors of . For small spins or large radial separations, the factor might be expanded using a binomial series with a modest amount of terms, providing an approximately separable source. It is also conceivable that a transformation might exist that would bring the source into a separable form. The benefits of the 1D-SSI method seen in the scalar case are compelling enough to justify a more thorough investigation of the gravitational Teukolsky source integration problem.
IV Generic Orbit SSF Regularization
IV.1 Mode-sum regularization review
Section II.3 provides a roadmap for calculating the retarded field, , including its decomposition in a spherical harmonic basis, and Sec. II.1 discusses using the gradient of that field and the singular field (with the vector components also expanded in the same basis) to yield the mode-sum regularized self-force
[TABLE]
This equation differs from (4) in making clear that individual -mode self-force components may differ in value in the limit as depending upon the direction of approach in . This notation aligns with that used in the EHS discussion (i.e., Eq. (62)) of mode functions. Using the spherical harmonic decomposition (65) of the retarded field, the -modes of three of the force components are
[TABLE]
The component222In the GSF case it is sufficient to regularize just three of the four force components because the final component is fixed by . In the SSF case the force has a tangential component along , leading to variation in mass Quinn and Wald (1997); Quinn (2000); Burko et al. (2002); Drasco et al. (2005); Poisson et al. (2011) and requiring calculation and regularization of . of the self-force is broken down into -modes, , only after the derivative is reprojected onto the basis.
To effect this change, we use the clever window function devised by Warburton Warburton (2015) (his Eq. 50). When multiplied with the field, affects neither the value of the field as nor its first derivative, yet produces a finite coupling between and (up to) four spherical harmonics . This allows the -modes of the -component to be reexpressed as
[TABLE]
where the are defined in terms of using the following condensed notation
[TABLE]
The coefficients , and a more detailed discussion of deriving Eqs. (130) and (131), are provided in Appendix A. Our expressions (130) and (131) are similar to ones found in Warburton (2015) with the exception of minor corrections.
To calculate from Eq. (3) we require an expansion of in terms of regularization parameters Barack and Ori (2000, 2003); Detweiler et al. (2003)
[TABLE]
where and the parameters , , and are all independent of . For each , the higher-order regularization terms (with coefficients ) have the property that the -dependent terms sum to zero Detweiler et al. (2003):
[TABLE]
As a consequence only the first two regularization parameters, and , are necessary to assure a convergent result and the regularized self-force can be calculated from just
[TABLE]
where we have defined for later convenience. While the sum in (134) gives a finite result, the higher-order terms drop off at a rate of . When the sum is approximated by being truncated at , there is a residual error that scales as . Due to computational costs, it is typically beneficial to truncate the SSF calculation at , which means that relying only upon the regularization parameters and will determine to just one or two digits of accuracy.
Including the higher-order parameters can improve the rate of convergence of the partial sums of Eq. (126), which are now written as
[TABLE]
Here there is a two-fold truncation, with determining the number of modes we calculate in the retarded field, , and setting the limit in the number of available higher-order regularization parameters. Eq. (135) converges at a rate of and therefore the SSF has an error that scales as . Unfortunately, only and are known analytically for generic orbits in Kerr Barack and Ori (2003) (although, terms up to are known for equatorial orbits in Kerr Heffernan et al. (2014)).
We overcome the lack of analytically known higher-order regularization parameters by fitting Detweiler et al. (2003) the high- contributions to the SSF to the assumed form in (132), similar to the means discussed in Sec. IVC of Warburton and Barack Warburton and Barack (2010). At high , the self-force contributions are primarily determined by the missing regularization parameters
[TABLE]
The number of regularization parameters that can be determined is limited by the precision of and . We take the last self-force -mode contributions, , and fit these values to regularization parameters by applying a least squares algorithm to Eq. (136). The value of is varied and a weighted average is taken as described in Warburton and Barack (2010). We also vary and use the standard deviation of the results to estimate the error produced by this fitting scheme. However, we do not use Eq. (47) in Warburton and Barack (2010), but instead reapply the fitted regularization parameters using Eq. (135) to improve the convergence of our SSF results. The estimated errors are also propagated to determine the accuracy of the SSF results. Errors due to fitting typically dominate over the error from terminating the -mode summation. The validity of these fits and their errors is further discussed in Sec. V.1, where we compare fitted conservative self-force data (for an inclined Schwarzschild orbit) to conservative self-force data that has been regularized with known higher-order parameters (for an equatorial Schwarzschild orbit).
IV.2 Conservative and dissipative self-force for generic orbits
As mentioned in Sec. II.1, the self-force can be decomposed into dissipative and conservative components, and , which have different physical effects on the orbital evolution Mino (2003); Diaz-Rivera et al. (2004); Barack (2009). Just as we defined the retarded force , we similarly define the advanced force from the advanced scalar field solution, along with its -mode contributions . Using the mode-sum scheme, the dissipative and conservative components to the self-force are constructed from symmetric and antisymmetric combinations of
[TABLE]
This decomposition is also beneficial for testing the numerical convergence of the self-force results: the dissipative component does not need to be regularized and will converge exponentially, while the conservative component requires regularization and will converge algebraically as discussed in Sec. IV.1.
As is well known Mino (2003); Hinderer and Flanagan (2008), the advanced and retarded forces may both be obtained from the retarded solution, being related at reflection point pairs in the orbital motion–points where the particle passes through the same radial and polar positions but with opposite radial and polar velocities, . Explicit calculations of the conservative and dissipative components of the self-force have been made by identifying these reflection points along restricted orbits (circular equatorial; eccentric equatorial; or inclined spherical) Barack (2009); Warburton and Barack (2011); Warburton (2015); Thornburg and Wardell (2017a).
For eccentric inclined orbits, these reflection points can be conveniently identified by mapping the particle’s motion to a two-torus, as shown in Fig. 4. In this figure we cover the torus using the coordinates and , related to the position in the polar () plane by Eq. (26). (Alternatively, some authors use the two angle variables Hinderer and Flanagan (2008); van de Meent (2016, 2018) to cover the torus.) The polar motion winds and wraps in this region, either a finite number of times for a resonant orbit or an infinite number of times for a nonresonant orbit. In the later case, the motion is ergodic and the motion will eventually pass all points arbitrarily closely. All of the field and self-force information can be projected onto the domain spanned by .
As an example, consider an orbit with geometric parameters and initial position set by taking and , where is measured in units of . The path of this orbit on the two-torus from to is traced out by the blue (solid) line in Fig. 4. For any point on this curve, its reflection point is identified by reflecting through the center of the plane at and (reflections can be made across any corner of the region equally well). The result of reflecting the entire blue (solid) curve is the red (dot-dashed) curve. This can be verified using Eqs. (26)-(29). Note that the red (dot-dashed) curve can also be described by an orbit moving backwards in time from to with the same geometric parameters as the blue (solid) line, but with opposite offset: but . This is in line with Eq. (2.46) in Mino (2003).
Therefore (up to a factor of ) the advanced force can be calculated by reflecting the retarded force data on the torus. Explicitly, the retarded and advanced forces are related by
[TABLE]
where and where there is no summation over . Eq. (139) can be extended to inclined spherical, eccentric equatorial, and resonant orbits as well, though the motions within the torus are severely restricted for these special orbits.
V Results
Our results are broken down into three categories:
- (a)
Eccentric inclined orbits in Schwarzschild spacetime; 2. (b)
Highly eccentric equatorial orbits about a rapidly rotating Kerr black hole, displaying quasinormal bursts; 3. (c)
Eccentric inclined (generic) orbits in Kerr spacetime.
V.1 Schwarzschild eccentric inclined orbits
We first examine eccentric inclined orbits in the Schwarzschild limit (). These models serve as a strong validation of the SSF code, since all elements of the field and self-force calculation are required, yet they can be compared to much simpler-to-compute eccentric equatorial models (i.e., ones with vastly fewer computed modes). The one-to-one correspondence results from spherical-symmetry of Schwarzschild spacetime, where two geodesics with the same eccentricities but different inclinations are related merely by a rotation.
In spherically symmetric spacetimes, the self-force for an eccentric inclined orbit can be compared to the force that is obtained through rotational transformation of the equatorial plane self-force . The transformation is
[TABLE]
where depends on the sign of ( when ).
The four SSF components for an orbit characterized by are plotted in Fig. 5. For equatorial orbits, the self-force is a periodic function of . This periodicity continues to be seen in Fig. 5 for the , , and components in the inclined model as these self-force components “loop” back onto themselves as the particle librates from to and then back to . This periodicity is evident in examining Eqs. (140), (141), and (143).
The behavior of is different. When the orbit is rotated out of the equatorial plane, the contribution is split between the rotated self-force components and . While differs from by a trigonometric factor, the projection of onto the new inclined basis depends on the longitudinal position of the particle. This causes to also depend upon (see Eq. (142)). The small body librates at different frequencies in and , which demonstrates why the inclined force component does not form a closed loop when plotted versus .
These inclined SSF results can be compared in quantitative detail, again via Eqs. (140)-(143), to results computed from an equivalent equatorial orbit . We refer to the self-force calculated directly using an inclined orbit as , while the force computed by rotating the equatorial orbit self-force remains being denoted by . The absolute residuals from comparing these orbits are plotted in Fig. 6. We also plot the errors and for both self-force calculations. The primary source of error comes from fitting the conservative component of the self-force. In Fig. 6 we see that the residual errors between the two calculations consistently fall below the errors that are estimated by our fitting procedure. This provides additional confidence in the validity of our error estimation, which is outlined in Sec. IV.1, and makes a strong case for having summed over all the required modes and correctly computed the regularization in the inclined model.
Additionally, we can compare specific numerical values of to previously and independently computed equatorial results published in Warburton and Barack (2011), by again using Eqs. (140)-(143) to transform the equatorial plane SSF. We compare both the conservative and dissipative parts of the self-force in Table 1. The fractional errors between the independently computed conservative parts typically fall below the estimated errors in the conservative parts themselves that owe to the high- fitting procedure. The dissipative part of our inclined SSF typically agrees with the transformed dissipative part from Warburton and Barack (2011) to 6 or more decimal places.
V.2 Highly eccentric orbit about a tapidly totating Kerr black hole
and quasinormal bursts in the waveform
Thornburg and Wardell Thornburg and Wardell (2014, 2016, 2017b, 2017a) were the first to demonstrate that, for highly eccentric orbits () about rapidly rotating black holes (), interesting “wiggles” arise in the scalar self-force. They further showed that these high frequency oscillations were attributable to excitation of a quasinormal mode (QNM), the least-damped mode, produced by periastron passage of the scalar-charged small body. Thornburg and Wardell observed these excitations for a number of orbital configurations. The most pronounced excitations were present in orbits with , though weak oscillations arise for the orbit (see Fig. 16 in Thornburg and Wardell (2017a)).
Thornburg and Wardell utilize a TD code, which can be well-suited for computing highly eccentric orbits. However, TD codes involve solving partial differential equations and have potential numerical issues with initial value transients, boundary conditions, and source modeling. Our code works in the frequency domain, where the numerical problem involves solving ordinary differential equations for large numbers of Fourier-harmonic modes. In general it is easier to attain higher accuracy with a FD code. However, a countering factor is that the required number of modes and computational demand in a FD code grows exponentially at high eccentricities. Accordingly, we have so far restricted ourselves to orbits with . On the positive side, a FD code only captures periodic behavior and is not subject to initial value transients. Given the many differences between the two approaches, a comparison between results seemed desirable.
Having said that, we have not made an exact comparison. We have so far not tried to make a very time consuming calculation with to duplicate one of the results in Thornburg and Wardell (2017a). At the same time, rather than replicating the results of Thornburg and Wardell, with , we decided to calculate the SSF and fluxes for the same orbital parameters but with a higher black hole spin: . The expectation was that we might see more pronounced ringing in the orbit if the QNM damping is lessened with a higher .
We also chose to model an orbit in the equatorial plane, which substantially offsets the computational cost of high eccentricity by restricting the mode spectrum to be bi-periodic and not tri-periodic. Additionally, higher-order regularization parameters are known for equatorial orbits Heffernan et al. (2014) and we were able to circumvent the fitting schemes discussed in Sec. IV.1 in this case, improving the convergence and reducing the estimated error.
Our FD SSF results for this model are plotted in Fig. 7. The closed loops in the force components are split out into conservative part, dissipative part, and total. We see the same oscillatory features in our self-force results as Thornburg and Wardell found, with the oscillations most prominent in the and self-force components. After the point charge’s periastron approach (), the ringing in the scalar field sweeps past the small body driving oscillations in the self-force, with the oscillations then decaying as the system approaches apastron. As expected, by increasing the black hole spin, we observe a more persistent ringing compared to that seen in the Thornburg and Wardell model.
V.2.1 Quasinormal bursts in the waveform and extracting multiple
quasinormal modes
As we mentioned in the Introduction, we decided to look at the waveform in this model to see if the excitations were present in the asymptotic field. While faint, there are indeed quasinormal bursts (QNBs) visible to most observers of the waveform. The waveform itself, highlighted in the Introduction with Fig. 1, appears devoid of ringing at any of three observer angles: , , and . However, high-pass filtering or emphasizing high frequencies, by taking two time derivatives of the waveform as shown in Fig. 2, makes the bursts visible. Figure 2 shows the second derivative measured by the observer at . Similar excitation is visible to an observer at , but the QNBs are not present for an observer at position (i.e., along the polar axis). As we show below, this is consistent with the ringing being due to (prograde) axial perturbations of the field in the Kerr geometry.
Rather than emphasizing high frequencies by taking time derivatives of the signal, one can instead apply a high-pass filter to attenuate the lower frequency “background.” We construct a high-pass Butterworth filter using Mathematica’s ButterworthFilterModel, ToDiscreteTimeModel, and RecurrenceFilter. We choose the filter’s parameters by inspecting the power spectrum of the waveform.
After applying the high-pass filter and observing the presence of QNBs, we attempted to extract a complex frequency for the excitation by (1) selecting a time window during which the excitation dominates the filtered signal and (2) then performing a least-squares fit of a burst template to the filtered data, as demonstrated in Fig. 8. The data was fitted to a real function of the form using Mathematica’s FindFit. Fitted complex frequencies have negative imaginary parts, consistent with damped bursts. The data in Fig. 8 was found to be best fit by the complex frequency (in units with ; henceforth assumed in this section).
We can compare this value to the spectrum of known QNM frequencies due to scalar perturbations of Kerr spacetimes published by Berti Berti et al. (2009). The QNMs depend on and are indexed by the spheroidal harmonic mode numbers and the overtone , where refers to the least-damped or “fundamental” overtone. Assuming but without assuming a value for , we find that the extracted complex frequency above most closely matches the QNM frequency for a spin of . In other words, by assuming that this complex frequency should be represented by a QNM, the extracted frequency accurately recovers the spin of the primary black hole to three digits. This result is consistent with those presented by Thornburg Thornburg and Wardell (2014, 2016, 2017b), who found that, across several orbital configurations and spin parameters, the QNM frequencies in his self-force data were best fit by the least-damped (smallest ) QNMs.
Surprisingly perhaps, our FD numerical results actually allow us to extract additional QNMs. To do so, we obtain the residuals between the high-frequency signal and its fit in Fig. 8 and apply the high-pass filter a second time to remove a remaining background (i.e., “flat-fielding” the signal). We fit and obtain the complex frequency of a second damped oscillation. By iterating this process, we managed to extract three additional QNM excitations in the filtered waveform. These are shown in Fig. 9. The numerical values of the frequencies of all extracted QNMs are presented in Table 2 and compared to the closest published QNMs for scalar perturbations of a Kerr spacetime with .
However, we can instead try to remain agnostic to the black hole spin and mode numbers and compare the extracted frequencies to all known QNM frequencies across Berti’s densely sampled set of Kerr spacetimes. Consulting Table 2, our second extracted frequency best fits a QNM in Berti’s table with frequency for . Our third extracted frequency best fits one with for and the fourth best fits Berti’s mode for . By simply looking for the best fit to known QNMs, we obtain multiple estimates of the black hole spin parameter. Multiple parameter estimates all yield values for the black hole spin that are surprisingly close to (with approximately three digits of agreement). If QNBs can be observed in highly eccentric EMRIs, it may well be possible to get repeated snapshot determinations of the mass and spin of the primary black hole. Furthermore, while the “orbital parts” of the EMRI waveform will evolve and move through the LISA passband, the frequencies of the QNB component of the waveform will remain invariant, as these depend upon the (essentially unchanging) primary mass and spin.
By reproducing Thornburg and Wardell’s “wiggles,” we affirm that these are integral components of the SSF. The finding of related QNBs in the scalar waveform suggests the strong likelihood that QNBs exist in the gravitational waveforms of (some) EMRIs. A gauge invariant signal of this type, from repeatedly “tickling” the primary black hole, might have important observational consequences in sufficiently high signal-to-noise ratio EMRIs.
V.3 Kerr inclined orbits
V.3.1 Spherical inclined orbits
We first examine inclined orbits in the Kerr background by calculating the SSF along spherical inclined orbits. Similar to other restricted orbits, spherical inclined orbits are biperiodic in their frequency spectrum, , rather than tri-periodic like eccentric inclined orbits. Additionally, while the number of summed radial-frequency modes in Eq. (64) rapidly grows with increasing eccentricities, the number of summed polar-frequency modes is not as dramatically affected by increasing the inclination. Calculating the radial mode functions is also one of the primary computational bottlenecks of our code. Altogether these factors significantly reduce computational costs, allowing us to compute the SSF along spherical orbits at large inclinations with high precision.
These orbits serve as a code test for us, since the SSF along spherical orbits was previously investigated by Warburton Warburton (2015). We reproduced the results from Warburton (2015) for the orbit with parameters . To match the conventions of Warburton (2015), the orbit is parametrized by the component of angular momentum instead of the inclination . The self-force data produced by our code are in good agreement with those of Warburton (2015). The conservative components agree to digits and dissipative components to 7 or more digits. Comparative SSF values are provided in Table 3.
V.3.2 Eccentric inclined orbits
The truly unique capability of our code is in being able to model the SSF on generic (bound) eccentric inclined orbits. We investigate in this paper the SSF on four different orbits of this type, with their characteristic parameters specified in Table 4. We refer to these orbits by their reference names: ‘base’, ‘large ’, ‘large ’, and ‘large ’. We use the orbit as our fiducial case and then vary either the orbital eccentricity, the orbital inclination, or the black hole spin to get a sense of how the self-force depends on these orbital and spin parameters. This also provides tests of our code’s ability to probe more challenging regions of parameter space. The ‘large ’ orbit is also used in Fig. 10 to demonstrate improved convergence of the mode-sum through incorporating additional numerically extracted regularization parameters.
While in restricted cases the self-force can be periodic, for generic orbits the self-force is instead biperiodic. As such, it is less practical to plot the self-force as a function of time or radial position as in Figs. 5 and 7. Instead, as long as the orbit is not resonant in and motion, we can map the self-force as contour levels on the torus spanned by the coordinates and , similar to the use of the torus in the discussion surrounding Fig. 4 of Sec. IV.2. The ergodic nature of the particle’s motion implies that the SSF is a smooth continuous field over and , with any given point eventually sampled by the motion (see also van de Meent (2018)). This representation of the SSF for the generic (nonresonant) orbits listed in Table 4 is shown in Fig. 11. (In these plots we use and as coordinates rather than angle variables as found in van de Meent (2018).)
For the orbits presented in Fig. 11, the largest variations in the scalar self-force occur in the radial direction, with the exception of the component. Consequently, despite the low eccentricities considered, , , and are most dependent on , i.e. the radial motion of the small body. We also see that the maxima and minima of each self-force component are shifted away from the turning points of the particle’s motion and the particle’s passage through the equatorial plane , as a result of conservative effects. These shifts are most easily recognized in .
Taking the “base” orbit as a fiducial result, we can also examine how the self-force changes as we vary the orbital parameters and or the spin parameter . For the ‘high ’ orbit, we increase the eccentricity from to . We see that the radial dependence of the self-force becomes further accentuated, due to the orbit’s increased eccentricity. Additionally, the maximum magnitude of the scalar self-force increases in every self-force component, most likely due to the particle’s smaller pericentric distance at the higher eccentricity.
For the high orbit, we increase the inclination from to . The dependence of the scalar self-force on the particle’s polar () motion becomes more pronounced, as the particle sweeps out a larger region above and below the equatorial plane. Additionally, the radial component of the scalar self-force shifts to become predominantly positive. A similar behavior is seen for inclined spherical orbits, where the average value of grows monotonically with inclination, as it ranges from to Warburton (2015). (Retrograde orbits are parametrized with in our code.)
For the high orbit, we increase the black hole spin parameter from to . We observe a stronger dependence of the scalar self-force on the polar position of the particle when is increased. Also, the radial component of the SSF becomes attractive () along the entire orbit in this case. This is consistent with previous work on circular equatorial orbits, where decreases with increasing Warburton and Barack (2010).
V.3.3 Flux balance
As a final self-consistency check, we analyze the balance between the asymptotic fluxes with the local dissipative self-force effects Mino (2003, 2005a, 2005b); Sago et al. (2006); Warburton and Barack (2011). The average work done on the particle by the SSF should be balanced by the rate of radiative energy loss. Likewise there should be a balance between the local torque on the particle due to the SSF and the angular momentum radiated away by the scalar field. The average local work and torque are given, respectively, by
[TABLE]
In practice, periodicity (or bi-periodicity) can be leveraged to compute Eqs. (144) and (145) with finite integrals over time or finite integrals over the two-torus. Note that only the dissipative component of the self-force contributes because both and are time-antisymmetric. Therefore the conservative pieces cancel when averaging.
The asymptotic energy and angular momentum fluxes can be calculated by analyzing the scalar field at and
[TABLE]
where , , an overdot represents a time derivative, and denotes a time () average. Also recall that . The flux balance formulas then take the form
[TABLE]
The fluxes and self-force are calculated independently from one another. Consequently, comparing our scalar self-force results with flux calculations provides a self-consistency check for our code. Flux balance comparisons are included in Table 5.
VI Summary
We considered a point scalar charge following generic bound geodesics in Kerr spacetime and have calculated the scalar self-force acting on it, as a model for the gravitational self-force problem. A Mathematica code was designed to perform these calculations in the frequency domain with arbitrary numerical precision (we are currently developing C code to accomplish the same goals with increased computational efficiency). Our numerical strategy includes novel features such as fast spectral source integration techniques that reduce expensive 2D source integrals to successive 1D Fourier sums. We apply the same techniques to integrate the geodesic equations of motion. The source calculation in the scalar case is sped up by orders of magnitude and argues for a thorough investigation of whether in the gravitational case the Teukolsky equation source can be similarly arranged to allow faster numerical integration.
The accuracy of our code was validated by comparing to prior calculations of and existing results on the SSF, such as for (1) eccentric equatorial orbits, (2) inclined spherical orbits, and (3) self-comparison between inclined eccentric Schwarzschild and equatorial eccentric Schwarzschild. In all cases we verify that we calculate the scalar field and self-force with accuracy.
In the process of computing the SSF on highly eccentric () equatorial orbits about a rapidly rotating () Kerr primary, we verified a result of Thornburg and Wardell Thornburg and Wardell (2017a)–the existence of “wiggles” in the self-force due to quasinormal-mode excitation of the primary following periastron passage. Their calculations were done with a time domain code while ours were done in the frequency domain. Given substantial differences in the methods, it is heartening to see the result confirmed.
Intriguingly, we further searched for and observed quasinormal bursts (shortened to QNBs earlier in the paper) in the asymptotic waveform. (This finding became a central highlight of the paper even though we have so far only computed it on equatorial orbits.) We found that the QNBs are a superposition of not just the least-damped QNM (as Thornburg and Wardell (2017a) had already discovered) but of the least-damped QNMs as well. While our calculations are of the scalar model problem, these QNBs are likely present in the gravitational waveform as well, which would provide a gauge-invariant indicator of the effect. If so, these faint repeated bursts offer a new opportunity in high signal-to-noise ratio EMRI observations to measure rotating black hole properties. In effect, each high , high EMRI waveform would have two components: a low frequency spectrum that evolves toward higher frequency as the inspiral (chirp) proceeds and a high frequency spectrum of superposed damped modes which remain fixed in frequency (though with evolving amplitudes and phases). It awaits future work to decide how practical measurement of QNBs might be in LISA observations given expected ranges on EMRI event rates.
Our results also focused on four different inclined eccentric orbits, with parameters given in Table 4, which represents the novel elements of our method and code. We displayed in Fig. 11 how the scalar self-force changes from one of these orbits to the next, by varying inclination, eccentricity, and black hole spin. Validations of the generic orbit SSF results included examining convergence rates of the conservative self-force and checking balance between local SSF work and torque done on the small body and asymptotic energy and angular momentum fluxes.
In future work we intend to apply the generic SSF code to study resonant orbits, directly measuring the size of jumps in the waveform that can be expected as a result of transient resonances and how those jumps vary with phase of the orbit upon entering the resonance Flanagan and Hinderer (2012); van de Meent (2014). We will also likely make a thorough survey of QNB strengths, including moving beyond equatorial orbits. Part of this work may focus on strategies for processing EMRI waveforms, e.g., matching templates or co-adding waveform segments, to try to draw QNBs up out of the detector noise.
Acknowledgements.
We thank Niels Warburton, Adrian Ottewill, Barry Wardell, Maarten van de Meent, Marc Casals, and Scott Hughes for helpful discussions. This work was supported in part by NSF Grants No. PHY-1506182 and No. PHY-1806447 and by the North Carolina Space Grant Graduate Research Fellowship. C.R.E. acknowledges support from the Bahnson Fund at the University of North Carolina at Chapel Hill.
Note added in proof.–After this paper was submitted a related paper by Thornburg, Wardell, and van de Meent was submitted that demonstrated the effect in the gravitational waveform Thornburg et al. (2019). Also, the authors were made aware of an earlier paper by O’Sullivan and Hughes O’Sullivan and Hughes (2016) where the quasinormal bursts were seen in the black hole horizon shear response.
Appendix A Regularizing the -component
As mentioned in Sec. IV.1, we use a window function discovered by Warburton Warburton (2015)
[TABLE]
This window function satisfies the necessary properties and as , ensuring that is unaffected by the transformation . Additionally Warburton’s window function cleverly avoids wide bandwidth coupling thanks to the compact relationship between and
[TABLE]
The coefficients are defined as
[TABLE]
where and are given in Barack and Sago (2010) as
[TABLE]
Under these considerations, efficient calculation of follows from the replacement
[TABLE]
[TABLE]
Refactoring Eq. (155), we recover Eq. (131)
[TABLE]
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Abbott, B. P. et al. (2016 a) (The LIGO Scientific Collaboration and VIRGO Collaboration) Abbott, B. P. et al. (The LIGO Scientific Collaboration and VIRGO Collaboration) (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. D 93 , 122003 (2016 a) . · doi ↗
- 2Abbott, B. P. et al. (2016 b) (The LIGO Scientific Collaboration and VIRGO Collaboration) Abbott, B. P. et al. (The LIGO Scientific Collaboration and VIRGO Collaboration) (LIGO Scientific Collaboration and Virgo Collaboration), Physical Review Letters 116 , 241103 (2016 b) . · doi ↗
- 3Abbott, B. P. et al. (2017 a) (The LIGO Scientific Collaboration and VIRGO Collaboration) Abbott, B. P. et al. (The LIGO Scientific Collaboration and VIRGO Collaboration) (LIGO Scientific and Virgo Collaboration), Physical Review Letters 118 , 221101 (2017 a) . · doi ↗
- 4Abbott, B. P. et al. (2017 b) (The LIGO Scientific Collaboration and VIRGO Collaboration) Abbott, B. P. et al. (The LIGO Scientific Collaboration and VIRGO Collaboration) (LIGO Scientific and Virgo Collaboration), The Astrophysical Journal Letters 851 , L 35 (2017 b) , ar Xiv:1711.05578 [astro-ph.HE] . · doi ↗
- 5Abbott, B. P. et al. (2017 c) (The LIGO Scientific Collaboration and VIRGO Collaboration) Abbott, B. P. et al. (The LIGO Scientific Collaboration and VIRGO Collaboration) (LIGO Scientific Collaboration and Virgo Collaboration), Physical Review Letters 119 , 161101 (2017 c) . · doi ↗
- 6Abbott, B. P. et al. (2017 d) (The LIGO Scientific Collaboration and VIRGO Collaboration) Abbott, B. P. et al. (The LIGO Scientific Collaboration and VIRGO Collaboration), Physical Review Letters 119 , 141101 (2017 d) , ar Xiv:1709.09660 [gr-qc] . · doi ↗
- 7Abbott, B. P. et al. (2016 c) (The LIGO Scientific Collaboration and VIRGO Collaboration) Abbott, B. P. et al. (The LIGO Scientific Collaboration and VIRGO Collaboration) (LIGO Scientific and Virgo Collaborations), Physical Review Letters 116 , 221101 (2016 c) . · doi ↗
- 8Abbott, B. P. et al. (2016 d) (The LIGO Scientific Collaboration and VIRGO Collaboration) Abbott, B. P. et al. (The LIGO Scientific Collaboration and VIRGO Collaboration) (LIGO Scientific Collaboration and Virgo Collaboration), The Astrophysical Journal Letters 818 , L 22 (2016 d) , ar Xiv:1602.03846 [astro-ph.HE] . · doi ↗