Exciting black hole modes via misaligned coalescences: II. The mode content of late-time coalescence waveforms

TL;DR

This paper investigates how the excitation of quasinormal modes in black hole mergers depends on source parameters like spin and orbital misalignment, suggesting mode analysis can reveal binary properties.

Contribution

It introduces a detailed analysis of mode excitation dependence on source parameters using a generalized Ori-Thorne algorithm and perturbation theory, advancing understanding of gravitational wave signatures.

Findings

Mode mixture is determined by black hole spin and orbital misalignment angles.

Measuring multiple ringdown modes can reveal binary properties.

Results may extend to less extreme mass ratios.

Abstract

Using inspiral and plunge trajectories we construct with a generalized Ori-Thorne algorithm, we use a time-domain black hole perturbation theory code to compute the corresponding gravitational waves. The last cycles of these waveforms are a superposition of Kerr quasinormal modes. In this paper, we examine how the modes' excitations vary as a function of source parameters, such as the larger black hole's spin and the geometry of the smaller body's inspiral and plunge. We find that the mixture of quasinormal modes that characterize the final gravitational waves from a coalescence is entirely determined by the spin $a$ of the larger black hole, an angle $I$ which characterizes the misalignment of the orbital plane from the black hole's spin axis, a second angle $\theta_{\rm fin}$ which describes the location at which the small body crosses the black hole's event horizon, and the direction sgn$(\dot\theta_{\rm fin})$ of the body's final motion. If these large-mass-ratio results hold at less extreme mass ratios, then measuring multiple ringdown modes of binary black hole coalescence gravitational waves may provide important information about the source's binary properties, such as the misalignment of the orbit's angular momentum with black hole spin. This may be particularly useful for large mass binaries, for which the early inspiral waves are out of the detectors' most sensitive band.

Figures

Tables

Table 1. Table 1: Mode amplitudes that we find compared with those presented by Taracchini et al. Taracchini  et al. ( 2014 ) for equatorial inspiral and plunge [cf. their Table II and Eq. (5)]. We calculate the mode amplitudes as described in Sec.  III.3 using waveform data provided by the authors of Ref.  Taracchini  et al. ( 2014 ) . In this table, we use the notation of Ref.  Taracchini  et al. ( 2014 ) , so that positive a 𝑎 a denotes prograde ( I = 0 ∘ 𝐼 superscript 0 I=0^{\circ} ), negative a 𝑎 a denotes retrograde ( I = 180 ∘ 𝐼 superscript 180 I=180^{\circ} ). Dash entries indicate that the mode, while present, could not be reliably fitted (using the method described in Ref.  Taracchini  et al. ( 2014 ) ). Blank entries indicate that the mode is not significantly excited and thus excluded from the model. At least for modest spins ( | a | ≤ 0.5 ​ M 𝑎 0.5 𝑀 |a|\leq 0.5M ), we find fairly good agreement. Because of systematic differences in the mode extraction algorithm, we do not expect perfect agreement; and, we expect the differences to be particularly marked for a 𝑎 a large. See text in Sec.  IV.3 for details of how to convert between mode amplitudes presented in Ref.  Taracchini  et al. ( 2014 ) A ℓ ​ m ​ n subscript 𝐴 ℓ 𝑚 𝑛 A_{\ell mn} , and mode amplitudes presented in this paper, 𝒜 ℓ ​ m ​ n , 𝒜 ℓ ​ m ​ n ′ subscript 𝒜 ℓ 𝑚 𝑛 subscript superscript 𝒜 ′ ℓ 𝑚 𝑛 \mathcal{A}_{\ell mn},\mathcal{A}^{\prime}_{\ell mn} .

$a/M$		$\left|A_{2-10}/A_{210}\right|$			$\left|A_{220}/A_{320}\right|$			$\left|A_{3-20}/A_{320}\right|$
	This paper		Ref. Taracchini et al. (2014)	This paper		Ref. Taracchini et al. (2014)	This paper		Ref. Taracchini et al. (2014)
$0.9$	0.0066		0.0025	2.40		2.11	0.044
$0.5$	0.011		0.010	0.34		0.34	0.0044
$0.0$	0.071		0.069	0.000			0.010		0.010
$-0.5$	0.20		0.21	0.10		0.13	0.096		0.093
$-0.9$	0.32		0.32	0.18		—	0.29		—

Equations

$$-\left[\frac{(r^2+a^2{)}^2}{\Delta }-a^2{\sin }^2\theta \right]{\partial }_{tt}\Psi -\frac{4Mar}{\Delta }{\partial }_{t\varphi }\Psi$$

$$-2s\left[r-\frac{M(r^2-a^2)}{\Delta }+ia\cos \theta \right]{\partial }_t\Psi$$

$$+{\Delta }^{-s}{\partial }_r\left({\Delta }^{s+1}{\partial }_r\Psi \right)+\frac{1}{\sin \theta }{\partial }_{\theta }\left(\sin \theta {\partial }_{\theta }\Psi \right)+$$

$$\left[\frac{1}{{\sin }^2\theta }-\frac{a^2}{\Delta }\right]{\partial }_{\varphi \varphi }\Psi +2s\left[\frac{a(r-M)}{\Delta }+\frac{i\cos \theta }{{\sin }^2\theta }\right]{\partial }_{\varphi }\Psi$$

$$-\left(s^2{\cot }^2\theta -s\right)\Psi =-4\pi \left(r^2+a^2{\cos }^2\theta \right)T,$$

$${\psi }_4=\frac{1}{2}\frac{d^2}{d{t}^2}\left(h_{+}-i{h}_{\times }\right).$$

$$T_{\alpha \beta }=\mu \frac{u_{\alpha }{u}_{\beta }}{\Sigma u^t\sin \theta }\delta \left[r-r(t)\right]\delta \left[\theta -\theta (t)\right]\delta \left[\varphi -\varphi (t)\right],$$

$$\Psi =\sum _{\ell =2}^{\infty }\sum _{m=-\ell }^{\ell }\sum _{n=0}^{\infty }{}_s{R}_{\ell mn}^{a\sigma }(r){}_s{S}_{\ell mn}^{a\sigma }(\theta ,\varphi )e^{-i\sigma t}.$$

$${\sigma }_{\ell mn}={\omega }_{\ell mn}-i/{\tau }_{\ell mn}.$$

$$h(t)$$

$$\text{tensor}\ast []h{}^{\mathrm{N}}(t)=\sum _{\ell ,m}\text{tensor}\ast []h{}_{\ell m}^{\mathrm{N}}(t){}_{-2}{Y}_{\ell m}(\theta ,\varphi ),$$

$${}_{-2}{S}_{kmn}^{a{\sigma }_{\ell mn}}(\theta ,\varphi )=\sum _{\ell }{\mu }_{m\ell kn}^{\ast }(a{\sigma }_{\ell mn}){}_{-2}{Y}_{\ell m}(\theta ,\varphi ),$$

$$\text{tensor}\ast []h{}_{\ell m}^{\mathrm{N}}(t)=\sum _{k=k_{\min }}^{\infty }\sum _{n=0}^{\infty }\left[a_{m\ell kn}(t){\mathcal{C}}_{kmn}+a_{-m\ell kn}^{\prime }(t){\mathcal{C}}_{k-mn}^{\prime }\right],$$

$${\mathcal{C}}_{kmn}\equiv {\mathcal{A}}_{kmn}{e}^{i{\varphi }_{kmn}},{\mathcal{C}}_{kmn}^{\prime }\equiv {\mathcal{A}}_{kmn}^{\prime }{e}^{i{\varphi }_{kmn}^{\prime }},$$

$$a_{m\ell kn}(t)$$

$$a_{m\ell kn}^{\prime }(t)$$

$${}_{-2}{Y}_{\ell m}(\pi -\theta ,\varphi {)}^{\ast }=(-1{)}^{\ell }{}_{-2}{Y}_{\ell -m}(\theta ,\varphi ).$$

$$\text{tensor}\ast []h{}_{\ell m}^{\mathrm{N}}(t)=\sum _{k=k_{\min }}^{\ell +K_{\ell }}\left[a_{m\ell k0}(t){\mathcal{C}}_{km0}+a_{-m\ell k0}^{\prime }(t){\mathcal{C}}_{k-m0}^{\prime }\right].$$

$$\text{tensor}\ast []h{}_{\ell m}^{RD0}(t)=\sum _{k=k_{\min }}^{\ell +K_{\ell }}\left[a_{m\ell k0}(t){\bar{\mathcal{C}}}_{km0}+a_{-m\ell k0}^{\prime }(t){\overline{{\mathcal{C}}^{\prime }}}_{k-m0}\right].$$

$$\delta \text{tensor}\ast []h{}_{\ell m}^{\mathrm{N}}(t)=\text{tensor}\ast []h{}_{\ell m}^{\mathrm{N}}(t)-\text{tensor}\ast []h{}_{\ell m}^{RD0}(t).$$

$$I=\pi /2-[\mathrm{sgn}(L_z)]{\theta }_m,$$

$$90^{\circ }-I\le$$

$$I-90^{\circ }\le$$

$${\mathcal{C}}_{kmn}$$

