On choosing the start time of binary black hole ringdown
Abstract
The final stage of a binary black hole merger is ringdown, in which the system is described by a Kerr black hole with quasinormal mode perturbations. It is far from straightforward to identify the time at which the ringdown begins. Yet determining this time is important for precision tests of the general theory of relativity that compare an observed signal with quasinormal mode descriptions of the ringdown, such as tests of the nohair theorem. We present an algorithmic method to analyze the choice of ringdown start time in the observed waveform. This method is based on determining how close the strong field is to a Kerr black hole (Kerrness). Using numerical relativity simulations, we characterize the Kerrness of the strongfield region close to the black hole using a set of local, gaugeinvariant geometric and algebraic conditions that measure local isometry to Kerr. We produce a map that associates each time in the gravitational waveform with a value of each of these Kerrness measures; this map is produced by following outgoing null characteristics from the strong and nearfield regions to the wave zone. We perform this analysis on a numerical relativity simulation with parameters consistent with GW150914 the first gravitational wave detection. We find that the choice of ringdown start time of after merger used in the GW150914 study Abbott et al. (2016a) to test general relativity corresponds to a high dimensionless perturbation amplitude of in the strongfield region. This suggests that in higher signaltonoise detections, one would need to start analyzing the signal at a later time for studies that depend on the validity of black hole perturbation theory.

