Comparison of post-Newtonian mode amplitudes with numerical relativity simulations of binary black holes

TL;DR

This study compares post-Newtonian and numerical relativity waveforms from binary black hole mergers, finding good agreement in many modes even during highly dynamical phases, which aids in improving gravitational waveform models.

Contribution

It demonstrates that leading order post-Newtonian amplitudes closely match numerical relativity results for many modes during strong-field black hole mergers.

Findings

Post-Newtonian amplitudes agree well with numerical relativity for many modes.

Modes with l = m and l = 2, m = 1 are least modified from post-Newtonian predictions.

Understanding these modes helps improve models of gravitational waveforms in strong gravity regimes.

Abstract

Gravitational waves from the coalescence of two black holes carry the signature of the strong field dynamics of binary black holes. In this work we have used numerical relativity simulations and post-Newtonian theory to investigate this dynamics. Post-Newtonian theory is a low-velocity expansion that assumes the companion bodies to be point-particles, while numerical relativity treats black holes as extended objects with horizons and fully captures their dynamics. There is a priori no reason for the waveforms computed using these disparate methods to agree with each other, especially at late times when the black holes move close to the speed of light. We find, remarkably, that the leading order amplitudes in post-Newtonian theory agree well with the full general relativity solution for a large set of spherical harmonic modes, even in the most dynamical part of the binary evolution, with only some modes showing distinctly different behavior than that found by numerical relativity simulations. In particular, modes with spherical harmonic indices l = m as well as l = 2, m = 1 are least modified from their dominant post-Newtonian behavior. Understanding the nature of these modes in terms of the post-Newtonian description will aid in formulating better models of the emitted waveforms in the strong field regime of the dynamics.

Figures

Tables

Table 1. Table 1: Correlation coefficients of the non-spinning fits to the SXS data for spherical harmonic modes (2,2), (2,1), (3,3), (3,2), (3,1), (4,4), (4,3), (4,2), (4,1), (5,5), (6,6), (7,7), and (8,8) at τ M = − 100 , − 10 , 0 , 10 𝜏 𝑀 100 10 0 10 \frac{\tau}{M}=-100,-10,0,10 .

Mode	$\tau =-100M$	$\tau =-10M$	$\tau =0M$	$\tau =10M$
(2,2)	0.999998	0.999779	0.999703	0.999798
(2,1)	0.999579	0.999854	0.999838	0.998717
(3,3)	0.999912	0.999787	0.999301	0.997751
(4,4)	0.999667	0.999299	0.998743	0.991130
(5,5)	0.999842	0.999204	0.997655	0.993545
(6,6)	0.999722	0.998421	0.994847	0.976633
(7,7)	0.999698	0.996487	0.982980	0.939223
(8,8)	0.999491	0.994833	0.965470	0.848470
(3,2)	0.999376	0.996502	0.972295	0.585903
(3,1)	0.997824	0.993977	0.981189	0.908449
(4,3)	0.999719	0.999168	0.998515	0.974125
(4,2)	0.998844	0.997680	0.995446	0.908354
(4,1)	0.984923	0.948784	0.863388	0.976998

Table 2. Table 2: Correlation coefficients of the aligned spinning fits to the SXS data for spherical harmonic modes (2,1), (3,2), (4,3), and (4,1) and mass ratios q = 1 , 1.5 , 2 , 3 , 7 𝑞 1 1.5 2 3 7 q=1,1.5,2,3,7 at τ M = − 100 , − 10 , 0 , 10 𝜏 𝑀 100 10 0 10 \frac{\tau}{M}=-100,-10,0,10 .

Mode	$q$	$\tau =-100M$	$\tau =-10M$	$\tau =0M$	$\tau =10M$
(2,1)	1	0.999325	0.999773	0.999988	0.999619
	1.5	0.999980	0.999931	0.999807	0.999291
	2	0.997846	0.995953	0.991763	0.984452
	3	0.997914	0.994403	0.991027	0.985538
	7	0.999152	0.998667	0.999101	0.998829
(3,2)	1	0.990234	0.981084	0.968417	0.897599
	2	0.980074	0.960144	0.914401	0.581929
	3	0.964765	0.952854	0.938971	0.784182
	7	0.993511	0.953335	0.056356	0.980378
(4,3)	1	0.998952	0.999922	0.998798	0.995412
	1.5	0.999844	0.999494	0.998659	0.874023
	2	0.961789	0.958457	0.960024	0.472553
	3	0.844141	0.843969	0.837170	0.222383
	7	0.987371	0.938732	0.759405	0.621659
(4,1)	1	0.991669	0.995539	0.985739	0.988915
	1.5	0.608697	0.623753	0.916280	0.954856
	2	0.684519	0.041862	0.220282	0.144261
	3	0.485186	0.233015	0.002791	0.334364
	7	0.926836	0.992870	0.981711	0.982951

Table 3. Table 3: List of 43 SXS simulations for initially non-spinning binary black holes, showing the SXS simulation ID, the mass ratio q 𝑞 q , and the numerical resolutions used for our analysis.

SXS Id	$q$	Resolutions
$2$	$1.00$	$4$, $5$, $6$
$180$	$1.00$	$2$, $3$, $4$
$198$	$1.20$	$3$, $4$, $5$
$7$	$1.50$	$4$, $5$
$8$	$1.50$	$4$, $5$
$194$	$1.52$	$2$, $3$
$169$	$2.00$	$3$, $4$, $5$
$184$	$2.00$	$2$, $3$, $4$
$201$	$2.32$	$1$, $2$, $3$
$259$	$2.50$	$3$, $4$, $5$
$191$	$2.51$	$2$, $3$
$30$	$3.00$	$3$, $4$, $5$
$168$	$3.00$	$3$, $4$, $5$
$183$	$3.00$	$2$, $3$, $4$
$200$	$3.27$	$1$, $2$, $3$
$193$	$3.50$	$2$, $3$
$294$	$3.50$	$3$, $4$
$182$	$4.00$	$2$, $3$, $4$
$190$	$4.50$	$2$, $3$
$54$	$5.00$	$3$, $4$, $5$
$56$	$5.00$	$3$, $4$, $5$
$107$	$5.00$	$3$, $4$, $5$

Table 4. Table 4: List of 121 SXS simulations for aligned-spin binary black holes, showing the SXS simulation ID, the mass ratio q 𝑞 q , the spins represented by χ 1 , 2 subscript 𝜒 1 2 \chi_{1,2} via χ → 1 , 2 = χ 1 , 2 ​ L ^ subscript → 𝜒 1 2 subscript 𝜒 1 2 ^ 𝐿 \vec{\chi}_{1,2}=\chi_{1,2}\hat{L} , and the numerical resolutions used for our analysis.