$${\mathcal{C}}_{kmn}^{\prime }$$

$${\theta }_{\mathrm{fin}}^{(2)}$$

$$\mathrm{sgn}({\dot{\theta }}_{\mathrm{fin}}^{(2)})$$

$${\mathcal{C}}_{kmn}^{\prime }(a,I,{\theta }_{\mathrm{fin}},\mathrm{sgn}({\dot{\theta }}_{\mathrm{fin}}))=$$

$${\mathcal{C}}_{kmn}^{\ast }(a,I,\pi -{\theta }_{\mathrm{fin}},-\mathrm{sgn}({\dot{\theta }}_{\mathrm{fin}})).$$

$${\mathbf{u}}^{(1)}={\mathbb{R}}^{-1}{\mathbf{u}}^{(2)}\mathbb{R}.$$

$$\text{tensor}\ast []h{}^{(i)}=\sum _{\ell ,m}\text{tensor}\ast []h{}_{\ell m}^{(i)}\text{tensor}\ast {[}_{-2}]Y{}_{\ell m}({\theta }_i,{\varphi }_i)(i=1,2),$$

$$\text{tensor}\ast []h{}_{\ell m}^{(1)}=\sum _{m^{\prime }}{\mathfrak{D}}_{m^{\prime }m}^{\ell }(\mathbb{R})\text{tensor}\ast []h{}_{\ell m^{\prime }}^{(2)},$$

$${\mathcal{C}}_{kmn}^{(1)}=\sum _{m^{\prime }}{\mathfrak{D}}_{m^{\prime }m}^k(\mathbb{R}){\mathcal{C}}_{k{m}^{\prime }n}^{(2)},$$

$$\left|\frac{A_{\ell -m0}}{A_{\ell m0}}\right|$$

$$\left|\frac{A_{{\ell }^{\prime }m0}}{A_{\ell m0}}\right|$$

Full text

I Introduction

In a companion paper (Ref. Apte and Hughes (2019), hereafter Paper I), we introduced a model describing the worldline followed by a small body that orbits a Kerr black hole and is driven by gravitational-wave (GW) emission to inspiral until it encounters a dynamical instability and plunges into the hole’s horizon. Our model describes the transition from inspiral to plunge for orbits that are misaligned from the larger black hole’s equatorial plane, generalizing earlier work by Ori and Thorne Ori and Thorne (2000) which did this for equatorial orbits. Our goal now is to use this model in order to study the GWs produced by such misaligned plunges.

Our motivation for this study is to understand how a coalescing binary’s final GW cycles depend on its orbital geometry as the binary enters its final plunge and merger. Especially for binaries with total masses greater than $50\,M_{\odot}$ or so (which are amply represented in the sample that LIGO and Virgo have discovered The LIGO Scientific Collaboration and the Virgo Collaboration (2018)), the early inspiral waves are emitted at low frequencies for which ground-based GW detectors have relatively poor sensitivity. The inspiral richly encodes information about the system’s masses and spins. However, if it is not in the detector’s band, we cannot measure these waves well, and we do not benefit from this rich encoding. By contrast, for systems with $M\gtrsim 50\,M_{\odot}$, the late merger and final ringdown waves are generated at frequencies which are nearly ideal for ground-based detectors. Especially as detectors’ mid- and high-frequency noise is reduced in future upgrades Punturo et al. (2010); Hild et al. (2011); Dwyer et al. (2015); Abbott et al. (2017); Adhikari and et al. (2017); Yu et al. (2018), we can expect these final merger cycles to be measured ever more precisely. These final cycles will also be important components of the waves that are measured by the space-based detector LISA Amaro-Seoane et al. (2017), which will measure processes involving black holes of tens of thousands to tens of millions of $M_{\odot}$.

Our particular goal here is to characterize how strongly different ringdown modes are excited as a function of the inspiral and plunge geometry, as well as on the spin of the binary’s larger black hole. Past work (e.g., Refs. Berti et al. (2016, 2018); Carullo et al. (2018) for recent examples) has examined the measurement of multiple ringdown modes. Such work has typically focused on the fact that, for a Kerr black hole, each mode’s frequency and damping time depends on the final merged remnant hole’s mass and spin in a unique way. Measuring two such modes and assuming the Kerr spacetime thus suffices, at least in principle, to measure the merged remnant’s mass and spin. Measuring more than two modes makes it possible to test the Kerr hypothesis.

To date, not as much attention has been given to what can be learned by measuring the amplitudes of ringdown modes (though see Ref. London (2018) for an important recent example to the contrary). The relative excitation of different modes depends on the geometry of the system as it approaches its final state. This geometry, in turn, depends upon the astrophysical history of the binary, with the spins and orbit expected to be substantially aligned for compact binaries that form in situ from stellar binaries, and with these angular momentum directions largely randomly oriented for binaries that form dynamically through multi-body interactions (e.g., in globular clusters); see Refs. Mandel and Farmer (2018); Abbott et al. (2016); Gerosa et al. (2018) for recent discussion and review.

For a coalescence in which the orbital angular momentum is nearly parallel to the large black hole’s spin, the $(\ell,m)=(2,2)$ ringdown mode is likely to be the most strongly excited. If the orbit is substantially misaligned from the larger hole’s spin axis, then modes with $(\ell,m)=(2,2)$, $(2,1)$, or $(2,0)$ might be excited by roughly the same amount. One must model all of the binary’s GW modes and multipoles to completely assess what can be learned from their measurement London et al. (2018). For the ringdown, a careful analysis is needed to understand how the excitation of these modes depends on the coalescence geometry.

Our goal is to begin developing such an analysis. We use black hole perturbation theory to provide an easily parameterized framework for studying how different modes are excited by binary black hole coalescence. Strictly speaking, our results thus only describe the limit in which one member of the binary is far more massive than the other. We expect, however, that insight from this limit will carry over to coalescences with general mass ratio at least qualitatively, and perhaps even provide good quantitative understanding for mass ratios larger than some threshold. Numerical relativity will be needed to explore late mode excitation by highly misaligned coalescences for general mass ratios.

We find that mode excitation varies in a predictable fashion as a function of certain parameters that describe the final coalescence geometry. Given the spin $a$ of the larger black hole, we find the relative excitation of black hole modes is determined by two angles — an angle $I$ describing the inclination of the binary’s orbit relative to the black hole’s equatorial plane, and an angle $\theta_{\rm fin}$ which describes the polar location at which the smaller body crosses the larger black hole’s event horizon — and its angular direction $\operatorname{sgn}(\dot{\theta}_{\rm fin})$ as it crosses. (Because we work in perturbation theory and linearize all deviations from Kerr in the small body’s mass $\mu$, the absolute excitation of ringdown modes also depends on this mass.)

This suggests that, if measurements by GW detectors can accurately measure multiple ringdown modes, it may be possible to use the relative amplitude of these modes to learn about the binary’s spin-orbit misalignment. This may be particularly valuable for ground-based measurements of high mass systems for which the inspiral waves (which encode spin-orbit misalignment through amplitude and frequency modulation of the waveform) are poorly measured. In other words, it may be possible to get information about the system’s spin-orbit misalignment from the late ringdown signal. At the very least, the ringdown may provide information that complements spin-orbit constraints obtained from earlier in the waveform, improving our ability to make inferences about the nature of a measured binary black hole.

We begin with a synopsis of how we compute GWs in Sec. II. We review the results of Paper I in Sec. II.1, describing how we build the worldline which the smaller member of the binary follows as it inspirals and then plunges into the larger black hole. As described in Sec. II.2, we then use this worldline to build the source term for the time-domain Teukolsky equation, which allows us to compute GWs produced by a small body following that worldline. The last several to several dozen cycles of the waveforms that we compute consist of ringing modes of the larger black hole. We characterize this mode content in Sec. III. We review the properties of Kerr black hole quasinormal modes in Sec. III.1, discuss important issues regarding the angular bases used to described these modes and our data in Sec. III.2, and then describe our algorithm for extracting the mode content of our waveforms in Sec. III.3. In most of the cases that we examine, the ringdown content we find is accurately described as entirely due to superpositions of the “fundamental” quasinormal black hole modes. In some cases, we believe that we may be able to discern the presence of the first overtone mode. Extracting this mode from our waveform data requires some care; our procedure for doing this extraction is described in Sec. V.3.

Section IV describes how we parameterize ringdown, a describes checks we instituted to make sure our results are reasonable. We begin in Sec. IV.1 by describing how the plunges which produce ringdown waveforms can be characterized using four parameters: the black hole spin $a$, an angle $I$ describing the inclination of the orbital plane, an angle $\theta_{\rm fin}$ that describes where the plunge terminates at the black hole, and the direction of the final plunge’s angular motion, $\operatorname{sgn}(\dot{\theta}_{\rm fin})$. As expected, we find that these parameters work well to describe the ringdown waveform. We next discuss in Sec. IV.2 certain symmetries which ringdown waveforms should inherit from the Kerr spacetime, and verify that our results respect these symmetries. Finally, we describe a comparison to previous results for the equatorial case in Sec. IV.3. As we describe, we do not expect perfect agreement because of differing methodologies, but we find good agreement over the regime where agreement is expected.