I Introduction
The quasinormal mode (QNM) spectrum seen during the ringdown of a perturbed black hole (BH) is determined by the Teukolsky equation; it carries the signature of the BH potential along with the BH horizon and asymptotic boundary conditions Teukolsky (1973); Press and Teukolsky (1973); Teukolsky and Press (1974). The recent detections of binary black hole (BBH) gravitational wave (GW) signals by LIGO (the Laser Interferometer GravitationalWave Observatory) Abbott et al. (2016b, c); Abbott et al. (2016d); Abbott et al. (2017a, b) allow us to begin to probe this QNM signature Abbott et al. (2016a). The QNM spectrum in a gravitationalwave observation allows us to perform tests of the nohair theorem. This theorem states that vacuum, asymptotically flat, stationary, axisymmetric, uncharged BHs are completely characterized by two parameters: the mass and the spin Mazur (2000); Misner et al. (1973); Dreyer et al. (2004); Gossan et al. (2012); Kamaretsos et al. (2012a). This allows us to constrain modified theories of gravity that violate the nohair theorem Berti et al. (2015); Yunes et al. (2016). Observing the QNM spectrum in GWs can be used to validate the BH uniqueness theorem. This theorem states that the exterior geometry of an vacuum, asymptotically flat, stationary, axisymmetric, uncharged BH must be Kerr Mazur (2000); Mars (1999).
However, testing the nohair and uniqueness theorems relies on observing GWs from the QNM perturbative regime (without additional transients remaining from the inspiral). This requires an appropriate choice of start time of this regime.^{1}^{1}1While conventions in the literature vary, in this paper, by ‘ringdown’, we explicitly mean the part of the postmerger gravitational waveform that can be described in terms of QNMs. Identifying this time in the signal is mathematically an illdefined problem, since QNMs form an incomplete and nonorthogonal basis Nollert (1996); Tufts and Kumaresan (1982). Hence, the conventions for choosing the start time of the ringdown have varied in the literature. Berti et al. Berti et al. (2007) and Baibhav et al. Baibhav et al. (2017) chose the start time based on maximizing the energy contained in the QNM. London et al. London et al. (2014) used after the peak of the dominant mode of (the NewmanPenrose scalar that encodes outgoing radiation) for fitting to NR waveforms.^{2}^{2}2 Since vacuum GR is a scaleinvariant theory, it is convenient to express distance and time in terms of source mass by setting . Explicitly, , where is the mass of the BH. Kamaretsos et al. Kamaretsos et al. (2012b) chose after the peak luminosity of the dominant mode of the waveform, while Thrane et al. Thrane et al. (2017) proposed a loudnessdependent start time. In the GW150914 testing general relativity (GR) paper Abbott et al. (2016a), different start times were used to perform the QNM analysis shown in Fig. 5 of that paper, and the results were consistent with GR when the start time was picked as (or later) after the merger.
None of these methods use information from the strong field to motivate the start times. The strong field refers to the region near the BHs (typically within a radius of few ), where the scale of the curvature is much smaller than the wavelength of a gravitational wave. In this paper, we develop an algorithmic method for validating choices of the start time of ringdown using strongfield features. Specifically, we measure the Kerrness, or closeness to Kerr, in the strongfield region of an NR simulation ringdown, and use null characteristics to map Kerrness onto the GW at asymptotic future null infinity, . We then demonstrate this method on a GW150914like system.
Determining Kerrness in the strongfield regime is nontrivial, since one needs a coordinateinvariant way of identifying a metric as Kerr. Necessary and sufficient conditions for a gaugeinvariant characterization of local isometry to a Kerr manifold were proposed by GarcíaParrado GómezLobo in GómezLobo (2016).^{3}^{3}3Throughout this text, isometry refers to the smooth mapping of manifolds equipped with metrics. We use this set of algebraic and geometric conditions to provide a numerical measure of Kerrness. Previous studies have used multipole moments of the BH apparent horizon Owen (2009), horizon spin measurement comparisons Scheel et al. (2009), or Petrov classification Baker and Campanelli (2000); Campanelli et al. (2009); Owen (2010) to characterize ringdown spacetimes. Our work is the first set of conditions that completely characterizes a spacetime as isometric to a Kerr manifold. We study the violation of these conditions postmerger in the strong field of a BBH simulation.
Connecting the strongfield region to the wave zone is a challenge, as the simulation gauge is different from the gauge in which GWs are observed. There is no straightforward way to transform between these gauges. Furthermore, establishing simultaneity between events is not possible in the GR framework, and thus there is no direct map between an event in the strongfield region and a point on the waveform. We therefore devise a scheme to approximately associate the two frames. The association used in this study is of a causeeffect nature: we follow the outgoing null characteristics from the strongfield region to the wave zone using a Cauchy Characteristic Extraction scheme (CCE) Handmer and Szilágyi (2015); Handmer et al. (2015, 2016), and associate events in the strong field to the wave zone. However, given that GR is a nonlinear theory, the source associated with a particular feature in the GW signal may not be well localized in the spacetime. Nevertheless, one would expect that the source dynamics that dominantly contribute to certain features in the waveform be localizable to a certain extent. Several such approximate localizations have been performed in linear perturbation theory Price et al. (2016); Cardoso et al. (2016).
This paper is organized as follows. Sec. II presents the theoretical methods used in this paper, and Sec. III discusses their implementation in NR simulations. Sec. IV then presents and discusses the results of applying these methods to an NR simulation with GW150914like parameters. We conclude in Sec. V. Figs. 15 and 22 are the flagship figures, presenting our major results. The the results are quantitatively summarized in Table 3.
Conventions
We work with the standard 3+1 decomposition of NR (cf. Baumgarte and Shapiro (2010) for an introduction). In this paper, refers to the spacetime metric, refers to the timelike unit normal vector, refers to the spatial metric on each slice, is the covariant derivative with respect to , and refers to the extrinsic curvature. We set and express all quantities in terms of , the sum of the Christodoulou Masses of the two BHs at the start of the simulation. Latin letters at the start of the alphabet, , refer to (4dimensional) spacetime indices, while Latin letters in the middle of the alphabet, are (3dimensional) spatial indices. We denote complex conjugation by an overbar (e.g. ). To avoid confusion among the multiple meanings of ‘data’ in this paper, we refer to the vacuum data on a spatial slice simply as ‘a slice’.^{4}^{4}4Vacuum data means that the spatial metric, , and the extrinsic curvature satisfy a set of constraint equations corresponding to the decomposition of the vacuum Einstein equations. Similarly, rather than being purely geometric, a ‘slicing’ in our case is a foliation equipped with a coordinate chart.
Ii Theory
ii.1 Characterizing strongfield Kerrness
First, we explain our method of measuring Kerrness in the strongfield region and develop a method to map it onto . Secs. II.1.1 and II.1.2 discuss theoretically characterizing Kerrness in the strongfield region, while Secs. II.2.1, II.2.2, and II.2.3 discuss mapping strongfield information onto the wave zone via null characteristics.
ii.1.1 Overview and historical background
Our overall goal in this section is to evaluate Kerrness: how close a numerical BH ringdown spacetime is to being locally isometric to the Kerr spacetime. In order to evaluate the Kerrness of a spacetime, we first need a set of theoretical conditions to evaluate whether a spacetime is isometric to Kerr. We can then turn these conditions into a set of measures, where deviation from zero indicates being farther from being locally isometric to Kerr. In a numerical simulation, one would evaluate these measures on spatial slices of a simulation. To characterize Kerrness in the strongfield region, one needs local quantifiers evaluated close to the BH, as opposed to looking at regions far away which are contaminated by gravitational radiation. Consequently, we seek a pointwise measure and do not use global measures on a slice such as those proposed in Backdahl and Valiente Kroon (2010a, b); Backdahl and Valiente Kroon (2012).
Uniquely characterizing a spacetime as Kerr has been historically challenging—until recently one could only classify spacetimes up to a Petrov type, which gives a weaker classification that admits several manifolds besides Kerr. The Petrov classification uses algebraic properties of the Weyl tensor based on the four principal null directions (PNDs), by solving the eigenbivector problem (cf. Stephani et al. (2009) for a review)
(1) 
where eigenbivectors have eigenvalues . The degeneracies of the PNDs give a unique algebraic classification of a spacetime. A spacetime with no repeated PNDs is fully general (Petrov Type I). A spacetime with at least one repeated PND is algebraically special. The Kerr metric belongs to a particular class of algebraically special spacetimes, the set of type D spacetimes, which have two double PNDs. A necessary condition for the manifold to be locally isometric to Kerr is to be type D.
Campanelli et al. Campanelli et al. (2009) used this approach to analyze a numerical BBH ringdown. They determined the degeneracies between the PNDs by solving the eigenbivector problem and measuring the difference between eigenvalues. Their analysis found that the spacetime first numerically settled to type II, which has only one double PND, and then to type D. Owen Owen (2010) later showed that this measure was sensitive to the choice of tetrad used to compute the Weyl scalars needed to solve the characteristic equation. He proposed a new measure, lesssensitive to the choice of tetrad, and showed that the spacetime settled to type D without first settling to type II.
A type D spacetime can then be shown to be locally isometric to Kerr through additional conditions. Kerr belongs to the KerrNUT subset of type D spacetimes. One needs to show that a spacetime is KerrNUT and then constrain the acceleration and the NUT parameters. We give more information on KerrNUT spacetimes and the various parameters in Appendix B. Ref. Campanelli et al. (2009) investigated the asymptotic behavior of the acceleration and the NUT parameter on a BBH simulation and showed they were constrained to be those of Kerr.
In this study, we do not solve the eigenbivector problem, but rather use a set of local algebraic and geometric conditions recently proposed by GarcíaParrado GómezLobo GómezLobo (2016) to show that a spacetime is locally isometric to Kerr. These conditions are formulated in a fully covariant way and thus avoid the complications in Campanelli et al. (2009) and Owen (2010) due to tetrad choice.
ii.1.2 Necessary and sufficient Kerrness conditions
To characterize a spatial Cauchy slice as isometric to Kerr, we first check if the slice is algebraically special. Next, we use two geometric conditions to check for the existence of Killing vectors (KVs) on the slice, and we impose two algebraic conditions to verify that the slice containing the KVs is type D. Then, we check the properties of the KVs and further classify the slice into the KerrNUT subfamily. Finally, imposing conditions on the acceleration and NUT parameters, we completely characterize the slice as locally isometric to Kerr. These conditions are summarized in Fig. 1.
All algebraic conditions are expressed in terms of electric and magnetic parts of the Weyl tensor,
(2)  
(3) 
where the left dual of the Weyl tensor is defined as . For a vacuum spacetime, these spatial tensors can be more readily evaluated on a slice as
(4)  
(5) 
where is the spatial Ricci tensor evaluated from . These can be combined into a complex quantity as
(6) 
In GómezLobo (2016), the algebraic condition for a slice to be locally algebraically special is given in Eq. 85 as
(7) 
where
This condition is equivalent to the speciality index in the Petrov classification literature (cf. Eq. 4.13 of Stephani et al. (2009)).
Recall that algebraic speciality corresponds to having one double PND, and hence is a weaker condition than being type D, which corresponds to having two double PNDs. A necessary algebraic condition for a slice to be type D is given in Theorem 4 of GómezLobo (2016) as
Type D 1  (8) 
which makes use of 4dimensional algebraic conditions proven in Ferrando et al. (2001) and orthogonally splits these onto the spatial slice. Here we have called the condition ‘Type D 1’ purely for bookkeeping purposes, in order to label each of the type D conditions.
The three sufficient conditions for a slice to be type D consist of two geometric conditions involving KVs and one algebraic condition which also includes the KV. As proven in Theorem 2 of GómezLobo (2016), a vacuum type D spacetime has a complex KV field which satisfies an algebraic condition
(9) 
where is derived from the Weyl tensor, and
(10) 
However, one must show that a KV field exists on the slice in the first place, and then that it satisfies the properties given in Eq. (9). The necessary and sufficient geometric conditions for a slice to contain a KV field are known as Killing Initial Data (KID), and for a vector , are given as
Type D 3  (11)  
Type D 4  (12)  
Satisfying these conditions guarantees that a KV field exists on the slice—note that these two conditions say nothing so far about type D.
We can then relate this KV field to the condition on the KV in a type D spacetime given in Eq. (9) by requiring
Type D 2  (13)  
where Eq. (13) is the orthogonal splitting of Eq. (9), and
(14)  
This procedure is shown in Theorem 6 of GómezLobo (2016).^{5}^{5}5The Type D 2 condition has a in the second term where GómezLobo (2016) has a . The sign error has been confirmed by the author of GómezLobo (2016). Similarly, The factor of in the definition of is not included in GómezLobo (2016), but is in the corresponding Mathematica notebook GómezLobo ().
Type D 3 and Type D 4 are independent geometric conditions that depend on the complex KV and show that the slice is KID. Type D 1 is a purely algebraic condition that informs us of the behavior of the PNDs. Type D 2 ties in the algebraic and geometric conditions, thereby completing the classification into type D. Speciality Index, meanwhile, is an independent algebraic condition.
In order to then show that an algebraically special, type D slice is locally isometric to Kerr, we must also show that it belongs to the KerrNUT subset of type D spacetimes. KerrNUT spacetimes have the symmetry property of two commuting KVs Stephani et al. (2009)  one spacelike and timelike, and thus if we impose this geometric condition on KV as defined above, we arrive at the condition given in Theorem 8 of GómezLobo (2016),^{6}^{6}6However, this has a typographical error (confirmed by the author GómezLobo ()), and should include , the complex conjugate, as given Eq. (15).
Kerr 1  (15) 
In order to further show that a slice is locally isometric to Kerr, we must place constraints on the parameters characterizing KerrNUT spacetimes. We summarize the parameters involved in Type D spacetimes in Appendix B. We require that , the NUT parameter, vanish, and , which is related to the acceleration of the BH, be greater than zero. These conditions are given in Theorem 8 of GómezLobo (2016) as
Kerr 2  (16) 
for the condition , where , and
Kerr 3  (17) 
for . However, the above expression only holds outside of the ergoregion GómezLobo () in Kerr. This condition is thus impractical to use in the this study, since it involves finding the ergoregion, and masking this region would introduce high levels of numerical error within a spectral code.
Thus, for a slice to be locally isometric to Kerr, it must satisfy all of the above conditions, which are summarized in Fig. 1. Since the vacuum spacetime at the start of a ringdown may be fully general, the left hand sides of the Kerrness conditions will not necessarily be zero on some slices. Instead, the Kerrness conditions turn into a set of Kerrness measures, where larger deviation from zero indicates a larger deviation from being isometric to Kerr.
ii.2 Connecting strongfield information to
ii.2.1 Motivation
Having characterized the Kerrness in the strongfield region, we connect this information to the GWs at . We develop a framework to map the evolution of the Kerrness measures computed during a postmerger simulation to the evolution of the postmerger waveform in the asymptotic frame. This provides a procedure to validate the choices of start time of ringdown when analyzing a gravitationalwave signal.
Just after the two BHs merge, the newly formed BH is expected to be highly distorted. The dynamics of the BH can be explained only via a full numerical simulation. At , where the GWs are observed, these strongfield dynamics are responsible for features in a small region close to the peak of the GW amplitude. Once the excitation amplitude in the strongfield region decays to a level when linear perturbation theory is valid the spacetime dynamics and the associated waveform is governed by the Teukolsky equation Teukolsky (1973); Press and Teukolsky (1973); Teukolsky and Press (1974). At , the waveform appears as a sum of exponentially damped sinusoids with a specific QNM frequency spectrum (with powerlaw tails that are usually very weak). Beyond this rough picture, the association of the specifics in the strongfield dynamics to the waveform is not well understood, especially during the merger and postmerger phases.
Understanding this association is crucial because several strongfield tests of GR rely on BH perturbation theory and thus, on identifying the perturbative regime in the waveform. These tests include the nohair theorem test, consistency tests of the QNM spectrum with the inspiral parameters, and the area theorem test. The start of ringdown in the GW is mathematically illdefined as damped sinusoids form an incomplete and nonorthogonal basis Nollert (1996); Tufts and Kumaresan (1982). Therefore, it is important that we validate the choices of start times in the data analysis of ringdown guided by the strongfield information, where the validity of perturbation theory can be better understood.
ii.2.2 Conceptual challenges
Mathematically, GR being a nonlinear theory does not allow for unambiguous localization of sources of GWs. However, to a certain extent, one expects that the dominant source of a particular feature in the wave zone be localizable to a relatively small region of the spacetime in the past light cone. For instance, studies like Price et al. (2016); Khanna and Price (2017) identify the dominant source for the peak of the waveform during the plunge of a test particle into a Schwarzschild BH with the particle crossing the lightring.^{7}^{7}7The lightring is the orbit of a massless particle around the BH, which corresponds to the peak of the BH potential located at in BoyerLindquist coordinates for a Schwarzschild BH. Furthermore, the last few cycles of the BBH GW signal are associated with the dynamics of a linearly perturbed BH Vishveshwara (1970). However, one needs to bear in mind that these studies are performed using linear perturbation theory where such localizations are better defined. For example, if one adds a massive particle instead of a test particle in the former case and makes the problem nonlinear, one would get some additional source contributions from selfforce, thus making the source localization trickier.
In the case of BBH postmerger, identifying specific events as a source of the features in the waveform cannot be done unambiguously owing to the nonlinear dynamics from merger. However, drawing intuition from analytical linear perturbation theory, we expect the region within the support of the analytical effective BH potential to contribute significantly to the waveform at . Thus, we argue that even in a nonlinear case, a small region in the spacetime around the BH containing the strongfield dynamics, can be associated as a dominant source of features in the GW.
Another challenge in performing this association is that the notion of simultaneity in GR is not absolute, which means that all spacelike slicings of the spacetime are equally valid. In numerical simulations however, we have to make a gauge choice. In our case this choice is made by the Cauchy evolution code. The spatial features corresponding to a particular timeslice are gauge dependent. We choose to monitor the Kerrness on a spatial coordinate 2sphere in the strongfield region, instead of computing a volume integral over the source region in a timeslice.^{8}^{8}8By doing so, the gauge effect is limited to uncertainty of picking the 2sphere, thereby avoiding contribution of gauge effects through the entire volume region.
Further, note that this association is only meant to be approximate, and is dependent on the radius of the 2sphere we monitor. Roughly speaking, the size of the 2sphere indicates the error bar in the association that originates from the gauge choice (specifically, error radius of extraction) .
ii.2.3 Forming a sourceeffect association via null characteristics
Given these challenges, we propose the following association scheme. By looking at evaluated on the simulation, we infer a 2sphere radius that lies within the strongfield region, containing and generating significant radiative fields. This 2sphere acts like an effective source for the GW seen at . We evaluate a surface integral of the Kerrness measures at each time slice during the ringdown on this 2sphere. Then, we connect the evolution of the Kerrness measures on this surface to the associated features in the GW by following the outgoing null characteristics emanating from this 2sphere. The details of this procedure are described below.
The GWs emanating from a source propagate to along outgoing null rays (since the spacetime is curved, a small portion of GWs also travel inside the light cone). We utilize this in constructing an association between strongfield information and the features on the GW. We associate a feature on the GW to a 2sphere in the strongfield region at a given time (in the simulation coordinates) if they lie on the same outgoing null hypersurface. This 2sphere can thus be interpreted as an effective source producing the point on the waveform. Hence, the choice of 2sphere should be close to the region generating GWs rather than farther out. Measuring Kerrness of such a surface would give an insight into validity of perturbation theory in the region that acts as a dominant source of the GWs.
A framework that is naturally suited for such connections is Cauchy Characteristic Extraction (CCE). CCE foliates the spacetime into a family of outgoing null hypersurfaces and formulates Einstein’s equations as an initialboundary value problem in a 2+2 characteristic decomposition. The mathematical details of this formalism can be found in Bishop et al. (1997); Handmer et al. (2016). CCE performs a characteristic evolution using the metric data on a timelike boundary of the Cauchy region (known as the worldtube) and propagates it to . At the radiation information is obtained as the Bondi news function Bondi et al. (1962). The GW strain can then be obtained from by a time integration,
(18) 
A key feature of this scheme is that each point at corresponds to a null hypersurface, which in turn corresponds to a particular (coordinate) time label on the world tube.
We can thus associate the average of the Kerrness on a 2sphere to spherical harmonic modes at . We illustrate this in Fig. 2. Here marks a specific timeslice (horizontal solid green line) in the Cauchy evolution region in a gauge chosen by the Cauchy code. The intersection of this timeslice with the worldtube boundary is a spatial (topological) 2sphere. The information on this 2sphere is propagated to along a null hypersurface labeled (solid purple line) as . The radiation feature carries the time stamp at , which, roughly speaking, arises from the 2sphere defined by the intersection of timeslice and the worldtube in the simulation and thus, we identify them to be associated.
Having established a framework to associate information on a 2sphere in the strongfield region to the waveform at , we now discuss the choice of the 2sphere used in this study. Motivated by analytical studies of test particles plunging into Schwarzschild BHs Price et al. (2016); Khanna and Price (2017), one might want to inspect the 2sphere associated with the peak of effective BH potential. However, locating it during the merger in a numerical simulation is nontrivial (if at all welldefined), and is beyond the scope of this paper. Furthermore, CCE cannot be performed from an arbitrarily small worldtube close to the horizon. This limitation arises because CCE is formulated in lightcone coordinates. In the regions very close to the horizon, lightcone coordinates can form caustics, leading to coordinate singularities. Because of these constraints, we choose the worldtube radius corresponding to the smallest coordinate 2sphere that is accessible to our procedure, but we visually verify that it contains strongfield dynamics by plotting in Figs. (b)b and (d)d.
ii.3 Inferring perturbation amplitudes via Kerrness
In order to give physical meaning to the values of the Kerrness measures outlined in Sec. II.1.2, we can compare their values (on a postmerger spacetime, for example) to those on a single BH with a known analytic perturbation. Specifically, we can compare the Kerrness measures during ringdown to those on a spheroidal QNM perturbed Kerr BH of the same final mass and spin, with varying dimensionless perturbation amplitude . This will provide a true physical comparison, as linearlyperturbed type D spacetimes are fully generic type I, and thus the Kerrness measures on the perturbed spacetime are expected to be nonzero Araneda and Dotti (2015). This comparison will allow us to infer the perturbation amplitude to which a particular coordinate time corresponds. We can then map this inferred amplitude onto the waveform using the methods in Sec. II.2.3.
ii.3.1 Kerrness measures on perturbed metrics
The perturbed metric is generated on a single slice for each by solving the Teukolsky equation and reconstructing the metric perturbation using a Hertzpotential formalism Yang et al. (2015); Lousto and Whiting (2002) (cf. Teukolsky (2015) for a general review). The resulting perturbation is then added to the background metric to give
(19) 
which is correct to linear order. The constraint equations for the metric are then solved to give a fully constraintsatisfying metric in KerrSchild coordinates using the extended conformal thinsandwich formalism (cf. Baumgarte and Shapiro (2010)).This introduces some nonlinear effects into the perturbed metric. Furthermore, the asymptotic radial behavior leads to blowup of the solution at large radii Ori (2003). Thus, before solving for , we multiply by an envelope of the form
(20) 
where is the radius of the outer horizon of the BH, is the width, and is the falloff of the envelope. Since the mapping of the Kerrness measures onto the waveform occurs at , as will be discussed in Sec. III.3, and the horizon typically has outer radius , we choose to give so as to minimally affect the perturbation at the extraction radius. Additionally, we choose in order to counteract the steep growth of the perturbation with radius. We plot the envelope in Fig. 4. In practice, the metric perturbation is generated using an extension of the code used in East et al. East et al. (2014), but with the QNM solution rather than an ingoing GW solution and using the full radial dependence.
Fig. 3 shows the behavior of the Kerrness measures averaged on a 2sphere of with on a BH of the same final mass and spin as the simulation outlined in Sec. III.1. The theoretical behavior of the Kerrness measures with perturbation amplitude is unknown GómezLobo (); Ionescu and Klainerman (2014), and thus this is the first (numerical) computation of the behavior. We first check that the measures converge to finite values with numerical resolution, thus representing real physical values. The Kerrness measures increase quadratically for small , then show higherorder effects for large . Type D 2 grows to (bestfit) quartic, Type D 3 and Kerr 1 become cubic, while Specialty and Type D 1 remain quadratic at , the largest amplitude for which we can solve for . In particular, the steep increase of the Type D 3 and Kerr 1 measures, which come from geometric conditions on KVs, indicates that at large enough perturbation amplitude, the slice fails to have even an approximate KV. Since the perturbation we are introducing is not axisymmetric, it makes sense that at large the slice loses this KV symmetry.
The linear perturbation regime corresponds to the region where the measures increase quadratically with , while the nonlinear regime approximately begins where one can see higherpower behavior. In this case, we see the transition from quadratic behavior around , suggesting that this is the approximate start of the nonlinear regime. In practice, one can normalize all of the values in this paper by . However, we do not do this for readability of the figures.
However, there are some sources of error in the analysis. The areal radius of the perturbed metric on a coordinate 2sphere of radius changes slightly with perturbation amplitude, changing by between and . Thus, a coordinateradius measure comparison does not happen on exactly the same 2sphere. Solving for changes the values of the mass and spin from the parameters used in creating . At the largest perturbation amplitude , the dimensionless spin changes by , while the mass changes by . We keep these errors in mind when computing the Kerrness values of the strongfield region in terms of and mapping them to the waveform for the binary case in Sec. IV.3.
ii.3.2 Mapping onto the waveform
A perturbation amplitude is associated with each timeslice of a postmerger spacetime in the strongfield region by the procedure described above. Since the procedure developed in II.2.3 allows us to associate simulation timeslices with the gravitational waveform at , we can map the perturbation amplitude associated with each timeslice to the corresponding parts of the waveform at . This gives an insight into deciding which portion of the waveform at can be modeled as being generated by linearly perturbed Kerr manifold, thus providing validation of start times chosen in data analysis that rely on perturbative description of Kerr.
ii.4 Measuring Kerrness on the horizon
In addition to local measures throughout a spatial slice discussed in Sec. II.1.2, Kerrness can also be evaluated on the postmerger apparent horizon (AH) (also known as the dynamical horizon in the literature). Owen describes a multipolar horizon analysis in Owen (2009), finding that the multipolar structure of a final BBH remnant was that of Kerr. Probing the multipolar structure also serves as a test of the nohair theorem Teukolsky (2015).
This formalism involves computing the mass multipole moments of the horizon as
(21) 
where is the scalar curvature of the horizon, is the metric volume element on the AH, and labels generalized (nonaxisymmetric) scalar spherical harmonics . These generalized spherical harmonics are computed from the eigenvalue problem
(22) 
where is the operator on the AH, and is its eigenvalue. In analogy with axisymmetric spherical harmonics , an effective is defined for the eigenvalues as
(23) 
where is the areal radius of the horizon. Since the values are timedependent, we refer to a given multipole by its final value.
As discussed in Owen (2009), the multipole moments that are zero on a Kerr BH either immediately vanish due to the symmetry of the dynamical horizon, or decay to zero from their excited values as the remnant BH settles to Kerr. The multipole moments that do not vanish on Kerr are functions of the mass and spin, and reach these values with increasing coordinate time. We use the code implemented and tested in Owen (2009) to compute the multipole moments. However, since the multipole moments are features of the horizon, we cannot map their behavior onto the waveform at . Moreover, CCE cannot be performed close to the horizon, as discussed in Sec. II.2.3. Nevertheless, we can compare the qualitative behavior of the multipole moments with those of the Kerrness measures as done in Secs. IV.1 and IV.2.
Iii Numerical implementation
iii.1 Description of simulation
We apply the methods outlined Sec. II to the numerical simulation presented in Fig. 1 of Abbott et al. (2016e), with similar parameters to GW150914, the first LIGO detection. The simulation is performed and the methods are applied using SpEC, the Spectral Einstein Code. The waveforms and parameters are available in SXS:BBH:0305 in the SXS Public Catalog SXS (). The simulation has initial mass ratio , and dimensionless spins and . The initial orbital frequency is . The final (postmerger) BH has dimensionless spin (within numerical error, as measured using the techniques in Scheel et al. (2009)) and mass . The inspiral proceeds for until the formation of a fullyresolved common AH. The visible part of the postmerger waveform on a linear scale has a temporal duration of .
Within a BBH simulation, the metric equations are evolved in a damped harmonic gauge Szilágyi et al. (2009); Lindblom and Szilágyi (2009), with excision boundaries just inside the apparent horizons Hemberger et al. (2013); Scheel et al. (2015), and minimallyreflective, constraintpreserving boundary conditions on the outer boundary Rinne et al. (2007). The spectral grid used during the inspiral of the simulation has an excised region for each BH. Once a common AH forms, the simulation proceeds for a few timesteps before switching to a new grid, in which there is one excision region for the new AH Hemberger et al. (2013). For this simulation, the gridswitch happens at . For more information on the code, see Lovelace et al. (2016).
iii.2 Implementation of Kerrness measures
We discuss the numerical implementation of the Kerrness measures outlined in Sec. II.1.2, and summarized in Fig. 1, on an NR BBH postmerger. Note that these measures will not be zero even on a numerical Kerr spacetime, due to the resolution of the simulation.
In order to quantify the Kerrness measures at each point, we convert the complex tensors into scalars by contracting them as
(24) 
where raising and lowering occurs using the spatial metric .^{9}^{9}9The Kerr 2 measure given in Eq. (16) requires that the imaginary part be zero, while the real part be . Hence, when evaluating Kerr 2, we measure the deviation of the imaginary part from zero, and the deviation of the real part from being positive (hence only including negative values). Throughout the rest of the paper, all of the measures will refer to their respective scalars generated using Eq. (24).
Because our simulations are performed using spectral methods, we expect errors to converge exponentially with increasing numerical resolution Press et al. (2007). In Fig. 5, we plot the Kerrness measures as a function of resolution for a single Kerr black hole; we see that the measures decay exponentially towards zero as expected.
SpEC solves a firstorder formulation of the Einstein equations, and therefore evolves both the spacetime metric and variables corresponding to its time and spatial derivatives Lindblom et al. (2006). The metric and first derivatives are available to the accuracy of the numerical simulation on each slice. Kerrness measures that require additional numerical derivatives, however, will have greater numerical noise and a higher numerical noise floor. The highest numerical order derivative needed to evaluate each measure is given in Fig. 1. Type D 4, which requires four numerical derivatives, is the noisiest measure and has a higher noise floor than the other measures, as shown in Fig. 5.
iii.3 Map from source to  implementation
In our study, we use a CCE implementation in SpEC (cf. Barkett and Scheel (2017), in prep). This implementation uses a no ingoing and outgoing radiation condition on the initial null hypersurface of the characteristic evolution. This means that the code treats the spacetime outside the worldtube as initially free of any gravitational radiation from the past.^{10}^{10}10During the Cauchy evolution, we perform the evolution with a boundary of and we do not neglect the backscatter from the region outside of the CCE extraction radius. Usually the CCE worldtube is placed at a large radius, and the CCE evolution begins at the start of the numerical simulation during early inspiral. However, here we begin CCE only at the merger portion of the Cauchy evolution, and in addition, we place the CCE worldtube at a very small radius. This means that extracted waveform does not contain contribution coming from the inspiral part of the dynamics.
By decreasing the radius of the extraction worldtube progressively by , we find the smallest radius of the worldtube that our procedure can be applied to occurs at a coordinate radius of . For a radius of , the CCE procedure can not be performed, presumably due to the formation of caustics. At , we get a very glitchy and unreliable extraction of the news function.
However, performing the CCE from such small radii gives rise to an additional complication. Since time stamps on the waveform at are induced by the simulation coordinates, the news function obtained is not necessarily in an inertial gauge. In a standard CCE scheme, a gauge transformation is applied to the news function in order to obtain it in an inertial gauge. To preserve the map between the time in simulation gauge and the time coordinate on the extracted news function, we do not perform this gauge transformation. We see the effect of the gauge transformation in the waveform at as a mixing of mode amplitudes. The effect is very small when the worldtube boundary for CCE is large i.e., lies in the weak field region. For instance, for a worldtube boundary of the effect of this transformation is negligible. To confirm this, we compute the overlap between the news extracted from with and without the gauge transformation. The overlap is defined as,
(25) 
where is the frequency domain Fouriertransformed news function, and denotes complex conjugation for ease of readability, and is the norm Apostolatos (1995).
We find that the mismatch, , is . This overlap computation uses only the merger and postmerger parts of the news function for the dominant () spinweighted spherical mode. However, for a worldtube radius of , there could be significant amplitude deviations between the waveforms in the simulationcoordinateinduced gauge and the inertial gauge. Because of technical difficulties in the code implementation, we could not apply the gauge transformation to an extraction from and quantify the difference. ^{11}^{11}11In the CCE implementation we use, the inertial coordinate at is setup by solving and quantity used in definition of the conformal factor in Handmer et al. (2016); Handmer (). The possible technical issues when performing extraction from small radius are  a) initializing the inertial coordinates are tied to the Cauchy evolution coordinates and b) the evolution scheme may not be accurate in the strong field regime. The interpolation scheme is not designed for extraction performed from very small radii.
Furthermore, before the noninertial to inertial gauge transformation, every point on at the same timestamp on the waveform corresponds to the same null hypersurface and therefore to the same simulation coordinate time. After the transformation, this is no longer true: the waveform seen at different sky directions with the same timestamp on the waveform corresponds to different null hypersurfaces and therefore different values of simulation coordinate time. This happens because the choice of the 2sphere is gaugedependent.Therefore, we omit the gauge transformation, as the aim in this paper is to connect the nearzone to the wave zone, requiring us to retain the timestamps.
Additionally, the initial noingoing radiation condition neglects gravitational radiation coming from the inspiral. This may be significant for extraction done at small radii, where the initial CCE null hypersurface connects the strongfield region close to merger to and may contain significant radiation from the inspiral. This could contribute towards the discrepancy between the and waveforms.
To assess this difference, we compare the news function obtained by extraction performed from with the extractions performed from the worldtubes of larger radii, all without the gauge transformation. The result of this is presented in Figure 6. We observe that all the extractions from radii greater than converge with radius, indicating that the effect of the gauge transformation is insignificant at these radii. Further, the extraction from has a significant amplitude discrepancy with the other extractions, particularly in its first cycle. Therefore, we would ideally wish to map the strongfield information computed on the 2sphere at a coordinate radius of on the news function that has been extracted from a larger radius like .
We do this mapping in two steps. First, we map the strongfield information computed on the 2sphere at a coordinate radius of onto the CCE performed from a worldtube of using the framework described above. Next, we note that the phase evolution of extraction from agrees with the extractions from larger radii.^{12}^{12}12The timederivative of the phase gives the instantaneous frequency of the gravitational radiation. We verify this in Fig. 7. Then we align the news function extracted from to the extraction from larger radii as shown in Fig. 6. The alignment is done such that the overlap between the CCE extracted news function from different world tube radii is maximized. The maximum normalized between the news function extracted from and is . Incidentally, this alignment is equivalent to aligning the real part of the news function at its global minima (or global maxima of the absolute value). Table 1 lists the time shifts that have been applied in order to align the news function extracted from a radius with extraction done at .
Worldtube radius  Alignment shift wrt  