Id	$q$	${\chi }_1$	${\chi }_2$	Res.
$4$	$1.00$	$-0.50$	$0.00$	$5$, $6$
$5$	$1.00$	$0.50$	$0.00$	$4$, $5$
$148$	$1.00$	$-0.44$	$-0.44$	$4$, $5$
$149$	$1.00$	$-0.20$	$-0.20$	$3$, $4$
$150$	$1.00$	$0.20$	$0.20$	$3$, $4$
$151$	$1.00$	$-0.60$	$-0.60$	$3$, $4$
$152$	$1.00$	$0.60$	$0.60$	$3$, $4$
$154$	$1.00$	$-0.80$	$-0.80$	$3$, $4$
$155$	$1.00$	$0.80$	$0.80$	$3$, $4$
$156$	$1.00$	$-0.95$	$-0.95$	$4$, $5$
$157$	$1.00$	$0.95$	$0.95$	$3$, $4$
$158$	$1.00$	$0.97$	$0.97$	$5$, $6$
$159$	$1.00$	$-0.90$	$-0.90$	$3$, $4$
$160$	$1.00$	$0.90$	$0.90$	$3$, $4$
$170$	$1.00$	$0.44$	$0.44$	$5$, $6$
$171$	$1.00$	$-0.44$	$-0.44$	$5$, $6$
$172$	$1.00$	$0.98$	$0.98$	$3$, $4$
$175$	$1.00$	$0.75$	$0.75$	$2$, $3$
$176$	$1.00$	$0.96$	$0.96$	$3$, $4$
$177$	$1.00$	$0.99$	$0.99$	$3$, $4$
$178$	$1.00$	$0.99$	$0.99$	$4$, $5$
$209$	$1.00$	$-0.90$	$-0.50$	$4$, $5$
$210$	$1.00$	$-0.90$	$0.00$	$4$, $5$
$211$	$1.00$	$-0.90$	$0.90$	$4$, $5$
$212$	$1.00$	$-0.80$	$-0.80$	$4$, $5$
$213$	$1.00$	$-0.80$	$0.80$	$4$, $5$
$214$	$1.00$	$-0.62$	$-0.25$	$4$, $5$
$215$	$1.00$	$-0.60$	$-0.60$	$4$, $5$
$216$	$1.00$	$-0.60$	$0.00$	$4$, $5$
$217$	$1.00$	$-0.60$	$0.60$	$4$, $5$
$218$	$1.00$	$-0.50$	$0.50$	$4$, $5$

Equations

$$h_{+}-i{h}_{\times }=\sum _{\ell =2}^{\infty }\sum _{m=-\ell }^{\ell }{h}_{\ell m}{}_{-2}{Y}_{\ell m}(\theta ,\varphi )\text{,}$$

$$h_{22}$$

$$h_{21}$$

$$h_{33}$$

$$h_{32}$$

$$h_{31}$$

$$h_{44}$$

$$h_{43}$$

$$h_{42}$$

$$h_{41}$$

$$h_{55}$$

$$h_{66}$$

$$h_{77}$$

$$h_{88}$$

$$A_{22}$$

$${\hat{A}}_{21}$$

$${\hat{A}}_{33}$$

$${\hat{A}}_{32}$$

$${\hat{A}}_{31}$$

$${\hat{A}}_{44}$$

$${\hat{A}}_{43}$$

$${\hat{A}}_{42}$$

$${\hat{A}}_{41}$$

$${\hat{A}}_{55}$$

$${\hat{A}}_{66}$$

$${\hat{A}}_{77}$$

$${\hat{A}}_{88}$$

$$h_{21}=C(v,\psi )\eta \left(\delta v-\frac{3}{2}\left({\mathbf{\chi }}_{\mathbf{a}}+\delta {\mathbf{\chi }}_{\mathbf{s}}\right)\cdot {\hat{\mathbf{L}}}_{\mathbf{N}}{v}^2\right)+O(v^3)\text{,}$$

$${\mathbf{\chi }}_{\mathbf{s}}$$

$${\mathbf{\chi }}_{\mathbf{a}}$$

$$h_{21}=C(v,\psi )\eta \left(\delta v-\frac{3}{2}{\chi }_{21}^{\text{eff}}{v}^2\right)+O(v^3)\text{,}$$

$${\hat{A}}_{21}={\gamma }_{21}{\hat{A}}_{21}^{\text{ns}}+{\beta }_{21}{\chi }_{21}^{\text{eff}}\text{.}$$

$${\hat{A}}_{\ell m}={\gamma }_{\ell m}{\hat{A}}_{\ell m}^{\text{ns}}(\eta )+{\beta }_{\ell m}{\chi }_{\ell m}^{\text{eff}}(\eta ,{\chi }_1,{\chi }_2)\text{,}$$

$${\chi }_{21}^{\text{eff}}$$

$${\chi }_{32}^{\text{eff}}$$

$${\chi }_{43}^{\text{eff}}={\chi }_{41}^{\text{eff}}$$

$$A_{\ell m}=\sqrt{\mathrm{Re}(h_{\ell m}{)}^2+\mathrm{Im}(h_{\ell m}{)}^2},$$

Full text

I Introduction and motivation

The Laser Interferometer Gravitational-Wave Observatory (LIGO) at two sites in the USA (Hanford, WA and Livingston, LA) and the Virgo detector in Pisa, Italy, have opened a new era in multi-messenger astronomy and fundamental physics via the discovery of binary black hole Abbott et al. (2016a, b, 2017a, 2017b, 2017c) and binary neutron star Abbott et al. (2017d) mergers. These discoveries have, for the first time, enabled tests of dynamical gravity in the strongly dissipative regime of the theory Abbott et al. (2016c, d, 2017a, 2017c), i.e. the period derivative of the binary $\dot{P}$ changes very rapidly during the time of observation (see Refs. Will (2006); Sathyaprakash and Schutz (2009); Yunes and Siemens (2013); Gair et al. (2013)). This is in contrast to the Hulse-Taylor binary Hulse and Taylor (1975) where the change in period $\dot{P}$ is essentially constant.

Radio measurements of the rate at which the orbital period decays in a binary neutron star allowed spectacular confirmation of the quadrupole formula Taylor et al. (1979); Taylor and Weisberg (1982); Lyne et al. (2004); Kramer et al. (2004). However, radio binary pulsars probe the weak field sector of the two-body dynamics111We note that the self-gravity of the neutron stars, which must be taken into account in the measurement of the various binary parameters, are large. Indeed, the compactness of neutron stars given by the dimensionless quantity ${\cal C}\equiv GM_{\rm NS}/c^{2}R_{\rm NS},$ where $M_{\rm NS}$ is the mass of the neutron star and $R_{\rm NS}$ its radius, is about ${\cal C}\sim 0.2.$ In this sense, the observations do probe the strong field regime of general relativity; however, the two-body dynamics is governed by weak fields., wherein the dimensionless gravitational potential $\phi$ of one of the bodies on the other is $\phi\ll 1,$ or, equivalently, the speed $v$ obeys $v/c\sim\sqrt{\phi}\ll 1.$ In contrast, gravitational wave observations of the merger make it possible to test general relativity when $\phi\sim 0.5$ (the largest it ever gets) and the system is strongly dissipative. Consequently, LIGO, Virgo, and other future ground-based gravitational wave detectors (KAGRA and LIGO-India) can test the validity of general relativity in an entirely new regime of the theory.