We show results in Sec. V. We begin with a catalog of modes with spheroidal index $\ell=2$ in Sec. V.1; additional modes are shown in Appendix D. This catalog indeed demonstrates that the excitation of each mode depends cleanly and predictably on the geometry of the final plunge, suggesting that the inverse problem — inferring the properties of the plunge geometry from a spectrum of measured modes — may be feasible. We describe in more detail interesting features we find in this catalog in Secs. V.2 and V.3, and discuss the impact of numerical errors in Appendix. B. In Sec. V.2, we report interesting universality behavior that appears to emerge, at least at shallow inclination angle, as we examine mode excitation. In brief, across a wide range of spins, excitation of the fundamental $(\ell,m)$ mode appears to follow a universal functional form that depends only on $\ell-m$ or $\ell+m$. In Sec. V.3, we present results which suggest that the first overtone mode can be discerned in our results for rapid spin and shallow inclination angle, and provide further evidence in Appendix C.

Our concluding discussion is given in Sec. VI. In addition to summarizing our findings, we suggest directions for future work that may improve our ability to model ringdown from misaligned binaries and to use these models to learn about coalescing black holes from their measured gravitational waves. A concise discussion of our key findings and their implications can be found in a companion Letter Hughes et al. (2019).

II Gravitational waves from large-mass-ratio inspiral and plunge

To model the GWs from a large-mass-ratio binary system, we first compute the worldline that a small body follows as it adiabatically evolves through a sequence of geodesic orbits, and then compute the transition to a plunging trajectory that carries it into the larger black hole. We use this worldline to build the source function to the time-domain Teukolsky equation Sundararajan et al. (2007, 2008), solve this equation Sundararajan et al. (2010); Zenginoğlu and Khanna (2011) to compute the GWs that arise from the small body’s motion on this worldline, and then characterize the quasinormal modes (QNMs) that terminate the waveform we find. In this section, we briefly summarize the key steps in this analysis in order to set up how we characterize the QNM content of these GWs.

II.1 Computing the inspiral, transition, and plunge worldline: A brief synopsis

Paper I discusses in detail how we construct the worldline describing how the small member of the binary inspirals and then plunges into its companion black hole. Here we recap the main points of Paper I, emphasizing aspects of that analysis that are important for this paper:

1. (1)

The small body begins on an inclined circular Kerr geodesic orbit, parameterized by initial radius $r_{0}$ and inclination angle $I$. All orbit properties vary smoothly over the domain $I=0^{\circ}$ (prograde equatorial) to $I=180^{\circ}$ (retrograde equatorial). The initial orbit is taken to adiabatically inspiral due to GW emission, shrinking in radius while keeping $I$ nearly constant Drasco and Hughes (2006). Inspiral describes the system until the small body comes near the separatrix separating stable from unstable orbits. 2. (2)

As the small body approaches the separatrix, inspiral accelerates until its inward motion is no longer adiabatic. We compute the transition connecting inspiral to the final plunge by following the principles developed by Ori and Thorne Ori and Thorne (2000), who showed how to do this for equatorial orbits. Our “generalized” Ori-Thorne procedure lifts this restriction, allowing us to compute the transition for arbitrarily inclined circular orbits. 3. (3)

At very late times, the small body’s motion is well approximated by a plunging geodesic (with constant parameters) that crosses the event horizon. In Boyer-Lindquist coordinates, this plunge terminates the small body’s motion at some $\theta_{\rm fin}$ on the horizon. This is an artifact of Boyer-Lindquist time, and reflects the fact that the event horizon is a surface of infinite redshift. As we discuss in Sec. II.2, this behavior ensures that the time-domain Teukolsky equation goes over to its homogeneous form at late times, which in turn ensures that the final GW cycles we compute are QNMs, coherently joined in phase to the preceding waveform.

Our worldline model requires us to make three ad hoc choices. The first is how to define the end of “inspiral” and the beginning of “transition.” The physics governing the transition picks out a range of times for this moment. We show in Paper I that the worldlines we develop vary very little over that range. The second choice is how to model the evolution of the orbit’s parameters $E$, $L_{z}$, and $Q$ during the transition. In Paper I, we refine the Ori-Thorne procedure to eliminate unphysical discontinuities present in their model, but note that there are many ways to implement this refinement. The two which we have investigated barely differ from one another. We suspect this would be the case for any reasonable model.

The third ad hoc choice describes when “transition” ends and “plunge” begins. As with the first choice, the physics of the transition picks out a range of times for this moment. Disconcertingly, we find that the worldlines we develop depend on this parameter in a non-negligible way (see Sec. V B of Paper I, especially Fig. 7). This could raise concerns that our conclusions will not be robust, but will depend upon how we make this choice.

Fortunately, we find that the ringdown modes we find are robust with respect to this choice, even though the plunge worldlines are not. As we have outlined in the Introduction, the ringdown modes we find depend on the spin $a$ of the larger black hole, the inclination angle $I$, and an angle $\theta_{\rm fin}$ which defines where the smaller body plunges into the black hole. As we vary the end of transition, the relation between $\theta_{\rm fin}$ and the worldline’s initial conditions can change by quite a bit. However, the dependence of the ringdown waves on $\theta_{\rm fin}$ does not depend on this choice. As long as we parameterize our modes using the parameter set $\bm{(}a,I,\theta_{\rm fin},\operatorname{sgn}(\dot{\theta}_{\rm fin})\bm{)}$, our conclusions are robust against how we choose to end the transition. See Appendix A for detailed discussion and results.

Although it is a relief that our physical conclusions are not impacted by how we choose these parameters, it is a fundamental flaw of the generalized Ori-Thorne model that these ad hoc choices exist. Work to improve this, or at least to better inform how these choices should be made, would be salubrious.

II.2 Solving the time-domain Teukolsky equation

Following the procedure of Paper I summarized in the previous section, we make the worldline which the small body follows as it inspirals and plunges into the large black hole. We then use this worldline to build the source term for the time-domain Teukolsky equation. This equation describes scalar, vector, and tensor field perturbations to the spacetime of a rotating black hole. In Boyer-Lindquist coordinates, it takes the form Teukolsky (1973)

[TABLE]

where $M$ is the mass of the black hole, $a$ is its angular momentum per unit mass, $\Delta=r^{2}-2Mr+a^{2}$, and $s$ is the “spin weight” of the field. For $s=-2$, this equation describes the radiative degrees of freedom of the gravitational field, and is related to the Weyl curvature scalar as $\Psi=(r-ia\cos\theta)^{4}\psi_{4}$. At future null infinity,

[TABLE]

With this quantity in hand, $h_{+}$ and $h_{\times}$ can be easily computed by a double time-integration.

The source $T$ in Eq. (II.2) is computed from the small body’s energy-momentum tensor,

[TABLE]

where $\Sigma=r^{2}+a^{2}\cos^{2}\theta$ and $u^{\alpha}$ denotes components of the small body’s 4-velocity along its worldline. To construct $T$, project $T_{\alpha\beta}$ onto certain legs of the Kinnersley tetrad, and then operate upon the resulting quantity with a second-order differential operator. See Ref. Teukolsky (1973) for detailed discussion.

We solve Eq. (II.2) in the time domain, dynamically providing information about the small body’s worldline to make the source $T$ for our inspiraling and plunging body; Fig. 1 shows a representative example of the waveform produced by this procedure. Details of our approach have been extensively described in past literature Sundararajan et al. (2007, 2008, 2010); Zenginoğlu and Khanna (2011), so we do not repeat this discussion here, modulo one remark that is significant for the purposes of this current work. Notice that the stress-energy tensor $T_{\alpha\beta}$ is inversely proportional to $u^{t}=dt/d\tau$. This factor converts between time $\tau$ along the worldline and time $t$ as measured by a distant observer. As the small body approaches the event horizon, $dt/d\tau\to\infty$. The source term thus “redshifts away,” smoothly converting the Teukolsky equation to its homogeneous form and connecting the GWs from the small body’s plunge phase to the Kerr black hole’s QNMs very naturally.

III Characterizing a waveform’s quasinormal mode content

The analysis presented in Sec. II yields the waveform produced by a small body that inspirals and plunges into a black hole; an example is shown in Fig. 1. The final cycles of these waveforms can be modeled as a linear superposition of QNMs. In this section, we describe how we characterize and extract the QNM content that describes the final cycles of inspiral and plunge waveforms like that shown in Fig. 1.

III.1 quasinormal modes

Teukolsky showed Teukolsky (1973) that Eq. (II.2) separates after decomposing in the frequency domain, yielding ordinary differential equations which govern the $\theta$ and $r$ dependence of the field $\Psi$. The solution for a given mode with (possibly complex) frequency $\sigma$ then becomes

[TABLE]

The functions ${{}_{s}}S^{a\sigma}_{\ell mn}(\theta,\phi)$ are known as spin-weighted spheroidal harmonics; when $a=0$, they reduce to the spin-weighted spherical harmonics ${{}_{s}}Y_{\ell m}(\theta,\phi)$ Goldberg et al. (1967). We discuss these functions and their properties in further detail below. See Ref. Berti et al. (2009) for a discussion of the radial function ${{}_{s}}R^{a\sigma}_{\ell mn}(r)$, as well as for additional details regarding the spin-weighted spheroidal harmonics.