=  .  
=  .  
=  .  
=  . 
Using this alignment we map the time stamps on the to those on . From this, we infer the mapping of strongfield information at on to the extraction done from , thus mapping the strongfield information onto the news function as seen in near inertial gauge.
We summarize our algorithm for mapping the strongfield information onto the news function:

Perform CCE from worldtube with radius of the 2sphere that lies in the strongfield region (whose evolution you wish to map on to the news function seen at ) without the final noninertial to inertial gauge transformation. The time stamps on this extracted news function are induced by the time coordinates in the simulation, thus providing a natural map between the evolution of the strongfield region and the wave zone.

Perform CCE from a large worldtube radius where the effect of the noninertial to inertial gauge transformation is negligible.

Align the news functions obtained in steps 1 and 2 such that the overlap between the waveform is maximized.

Use this alignment to map the time stamps of the news function extracted in step 1 to that in step 2. The 2sphere chosen in step 1 at the timeslice marked with the simulation time coordinate can be associated as the dominant source of the feature at with the same time stamp.
Iv Results
We now present the results of performing the analysis outlined in Secs. II and III on the GW150914like simulation detailed in Sec. III.1. Sec. IV.1 presents the behavior of the multipole moments of the AH, which provides a comparison for the Kerrness measures on the simulation volume. Sec. IV.2 discusses the results of evaluating the Kerrness measures on the postmerger spacetime and mapping them onto the waveform at , presenting them in terms of the percentage decrease from their peak values. Sec. IV.3 presents the results of comparing the Kerrness measures on the postmerger spacetime to values on perturbed data, in order to infer the perturbation amplitude in the strongfield region, and mapping them onto the waveform, presenting them in terms of the inferred perturbation amplitude . The percentage decrease from the peak value and can then be used to estimate the overall level of Kerrness and validate choices for the start time of ringdown. Finally, in Sec. IV.4, we discuss the implications of these results on analyzing ringdown in GW data, and in Sec. IV.5 we compare our results to the ringdown start times chosen in the GW150914 testing GR study Abbott et al. (2016a).
iv.1 Horizon behavior and multipolar analysis on BBH ringdown
As a first measure of Kerrness, we apply the horizon multipolar analysis outlined in Owen (2009) and summarized in Sec. II.4 to the simulation described in Sec. III.1. Fig. 8 presents the behavior of the AH. The areal mass of the AH, given by where is the proper area of the AH, sharply settles to a final value. The minimum and maximum radii are initially noisy, as the AH is initially peanut shaped, but they decrease exponentially with coordinate time, showing a settling of the AH to the final state. However, the radii are coordinatedependent measures, and thus to check if the BH settles to Kerr it is more instructive to look at the AH multipole moments.
Fig. 8 shows the behavior of the initially nonvanishing quadrupole and hexadecupole moments, labeled by their corresponding at the final time, as given in Eq. (23). The quadrupole moments correspond to and the hexadecapole moments correspond to . The multipole moments behave as expected for a generic simulation remnant settling to a Kerr BH. As explained in Owen (2009), two of the five quadrupole moments immediately vanish by reflection symmetry, while two others exponentially go to zero (eventually hitting a numerical noise floor) as the final remnant settles to Kerr. Four of the nine possible hexadecupole moments immediately vanish from reflection symmetry, while four go exponentially to zero as the remnant settles to Kerr. Note that the and moments vanish on Kerr due to symmetry. As in Owen (2009), one quadrupole moment () and one hexadecupole moment (), both corresponding to , do not vanish, but rather attain a constant value in line with that of a Kerr BH of the same final mass and spin.
The multipolar behavior thus confirms that the final state of the AH is that of a Kerr BH. This serves as an independent test of Kerrness, and thus one would expect the Kerrness measures presented in Sec. II.1.2 to also show the strongfield region exponentially settling to Kerr. This also serves as numerical evidence for BH uniqueness, as the final remnant of a BBH merger is indeed Kerr, as also discussed in Owen (2009). Similarly, since the final multipolar structure can be described completely by the mass and spin, this serves as numerical validation of the nohair theorem.
iv.2 Measuring and mapping Kerrness onto the waveform
The goal in this section is to validate choices of the start time of ringdown using Kerrness measures on the GW150914like system described in Sec. III.1. We now present the results of evaluating the Kerrness measures outlined in Secs. II.1.2 and III.2 (and summarized in Fig. 1) in the strongfield region and mapping them onto the waveform at using the procedure given in III.3. These measures are evaluated pointwise on each slice, and we map the value on a 2sphere at a radius of onto the news function. Recall that larger values of the Kerrness measures indicate greater deviation from being locally isometric to Kerr.
Fig. 9 shows the Kerrness measures averaged at various coordinate radii on each slice of the postmerger spacetime, presented as a function of coordinate time. All of the measures decay exponentially toward zero, showing that the spacetime approaches an isometry to Kerr. This confirms the results of the multipolar analysis in Sec. IV.1. Additionally, this serves as a numerical verification of BH uniqueness, as the final state of a BBH merger is isometric to Kerr. The behavior of the measures at large radii (such as in this case) is especially interesting to the question of BH uniqueness, which is particularly concerned with the domain of outer communication Ionescu and Klainerman (2014).
Fig. 10 shows the behavior of the Speciality Index, an algebraic measure (Type D 1) and a geometric measure (Kerr 1) in the volume, as a function of increasing coordinate time. We see a distinct quadrupolar pattern in all our measures (the equatorial plane has a modal pattern that corresponds to ), consistent with the dominant mode of gravitational radiation. Furthermore, the Speciality Index and Type D 1 measures, which determine properties of the PNDs, settle first further from the BH, while the geometric Kerr 1 measure, which is determined by properties of the KV, first settles closer to the BH.