I.1 Modeling binary black hole dynamics

The dynamics of a binary black hole consists of three phases: inspiral, merger, and ringdown. Inspiral refers to the early phase of the binary evolution when the effect of radiation reaction on the orbital motion is small. The slow-motion, weak-field dynamics of this phase, when the two black holes are far apart, is well-modeled by post-Newtonian (PN) theory (see Ref. Blanchet (2006) for a review) where all the observables are expressed as a power series in $v/c$. The strong field dynamics close to the merger and the dynamics of the highly deformed remnant black hole can only be modeled using numerical relativity, where one solves Einstein’s equation for the two-body problem using numerical techniques (see Ref. Pretorius (2007) for a review). The ringdown phase of the dynamics occurs when the remnant black hole has become less deformed and can be well-approximated as a perturbation of a Kerr black hole and modeled using black hole perturbation theory (see Ref. Sasaki and Tagoshi (2003) for a review).

The waveform emitted by an inspiralling compact binary predominantly consists of the quadrupole mode. It was pointed out that controlling the evolution of the orbital phase of the dominant mode was far more important Cutler et al. (1993) for the detection problem than controlling the correction to its amplitude or the inclusion of higher order modes that contain wave frequencies other than twice the orbital frequency of the quadrupole mode. In this so-called restricted post-Newtonian approximation one neglects the correction to the amplitude of the waveform arising from higher order multipoles. However, higher modes are critical for an unbiased estimation of both the intrinsic parameters of a binary (e.g. companion masses and spins) but also the orientation of the binary relative to a detector and its position on the sky Van Den Broeck (2006); Van Den Broeck and Sengupta (2007); Arun et al. (2007a, b, 2009a). Waveforms based on this new information could be useful to test general relativity in the high-curvature regime of the theory, which is one of the principal goals of gravitational wave astronomy.

Understanding the structure of the multipole modes and their dependence on the intrinsic parameters of the binary will be critical in building more refined waveform models and using gravitational wave observations to test predictions of general relativity in dynamical spacetimes.

I.2 Imprints of progenitors on the black hole ringdown spectrum: Past

studies

Previous studies Kamaretsos et al. (2012a, b) found that the properties of the progenitor system such as the symmetric mass ratio $\eta=m_{1}m_{2}/(m_{1}+m_{2})^{2}$ (where $m_{1,2}$ denote the binary component masses) or mass ratio-weighted combinations of component spins (referred to as “effective spin” parameter) leave their unique imprints on various ringdown modes of the remnant black hole. More specifically, Ref. Kamaretsos et al. (2012a) had found that for mergers of nonspinning black holes, the amplitudes of the four strongest modes, (2,2), (2,1), (3,3), (4,4), continue to depend on mass ratio even $15M$ after the luminosity of the (2,2) mode peaks. In a follow-up work, Ref. (Kamaretsos et al., 2012b) further investigated binaries whose component spins are aligned with the orbital angular momentum. The authors found that the relative amplitude of the (2,1) mode $10M$ after the peak luminosity could be captured by a fitting formula with two variables: $\eta$ and an effective spin parameter.

While the fitting function for the (2,2) mode in Ref. Kamaretsos et al. (2012a) was based on physical intuition gained from PN theory, the ansatz for the other modes were based on fits to the numerical data. Similarly, Kamaretsos et al. (2012b) uses a functional form for the (2,1) mode inspired by PN theory including not only symmetric mass ratio, but also spin dependencies. Achieving a good fit required an effective spin combination that was slightly different from the one found in PN approximation (Mishra et al., 2016). These results pointed to the interesting possibility of inferring the properties of the progenitor black holes just from the late ringdown signal.

Following a different approach London et al. (2014), London et al. studied the $\eta$-dependence of the higher modes of the post-merger amplitudes. Their study is based on fitting the amplitudes of higher modes from numerical relativity simulations to high-order polynomials in the symmetric mass ratio. Such fits are useful in building analytical models of the post-merger waveforms. Indeed, in a more recent study  London (2018), London et al. developed a new phenomenological waveform model that includes higher modes.

I.3 Present work

In the present work we use a combination of numerical relativity simulations and PN theory to study the evolution of different modes of gravitational waves as a function of time, mass ratio and mode-dependent “effective spin” parameters (see Eq. 10). The phase evolution of each mode, being a multiple of the orbital phase, is essentially the same for all modes and has been amply treated in the literature; we, therefore, restrict our study entirely to the mode amplitudes. The two gravitational wave polarisations $h_{+}$ and $h_{\times}$ from an inspiralling binary are, in principle, composed of infinitely many modes as exemplified by the relation Thorne (1980)

[TABLE]

where ${}_{-2}Y_{\ell m}$ are the $-2$ spin-weighted spherical harmonics, $(\theta,\phi)$ define the direction of propagation of the wave, and $h_{\ell m}$ are the spherical harmonic wave modes. Although the quadrupole $(\ell,m)=(2,2)$ is the dominant mode, higher order modes can have comparable, albeit smaller, amplitudes relative to the quadrupole when the component masses are very different or the compact objects have significant spin. In the inspiral regime, using a cocktail of approximation schemes, PN theory provides an effective framework to relate the radiative multipoles observed at infinity to the source multipoles Blanchet (2006); Blanchet et al. (1995), thereby expressing the observed gravitational waveform in terms of the multipole moments of the source.

Kelly and Baker Kelly and Baker (2013) investigated the effects of mode mixing between the spherical harmonics, used in numerical relativity and PN theory, and spheroidal harmonics. The latter capture the axial symmetry of the Kerr spacetime and hence are a more suitable basis to describe perturbations of the Kerr metric during the ringdown. In particular, they showed that the $(3,2)$ spherical harmonic mode has significant contributions from different spheroidal harmonic modes, which is referred to as mode mixing. A spheroidal harmonic decomposition renders the modes to fall off more smoothly as a function of time, thus allowing a more simplified modeling of the waveform. While this is true, we provide an alternative interpretation of mode mixing as arising due to the failure of the point particle description of PN theory close to the formation of a common horizon.