For certain frequencies, the modes (4) satisfy physical boundary conditions: they describe radiation that is purely ingoing on the black hole’s event horizon, and purely outgoing at null infinity. Such solutions are the black hole’s QNMs. Frequencies for which such solutions hold are written $\sigma_{\ell mn}$, each labeled by the mode’s spheroidal harmonic indices plus an overtone index $n=0,1,\ldots$ Press and Teukolsky (1973); Berti et al. (2009). These frequencies are complex, and can be written

[TABLE]

Given mode indices $(\ell,m,n)$ and assuming that the endstate of the merged system is a Kerr black hole, these frequencies depend only on the black hole’s mass $M$ and spin parameter $a$. Code for computing $\omega_{\ell mn}$ and $\tau_{\ell mn}$ for ringing modes of Kerr black holes, as well as tables describing the results, are provided by Berti Berti et al. (2009, 2006).

The Teukolsky equation’s symmetry properties dictate that for each mode $(\ell,m,n)$, the separated angular equation is actually satisfied by two eigenfunctions: ${}_{s}S^{a\sigma}_{\ell mn}(\theta,\phi)$ and ${}_{s}S^{a\sigma}_{\ell-mn}(\pi-\theta,\phi)^{*}$, where ∗ denotes complex conjugation. A given QNM is thus actually characterized by pairs of frequencies $\omega_{\ell mn}$ and $-\omega^{*}_{\ell-mn}$ Berti et al. (2006). Specializing to spin weight $s=-2$ and bearing this symmetry in mind, we follow Ref. Berti et al. (2006) [cf. their Sec. II A, especially Eq. (2.9) and nearby text] and write the gravitational waveform for a mode as

[TABLE]

(We explain why we have shifted the spheroidal index from $\ell\to k$ in the next section.) Here ($\mathcal{A}_{kmn},\mathcal{A}^{\prime}_{kmn})$ and $(\phi_{kmn},\phi^{\prime}_{kmn})$ are the mode amplitude magnitudes and phases111Our notation, which associates the QNM amplitudes labeled by $(k,m)$ with $\mathcal{A}_{kmn},\mathcal{A}^{\prime}_{kmn},\phi_{kmn},\phi^{\prime}_{kmn}$, is aligned with Refs. Berti et al. (2006); Buonanno et al. (2007), but differs from Ref. Taracchini et al. (2014)[cf.Eq.(5)]. See Sec. IV.3 for discussion on converting between these different conventions. and $t_{0}$ marks the time at which quasinormal ringing begins. Each mode’s absolute excitation depends in addition on the small body’s mass $\mu$ and the distance $D$ to the binary: ${\cal A}_{\ell mn}\propto\mu/D$. In the results we present, we set this factor to 1; one can then scale the amplitudes by multiplying by $\mu/D$.

We will imagine a waveform that is ringdown dominated and thus is accurately described by Eq. (LABEL:eq:spheroidaldecomp) for $t\geq t_{0}$. However, for perturbations sourced by a plunging body, it is not possible to absolutely determine the start of the ringdown Dorband et al. (2006); Berti and Cardoso (2006). Previous studies considering comparable mass binaries have associated $t_{0}$ with various phenomenological indicators such as the peak of the gravitational-wave amplitude Buonanno et al. (2007); London et al. (2014) or the peak orbital frequency Taracchini et al. (2014). Such associations are not well adapted to inclined, extreme mass ratio orbits since the radiation and orbital frequency do not always exhibit a clear peak. In the analysis we describe below, we will take $t_{0}$ to be the time at which the small body reaches a radius equal to that of the prograde equatorial photon orbit (see Refs. Khanna and Price (2017); Price et al. (2016) for further discussion). Although the mode amplitudes $\mathcal{A}_{kmn},\mathcal{A}^{\prime}_{kmn}$ and phases $\mathcal{\phi}_{kmn},\mathcal{\phi}^{\prime}_{kmn}$ that we find then depend on $t_{0}$ (see Ref. Dorband et al. (2006) for further discussion on the time-shift problem), our discussion of how the mode excitation depends on $\bm{(}I,a,\theta_{\rm fin},\operatorname{sgn}(\dot{\theta}_{\rm fin})\bm{)}$ will be independent of $t_{0}$.

III.2 Spherical and spheroidal expansions

The numerical code which we use to solve Eq. (II.2) Sundararajan et al. (2007, 2008, 2010); Zenginoğlu and Khanna (2011) decomposes the gravitational radiation as

[TABLE]

where $h_{\ell m}^{\rm N}(t)=h_{\ell m,+}^{\rm N}(t)-ih_{\ell m,\times}^{\rm N}(t)$ is the $(\ell,m)$ spherical multipole component. The superscript N emphasizes that each component is output from our numerical code. Since we model the ringdown in the spheroidal basis, we must take into account spherical-spheroidal mode mixing Berti and Klein (2014); Cook and Zalutskiy (2014); London and Fauchon-Jones (2018). This occurs because basis functions with the same degree $m$ overlap:

[TABLE]

where $k$ is the spheroidal harmonic order and $\ell$ is the spherical harmonic order. In the Schwarzschild limit, $\mu_{m\ell kn}$ collapses to $\delta_{k\ell}$. Our definition of the overlap coefficient coincides with that used in Ref. Berti and Klein (2014), cf. their Eq. (5). To avoid confusion, we use $k$ to label spheroidal harmonic decomposed QNMs, and label decompositions based on spherical harmonics with $\ell$. Our code for computing the coefficients $\mu_{m\ell kn}$ is based on the algorithm described in Appendix A of Ref. Hughes (2000) (see also Ref. Hughes (2003)). It has been validated to numerical precision using comparison data kindly provided by Cook, and agrees with Berti’s online data tables (modulo an unimportant overall phase).

Equating the left-hand side of Eq. (7) to the left-hand side of (LABEL:eq:spheroidaldecomp), multiplying both sides by ${}_{-2}Y^{*}_{\ell m}(\theta,\phi)$, integrating over the sphere and using Eq. (8), we find

[TABLE]

where $k_{\rm min}={\rm max}(2,|m|)$. In practice, the sums are truncated at some finite maximum index; we discuss this truncation in detail in the next section. The mode amplitudes and phases have been absorbed here into complex amplitudes,

[TABLE]

and the time dependent coefficients are given by

[TABLE]

To derive these results, we have used the fact that

[TABLE]

It is important to note that the spherical multipole $h^{\rm N}_{\ell m}$ contains, in general, contributions from spheroidal modes of both $m$ and $-m$.

III.3 Mode extraction

Our work to characterize mode amplitudes can be considered an extension of Taracchini et al. Taracchini et al. (2014), whose analysis was limited to equatorial trajectories ($I=0^{\circ}$ and $180^{\circ}$). Initially, we extended their algorithm to include higher-order modes. However, with this algorithm the calculated mode amplitudes, especially sub-dominant modes, are rather sensitive to the choice of fitting interval. This makes it essentially impossible to draw robust conclusions about the waveform’s mode content. Other studies in the literature that rely on fits to damped sinusoids in either the time or frequency domain have dealt with this ambiguity by choosing the fitting interval with lowest error (Dorband et al. (2006); Kelly and Baker (2013); London et al. (2014); Buonanno et al. (2007)). This choice leads to discontinuities in the mode amplitudes as a function of $\theta_{\rm fin}$, which is a general symptom of the fact that the fit error (as a function of chose fit interval) suffers from multiple local minima. Oscillations in fit error, similar to those discussed in Ref. Berti et al. (2007), lead to jumps in the calculated mode amplitude as $\theta_{\rm fin}$ is incremented.

Our aim also differs from that of numerical relativity analyses. Since we are working within black hole perturbation theory, we know the BH spin and mass a priori and can assume that the ringdown comprises those Kerr QNM frequencies. In addition, numerical relativity analyses have primarily focused on nonprecessing comaprable mass systems, where the effect of higher-order mode mixing [specifically beyond the (2,2) and (3,2) QNMs] can be neglected. It remains to be determined whether the same trend holds for highly precessing systems with larger mass ratios. Instead of neglecting mode mixing, our algorithm specifically addresses how to disentangle mode mixing by factoring in the higher-order spherical waveform modes, for example in utilizing $\psi^{\rm N}_{22},\psi^{\rm N}_{32}$ and $\psi^{\rm N}_{42}$ simultaneously. An exception is Ref. London et al. (2014), which adopts an algorithm to measure higher-order modes and overtones. However, this algorithm is based on a fit to a frequency domain waveform, which leads to mode amplitudes that still depend on the choice of fitting interval.

Before describing our mode extraction algorithm in detail, we introduce some simplifications to the general framework we have described. First, we will typically neglect overtone modes, $n\geq 1$. The damping times $\tau_{kmn}$ tend to decrease fairly rapidly with $n$ Leaver (1985). The fundamental mode $n=0$ thus dominates the ringdown, especially for “late” times $t-t_{0}\gtrsim 25M$. We thus mostly focus on $n=0$, and drop the sum over $n$ (though see the following section for a discussion about including the first overtone).