The Kerr 2 measure, which constrains the NUT parameter, is effectively constant throughout the ringdown, as shown in Fig. 11. Since the NUT parameter is one of the hairs of a generic type D manifold, Fig. 11 confirms that a NUT charge is not generated during a BBH merger. We thus do not include it further in our analysis.
Of these measures, two are algebraic constraints—Type D 1 and Type D 2—and three are geometric constraints on the KV, Type D 3, Type D 4, and Kerr 1. In Fig. 9 we see that all the algebraic measures decay in a similar fashion and all the geometric measures decay similarly. Type D 4, which requires 4 numerical derivatives, is visibly noisier than the other measures. This measure checks if the vector identified as satisfies the Killing equation and is crucial for a rigorous mathematical characterization of Kerr manifold. However, all geometric measures depend on the same Killing vector and we observe that Type D 4 has a similar decay property as Type D 3 and Kerr 1. Thus, we do not include the noisier Type D 4 in our analysis, rather treating Type D 3 as a proxy for both.
Each measure at each radius in Fig. 9 eventually reaches a floor. This is confirmed to be a numerical noise floor in Fig. 12, where the floor is shown to exponentially converge to zero with numerical resolution. The radial behavior of the Kerrness measures stems from the radial behavior of the Weyl tensor and the metric quantities. For example, for a stationary background, and , and thus Speciality Index given in Eq. (7) should be , which we indeed observe.
The analysis outlined in Sec. III.3 requires the Kerrness measures to be extracted at in order to map them to the news function. Fig. 10 shows that the Kerrness measures have strong support at , thus justifying the choice of radius as being in the near field.^{13}^{13}13The measures at in Fig. 9 behave similarly to those at indicating that also behaves like the near field region, but unfortunately we have not been able to perform CCE from this small a radius.
The Kerrness measures quantify the violation of the conditions for a manifold to be isometric to Kerr and therefore, they need not have the same dimensions and sensitivities. Thus, one cannot compare the absolute magnitudes of these measures with each other and directly translate their value into statements on validity of start time of perturbative regime. In order to normalize and combine them into an overall measure of Kerrness, we use the concomitant percentage decrease from their peak values.
We present the percentage decrease of each of these measures from their peak values mapped on to the news function in Fig. 13 and Fig. 14. In the bottom panels of these figures, the news function is plotted as a function of time. On the same time axis, the top panel depicts the corresponding evolution of the Kerrness measure in the strongfield region. The waveform feature in the bottom panel at a particular time coordinate is associated to the timeslice carrying the same time label, via sourceeffect association outlined in Sec. II.2.3. In the bottom panel, the Kerrness value at this time characterizes the deviation from Kerr.