This formation marks a rather intriguing transition in the binary black hole dynamics: from the perturbative dynamics of the two black holes to the perturbative dynamics of the remnant black hole via this highly non-perturbative merger. This transition is captured in the full general relativity solution, provided by numerical relativity simulations that by design track the dynamics of black hole horizons. Hence, given the availability of numerical relativity catalogs for binary black hole mergers, it is interesting to ask how the information about the two black holes, encoded in PN expressions, propagates from the inspiral to the merger and ringdown phases and whether with just observing the latter two, one can infer the properties of the binary components. Here, extending the works of Kamaretsos et al. (2012a, b), we compare several spherical harmonic mode amplitudes from SXS numerical simulations SXS with the leading terms in the corresponding PN expressions; allowing for one free parameter in the nonspinning case and two free parameters in the spinning case. Our aim is to search for those modes that are fitted very well by the aforementioned PN-based fits, and hence retain information about the progenitor system.

We find that the most dominant mode amplitudes hardly change their dependence on the symmetric mass ratio, given from PN theory, throughout the evolution of the binary. Further, we find that the signature of the strong field regime is encoded in a small number of modes that are sub-dominant, with their amplitude being less than 10% of that of the quadrupole.

The paper is organized as follows: Section II describes the SXS numerical simulations we employ for the study. The fitting model we use, based on PN expressions for leading order spherical harmonic modes of the waveform, is explained in Sec. III. Our results on the PN signatures in the spherical harmonic modes of numerical relativity are described in Sec. IV and the implications of these results for modelling waveforms from binary black holes are discussed in Sec. V. Appendix A provides an alternative representation of some rmesults in Section IV. Some of the technical details of the simulations are elaborated in appendix B.

II Numerical simulations

This study utilizes publicly available binary black hole gravitational waveforms from the SXS collaboration SXS . The concrete simulations used are listed in Appendix B. Specifically, Table 3 lists the 43 non-spinning simulations that were used, while Table 4 lists the 121 aligned-spin simulations. The simulations were originally presented as follows:

• •

The first SXS waveform catalog Mroue et al. (2013) ($1\leq$ SXS id $\leq 174$).

• •

Simulations for developing techniques for very high black hole spins Lovelace et al. (2015); Scheel et al. (2015) ($175\leq$ SXS id $\leq 178$).

• •

Simulations for a waveform surrogate model for non-spinning binary black hole systems Blackman et al. (2015) $180\leq$ SXS id $\leq 201$).

• •

Binary black hole simulations at mass-ratio 7 with particularly many inspiral cycles Kumar et al. (2015) ($202\leq$ SXS id $\leq 207$).

• •

A study of aligned spin binary black hole systems Chu et al. (2016); Kumar et al. (2016) ($209\leq$ SXS id $\leq 304$).

The simulations were computed with the Spectral Einstein Code (SpEC) SXS Collaboration , a multi-domain pseudo-spectral code designed to solve elliptic and hyperbolic partial differential equations, in particular the Einstein equations. SpEC computes initial data with the extended conformal thin sandwich method York (1999); Pfeiffer and York (2003) utilizing quasi-equilibrium black hole excision boundary conditions Cook and Pfeiffer (2004); Caudill et al. (2006); Lovelace et al. (2008) and iterative eccentricity reduction Buonanno et al. (2011) to achieve quasi-circular inspirals. SpEC evolves the Generalized Harmonic form of Einstein’s equations Friedrich (1985); Pretorius (2005) in first order form Lindblom et al. (2006) with constraint damping Gundlach et al. (2005); Pretorius (2005); Lindblom et al. (2006) and constraint preserving boundary conditions Lindblom et al. (2006); Rinne (2006); Rinne et al. (2007). The code uses black hole excision Scheel et al. (2009); Szilagyi et al. (2009); Hemberger et al. (2013), coupled with a dual-frame approach to have the computational grid track the motion of the black hole horizons Scheel et al. (2006). The gravitational wave data used in our study was extracted with Regge-Wheeler-Zerilli wave-extraction Sarbach and Tiglio (2001); Regge and Wheeler (1957); Zerilli (1970), and corrected for time-dilation effects at the extraction radius Boyle and Mroue (2009); Taylor et al. (2013) and for mode-mixing arising from small residual motion of the center of mass Boyle (2016). More technical details are given in the original publication presenting the simulations Mroue et al. (2013); Lovelace et al. (2015); Scheel et al. (2015); Blackman et al. (2015); Kumar et al. (2015); Chu et al. (2016); Kumar et al. (2016).

SpEC simulations are generally run at multiple numerical resolutions, in order to be able to assess numerical convergence and numerical truncation error. Indeed, we have restricted the present study only to simulations that are available at multiple resolutions. The last column in Tables 3 and 4 lists the resolutions of each simulations that were used. For each simulation, the accuracy increases with a larger numerical value in this column. However, because of improvements to SpEC’s numerical algorithms in the course of time, it is not possible to assign an absolute meaning to these resolution values. Using the different numerial resolutions, we compute an error bar for every numerical value extracted from the numerical relativity data based on the difference in this value when extracted from the numerical relativity data of different resolution.

Visual inspection of the $(\ell,m)$ modes indicate that the $(2,2)$ mode is well-behaved for the simulations considered here. The leading sub-dominant modes with $m\leq 4$, specifically (2,1), (3,3), (3,2), (3,1), (4,4), (4,3), (4,2), (4,1), are also generally well-behaved with only rare visible unphysical features, like for instance unexpected extraneous oscillations during ringdown. Further, the numerical errors of these modes, see Figs. 2 and 4, indicate good numerical convergence of the considered runs.

Modes with higher frequency, (5,5), (6,6), (7,7), (8,8), unfortunately, appear often compromised during merger and ringdown. The most common symptom is that these modes reach their maximum a few $M$ earlier than expected, and do not exhibit a clear exponential decay thereafter. These symptoms are consistent with an insufficiently fine radial grid, on which the short-wavelength high-frequency merger waves would not be resolved well enough as they propagate to the extraction spheres, and are thus unphysically damped away. Moreover, extrapolation to infinite extraction radius appears to magnify non-physical features in these high-frequency modes, in about half of the simulations considered. Gravitational wave extrapolation is most important for the early inspiral, where the wavelength is long Boyle et al. (2007); Boyle and Mroue (2009) and is less important for the merger portion considered here. To mitigate impact on the high-$m$ modes –(5,5) and above–, we therefore decided in the present study to utilize the gravitational waveforms extracted at the largest available extraction radius.