Second, for the situations that we consider, the coefficients $\mu_{m\ell kn}$ tend to peak at $k=\ell$, falling off rapidly in magnitude away from this peak. As such, we truncate the sum over $k$ at $k_{\rm max}=\ell+K_{\ell}$, where $K_{\ell}$ is found by determining how many spheroidal modes must be included to accurately fit the spherical mode numerical data $h^{\rm N}_{\ell m}(t)$; we elaborate on this below. Equation (9) becomes

[TABLE]

With this framework in hand, our goal now is to calculate the mode amplitudes $\mathcal{C}_{km0}$ and $\mathcal{C^{\prime}}_{km0}$ given $\tensor{*}[]{h}{{}_{\ell m}^{\rm N}}(t)$, for which we have developed the following algorithm:

• •

Consider a set of ${\cal N}$ spherical multipoles $(\ell_{i},m)$. Choose $k_{\rm min}=\ell_{1}$ where $\ell_{1}={\rm max}(2,|m|)$, $\ell_{2}=\ell_{1}+1$, $\ldots$, $\ell_{\cal N}=\ell_{1}+{\cal N}-1$. Choose the factor $K_{\ell}$ appearing in Eq. (13) such that $K_{\ell_{1}}={\cal N}-1$, $K_{\ell_{2}}={\cal N}-2$, $\ldots$, $K_{\ell_{\cal N}}=0$.

• •

Consider a moment $t=t_{j}$; evaluate Eq. (13) and its time derivative at this moment. This yields $2{\cal N}$ linear equations which can be inverted for the $2{\cal N}$ unknown mode amplitudes ${\cal C}_{km0}$ and ${\cal C}^{\prime}_{k-m0}$.

• •

Check consistency of the algorithm by computing the mode amplitudes at multiple times. If the ringdown model is consistent with the data, then we find that each ${\cal C}_{km0}$ and ${\cal C}^{\prime}_{k-m0}$ settles down to a constant during the time period in which the waveform is QNM dominated. If these amplitudes do not stabilize to a constant, then the model is not adequate. It may be that the model needs to include more spheroidal modes (in which case we increase ${\cal N}$ and repeat the algorithm), or that the radiation is simply not QNM dominated.

• •

When the model is consistent with data, we denote by an overbar the values to which the amplitudes settle down: $\bar{\mathcal{C}}_{km0}$, $\bar{\mathcal{C^{\prime}}}_{k-m0}$. In practice, we determine $\bar{\mathcal{C}}_{km0}$ and $\bar{\mathcal{C^{\prime}}}_{k-m0}$ by calculating a moving average with fixed size $\Delta t$. We then associate $\bar{\mathcal{C}}_{km0}$ and $\bar{\mathcal{C^{\prime}}}_{k-m0}$ with the average in the interval with the least variance. The ringdown model becomes

[TABLE]

The “0” in the superscript on the right-hand side of Eq. (14) labels the fact that this model is based on the $n=0$ fundamental ringdown modes.

We illustrate this algorithm with an example. The left panel of Fig. 2 shows the final several numerical GW cycles $h^{\rm N}_{22}$ arising from a small body inspiraling and plunging into a black hole with $a=0.9M$. The orbit is inclined at $I=20^{\circ}$, and crosses the horizon at a final polar angle $\theta_{\rm fin}=72.1^{\circ}$ with angular velocity $\dot{\theta}_{\rm fin}<0$. We first try to fit these waves using a model with ${\cal N}=1$ (i.e., using only the $\ell=2$ modes). The middle top panel of Fig. 2 shows the coefficients ${\cal A}_{220}$ and ${\cal A}^{\prime}_{2-20}$ we find in this case. We find no span of time at which ${\cal A}_{220}$ and ${\cal A}^{\prime}_{2-20}$ settle down to constants. The choice ${\cal N}=1$ poorly describes these data.

Consider next fits using a model with ${\cal N}=2$ (i.e., now using the $\ell=2$ and $\ell=3$ modes). As the middle bottom panel of Fig. 2 shows, the fit is improved, but we still see oscillations in the extracted amplitudes. This indicates that there is still room for improvement in this model. Finally, the top right-panel shows the fit for ${\cal N}=3$ (now including modes $\ell=2$, $\ell=3$, and $\ell=4$). Here at last we find that the amplitudes have settled down to a nearly constant level, at least over the time interval $25M\lesssim t-t_{0}\lesssim 120M$. At earlier times, the signal is not yet QNM dominated; at later times, the modes have decayed away, and our fit becomes noise dominated. The residual $\delta h^{\rm N}_{22}\equiv h^{\rm N}_{22}-h^{\rm RD0}_{22}$ shown in the lower-right panel illustrates that the fit describes the numerical data well over this time interval.

III.4 Overtones

The ringdown model $\tensor{*}[]{h}{{}_{\ell m}^{\rm RD0}}(t)$ only includes fundamental QNMs. As we will show in Sec. IV, this model is consistent with data during the late ringdown. However, in the early ringdown, overtone modes (which tend to decay much more quickly than the fundamental) may be present, and a “fundamentals-only” model will not capture their contributions to the radiation that arises from the small body’s final plunge. For larger spins (for which the rapid decay of the overtones is not so rapid), it may be feasible to isolate their contribution to the early ringdown. Past work has similarly proposed techniques to isolate overtones from numerical relativity models of binary black hole coalescence London et al. (2014), where the initial “perturbation” is very large and overtone excitation is important.

Suppose we have computed a fundamental-only ringdown model $h^{\rm RD0}_{\ell m}(t)$. Define the residuals from this model as

[TABLE]

In principle, we could repeat the algorithm described in Sec. III.3, but using $\delta h^{\rm N}_{\ell m}$ rather than $h^{\rm N}_{\ell m}$, and fitting with modes which have $n=1$. In practice, because the overtones are quite short-lived, we do not find a time domain over which their amplitudes settle down to a constant level. In Sec. V.3, we describe a modification to our mode extraction algorithm which we use to account for this difficulty and to estimate overtone amplitudes.

One could imagine iterating further, yielding ever higher-order fits to the ringdown overtones. In practice, we expect that this method will be greatly limited by numerical accuracy, and that it is likely to be quite challenging to find all but perhaps the $n=1$ overtones. Even in that case, overtones are most likely to be discernible only if the black hole’s spin is quite rapid (so that the overtone’s decay is relatively slow). Evidence that we may be finding the first overtone with $a=0.99M$ is presented in Sec. V.3.

IV Parameterization, checks, and comparisons with past work

Before discussing the results we find, we describe how we parameterize our data, examine symmetries that our results should respect, and check that our results behave as expected in the Schwarzschild limit, where spherical symmetry implies certain relations among the different amplitudes. We also verify that our results agree with past work in the equatorial limit. We then explore the inclined Kerr case in Sec. V.

IV.1 A clean and complete parameterization

The ringing cycles that we wish to study are sourced by the final moments of the small body’s worldline. To describe those waves, we need a parameterization that completely characterizes the small body’s final motion on its worldline as it plunges into the black hole.

Begin by considering the circular and inclined geodesic orbits on which the small body initially moves. Circular geodesics of Kerr black holes are generally described using the orbit’s radius $r$ and some angle describing the orbit’s tilt from the equatorial plane. We use the angle $I$, defined (in radians) by

[TABLE]

where $\theta_{m}$ is the minimum value of the Boyer-Lindquist angle $\theta$ reached on an orbit. $I$ smoothly varies from $0^{\circ}$ for prograde equatorial orbits to $180^{\circ}$ for retrograde equatorial orbits. For Schwarzschild, $I$ corresponds exactly to the angle at which the orbit is inclined from the equatorial plane. Although not quite amenable to such a simple interpretation for general spin, it provides a useful notion of orbit tilt for Kerr as well.

An orbit with inclination $I$ oscillates in $\theta$ in the range

[TABLE]

It turns out that $I$ remains nearly constant during the small body’s inspiral and plunge Hughes (2000). This means that over the small body’s worldline, its polar position oscillates over the range given by Eq. (17) until the moment that it crosses the larger black hole’s event horizon. In Boyer-Lindquist coordinates, the small body freezes at the angle $\theta_{\rm fin}$ where it enters the event horizon.

As discussed at length in Paper I, two worldlines with the same initial values of $r$ and $I$ can begin with different values of $\theta$. Although these worldlines will follow the same trajectories in $(r,I)$ during inspiral and plunge, their worldlines will follow different trajectories in $\theta$, eventually crossing the event horizon at different final polar angles $\theta_{\rm fin}$. An example of this is shown in Fig. 5 of Paper I.

Given $I$, there are in fact two worldlines which terminate at each value of $\theta_{\rm fin}$ allowed by the range (17): one with $\dot{\theta}>0$ during the final plunge, and one with $\dot{\theta}<0$ in these moments. This is illustrated in Fig. 3. Although the final state in both cases is identical, how the system reaches that final state is quite different. The QNM signature of these two situations differs accordingly as well.