In these figures, we delineate 6 lines marking the percentage decrease from the peak value of each of the Kerrness measures as a function of time—both in the strongfield region and on the news function at . As stated before, these measures have different decay properties and so do not decay to a particular percentage of their peak value at the same time. The difference between the time at which measures decay to a particular percent is tabulated in Table 2.
of peak value  Spread in time  Combined % time  
100  12.  3847.  5  
50  9.  8  3857.  
30  9.  3860.  7  
10  8.  3  3867.  7  
5  8.  7  3871.  9  
1  6.  1  3881.  3 
We present the combined percentage decrease from the peak value on the news function in Fig. 15. The shaded bands correspond to spread in percentage decay on the news function. The widths of these bands are given in Table 2. The solid line at the end of each band marks the time when all these measures have decayed to the indicated percentages and this can be used to conservatively choose the start time.
Furthermore, in this figure we do not include the Specialty Index. The Specialty Index is an independent measure that quantifies if the manifold is algebraically special. Since this is the weakest condition in our Kerrness characterization scheme, we see that it gets satisfied earliest on the postmerger simulation from Fig. 14. The of peak line which occurs unexpectedly late arises because of numerical reasons. We assert this by looking at the nearly flat nature of Specialty Index curves in Fig. 9 at late times, very close to the numerical noise floor.
We observe that all measures decay to of their peak value within half a cycle from the peak of the news function. Further, in approximately one cycle, all the measures are reduced to of their peak values. The spread in each of the bands is about when we include all the Kerrness measures in computing the band, and this shrinks to when we exclude Specialty Index.
We combine the measures with equal weights, thereby presenting a conservative result. Furthermore, we have repeated our analysis with larger worldtube radii and confirmed that our results do not change significantly; specifically, for and the results change only by .
iv.3 Estimating and mapping the perturbation amplitude onto the waveform




In order to provide a physical understanding for the values of the measures in the strongfield region shown in Figs. 9 and 10, we can compare the values to those on an initial slice of a perturbed Kerr BH with the same final mass and spin as the BBH simulation, as outlined in Sec. II.3.1. We can then map the inferred strongfield perturbation amplitude onto the waveform using the procedure outlined in Secs. II.2.3 and III.3. This procedure involves the following steps:

Generate perturbed Kerr manifolds for a range of amplitudes .

Compute the Kerrness measures at on these slices.

Compute the Kerrness measures at on the postmerger BBH simulation (verifying that the gaugeinvariant areal radii of the coordinate 2spheres are approximately equal for the singleBH and the BBH case).

Identify the coordinate time in the postmerger BBH simulation at which the Kerrness measures at agree with those on the perturbed Kerr slice for a given — this gives a crossing time for this .

Use this crossing time to map the inferred onto the waveform.
Fig. 17 shows the inferred for the BBH ringdown simulation as a function of coordinate time in the simulation. The gaugeinvariant areal radii at on the BBH simulation slices and on the metric perturbation are within . The values of the Kerrness measures on the perturbed data vary quadratically with , as shown in Fig. 3. At higher values of , they obtain higherpower dependence, as discussed in Sec. II.3.1. Each Kerrness measure decays through various as the simulation progresses. Type D 1 and Type D 2, the two algebraic conditions, have comparable crossing times for a given , while the two geometric KV conditions, Type D 3 and Kerr 1, also have comparable crossing times. Speciality Index crosses around before the other measures, in part because it is a weaker condition that the others. Each crossing time has an intrinsic spread due to sampling, and not all measures cross each due to numerical noise floors, leading to spreads in crossing time.
In Fig. 16, we visually check the spacetime features by comparing corresponding to and on the perturbed Kerr metric with the corresponding timeslice during the postmerger simulation. The crossing time spread for a particular arises because of the imperfect mapping between an analytically perturbed Kerr BH and the postmerger spacetime. Therefore, unlike in an ideal mapping, the combined crossing times will have a spread. In particular, the difference in the features between the postmerger and the perturbed Kerr slice indicates a difference in symmetry and explains the larger spread in the crossing time between the KVdependent measures. We see that the spread in the combined crossing times using only algebraic measures is much smaller than when we include the geometric measures.
We next map the inferred perturbation amplitude to the news function, using a procedure similar to the one in the previous section, and present the result in Fig. 18. The top panel of the figure indicates the crossing time for the Specialty Index, the middle panel for the algebraic measures, and the bottom panel shows that for geometric measures. The spread in the crossing time for the algebraic measures decreases from at the start, to our sampling rate, . This occurs because at the very start of postmerger, the system is not yet in a perturbative regime and therefore, our mapping contains a larger error. Geometric measures are more drastically affected by the imperfections in the mapping, indicating the differences in the symmetries of the two systems. On including the geometric measures, the crossing time spreads to . We confirm that the spread of the crossing times calculated using the algebraic measures is always contained within the spread of crossing times calculated using the geometric measures.
As the signal decays from the peak to a barely visible amplitude on a linear scale ( cycles) at , the corresponding perturbation in the strongfield region decreases by an order of magnitude. The peak of the news function corresponds to a perturbation amplitude of . Further, it takes about 2 cycles in the wave zone for the perturbation amplitude to decay to half its peak value. Also, by the time the perturbation amplitude decays by an order of magnitude, there is hardly any power left in the signal.
iv.4 Implication of the start time on data analysis
iv.4.1 From news to h
In order to compare the Kerrness measures on the GW to the loss in signaltonoise ratio (SNR) at the times used in Abbott et al. (2016a), we must first calculate the strain from the news function, and then calculate the merger time. As outlined in Sec. III.3, can be calculated by integrating the CCE news function. One can also independently calculate using the ReggeWheelerZerilli (RWZ) (cf. Rinne et al. (2009) for details on the method) Regge and Wheeler (1957); Zerilli (1970a, b); Moncrief (1974) method and then extrapolating it in powers of the extraction radius (cf. Boyle and Mroue (2009) for details). The RWZ method and extrapolation have been implemented and tested in SpEC Boyle and Mroue (2009); Taylor et al. (2013), and the strain calculated by this method was presented in the GW150914 detection paper Abbott et al. (2016e). This method, however, has a different retarded time axis Boyle and Mroue (2009) than the CCE news function. Thus, we differentiate the RWZ strain to get a news function, and shift it to align in phase with the CCE news function. We check the CCE results by comparing the output of the two methods, presenting the results in Fig. 19.
In the GW150914 testing GR study Abbott et al. (2016a), is defined as the point at which the quadrature sum of the and polarizations of the mostprobable, or maximum a posteriori (MAP) waveform, produced by EffectiveOneBody (SEOBNRv4) template Pürrer (2016) is maximal. For this study, we use the spinweighted spherical harmonic mode of the MAP waveform, as this is the leastdamped QNM. In this study, rather than using the EOBNR waveform, we calculate based on the time of maximum amplitude of the timeshifted RWZ strain, as
(26) 
where
(27) 
We find .
iv.4.2 Start time and the SNR
No. of cycles  SNR  Kerrness  