The impact on our analysis of imperfections in the underlying numerical data can be judged in two ways: First, Figs. 2, 4, and 5 show error bars for each data-point, obtained from the difference in value from numerical simulations of different numerical resolution. These error bars are generally small compared to the physical effects being analysed. Second, our analysis should produce results that are slowly and smoothly varying with change of the underlying binary black hole parameters like mass-ratio or black hole spin. In particular, simulations at nearby parameter points should yield similar answers, and indeed they do, even if the simulations come from very different epochs. The results obtained here (e.g. in Fig. 2) vary smoothly with parameters, with the scatter being consistent with the error bars. As such, we believe the underlying numerical data to be reliable for our purposes, except perhaps, for the analyses of $(5,5)$ through $(8,8)$ in the regime after the $(2,2)$ mode reached peak amplitude.

III Leading order post-Newtonian approximations of gravitational wave mode

amplitudes

The ‘plus’ and ‘cross’ gravitational wave polarizations can be decomposed in terms of spherical harmonics as shown in Eq. (1). The complete PN expressions for various spherical harmonic modes $h_{\ell m}$, given the currently available accuracies of the multipole moments, are reported in Refs. Blanchet et al. (2008) and Arun et al. (2009b); Mishra et al. (2016), for nonspinning binaries moving in circular orbits and for systems whose spins are aligned or anti-aligned with respect to the orbital angular momentum, respectively.

As we are going to crucially exploit the leading order dependencies of these modes, we list them below for convenience. Note that in these expressions $v$ is the PN velocity parameter and $\eta$ the symmetric mass ratio defined as $\eta=\frac{m_{1}m_{2}}{M^{2}}$, with binary component masses $m_{1},m_{2}$ and total mass $M$.

The structure of the various modes in PN theory reads as Blanchet et al. (2008)

[TABLE]

where $\mathcal{C}_{\ell m}$ are complex constants, $v$ is the PN velocity parameter which captures the time dependency of the wave modes, $\psi$ is the PN phase variable, and $\delta=\frac{m_{1}-m_{2}}{m_{1}+m_{2}}$ is a mass asymmetry parameter which can be rewritten as $\delta=\sqrt{1-4\,\eta}$ for $m_{1}>m_{2}$. It vanishes for equal mass binaries.

Based on the structure of the expressions in Eqs. (2) we introduce the leading order PN approximations which capture the leading order $\eta$ and spin dependencies for fixed $v$—i.e. at a fixed time—and thus allow us to examine the numerical relativity waveforms for PN signature or rather deviations from it. Our goal is somewhat diffent from the usual approach in the literature as we are aiming to study the behavior of the mode amplitudes in terms of the intrinsic parameters of the binary system and not as a function of time.

III.1 Nonspinning binaries

In order to gain insight into the behavior of the amplitudes of the nonspinning modes, we choose the following fitting functions $A_{\ell m}=|h_{\ell m}|$ which capture the leading order dependencies of the PN expressions (2) on the mass ratio parameters $\eta$ and $\delta$:

[TABLE]

where $\alpha_{\ell m}$ are the scaling factors that we fit for. The hatted amplitudes $\hat{A}_{\ell m}={A_{\ell m}}/{A_{22}}$ have been normalized with respect to the (2,2) mode to cancel the overall $\eta$-factor present in every mode.

III.2 Aligned spin binaries

In PN theory, spin effects are sub-dominant and are not present at leading order for any mode Mishra et al. (2016). Current-multipole modes which obey $\ell+m=\text{odd}$ contain spin-dependent terms at 0.5 PN order above the leading term and thus are more likely to exhibit spin effects Kamaretsos et al. (2012b). We focus on the four current-multipole modes with $\ell\leq 4$, $(2,1)$, $(3,2)$, $(4,3)$, and $(4,1)$. The PN expression for the (2,1) mode to the next-to-leading order in $v$ is given by

[TABLE]

where $C$ is a function of the orbital velocity $v$ and the PN phase variable $\psi$, $\bm{L_{N}}$ is the orbital angular momentum, and

[TABLE]

denote, respectively, the symmetric and antisymmetric spin combinations of the initial black hole spins $\bm{\chi_{1}}$ and $\bm{\chi_{2}}$. Since we assume that spins and the orbital angular momentum are aligned, we can write instead

[TABLE]

where $\chi^{\text{eff}}_{21}=\chi_{a}+\delta\,\chi_{s}$, with $\chi_{a,s}=\bm{\chi}_{a,s}\cdot\bm{\hat{L}_{N}}$ being the projection of the symmetric/antisymmetric spin vectors along the orbital angular momentum. This form motivates the fitting ansatz, with non-spinning $\hat{A}_{21}^{\text{ns}}=\delta$ from (3), which reads as

[TABLE]

This can be generalized to arbitrary $\ell m$ as

[TABLE]

with different effective spin parameters for different modes defined by the linear combination of the spin parameters in the PN expressions for those modes. The functional forms of the effective spin parameters for the different modes are motivated by Eqs. (12) of Ref. Arun et al. (2009b) and are given by

[TABLE]

Equation (9) has two fit parameters $\beta_{\ell m},\gamma_{\ell m}$ whereas Eqs. (3) only need one. The additional parameter is aimed to capture the extra degrees of freedom due to spins and account for the fact that the nonspinning and spinning effects enter at different PN orders.

IV Post-Newtonian signature in numerical relativity waveform amplitudes

IV.1 Numerical wave mode amplitudes

From each of the numerical relativity simulations summarized in Sec. II and detailed in Appendix B, we extract the time series of real and imaginary components of the spherical harmonic wave modes $h_{\ell m}$ (see Eq. (1)) for the 13 modes with $(\ell,m)\in\{$(2,2), (2,1), (3,3), (3,2), (3,1), (4,4), (4,3), (4,2), (4,1), (5,5), (6,6), (7,7), (8,8)$\}$. The remaining numerical modes with $\ell\geq 5,m\neq\ell$ are excluded from this study due to their relatively small amplitudes and large numerical errors. Fig. 1 shows the real and imaginary components of six of the 13 wave modes together with their amplitude for an example SXS run: SXS:BBH:0169, mass ratio $q=2$, and non-spinning.

Our goal is to examine the behavior of various spherical harmonic modes for PN signature and deviations from it around the time of merger. We concentrate this study on the evolution of the real amplitude

[TABLE]

of the spherical harmonic wave modes $h_{\ell m}$, where both $A_{\ell m}$ and $h_{\ell m}$ are functions of time $\tau$, as well as binary parameters $M,\eta,\vec{\chi}_{1},\vec{\chi}_{2}$. The time variable $\tau=t-t_{\text{max}}^{22}$ has been shifted such that the peak amplitude of the (2,2) mode is located at $\tau=0$. In the case of non-spinning binary black holes with quasi-circular orbits the parameter tuple $\theta=(M,\eta,\vec{\chi}_{1},\vec{\chi}_{2})$ reduces to the two mass parameters, the total mass $M$, and the symmetric mass ratio $\eta$. If the black holes are aligned spinning, $\theta=(M,\eta,\chi_{1},\chi_{2})$ also contains the spin magnitudes.