We thus find that a waveform’s QNM content depends on 4 parameters which characterize its final plunging behavior: the black hole spin $a$, the orbital inclination $I$, the final polar angle $\theta_{\rm fin}$, and the sign of angular velocity in the plunge’s final moments, $\operatorname{sgn}(\dot{\theta}_{\rm fin})$. As we show in Sec. V, these parameters completely describe the mode amplitudes $\mathcal{C}_{kmn}$ and $\mathcal{C}^{\prime}_{kmn}$:

[TABLE]

IV.2 Mode excitation symmetry characteristics

The properties of black hole spacetimes imply that certain symmetries should exist in the mode excitation. First consider two plunges that are on worldlines $u^{(1)}_{\alpha}$, $u^{(2)}_{\alpha}$ with the same orbital inclination about a given black hole (so that they share values of $a$ and $I$). Imagine that these plunges approach the horizon such that

[TABLE]

The Kerr spacetime is reflection symmetric about the equatorial plane ($\theta\rightarrow\pi-\theta$), so these infalling bodies have identical worldlines modulo a reflection which inverts in $\cos\theta$ (and possibly modulo a rotation $\phi\rightarrow\phi+\delta\phi$, which is also a continuous symmetry of Kerr).

As the worldlines are related by a reflection, the radiation they produce, which is decomposed onto the harmonics $\tensor{*}[_{-2}]{Y}{{}_{\ell m}}(\theta,\phi)$, should also be related by a reflection. The spin-weighted spherical harmonics transform under reflection by Eq. (12). This implies that the primed and unprimed mode amplitudes of the two worldlines are related by the following reflection symmetry:

[TABLE]

For each QNM, it suffices to only specify either $\mathcal{C}_{kmn}\bm{(}a,I,\theta_{\rm fin},\operatorname{sgn}({\dot{\theta}_{\rm fin}})\bm{)}$ or $\mathcal{C}^{\prime}_{kmn}\bm{(}a,I,\theta_{\rm fin},\operatorname{sgn}({\dot{\theta}_{\rm fin}})\bm{)}$. In all of our calculations, we find that the mode amplitudes we find are consistent with the symmetry (21). Figure 4 shows this explicitly for one example (which we discuss in much greater depth in Sec. V). Moving forward, we will generally only show one of $\mathcal{C}_{kmn}$ or $\mathcal{C}^{\prime}_{kmn}$ in our results.

Second, consider the Schwarzschild limit, $a=0$. Because of this spacetime’s spherical symmetry, there is no unique notion of an “equatorial plane,” so any two worldlines $u^{(1)}_{\alpha}$ and $u^{(2)}_{\alpha}$ that begin from circular orbits can be related to each other by a rotation. Consider the rotation generated by the quaternion $\mathbb{R}$ which relates the spatial components of the worldlines $\mathbf{u}^{(1)}$ and $\mathbf{u}^{(2)}$ as

[TABLE]

Now, decompose the radiation resulting from $u^{(i)}_{\alpha}$ onto the spin-weighted spherical harmonics,

[TABLE]

where $(\theta_{i},\phi_{i})$ are the coordinates in which the angular components of $\mathbf{u}^{(i)}$ are expressed. Then the radiation multipoles $h^{(1)}_{\ell m}$ and $h^{(2)}_{\ell m}$ can be related to each other by the basis transformation formula,

[TABLE]

where $\mathfrak{D}^{\ell}_{m^{\prime}m}(\mathbb{R})$ is the Wigner rotation matrix, which takes on a simple form in terms of $\mathbb{R}$ Boyle (2016). Using Eq. (24) and the fact that $\tensor{*}[_{-2}]{S}{{}^{a\sigma}_{kmn}}=\tensor{*}[_{-2}]{Y}{{}_{km}}$ for $a=0$, we can relate the mode amplitudes between any two worldlines as

[TABLE]

provided that their trajectories satisfy Eq. (22).

We use Eq. (25) to test our mode extraction algorithm as follows. Let $u^{(1)}_{\alpha}$ be the worldline of a small body on an equatorial trajectory ($I^{(1)}=0$), and let $u^{(2)}_{\alpha}$ be the worldline of a small body on an inclined trajectory ($I^{(2)}>0$). By Eq. (25), the mode amplitudes calculated from the $u^{(2)}_{\alpha}$ plunge should transform to amplitudes computed from the equatorial $u^{(2)}_{\alpha}$ plunge.

Figure 5 shows our results for $I^{(2)}=20^{\circ}$ and $I^{(2)}=60^{\circ}$ using several different values of $\theta^{(2)}_{\rm fin}$. The black circle at $(\theta_{\rm fin},{\cal A}_{220})=(90^{\circ},2.3)$ shows the mode amplitude we find for the equatorial worldline $u^{(1)}_{\alpha}$. The blue and red dots show the amplitudes ${\cal A}_{220}$ we find for the various inclined worldlines $u^{(2)}_{\alpha}$ we consider. The black dashes show the amplitude that we find by applying Eq. (25) to rotate the inclined case into the equatorial plane. At high inclination ($I^{(2)}=60^{\circ}$), the rotated mode amplitudes agree with the equatorial amplitude to within $8.5\%$. As we’ll discuss in Secs. V.1 and B, numerical noise appears to limit our accuracy for high orbit inclination. At lower inclination ($I^{(2)}=20^{\circ}$) the numerical error is reduced, and we find that all rotated mode amplitudes agree with the equatorial amplitudes to within $2.0\%$. In Appendix B, we further discuss the impact of numerical errors on our results. Code for computing the Wigner rotation matrix elements in terms of quaternions is provided by M. Boyle Boyle (2018).

IV.3 Comparison with past equatorial results

As a final check, we compare our results for equatorial inspiral and plunge ($I=0^{\circ}$ and $I=180^{\circ}$) with those given by Taracchini et al., Ref. Taracchini et al. (2014). We do not expect these two analyses to agree perfectly. This is in part due to differences in the trajectory and waveform calculations. Also, the method used in Ref. Taracchini et al. (2014) to extract QNMs from late waveforms is quite different from that which we developed. In particular, the number of spheroidal modes that are included in the fit in Ref. Taracchini et al. (2014) is typically fewer than the number of modes that we use here. Since mode mixing is strongest at large spin, we expect the largest systematic differences at large $a$. Given these considerations, we limit our comparison by analyzing the same waveform data from Ref. Taracchini et al. (2014), which was provided by the authors (two of whom are authors of this paper). We also only model the modes which were considered in Ref. Taracchini et al. (2014). To properly compare with Ref. Taracchini et al. (2014), we also have to adjust our notation slightly by relating the amplitudes defined here [appearing in Eq. (LABEL:eq:spheroidaldecomp)] to those defined in Ref. Taracchini et al. (2014) [cf. their Eq. (5)]. For prograde orbits $a\geq 0$ ($I=0^{\circ}$),

[TABLE]

For retrograde orbits, the waveform multipoles in Ref. Taracchini et al. (2014) (which we refer to as $h^{\rm N}_{\ell m}$) were computed in a coordinate system where $L_{z}$ is positive. Although this is opposite of our convention, the computed modes in this study can be related to those modes computed in Ref. Taracchini et al. (2014) by a $180^{\circ}$ coordinate rotation [see Eq. (25)]. Thus for $a<0$,

[TABLE]

We also adapt the conventions used in Ref. Taracchini et al. (2014) for choosing $t_{0}$.

Table 1 shows the result of this comparison. We indeed find decent agreement with their mode amplitudes, at least for $a\leq 0.5M$. The lack of perfect agreement, and disagreement at $a=0.9M$, is to be expected given the rather different methods of calculating the mode excitation. For instance, our analysis suggests that the $(k,m,n)=(3,1,0)$ QNM is excited, but this mode is not modeled in Ref. Taracchini et al. (2014). This leads to larger systematic differences between the two methods at high spin, $a=0.9M$.

V Results

We now present detailed results describing the mode excitation we find for inclined plunges of rotating black holes. Figure 6 shows the range of cases that we examine. Each point in this figure represents a particular choice of $a$ and $I$ that we study. For each non-equatorial case ($I\neq 0^{\circ}$ and $I\neq 180^{\circ}$), we examine 36 different values of $\theta_{\rm fin}$, 18 for each sign of $\dot{\theta}_{\rm fin}$. For the equatorial cases, we always have $\theta_{\rm fin}=90^{\circ}$, and $\dot{\theta}_{\rm fin}=0$.

V.1 A mode excitation catalog

We begin by presenting a catalog of spheroidal modes with $k=2$ for a range of black hole spins, $a\in(0.1M,0.5M,0.9M$), and for four values of orbital inclination, $I\in(20^{\circ},60^{\circ},120^{\circ},160^{\circ})$. The inclinations $20^{\circ}$ and $60^{\circ}$ represent low and high prograde values (i.e., inclinations for which $L_{z}$ is positive); the inclinations $160^{\circ}$ and $120^{\circ}$ represent low and high retrograde values (with $L_{z}<0$).