Peak  60  100  7.5  
cycle  30  40  40  50  7.5  
1  cycle  20  25  35  30  5 
1  cycles  10  20  8  12  3.5 
2  cycles  10  7  5  2  2.5 
2  cycles  10  1  12 
3  cycles  5  1  0.5  1 
While picking too early a start time for an analysis that relies on being in ringdown gives inaccurate and biased results, picking a start time too late leads to a large statistical error. Since the amplitude of the signal decays exponentially with time, the SNR in ringdown decreases as exponentialsquared with the start time. Consequently, the spread in the posteriors during estimation of ringdown parameters, which goes inversely proportional to matchfiltered SNR, increases and gives rise to large statistical uncertainties. Therefore, one must chose an optimal middle ground considering both these factors.
In the top panel of Fig. 20, we show the percentage decrease in matchfiltered SNR with the start time of the ringdown. A matchfiltered SNR is a noiseweighted scalar product between the signal and the template and is defined as
(28) 
where denotes complex conjugation for ease of readability. Here, corresponds to a ringdown waveform that is tapered at and acts as a signal. We filter this against the template, , which is tapered with varying start time. Further, corresponds to power spectral density (PSD) of aLIGO at design sensitivity Aasi et al. (2015). However, since we present our results in terms of ratios, our analysis remains valid to any detector noise curve. Then, a Fourier transform is taken to evaluate Eq. (28). Here we use only the spinweighted spherical harmonic mode of the RWZ strain waveform computed in Sec. IV.4.1. The system is considered to be optimally oriented with respect to the detector for this calculation.
The tapering is done with a tanh window function defined as
(29) 
is the time around which the tapering is centered and it is set to the start time of the perturbative regime. is set to 10 in making the top panel of Fig. 20. Furthermore, we confirm that our results do not change significantly with the tuning parameter using .
We then present percentage decrease of SNR in the top panel of Fig. 20 by defining for start time at . Further, on this same plot we also indicate the amplitude of perturbation in the strongfield region (as calculated using the algebraic measures) at the start time, giving an insight into how statistical and systematic errors vary with the choice of start time.
The bottom panel of Fig. 20 presents the total energy radiated through the mergerringdown as a function of time. This indicates the strength of GW signal and is calculated by integrating Ruiz et al. (2008)
(30) 
Furthermore, on the same plot we mark the percentage decrease of the Kerrness measures from their peak values, providing a comparison between the strength of the signal left for performing the analysis versus Kerrness evaluated in the strongfield region.
To impress the sharp tradeoff in systematic and statistical uncertainties in choosing the start time of the ringdown and, to provide an intuition of implication of Fig. 20, we present the spread in estimation of dominant QNM frequency, in Fig. 21. For this, we calculate the spread using the Fisher information matrix formalism similar to that in Eq. 4.1a of Berti et al. (2007), for a GW150914like system. In particular, we set to 253.7 Hz and the quality factor, to 3.2. However, we emphasize that this is a rough estimate intended only to provide intuition and, we plan to follow this up by a rigorous Bayesian parameter estimation in the future.
We present the interplay between the systematic and statistical uncertainty concisely in Table 3. Furthermore, we find that by the time the news function peaks, the SNR has already dropped down to . However, at this time the algebraic Kerrness measures are at their peak value. We also observe that by about a cycle of news function, there is less than 20 percent SNR left in the signal. Therefore, there seems to be a sharp tradeoff between the systematic modeling error and statistical uncertainties.
iv.5 Comparison with GW150914 testing GR paper
The test of consistency of ringdown of the GW150914 signal with the analytically predicted QNM frequency is given in Fig. 5 of Abbott et al. (2016a). The analysis chooses various start times of ringdown, namely . At a start time of (or later), parameter estimation of the dominant QNM in ringdown is consistent with predictions using initial masses and spins.
A time for the system corresponds to from . In our analysis, (cf. Eq. (26)), while the peak of the news function occurs at . Thus, corresponds to after the peak of the news function. In this region, as shown in Fig. 22, the perturbation amplitude is . Our analysis indicates that this corresponds to a relatively large deviation from Kerr. Recall that Fig. 3 suggests that is the approximate start of the nonlinear regime.
With a start time of , the SNR was about and the spread in the estimate of QNM frequency was roughly Abbott et al. (2016a). Because of this low SNR and high spread, the GW150914 testing GR analysis may not have been sensitive to the large nonKerrness we see close to the BH. However, in the case of higher SNR, where the analysis is sensitive to the systematics of the ringdown model, our study suggests picking a later start time.
Our analysis uses geometric and algebraic conditions to identify isometry to Kerr. However, these conditions do not directly measure the deviation of the curvature BH potential from that of Kerr. Since the QNM are intrinsic tests of the BH potential along with the boundary conditions, deviation of QNM frequencies will depend on details of the BH potential, and thus are not directly quantified in our measures. Additionally, the parameters used in this study correspond to SXS:BBH:0305 waveform used in the GW150914 detection paper Abbott et al. (2016b), which are slightly different from those of the MAP waveform used in the testing GR paper.
V Conclusion
In this study, we present a method for validating choices of the time at which a BBH GW signal can be considered to enter the ringdown stage. This is done by computing algebraic and geometric measures of Kerrness in the strongfield region of an NR simulation of a BBH ringdown, and then associating each point on the asymptoticframe waveform with a particular value of these Kerrness measures. Thus, for each point on the asymptoticframe waveform there is an estimate for how close the BH spacetime is to Kerr spacetime. This is the first time this set of measures, proposed in GómezLobo (2016), are evaluated in the strongfield region. This is also the first time measures of Kerrness in the strongfield region is mapped onto the waveform. We outline this method in Secs. II and III, and implement this analysis in Sec. IV on a GW150914like NR simulation.
We observe that after merger, the Kerrness measures of a BBH ringdown simulation decrease exponentially with coordinate time in the simulation, eventually settling to a numerical noise floor, as shown in Fig. 9. This decay is consistent with measuring Kerrness using multipole moments of the apparent horizon, as in Fig. 8 and Owen (2009). In all cases, the measures on the final slice of the simulation indicate that the final remnant is a Kerr BH, thus providing numerical consistency with the BH uniqueness theorem. Moreover, we find that the final state in the multipolar analysis depends just on mass and spin, which serves as a confirmation of the nohair theorem in the strongfield region. Additionally, as shown in Fig. 10, the Kerrness measures have a quadrupolar (with ) structure consistent with the dominant gravitational radiation. The geometric measures, which rely on the existence of a Killing vector field, first decay to zero close to the horizon, then later they decay at larger radii as gravitational radiation propagates out. On the other hand, algebraic measures, which depend on principal null directions, first decay to zero at larger radii, before decaying near the BH. We also see that the NUT parameter remains zero during merger and ringdown, as shown in Fig. 11.
These gaugeindependent Kerrness measures are crucial to the nonlinear stability analysis of Kerr, as they quantify the deviation from being isometric to Kerr. The analytical behavior of these measures with perturbation amplitude is unknown Ionescu and Klainerman (2014); GómezLobo (). Through this study we provide insights into their numerical behavior in Fig. 3. We find that all of these measures scale quadratically with for low amplitude perturbations, but acquire higherorder nonlinearities for larger perturbation amplitudes. Furthermore, in Figs. 9 and 10, we provide the radial behavior of these measures, up to a large radius of . For a BBH simulation, we track these measures starting from merger, where linear perturbation theory is not expected to hold. Despite the large initial excitation, our study shows that the BBH ringdown simulation evolves to a final Kerr state, providing a numerical validation of the nonlinear stability of Kerr.
To connect the Kerrness measures in the strongfield region to the asymptotic waveform at , we use CCE, which evolves Einstein’s equations on a foliation of outgoing null hypersurfaces. A null characteristic evolution can be done only in a region free from caustics. We demonstrate that CCE results using a worldtube at are consistent with those done from larger radii. This implies that during ringdown, caustics only exist very close to the BH. Furthermore, we show that the map between the strongfield region and the wave zone can be extended all the way in to .
Although caustics do not form, we see in Figs. 10, (b)b, and (d)d strong features in the curvature quantity in the region enclosed by . This implies that our extraction radius choice of lies within the strongfield and within the support of the BH potential.
In Fig. 13, we label each point of the BBH ringdown waveform with the percentage decrease of the Kerrness measures in the strongfield region relative to their maximum values. In order to give a physical interpretation of the values of the Kerrness measures, we compare them throughout the postmerger spacetime to those evaluated on a QNM perturbed Kerr BH of the same final mass and spin. From this we infer the amplitude of BH perturbation during ringdown and map onto a particular point in the BBH ringdown waveform; this is marked on the BBH ringdown waveform in Fig. 18.
As the BBH simulation proceeds after merger, the strongfield region starts looking like Kerr, indicating the validity of perturbative analysis. However, as time progresses, the amplitude of the ringdown decreases, leading to a rapid decay in SNR in a GW detection. We find that by the time the Kerrness measures decrease to of their peak values, there is only about SNR left in the signal. In terms of perturbation amplitude close to the BH, this maps to an amplitude between . This occurs after cycles of the news function, which corresponds to after . Additionally, we find that the start time of ringdown used in Abbott et al. (2016a), , corresponds to an amplitude of . In future detections with higher SNR, where the statistical noise is significantly smaller, one may need to choose a later start time to perform precision tests of GR such as nohair theorem tests.
A future extension to this project would be to investigate methods that allow us to perform similar sourceasymptotic frame associations for smaller radii. For instance, the light ring would be an interesting region to monitor as it is crucial to the QNM structure. This can perhaps be done numerically through raytracing methods such as those used in Bohn et al. (2016) and Bohn et al. (2015), to understand the evolution of the peak of the BH potential (if it forms). Another possible technique could be to try performing CCE from smaller radii after the high amplitude of the initial excitation has reduced. Additionally, being able to perform this association at smaller radii would allow one to understand the propagation of perturbations very close to the BH horizon onto the waveform; these are expected to be redshifted and appear on the waveform with a large time delay.
Another avenue of future work would be investigating the effects of implementing a more realistic condition on the initial null hypersurface by relaxing the noingoingwaves condition used in performing CCE. In addition, we can study the tradeoff involved in choosing an earlier ringdown time, which will decrease the spread in recovered ringdown parameter posterior distributions and increase the systematic errors that arise because of deviations from Kerr in the strongfield region.
The methods used in this paper can be applied to future BBH detections in order to guide the choice of the start time of ringdown. For the sake of quick reference to the procedure described in this paper, we concisely enumerated the steps in Appendix A. Note that the results of this paper approximately holds for any equal mass system with an appropriate mass rescaling (cf. footnote 2). Our method would better allow one to perform precision tests of GR that depend on being in perturbative regime, such as tests of the nohair theorem and area theorem. With this procedure, we provide an algorithmic way to check whether an unexpected deviation in a QNM analysis is due to not being in the perturbative regime, rather than due to a violation of GR or corresponding theorems. For future detections, we plan to repeat this analysis using an NR simulation with the MAP waveform parameters.
Acknowledgements.
We would like to thank Alfonso GarcíaParrado GómezLobo, Kevin Barkett, Mike Boyle, Yanbei Chen, Joshua Goldberg, Casey Handmer, Daniel Hemberger, Maximilliano Isi, Badri Krishnan, Nicholas Meyer, Harald Pfeiffer, Leo Stein and Vijay Varma for many valuable conversations. In particular, we would like to thank Mike Boyle, Casey Handmer, Harald Pfeiffer, and Leo Stein for careful reading of this manuscript. We would like to thank William East for helping to generate the QNM metrics, and Geoffrey Lovelace for supplying the BBH simulation data used in this study. This work was supported in part by the Sherman Fairchild Foundation, the Brinson Foundation, NSF grants PHY–1404569 and AST–1333520 at Caltech, and NSF grant PHY–1606654 at Cornell University and PHY–1404395, PHY1707954 and PHY1352511 at Syracuse University. MO gratefully acknowledges the support of the Dominic Orr Graduate Fellowship at Caltech. We used SpEC (Spectral Einstein Code) to perform the simulations and analysis spe (). Computations were performed on the Zwicky and Wheeler clusters at Caltech, which are supported by the Sherman Fairchild Foundation and by NSF award PHY–0960291. The BBH simulation was performed on the ORCA cluster at California State University, Fullerton (CSUF), supported by CSUF, NSF grant No. PHY–142987, and the Research Corporation for Science Advancement.Appendix A Quick outline of our procedure
In this appendix, we concisely provide an outline of the algorithmic procedure developed in this paper. For future BBH detections, we propose to follow the steps enumerated below