IV.2 Results for non-spinning binaries

Figure 2 contains the condensed results of our study of PN signature in the waveform amplitudes of initially non-spinning binary black holes around the time of merger. The individual symbols in these plots are the amplitudes of the various gravitational wave modes extracted from the SXS simulations. Each symbol carries an error bar, often too small to be visible, that is derived as the difference in the extracted amplitude for at least two different numerical resolutions. The solid lines in Fig. 2 represent our leading order PN approximations in Eq. (3) fitted to the numerical data for various wave mode amplitudes as a function of the symmetric mass ratio. The temporal evolution is presented via snapshots at four different times ${\tau}/{M}=-100,-10,0,10$, corresponding to the columns in Fig. 2. The three rows group the different modes by numerical strength. The first row contains the data and fits for the four strongest modes (2,2), (2,1), (3,3), and (4,4), the second row shows the remaining $\ell=m$ modes, and the sub-dominant modes with $\ell=3,4$ are bundled in the last row. The larger subfigures plot relative amplitudes $\hat{A}_{\ell m}(\tau)={A_{\ell m}(\tau)}/{A_{22}(\tau)}$, $\ell m\neq 22$, with only $A_{22}$ being shown as an absolute amplitude. The error bars correspond to twice the numerical errors $\sigma_{\ell m}$ shown in the smaller subplots.

The amplitude of the (2,2) mode behaves as expected and increases towards its maximum at $\tau=0$. Due to the suppression of the two next strongest modes, (3,3) and (2,1), for equal mass binaries, the (2,2) mode stays most significant in the realm between $q=1$ to $q=2$ which is where all detections by the LIGO Virgo Collaboration were made Abbott et al. (2016d, 2017a, 2017c, 2017b). The situation for the other modes paints a more interesting picture for low symmetric mass ratios $\eta<0.15$ where their amplitudes increase more quickly relative to the (2,2) mode. This tendency shows the importance of the inclusion of higher modes for medium to extreme mass ratio binary coalescences.

The first column in Fig. 2 shows the comparison of the leading order PN approximations to the numerical data at time $\tau=-100M$. The approximations work beautifully and confirm the expectation that PN theory describes the functional dependence of the gravitational wave amplitudes on $\eta$ very well during the inspiral. The level of agreement between the data and the fits is quantified by the correlation coefficients222The correlation coefficient $C$ between the data vector $\bm{d}$ and the appropriate fit vector $\bm{f}$, with averages $\bar{d},\bar{f}$, is defined as $C=\frac{(\bm{f}-\bar{f})\cdot(\bm{d}-\bar{d})}{\sqrt{(\bm{f}-\bar{f})^{2}(\bm{d}-\bar{d})^{2}}}$. in Table 1. The situation stays very similar close to merger at $\tau=-10M$, even though the amplitude of the (3,2) mode is starting to show deviations from the PN inspired fitting. The fits for the remaining modes capture the data extremely well despite the common belief that PN theory should fail in this regime due to the increase in the orbital velocity parameter $v$.

The picture becomes truly exciting at $\tau=0$, after a common horizon has already formed. The amplitudes of the four dominant modes are fitted exceptionally well by the leading order PN approximations, hence giving us an insight into how little these amplitudes are affected by the dynamics during the coalescence of the binary system. The PN-inspired fits to the $\ell=m$ modes in the second row of Fig. 2 are still remarkably well captured, especially for $\ell=5,6$. The amplitudes of the modes with $\ell=3,4$, $m<\ell$ exhibit a different behavior: their numerical amplitudes deviate strongly from the PN-inspired fitting and thus indicating that the merger process affects the dynamics of these mode amplitudes more than the four dominant modes or the modes with $\ell=m$.

Finally, the last column of Fig. 2 contains the data and fits during the early ringdown at $\tau=10M$. The deviations from the leading order PN approximations have increased, compared to time $\tau=0$ which is reflected in the correlation coefficients in Table 1. The amplitude data for $\tau=10M$ in the second row of Fig. 2 appears to be well captured by the leading order PN approximations, but it exhibits large numerical errors that make a quantitative evaluation of the approximations impractical, see Table 1. The dominant modes show a very intriguing outcome. Their amplitudes seem to maintain the PN signature from earlier times fairly well. This reproduces the earlier findings (Kamaretsos et al., 2012a) that found $\eta$-dependences in these amplitudes during the ringdown. Our analysis goes beyond that and shows that this dependence is still mostly of PN signature $10M$ after the merger.

In summary, we can say that the four dominant modes with large amplitudes and the wave modes with $\ell=m$ maintain the PN signature of the inspiral phase exceptionally well in their amplitudes, from the inspiral throughout the merger into the ringdown, while the spherical harmonic wave modes with $\ell\neq m$ deviate from this PN-like behavior as the evolution of the binary approaches the merger. Thus, these present interesting candidates for binary black hole merger studies and strong field tests of general relativity, with the (3,2) mode being especially intriguing as it is the strongest of these modes and hence the most significant for future detections.

IV.3 Results for aligned spins

Let us now discuss the aligned-spin simulations. We extract amplitudes $A_{22}$ and $\hat{A}_{\ell m}$ as above for all simulations with aligned spins. Fig. 3 plots the amplitudes of all 13 modes for each aligned spinning SXS simulations in a distinct data point. For a given value of $\eta$, there are generally multiple simulations with different spins; these simulations lead to different amplitudes, resulting in the vertical scatter of data-points at each $\eta$. Most modes show a small amount of vertical scatter whereas the (2,1) and (3,2) mode exhibit significant variation amongst the different simulations at a given symmetric mass ratio. This large amount of spread in the amplitudes of the aforementioned modes is a consequence of a strong dependence on the omitted two parameters, the spin magnitudes $\chi_{1}$ and $\chi_{2}$ of the component black holes.

Not only do the (2,1) and (3,2) modes show the most pronounced scatter in Fig. 3, but they are also among the current-multipole modes with $\ell+m=\text{odd}$, which exhibit spin effects at low PN order (cf. Sec. III.2). We will therefore now investigate the spin dependence of the relevant four modes (2,1), (3,2), (4,3), and (4,1) in more detail. The results and corresponding correlation coefficients for the amplitude data and the PN-inspired fits for the (2,1), (3,2), (4,3), and (4,1) modes are presented in Fig. 4 and Table 2, respectively. While Fig. 4 spreads the temporal evolution of the mode amplitudes for given mass ratios $q$ in columns, we give another representation of the same information in Fig. 5 in Appendix A, where the roles of mass ratio $q$ and time $\tau$ are flipped. This gives a better intuition of how each mode amplitude evolves for fixed mass ratio throughout the merger.