Figures 7 – 10 show the magnitude of the fundamental mode excitations ${\cal A}_{2m0}$ that we find in these 12 different cases. In each plot, the top row shows results for $a=0.1M$, the middle for $a=0.5M$, and the bottom row for $a=0.9M$. Going from left to right, the columns present data for $m=-2$ through $m=2$. Each point in the panels shows ${\cal A}_{2m0}$ for a different value of $\theta_{\rm fin}$; red dots are for $\dot{\theta}_{\rm fin}>0$, blue are for $\dot{\theta}_{\rm fin}<0$. Additional plots, presenting phases and additional values of $k$, are shown in Appendix D.

Several modes for the high inclination cases we examine ($I=60^{\circ}$ and $I=120^{\circ}$) appear to be affected by numerical noise. We believe this is because in these cases the plunging body rapidly crosses multiple angular grid zones; the same effect led to the relatively large errors we find for the $I=60^{\circ}$ Schwarzschild rotation test described in Sec. IV. The shallow inclination cases ($I=20^{\circ}$ and $160^{\circ}$) are substantially less affected by noise. In Appendix B, we estimate this numerical error by comparing results with a set of higher resolution waveforms — calculated with a higher resolution integration of Teukolsky’s equation. The least reliable mode extractions appear to be ${\cal A}_{k20},{\cal A}^{\prime}_{k-20}$ for $I=60^{\circ}$ and $I=120^{\circ}$. Code enhancements to improve this behavior are under study right now.

The key result we wish to illustrate is that each spheroidal fundamental mode $(k,m,0)$ is excited in a way that depends uniquely and predictably on the parameters $[a,I,\theta_{\rm fin},\operatorname{sgn}(\dot{\theta}_{\rm fin})]$ characterizing its final plunge. Figure 11 shows another view of this, illustrating how mode excitation varies as a function of the angles $(I,\theta_{\rm fin})$ at $a=0.5M$. This figure illustrates how the spectral content of ringdown modes varies, in some cases significantly, as a function of orbit inclination: some modes, such as $(k,m)=(2,0)$ are absent or weak at small inclination, but are very strong for $I$ large; others, such as $(k,m)=(2,1)$ are present at all inclinations, but show large changes in how they depend on $\theta_{\rm fin}$ as $I$ increases.

At least in principle, the clean mapping between source geometry and mode excitation suggests that the inverse problem may be feasible: inferring the properties of the source geometry given knowledge of the excitation of multiple QNMs. We discuss this further in our conclusions.

V.2 Universal mode excitation for shallow inclination?

In assembling this catalog, we have found intriguing trends in how certain modes are excited. A particularly interesting one occurs at low inclination: for $I=20^{\circ}$, we find that modes with the same $k-m$ are excited in largely the same manner, showing a nearly universal functional dependence on $\theta_{\rm fin}$; the same behavior is seen for $I=160^{\circ}$ for modes with the same $k+m$. This behavior is only weakly dependent on spin in the range that we have investigated.

This is illustrated in Fig. 12. Here, we overlay normalized mode amplitudes for several different values of $k$, grouping mode amplitudes with the same black hole spin and the same $k-m$ or $k+m$ values. Amplitudes with given $(a,I,k,m)$ are normalized by averaging over $\cos\theta_{\rm fin}$:

[TABLE]

The range $\theta_{\rm max/min}$ is given by Eq. (17); for $I=20^{\circ}$, $\theta_{\rm max}=110^{\circ}$ and $\theta_{\rm min}=70^{\circ}$. The integral in Eq. (28) is numerically evaluated by interpolating between mode amplitudes calculated at various $\theta_{\rm fin}$.

As Fig. 12 shows, for each value of spin and each value of $k-m$ or $k+m$, a nearly universal functional dependence emerges: In each panel, the normalized amplitude traces out a figure that is nearly the same at that spin for all spheroidal mode indices that have a given value of $k-m$ (for $I=20^{\circ}$) or $k+m$ for ($I=160^{\circ}$). It is noteworthy that the dependence on spin is quite weak: the results at $a=0.1M$ and $a=0.5M$ are nearly identical, and are not dramatically different from the results at $a=0.9M$.

These trends break down at higher inclinations $I=60^{\circ},120^{\circ}$; presumably there is some maximum angle at which this nearly universal excitation form holds (although it should be noted that the highly inclined results are significantly more polluted by numerical noise). As we discuss in the Conclusions, further work may show that it will be possible to exploit this behavior to better understand mode excitation for misaligned plunges.

V.3 Overtones

Another intriguing feature which we noticed in assembling our catalog of mode excitation is, possibly, the presence of overtone modes (i.e., modes with $n\geq 1$) excited by the final plunge and merger. Here we briefly investigate this idea, and present evidence that we may be finding the first overtone in some cases.

Since the overtones are short-lived, amplitudes calculated using the algorithm described in Sec. III.3 do not stabilize to constant values as in the example in Fig. 2. Therefore we conduct a Levenberg-Marquardt nonlinear least squares fit Marquardt (1963) to model the residuals $\delta h^{\rm N}_{22}$ with overtones. We expect this fit to work best for rapidly rotating black holes as excited overtones will decay more slowly. To test this, we examined an equatorial $I=0^{\circ}$ inspiral and plunge for $a=0.99M$ and modeled the early ringdown waves with a single complex frequency $\sigma_{221}$ and amplitude $\mathcal{C}_{221}$. While $t_{0}$ is still defined as the time at which the small body crosses the light ring, we still need to choose a time window over which to fit the model. To choose this the window we examine the following notion of relative fit error:

[TABLE]

where $t_{i}$ is sampled from $t_{\rm min}<t_{i}<t_{0}+70M$. A similar equation defines the relative fit error $\langle\varepsilon_{331}\rangle$.

The right panel of Fig. 13 shows the relative fit error $\langle\varepsilon_{221}\rangle$ for the $I=0^{\circ}$, $a=0.99M$ test case. Not surprisingly, the error is large for $t_{\rm min}-t_{0}$ near zero, since at these times the radiation is not yet fully described by QNMs. The error decreases to a minimum when $t_{\rm min}-t_{0}\gtrsim 30M$; concomitant to this, the fitted amplitude ${\cal A}_{221}$ settles down to a constant. The left panel of Fig. 13 shows the resulting waveform fits, where the top left panel shows the spherical $(2,2)$ mode waveform, along with the fundamental ringdown model constructed with the algorithm discussed in Sec. III.3. The bottom left panel of this figure shows the overtone model.

Encouraged by this test case, we applied this method to the cases in our ringdown catalog. For $I=20^{\circ}$, we construct the overtone model $h^{{\rm RD}1}_{22}$ with two frequencies $\sigma_{221},\sigma_{321}$ and amplitudes $\mathcal{C}_{221},\mathcal{C}_{321}$. The most promising results are for $n=1$ overtones at low inclination and high spin, $(a,I)=(0.9M,20^{\circ})$. In this case, we find that the residuals $\delta h_{22}^{\rm N}$ are well fit by a superposition of one or two overtones. Interestingly, the $n=1$ amplitudes we find are similar in form to the $n=0$ mode excitation. In particular, the modes ${\cal A}_{221}$ and ${\cal A}_{331}$ depend on $\cos\theta_{\rm fin}$ in a manner that is reminiscent of the universal form we discussed in Sec. V.2.

This iterative procedure of calculating overtone amplitudes relies on the assumption that the initial (fundamentals-only) model $\tensor{*}[]{h}{{}_{\ell m}^{\rm RD0}}$ is not biased by contaminating overtones. As a check on this assumption, in Appendix C, we calculate the systematic error incurred when a single overtone is present using toy waveform data. The small systematic errors in the recovered $n=0$ and $n=1$ amplitudes provide further evidence of present overtones plotted in Fig. 14.

VI Conclusion

We have calculated the black hole QNM excitation resulting from a plunging small body whose trajectory is calculated using a generalized Ori-Thorne algorithm. Our method to extract the $n=0$ mode amplitudes does not involve choosing a fitting region, but instead determines when the ringdown model is self-consistent. We find that the mode amplitudes are cleanly parameterized in terms of the black hole spin $a$, the orbital inclination $I$, the small body’s final polar angle $\theta_{\rm fin}$ and its final angular direction $\operatorname{sgn}(\dot{\theta}_{\rm fin})$. We have tabulated the results of our analysis and will provide a Mathematica notebook that can be used to plot mode amplitudes as a function of different combinations of $\bm{(}a,I,\theta_{\rm fin},\operatorname{sgn}(\dot{\theta}_{\rm fin})\bm{)}$ mat . Importantly, we find that this parameterization removes the influence of the final plunge on ad hoc characteristics of the inspiral and plunge model; this is discussed in detail in Appendix A. As long as we use the parameter set $\bm{(}a,I,\theta_{\rm fin},\operatorname{sgn}(\dot{\theta}_{\rm fin})\bm{)}$ to characterize the data, our conclusions about mode excitation appear to be robust.

At least in the large mass ratio limit, our results indicate that there is a clean map between ringdown mode excitation and the properties of the merging binary at plunge. This supports our motivating idea that, by measuring multiple ringdown modes, it may be possible to learn about a binary’s characteristics. By measuring a set of mode amplitudes, a ringdown measurement may provide useful information about the orbit inclination $I$, similar to how Refs. Kamaretsos et al. (2012a, b) demonstrated that ringdown preserves the memory of a binary’s mass ratio and aligned spin. For large-mass systems (which radiate relatively few inspiral cycles in band), and especially as the high- and mid-frequency sensitivity of ground-based detectors is improved, this could significantly increase what gravitational-wave observations can learn about spin-orbit alignment, a property that is particularly important for constraining the formation history of binary black holes. We refer to the reader to Ref. Mandel and Farmer (2018) for overview discussion, Abbott et al. (2016) and references therein for discussion of these observables and source astrophysics, and Gerosa et al. (2018) for recent discussion.