Performing an NR simulation with waveform parameters inferred from parameter estimation, and saving the metric data,

Generating worldtube data for various extraction radii and performing CCE from the innermost possible radius,

Evaluating the Kerrness measures on the metric data at this radius for BBH ringdown,

Evaluating the Kerrness measures on QNM perturbed data with the same final mass and spin, and inferring corresponding perturbation amplitude from the Kerrness values,

Mapping the Kerrness measures and inferred perturbation amplitudes to the waveform via nullcharacteristics,

Using these results to validate choices for the start time of ringdown in detector data analysis.
Appendix B KerrNUT parameters
In this appendix, we provide a review of the parameters of the KerrNUT solution. The Kerr family of vacuum solutions is unique when one imposes axisymmetry, stationarity and regularity on the BH horizon along with asymptotic flatness. However, if one allows for generalization by relaxing the asymptotic flatness condition, one arrives at a family of solutions called KerrNUT. This solution is a part of the broader family of EinsteinMaxwell type D solutions. This generalized family of spacetimes is parameterized by 6 parameters (potentially 7 if one includes the cosmological constant ). In Table 4, we summarize the parameters, as well as their physical meaning and symbols used in various texts.
The general EinsteinMaxwell Type D solution (including cosmological constant ) has the form given in Eq. 21.11 of Stephani et al. (2009), with parameters , , , , , and . refers to the mass parameter (closely related to the mass of the BH), is related to the angular momentum parameter (closely related to the spin of the BH), is related to the acceleration , is the electric charge, is the magnetic charge, and is known as the NUT parameter. As outlined in Plebanski and Demianski (1976), the mass and the NUT parameter form a complex quantity, as do the angular momentum and the acceleration, similarly to the electric and magnetic charges. In Plebanski and Demianski (1976), and do not appear in the curvature quantities, and are called kinematical parameters, while the others are dynamical parameters.
As shown in Table 21.1 of Stephani et al. (2009), setting all of the parameters to zero except for , (and hence and ), and to zero yields the KerrNewman solution, while also setting yields the ReissnerNordstrom solution. KerrTaubNUT metrics, meanwhile, are parametrized by mass, spin, and , with , and are thought to be unphysical Aliev et al. (2008). The vacuum BBH case considered in this study, meanwhile, sets and , since there are no electric or magnetic charges at the start of the simulation, and no sourcing of them during the simulation.
An accelerating and rotating BH with a NUT charge will have nonzero , , , and , with . A Kerr solution with a NUT charge will then have . An accelerating and rotating BH, meanwhile, will have . Finally, the Kerr solution has both and . An illustration of this is provided in Fig. 1 of Griffiths and Podolsky (2005). The condition gives the Kerr 2 condition considered in this paper, given in Eq. (16).
After setting , the parameters , and are related to the mass and spin of a BH are as follows,
(31) 
Since, and for a Kerr BH, the condition that gives , which corresponds to the Kerr 3 condition given in Eq. (17).
Stephani Stephani et al. (2009) 
GarcíaParrado GómezLobo (2016) 
Plebanski Plebanski and Demianski (1976) 
Griffiths Griffiths and Podolsky (2005) 