The (2,1) modes appears to behave similarly to the non-spinning case and maintain the PN signature of the inspiral throughout the merger into the early ringdown. The leading order PN approximations capture the numerical relativity data exceptionally well for mass ratios $q=1,1.5,7$. Mass ratios $q=2,3$ show a larger scatter around the linear fit line which appears to be a result of the larger variation of initial spin combinations of the SXS waveforms for lower mass ratios, see Table 4 (in particular, the SXS simulations explore variations of anti-symmetric spin $\chi_{1}-\chi_{2}$ much more exhaustively for $q=1,2,3$). This scatter is not visible in the data for mass ratio $q=1$ due to the suppression of asymmetries for equal mass binaries.

The results are similar for the (4,3) and (4,1) modes. Mass ratio $q=1$ is again extremely well captured by the leading order PN approximations. While mass ratios $q=1.5,7$ appear to be fitted well in Fig. 4, the correlation coefficients show that the PN inspired fit performs less admirably than for the (2,1) mode. The scatter around the approximation lines for $q=2,3$ is also much more prevalent, especially at times $\tau/M=0,10$.

Our analysis of the (3,2) mode does not include mass ratio $q=1.5$: Two of the three distinct initial spin combinations reduce to the same effective spin $\chi^{\text{eff}}_{32}=\eta\,\chi_{s}$ due to its symmetry and thus make a linear fit to two points a moot exercise. The analysis of the (3,2) mode shows that it takes a special role amongst the four studied modes when the system includes spins. The (3,2) mode amplitude is the only to show major deviations from the PN signature for equal mass binary waveforms, $q=1$, and even during the inspiral at $\tau=-100M$. The picture for mass ratios $q=2,3,7$ is the same as what we observed for the (4,3) and (4,1) modes.

In summary, the (2,1) mode seems to do as well as we saw from the non-spinning scenario, while the three weaker modes exhibit various different effects. All modes, even (2,1), showed some level of scatter for mass ratios $q=2,3$ that probably is a result of variations in the initial spin data, which did not show up for equal mass binaries for the three modes with odd $m$, (2,1), (4,3), and (4,1). Hence, it appears that the addition of aligned spin has a stabilizing effect against deviations from a PN signature for these modes, if the mass ratio is $q=1$. The (3,2) mode takes a distinct role as its amplitude for mass ratio $q=1$ shows deviations from the leading order PN approximation already during the inspiral. Hence, it again stands out as the mode of interest in strong field tests of general relativity.

V Implications of the results for modelling waveforms from binary black holes

In this Section we will discuss the behavior of the various modes as a function of time and the system’s mass ratio $\eta$ and mode-dependent ‘effective spin’ (which is different for different modes).

Time evolution

Our study has shown that for most part of the adiabatic evolution when $\dot{\omega}/\omega^{2}\ll 1,$ where $\omega$ is the orbital frequency, the mass-ratio dependence of the amplitudes of the various spherical harmonic modes are as predicted by PN theory to lowest order. The overall multipole structure is set in when the two black holes are well-separated; it is difficult to deform this multipolar structure because spacetime has a large bulk modulus. Remarkably, most $\ell=m$ mode amplitudes continue to agree with the PN prediction well after the common horizon has formed. This includes the dominant $\ell=m=2,3,4$ modes as well as the weaker $\ell=m>4$ modes (cf. Fig. 2, first and second row). The strong field dynamics does affect the $\ell\neq m$ modes (except the strong $(2,1)$ mode), especially when the two bodies get closer together as evidenced by the change in the weaker $\ell\neq m$ modes (cf. Fig. 2, last row). Our analysis shows that as we get close to the merger phase, say $\tau\sim-10M$ (this is about when the common apparent horizon forms) the moments begin to deform from their PN behavior.

$\eta$-dependence

Figures 2 and 4 show the behavior of the mode amplitudes as a fucntion of symmetric mass ratio $\eta$ at different epochs and as a function of ‘effective spin’ for different mass ratios $q$ and epochs, respectively. For nonspinning systems, the $\ell=m$ modes are in pretty good agreement with the leading order PN behavior as a function of $\eta$ (see Eqs. 3). This is true both at earlier times $\tau\sim 100\,M$ when PN equations are expected to provide a good description of the mode amplitudes, as well as at epochs when the PN equations are believed not to be accurate. In fact, even at the onset of merger at $\tau\simeq-10\,M$ and beyond $\tau=0$ when the black hole begins to settle down (i.e. $\tau\sim 10\,M$) $\ell=m$ modes show little departure from the PN behavior.

However, the weaker $\ell\neq m$ modes are altered significantly already at the onset of the merger ($\tau\sim-10\,M$), especially for comparable mass binaries (i.e. $\eta\simeq 1/4$). One exception to this rule is the $\ell=2,$ $m=1$ mode. This mode is the strongest sub-dominant mode after $\ell=m=3$ (see Fig. 2) and is not easily modified by the strong field dynamics. The amplitude of the other $\ell\neq m$ modes are at the level of $\lesssim 8\%$ (for highly asymmetric systems) of the $(2,2)$ mode amplitude, while the $\ell=2,$ $m=1$ mode could be as large as 30% of the overall amplitude.

Spin-dependence

Figures 4 and 5 present the mode amplitudes of the four $\ell+m=\text{odd}$ modes (2,1), (3,2), (4,3), and (4,1) as a function of their respective effective spins (see Eq. 10) for different mass ratios and epochs. Again, the PN approximation does remarkably well at capturing the behavior of the numerical data for the strong (2,1) mode, for all epochs and mass ratios. However, the three weaker modes whose amplitudes show deviations from the PN behavior for non-spinning systems as early as $\tau\sim-10\,M$, agree with the PN approximation to some extent when the system includes aligned spins. Hence, it appears that the addition of spin to the system has a stabilizing effect on the PN signature.

That being said, the situation is much more complicated than for non-spinning systems as the quality of the agreement with the PN signature depends not only on the mode and epoch, but also the mass ratio and the sample spread of initial black hole spins. The agreement is good for mass ratios $q=1.5,7$ for which the data shows the linear behavior in the respective effective spin combination for all modes. Mass ratios $q=2,3$ were sampled with a much larger distribution in the initial spins (cf. Table 4), resulting in an envelope of data points around the linear PN approximation in $\chi^{\text{eff}}_{\ell m}$. These envelopes widen during later epochs and for weaker modes. The found spread in the data hints that the relationship between initial black hole spins and mode amplitudes during the merger and ringdown cannot be captured in one effective spin combination. Finally, the equal mass systems were sampled with a similarly large spread in the intial spins, but do not show the envelope characteristics of the higher mass ratios which is a result of the symmetry in the system. The amplitudes of the three odd-$m$ modes (2,1), (4,3), and (4,1) are wonderfully captured by the linear PN approximations at all epochs.