Much work must be done to see whether these measurements can be done in practice. First, we must determine to what extent inferences based on black hole perturbation theory can be trusted. In the large-mass-ratio limit, the background spacetime has a well-defined spin which remains constant during the entire coalescence. When the mass ratio is not so large, the spin of the merged remnant is dominated by the binary’s orbital angular momentum at plunge, and information about the members’ spins will be washed away. In addition, studies must be done to determine how accurately detectors can measure the amplitude of ringdown modes. Although there has been much work studying how well the modes’ frequencies and damping times can be measured, less attention has been paid to date on the modes’ amplitudes (though see Ref. Carullo et al. (2018)). We hope this work will motivate more research in this vein.

Our work has also uncovered interesting properties of mode excitation for misaligned coalescence. An intriguing behavior we discussed in Sec. V.2 is that, at least for shallow inclination angle, the excitation of modes has a nearly universal form: the dependence of the fundamental mode excitation ${\cal A}_{km0}$ on $\cos\theta_{\rm fin}$ and $\operatorname{sgn}(\dot{\theta}_{\rm fin})$ is nearly identical across $k$ and $m$ for fixed $k-m$ (for prograde orbits) or $k+m$ (for retrograde orbits), for a wide range of spin $a$. It may be possible to exploit this relation to make a “quick and dirty” assessment of mode excitation across a wide range of mode values and black hole spins. If this universality includes the Schwarzschild limit, one might even be able to use Eq. (25) in concert with a catalog of Schwarzschild equatorial mode excitations to rapidly estimate the QNM strength for a wide range of physically relevant parameters.

Finally, we see indications that it may be possible to extract information about overtone modes, at least when the black hole spin is large and the plunge inclination is shallow. This is a somewhat delicate operation, so we are cautious about claims in this regime. However, if it can be shown that there is a regime in which the overtones can be reliably understood, this may make possible additional consistency tests and ways of testing the nature of the Kerr metric using ringdown gravitational waves.

Acknowledgments

Our work on this problem was supported at Massachusetts Institute of Technology (MIT) by National Science Foundation Grants No. PHY-1403261 and No. PHY-1707549; H. L. was in addition supported by an MIT Dean of Science Graduate Fellowship. G. K. acknowledges research support from NSF Grant No. PHY-1701284 and Office of Naval Research/Defense University Research Instrumentation Program Grant No. N00014181255. We thank Lionel London for helpful discussions and feedback on this work, Emanuele Berti for helpful discussions and for providing tools for computing ringdown mode properties, Michael Boyle for providing tools which enabled our rotation tests, and Gregory Cook for providing data (based on Ref. Cook and Zalutskiy (2014)) that allowed us to validate our spherical-spheroidal mode mixing coefficients.

Appendix A Robustness of mode excitation to the plunge worldline prescription

As we strongly emphasized in Paper I and re-emphasized in Sec. II.1, the generalized Ori-Thorne algorithm we use to construct the worldline followed by our inspiraling and plunging small body requires us to make three ad hoc choices. Two of these (a time parameter $L_{i}$ which defines when we end “inspiral” and begin “transition,” plus the model we use to smooth the behavior of orbital constants during the transition) turn out to have very small or even negligible impact on the worldline. The impact on the worldline of the third choice, the time $L_{f}$ at which “transition” ends and begin “plunge” begins, is not negligible. Worldlines which start with the same initial conditions but use different choices for $L_{f}$ will end with different values of the final polar angle variable $\theta_{\rm fin}$.

Figure 15 illustrates that the impact of this ad hoc choice on QNM excitation is nugatory. Although changing the parameter $L_{f}$ for fixed initial conditions changes $\theta_{\rm fin}$, it does not affect the manner in which QNM excitation depends on $L_{f}$. In Fig. 15, we show mode excitation for $(a,I)=(0.5M,60^{\circ})$ as a function of $\cos\theta_{\rm fin}$ for the $(k,m,n)=(2,2,0)$ spheroidal mode. Data represented dots, squares, triangles, and diamonds were computed using $L_{f}=2.5$, $2.59$, $2.83$, and $3.07$, respectively. The mode excitation shows no dependence on $L_{f}$, only on the value $\theta_{\rm fin}$ at which the plunge terminates.

We find the same effect for all other modes that we have examined. Changing $L_{f}$ changes the relation between a worldline’s initial conditions and $\theta_{\rm fin}$, but does not affect the fundamental manner in which mode excitation depends on a plunge’s final geometry. Although, as we have emphasized elsewhere, it would extremely useful to eliminate the need for the ad hoc parameters introduced by the generalized Ori-Thorne model, we nonetheless can make robust assessments of how mode excitation behaves as a function of parameters which characterize the final geometry of a coalescing binary.

Appendix B Resolution effects on mode amplitudes

As discussed in Secs. IV.2 and V.1, numerical errors dominate the accuracy of our mode extraction at high orbital inclination. Resolution effects are apparent when, for example, plotting the mode amplitudes and phases as a function of the plunge angle, $\cos(\theta_{\rm fin})$. In Figs. 23 ($I=60^{\circ}$) and 24 ($I=120^{\circ}$), there are small jumps in the $m=2$ mode amplitudes. In contrast, the mode amplitudes for low orbital inclination (Figs. 21 and 22) do not exhibit such jumps.

To further quantify our numerical error, we calculated the mode amplitudes from a set of higher resolution waveforms (Fig. 17). With increased resolution, the low inclination results do not change much, indicating near convergence. However, the high inclination results are affected, indicating larger numerical errors.

Appendix C Estimating systematic errors from a toy model with overtones

At high spin, $a=0.9M,0.99M$, when overtones may be resolved, we employ an iterative process to extract the overtone amplitudes from the numerical data. Initially, we construct a ringdown model with only fundamental modes, $\tensor{*}[]{h}{{}_{\ell m}^{\rm RD0}}$, as described in Sec. III.3. Next, we model any remaining contributions to the ringdown with overtones (see Sec. III.4). However, this iterative process relies on the assumption that the overtones will not bias the initial fundamentals-only model. If this assumption does not hold, then $\tensor{*}[]{h}{{}_{\ell m}^{\rm RD0}}$ may be biased, which can lead to systematic errors in extracting the overtone amplitudes.

To estimate the magnitude of such systematic errors, we analyzed toy waveform data — consisting of only QNMs — with known (injected) mode amplitudes. Since we would like the toy waveform data to resemble our real waveform data, we fix the injected $n=0$ mode amplitudes to some representative values. Here we have adopted the $n=0$ mode amplitudes for $a=0.9M$, $I=0^{\circ}$, which were calculated using the fundamentals-only model. Before $n=1$ overtones are added, the leading order spherical mode $\tensor{*}[]{h}{{}_{22}^{\rm N}}$ comprises QNM contributions with magnitude $|a_{m2k0}(t)\ \mathcal{C}_{k20}|$ (for QNMs with $m\geq 0$) and $|a^{\prime}_{-22k0}(t)\ \mathcal{C}^{\prime}_{k-20}|$ (for QNMs with $m\leq 0$). We then add the $(k,m,n)=(2,2,1)$ overtone, and generate a 11 different waveforms each with various overtone amplitudes in the range $10^{-2}\leq\mathcal{A}_{221}/\mathcal{A}_{220}\leq 10^{2}$, while fixing $\mathcal{A}_{220}$ to $2.35$ (Fig. 18).

Our results indicate that the calculated $n=0$ mode amplitudes are not significantly biased, and can be recovered well within $0.01$ even when the overtone amplitude is large (Fig. 19). The $n=1$ mode amplitude is also faithfully recovered with the iterative procedure, incurring either a relative error under 1% (for a loud overtone) or an absolute error under $0.01$ (for a weak overtone) (Fig. 20). This is primarily because the algorithm extracts the amplitudes at times when the fundamentals-only model is most self consistent with the data; as the injected overtone amplitude increases, the mode extraction algorithm compensates by adjusting the fitting interval to later times in the ringdown (see the caption of Fig. 18). Bias is also minimized by the averaging procedure (i.e. associating $\bar{\mathcal{C}}_{km0}$ and $\bar{\mathcal{C^{\prime}}}_{km0}$ with the mode amplitudes), even when overtone contribution $|a_{2221}(t)\ \mathcal{C}_{221}|$ is exceeds than that of subdominant higher-order modes.

Appendix D Further mode excitation catalog entries

Here we present additional catalog entries describing QNM excitation. In Fig. 16 we show the mode excitation phase $\phi_{2m0}$ for the spins $a=0.1M$, $0.5M$, and $0.9M$ for $I=60^{\circ}$. In Figs. 21–24 and we show the amplitude magnitude ${\cal A}_{km0}$ for $k=3$.