Cosmological constant  
Mass parameter  
NUT parameter  
Angular momentum parameter  
Acceleration parameter  
Electric charge  
Magnetic charge 
References
 Abbott et al. (2016a) B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev. Lett. 116, 221101 (2016a), eprint 1602.03841.
 Teukolsky (1973) S. A. Teukolsky, Astrophys. J. 185, 635 (1973).
 Press and Teukolsky (1973) W. H. Press and S. A. Teukolsky, Astrophys. J. 185, 649 (1973).
 Teukolsky and Press (1974) S. A. Teukolsky and W. H. Press, Astrophys. J. 193, 443 (1974).
 Abbott et al. (2016b) B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev. Lett. 116, 061102 (2016b), eprint 1602.03837.
 Abbott et al. (2016c) B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev. Lett. 116, 241103 (2016c), eprint 1606.04855.
 Abbott et al. (2016d) B. P. Abbott et al. (Virgo, LIGO Scientific), Phys. Rev. X6, 041015 (2016d), eprint 1606.04856.
 Abbott et al. (2017a) B. P. Abbott et al. (VIRGO, LIGO Scientific), Phys. Rev. Lett. 118, 221101 (2017a), eprint 1706.01812.
 Abbott et al. (2017b) B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, V. B. Adya, et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 119, 141101 (2017b), URL https://link.aps.org/doi/10.1103/PhysRevLett.119.141101.
 Mazur (2000) P. O. Mazur (2000), eprint hepth/0101012.
 Misner et al. (1973) C. Misner, K. Thorne, and J. Wheeler, Gravitation, no. pt. 3 in Gravitation (W. H. Freeman, 1973), ISBN 9780716703440, URL https://books.google.com/books?id=w4Gigq3tY1kC.
 Dreyer et al. (2004) O. Dreyer, B. Kelly, B. Krishnan, L. S. Finn, D. Garrison, and R. LopezAleman, Classical and Quantum Gravity 21, 787 (2004), eprint grqc/0309007.
 Gossan et al. (2012) S. Gossan, J. Veitch, and B. S. Sathyaprakash, Phys. Rev. D 85, 124056 (2012), eprint 1111.5819.
 Kamaretsos et al. (2012a) I. Kamaretsos, M. Hannam, S. Husa, and B. S. Sathyaprakash, Phys. Rev. D 85, 024018 (2012a), eprint 1107.0854.
 Berti et al. (2015) E. Berti et al., Class. Quant. Grav. 32, 243001 (2015), eprint 1501.07274.
 Yunes et al. (2016) N. Yunes, K. Yagi, and F. Pretorius, Phys. Rev. D94, 084002 (2016), eprint 1603.08955.
 Mars (1999) M. Mars, Class. Quant. Grav. 16, 2507 (1999), eprint grqc/9904070.
 Nollert (1996) H.P. Nollert, Phys. Rev. D 53, 4397 (1996), URL https://link.aps.org/doi/10.1103/PhysRevD.53.4397.
 Tufts and Kumaresan (1982) D. W. Tufts and R. Kumaresan, Proceedings of the IEEE 70, 975 (1982), ISSN 00189219.
 Berti et al. (2007) E. Berti, V. Cardoso, J. A. Gonzalez, U. Sperhake, M. Hannam, S. Husa, and B. Brügmann, Phys. Rev. D 76, 064034 (2007), eprint grqc/0703053.
 Baibhav et al. (2017) V. Baibhav, E. Berti, V. Cardoso, and G. Khanna (2017), eprint 1710.02156.
 London et al. (2014) L. London, D. Shoemaker, and J. Healy, Phys. Rev. D 90, 124032 (2014), URL https://link.aps.org/doi/10.1103/PhysRevD.90.124032.
 Kamaretsos et al. (2012b) I. Kamaretsos, M. Hannam, and B. S. Sathyaprakash, Physical Review Letters 109, 141102 (2012b), eprint 1207.0399.
 Thrane et al. (2017) E. Thrane, P. Lasky, and Y. Levin, ArXiv eprints (2017), eprint 1706.05152.
 GómezLobo (2016) A. G.P. GómezLobo, Class. Quant. Grav. 33, 175005 (2016), eprint 1602.08075.
 Owen (2009) R. Owen, Phys. Rev. D80, 084012 (2009), eprint 0907.0280.
 Scheel et al. (2009) M. A. Scheel, M. Boyle, T. Chu, L. E. Kidder, K. D. Matthews, and H. P. Pfeiffer, Phys. Rev. D79, 024003 (2009), eprint 0810.1767.
 Baker and Campanelli (2000) J. Baker and M. Campanelli, Phys. Rev. D 62, 127501 (2000), eprint grqc/0003031.
 Campanelli et al. (2009) M. Campanelli, C. O. Lousto, and Y. Zlochower, Phys. Rev. D79, 084012 (2009), eprint 0811.3006.
 Owen (2010) R. Owen, Phys. Rev. D81, 124042 (2010), eprint 1004.3768.
 Handmer and Szilágyi (2015) C. J. Handmer and B. Szilágyi, Class. Quant. Grav. 32, 025008 (2015), eprint 1406.7029.
 Handmer et al. (2015) C. J. Handmer, B. Szilágyi, and J. Winicour, Class. Quant. Grav. 32, 235018 (2015), eprint 1502.06987.
 Handmer et al. (2016) C. J. Handmer, B. Szilágyi, and J. Winicour, Classical and Quantum Gravity 33, 225007 (2016), eprint 1605.04332.
 Price et al. (2016) R. H. Price, S. Nampalliwar, and G. Khanna, Phys. Rev. D 93, 044060 (2016), eprint 1508.04797.
 Cardoso et al. (2016) V. Cardoso, E. Franzin, and P. Pani, Physical Review Letters 116, 171101 (2016), eprint 1602.07309.
 Baumgarte and Shapiro (2010) T. W. Baumgarte and S. L. Shapiro, Numerical Relativity: Solving Einstein’s Equations on the Computer (Cambridge University Press, New York, 2010).
 Backdahl and Valiente Kroon (2010a) T. Backdahl and J. A. Valiente Kroon, Phys. Rev. Lett. 104, 231102 (2010a), eprint 1001.4492.
 Backdahl and Valiente Kroon (2010b) T. Backdahl and J. A. Valiente Kroon, Annales Henri Poincare 11, 1225 (2010b), eprint 1005.0743.
 Backdahl and Valiente Kroon (2012) T. Backdahl and J. A. Valiente Kroon, J. Math. Phys. 53, 042503 (2012), eprint 1111.6019.
 Stephani et al. (2009) H. Stephani, D. Kramer, M. MacCallum, C. Hoenselaers, and E. Herlt, Exact Solutions of Einstein’s Field Equations, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2009), ISBN 9781139435024, URL https://books.google.com/books?id=SiWXP8FjTFEC.
 Ferrando et al. (2001) J. J. Ferrando, J. A. Morales, and J. A. SÃÂ¡ez, Classical and Quantum Gravity 18, 4939 (2001), URL http://stacks.iop.org/02649381/18/i=22/a=315.
 (42) A. G.P. GómezLobo, Private communication.
 Khanna and Price (2017) G. Khanna and R. H. Price, Phys. Rev. D 95, 081501 (2017), eprint 1609.00083.
 Vishveshwara (1970) C. V. Vishveshwara, Nature 227, 936 (1970).
 Bishop et al. (1997) N. Bishop, R. Gomez, L. Lehner, and J. Winicour, ArXiv General Relativity and Quantum Cosmology eprints (1997), eprint grqc/9705033.
 Bondi et al. (1962) H. Bondi, M. G. J. van der Burg, and A. W. K. Metzner, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences 269, 21 (1962), ISSN 00804630, URL http://www.jstor.org/stable/2414436.
 Araneda and Dotti (2015) B. Araneda and G. Dotti, Class. Quant. Grav. 32, 195013 (2015), eprint 1502.07153.
 Yang et al. (2015) H. Yang, A. Zimmerman, and L. Lehner, Phys. Rev. Lett. 114, 081101 (2015), eprint 1402.4859.
 Lousto and Whiting (2002) C. O. Lousto and B. F. Whiting, Phys. Rev. D66, 024026 (2002), eprint grqc/0203061.
 Teukolsky (2015) S. A. Teukolsky, Class. Quant. Grav. 32, 124006 (2015), eprint 1410.2130.
 Ori (2003) A. Ori, Phys. Rev. D67, 124010 (2003), eprint grqc/0207045.
 East et al. (2014) W. E. East, F. M. RamazanoÄlu, and F. Pretorius, Phys. Rev. D89, 061503 (2014), eprint 1312.4529.
 Ionescu and Klainerman (2014) A. D. Ionescu and S. Klainerman (2014), eprint 1412.5679.
 Abbott et al. (2016e) B. P. Abbott, R. Abbott, T. D. Abbott, M. R. Abernathy, F. Acernese, K. Ackley, C. Adams, T. Adams, P. Addesso, R. X. Adhikari, et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 061102 (2016e), URL http://link.aps.org/doi/10.1103/PhysRevLett.116.061102.
 (55) http://www.blackholes.org/waveforms.
 Szilágyi et al. (2009) B. Szilágyi, L. Lindblom, and M. A. Scheel, Phys. Rev. D 80, 124010 (2009), eprint 0909.3557.
 Lindblom and Szilágyi (2009) L. Lindblom and B. Szilágyi, Phys. Rev. D 80, 084019 (2009), eprint arXiv:0904.4873.
 Hemberger et al. (2013) D. A. Hemberger, M. A. Scheel, L. E. Kidder, B. Szilágyi, G. Lovelace, N. W. Taylor, and S. A. Teukolsky, Class. Quant. Grav. 30, 115001 (2013), eprint 1211.6079, URL http://stacks.iop.org/02649381/30/i=11/a=115001.
 Scheel et al. (2015) M. A. Scheel, M. Giesler, G. Hemberger, D. A. andÂ· Lovelace, K. Kuper, M. Boyle, and L. E. Szilágyi, B. andÂ· Kidder, Class. Quant. Grav. 32, 105009 (2015), eprint 1412.1803.
 Rinne et al. (2007) O. Rinne, L. Lindblom, and M. A. Scheel, Class. Quant. Grav. 24, 4053 (2007), URL http://stacks.iop.org/02649381/24/4053.
 Lovelace et al. (2016) G. Lovelace et al., Class. Quant. Grav. 33, 244002 (2016), eprint 1607.05377.
 Press et al. (2007) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes 3rd Edition: The Art of Scientific Computing (Cambridge University Press, New York, NY, USA, 2007), 3rd ed., ISBN 0521880688, 9780521880688.
 Lindblom et al. (2006) L. Lindblom, M. A. Scheel, L. E. Kidder, R. Owen, and O. Rinne, Class. Quant. Grav. 23, S447 (2006), eprint grqc/0512093.
 Barkett and Scheel (2017) K. Barkett and M. Scheel (2017), forthcoming.
 Apostolatos (1995) T. A. Apostolatos, Phys. Rev. D 52, 605 (1995), URL https://link.aps.org/doi/10.1103/PhysRevD.52.605.
 (66) C. Handmer, Private communication.
 Rinne et al. (2009) O. Rinne, L. T. Buchman, M. A. Scheel, and H. P. Pfeiffer, Classical and Quantum Gravity 26, 075009 (2009), URL http://stacks.iop.org/02649381/26/i=7/a=075009.
 Regge and Wheeler (1957) T. Regge and J. A. Wheeler, Phys. Rev. 108, 1063 (1957), URL https://link.aps.org/doi/10.1103/PhysRev.108.1063.
 Zerilli (1970a) F. J. Zerilli, Phys. Rev. Lett. 24, 737 (1970a), URL https://link.aps.org/doi/10.1103/PhysRevLett.24.737.
 Zerilli (1970b) F. J. Zerilli, Phys. Rev. D 2, 2141 (1970b), URL https://link.aps.org/doi/10.1103/PhysRevD.2.2141.
 Moncrief (1974) V. Moncrief, Annals Phys. 88, 323 (1974).
 Boyle and Mroue (2009) M. Boyle and A. H. Mroue, Phys. Rev. D80, 124045 (2009), eprint 0905.3177.
 Taylor et al. (2013) N. W. Taylor, M. Boyle, C. Reisswig, M. A. Scheel, T. Chu, L. E. Kidder, and B. Szilágyi, Phys. Rev. D88, 124010 (2013), eprint 1309.3605.
 Pürrer (2016) M. Pürrer, Phys. Rev. D93, 064041 (2016), eprint 1512.02248.
 Aasi et al. (2015) J. Aasi et al. (LIGO Scientific), Class. Quant. Grav. 32, 074001 (2015), eprint 1411.4547.
 Ruiz et al. (2008) M. Ruiz, R. Takahashi, M. Alcubierre, and D. Nunez, Gen. Rel. Grav. 40, 2467 (2008), eprint 0707.4654.
 Berti et al. (2007) E. Berti, J. Cardoso, V. Cardoso, and M. Cavaglia, Phys. Rev. D76, 104044 (2007), eprint 0707.1202.
 Bohn et al. (2016) A. Bohn, L. E. Kidder, and S. A. Teukolsky, Phys. Rev. D94, 064008 (2016), eprint 1606.00437.
 Bohn et al. (2015) A. Bohn, W. Throwe, F. Hébert, K. Henriksson, D. Bunandar, M. A. Scheel, and N. W. Taylor, Class. Quant. Grav. 32, 065002 (2015), eprint 1410.7775.
 (80) Spectral einstein code, https://www.blackholes.org/code/SpEC.html.
 Plebanski and Demianski (1976) J. Plebanski and M. Demianski, Annals of Physics 98, 98 (1976), ISSN 00034916, URL http://www.sciencedirect.com/science/article/pii/0003491676902402.
 Aliev et al. (2008) A. N. Aliev, H. Cebeci, and T. Dereli, Phys. Rev. D 77, 124022 (2008), URL https://link.aps.org/doi/10.1103/PhysRevD.77.124022.
 Griffiths and Podolsky (2005) J. B. Griffiths and J. Podolsky, Class. Quant. Grav. 22, 3467 (2005), eprint grqc/0507021.