The major exception is presented by the (3,2) mode’s amplitude for equal mass systems which shows a curious, but definite nonlinear dependence on its effective spin at all studied epochs. This behavior is curious for two reasons: The PN signature is the strongest for the other three modes at mass ratio $q=1$ and all epochs (cf. Fig. 5, first column). Further, it is the only case (i.e. the only mode for both the aligned and non-spinning simulations) where the PN approximation seems to already fail at $\tau=-100\,M$. It shows that the (3,2) modes takes a special place amongst the $\ell+m=\text{odd}$ modes with $\ell\leq 4$. This is captured by its even azimuthal number $m=2$ whereas the other three are odd $m$ modes. Thus, the (3,2) mode is the most interesting mode amongst all the weaker $\ell\neq m$ modes to study deviations from the PN signature: It is the strongest of these modes and thus the easiest to detect, it does not vanish for non-spinning, equal mass systems, and it is affected by spin effects where it can capture departures from PN theory well into the inspiral-regime.

VI Conclusions

In this paper we have provided a comparison of the amplitudes of spherical harmonic modes of gravitational waves from merging binary black holes computed using the leading order PN approximation with those obtained from numerical relativity simulations.

The post-Newtonian approximation is based on the point-particle description of the two-body problem in general relativity. It is a good approximation when the two bodies are far from each other (i.e., their distance of separation $r$ is far greater than the scale of the horizon $R_{s}\sim 2GM/c^{2}$ of the component masses), but expected to breakdown when the two bodies are close to coalescence $r\sim\mbox{few}\times 2GM/c^{2}.$ While the post-Newtonian approximation is now known to a high order in the expansion parameter $v/c=\sqrt{GM/c^{2}r},$ it is not expected to capture the strong field dynamics of the theory close to merger.

Numerical relativity simulations, on the contrary, are exact solutions to Einstein’s equations for the two-body problem. They capture the strong field dynamics, including the dynamics of the common horizon and how that horizon approaches the final Kerr state. While these simulations can, in principle, resolve the full spectrum of modes emitted by the binary, in practice finite resolution and numerical accuracy limit the number of modes that can be extracted reliably to the quadrupole, octupole, and hexadecapole modes, corresponding to spherical harmonic index of $\ell=2,3,$and $4,$ respectively.

The chief finding of our study is that the dependencies of these dominant mode amplitudes on the symmetric mass ratio and the binary’s spins, computed in the leading order post-Newtonian approximation, agree remarkably well with those extracted from numerical relativity simulations, deep into the regime where the approximation should not have worked. In particular, the quadrupole modes $(2,2)$ and $(2,1)$, extracted from numerical relativity simulations, show little departure from the leading order post-Newtonian expression throughout the inspiral and merger. This is also true for the $(3,3)$ and the $(4,4)$ modes. This implies that the dominant multipole structure of the system remains frozen as determined by the point-particle approximation. All the modes begin to show significant departure from post-Newtonian description in the quasi-normal mode regime, $\sim 10GM/c^{3}$ after the waveform reaches its peak amplitude.

The weaker modes with $\ell\!=\!3,4$, $m\!\neq\!\ell$ modes also agree with the leading order post-Newtonian expressions when the system is $\sim 100GM/c^{3}$ away from coalescence, but begin to show significant departure from the leading order post-Newtonian behavior well-before the epoch when the waveform reaches its peak amplitude. In other words, the point-particle approximation of post-Newtonian theory to the leading order is no longer adequate in describing the behavior of the amplitude of these modes. It is for this reason that we conclude that these weaker modes are affected far more by the strong field regime of the binary evolution than the stronger modes $(2,2),(2,1),(3,3)$ and $(4,4)$.

It is well known that the $(3,2)$ spherical harmonic mode is a mixture of several spheroidal harmonic modes, which causes it to decay non-monotonically in the ringdown regime of the signal Kelly and Baker (2013); London et al. (2014); London (2018). While this is true, the new insight from our study is that we can exploit the leading order post-Newtonian expressions in any analytical modeling of the mode amplitudes. Furthermore, we believe that understanding the multipole structure of the common horizon could provide further insight into why certain modes are affected far more by the strong field dynamics than others.

Acknowledgements.

We thank Abhay Ashtekar, Michael Boyle, Mark Hannam, Lionel London, Sean McWilliams and Leo Stein for helpful comments on the manuscript. KGA and BSS acknowledge the support by the Indo-US Science and Technology Forum through the Indo-US Centre for the Exploration of Extreme Gravity, grant IUSSTF/JC-029/2016. SB and BSS are supported in part by NSF grants PHY-1836779, AST-1716394 and AST-1708146. KGA is partially support by a grant from Infosys Foundation. KGA also acknowledge partial support by the grant EMR/2016/005594. Computing resources for this project were provided by The Pennsylvania State University. This document has LIGO preprint number LIGO-P1800367.

Appendix A Alternative representation of the results for aligned spin fits

Fig. 5 presents the same information as Fig. 4, but with the roles of the time $\tau$ and mass ratio $q$ flipped in the figure. This presentation allows a more streamlined look at how each mode behaves as function of time for a given mass ratio, thus making very evident, how strongly the (3,2) modes deviates from the leading order post-Newtonian approximation for mass ratio $q=1$.

Appendix B Numerical relativity simulations from the SXS project

For each numerical resolution, the SXS waveform catalog provides a metadata file with information about the specifics of the run as well as the gravitational waveforms decomposed into spherical harmonics for both the Newman-Penrose scalar $\Psi_{4}$ and the gravitational wave strain $h$. Our analysis focuses on the latter and was conducted with the data contained in the files ‘rhOverM_Asymptotic_GeometricUnits_CoM.h5’, which provide the spherical harmonic modes of $h$, at the outermost extraction radius, and extrapolated to asymptotic null infinity. Furthermore, the data in these files are corrected for mode mixing that can arise if initial transients during start of the evolution induce a motion of the center of mass of the binary black hole Boyle (2016); Ossokine et al. (2015). The retarded time-coordinate is corrected for gravitational redshift effects Boyle and Mroue (2009).

These HDF5 files structure the data into four groups containing the same 77 datasets, but for different extrapolation orders $N=2,3,4$, as well as the outermost extraction radius. The datasets store the simulation output as a time series of real and imaginary components of the coefficients $h_{\ell m}$ in the expansion in spherical harmonics with spin-weight $s=-2$ (1) of the gravitational wave strain $h$ for all 77 modes with $\ell=2,\dots,8$, $m=-\ell,\dots,\ell$. In order to put errors on the numerical data we restricted our analysis to the 43 non-spinning and 121 aligned spinning simulations that are provided at a minimum of two different resolution levels, see Tables 3 and 4. Further, we restrict our analysis to the outermost extraction radius which yields the most accurate numerical results for the merger and ringdown regimes.