# Mergers of nonspinning black-hole binaries: Gravitational radiation characteristics

###### Abstract

We present a detailed descriptive analysis of the gravitational radiation from black-hole binary mergers of nonspinning black holes, based on numerical simulations of systems varying from equal-mass to a 6:1 mass ratio. Our primary goal is to present relatively complete information about the waveforms, including all the leading multipolar components, to interested researchers. In our analysis, we pursue the simplest physical description of the dominant features in the radiation, providing an interpretation of the waveforms in terms of an implicit rotating source. This interpretation applies uniformly to the full wave train, from inspiral through ringdown. We emphasize strong relationships among the modes that persist through the full wave train. Exploring the structure of the waveforms in more detail, we conduct detailed analytic fitting of the late-time frequency evolution, identifying a key quantitative feature shared by the modes among all mass ratios. We identify relationships, with a simple interpretation in terms of the implicit rotating source, among the evolution of frequency and amplitude, which hold for the late-time radiation. These detailed relationships provide sufficient information about the late-time radiation to yield a predictive model for the late-time waveforms, an alternative to the common practice of modeling by a sum of quasinormal mode overtones. We demonstrate an application of this in a new effective-one-body-based analytic waveform model.

###### pacs:

04.25.Dm, 04.30.Db, 04.70.Bw, 04.80.Nn 95.30.Sf, 95.55.Ym 97.60.Lf## I Introduction

The final merger of two black holes (BHs) having comparable masses will produce an intense burst of gravitational radiation, and is expected to be one of the strongest sources in the gravitational-wave sky. Mergers of stellar black holes are key targets for ground-based detectors such as LIGO, VIRGO, and GEO600, and knowledge of the merger waveforms is an important component of improving the detectability of such systems. The space-based LISA detector will observe mergers of massive black holes at high signal-to-noise ratios, allowing tests of general relativity in the strong-field, dynamical regime.

Today, numerical relativity (NR) studies are beginning to progress toward a full description of black-hole binary merger systems. For noneccentric inspirals, this space is spanned by seven parameters: the symmetric mass-ratio , and the six combined components of the black holes’ spin vectors. Considerable study has been focused on the fiducial center point of this parameter space, the case of equal-mass nonspinning black-hole mergers. After the series of breakthroughs that ushered in an era of rapid progress in the field Brügmann et al. (2004); Pretorius (2005); Campanelli et al. (2006); Baker et al. (2006a), several investigations assessing the accuracy of the available equal mass waveforms and applying them to data analysis were conducted Baker et al. (2007b); Buonanno et al. (2007a); Baker et al. (2007a, b); Boyle et al. (2007).

In this paper, we undertake a descriptive study of the waveforms generated in the late inspiral and merger of black-hole binaries for the subspace of nonspinning black holes, parametrized only by . Our study is based on a series of numerical simulations, discussed in Sec. III, covering at least the last orbits of nonspinning black-hole binary mergers with mass ratios extending to 6:1 (). Several of the simulations presented here have already been applied in a recent paper, focusing on the development of a faithful analytic waveform model Buonanno et al. (2007b). Here we provide details of these and additional simulations, together with considerable analysis, focused on providing a qualitative and quantitative picture of how the waveforms from nonspinning black-hole mergers depend on . Nonspinning black-hole binary merger waveforms were previously examined in Ref. Berti et al. (2007), but our analysis is novel and complementary to that work. Our descriptive presentation puts emphasis on the relationships between waveforms from the different mass-ratio cases and different harmonic modes, with references to Ref. Berti et al. (2007) where related observations have been made. Our approach to describing the inspiral-merger-ringdown transition is particularly distinct, founded in a uniform approach that describes all stages of this process in similar terms, and ultimately suggesting a complementary physical picture.

Black-hole-binary merger waveforms have been noted for their “simplicity.” For the nonspinning systems the simple physics of the coalescence is exposed by a spherical harmonic decomposition of the waveforms. In Sec. IV we walk readers through the basic features of the radiation, characterizing amplitude and phase evolution of the multipolar components, and discussing relationships among the simulations representing different mass ratios, and among the multipolar components of each simulation. As we analyze the waveforms we develop a conceptual interpretation of the basic waveform features. In this interpretation we consider the structure of an implicit rotating source, which could have generated the measured radiation through its rotational motion. This allows a uniform interpretation that applies throughout the coalescence process: inspiral, merger and ringdown.

In Sec. V, we examine the strong final burst of radiation beginning before the formation of a common horizon. We quantitatively describe the phasing in terms of an analytic model, based on a continuous, monotonically increasing frequency. We find, in particular, that the peak rate of change in frequency, appropriately scaled, is the same across all modes and mass ratios. We also identify relationships among the mode amplitudes and phases, which are connected to an approximately linear relationship between angular momentum and frequency: . We interpret these relationships in terms of the implicit source.

Finally, In Sec. VI, we demonstrate the utility of what we have learned in our waveform characterization by applying some of the quantitative features we have uncovered in a new variation on the analytic waveform model in Buonanno et al. (2007b), which was based on the Effective-One-Body (EOB) resummation of the Post-Newtonian(PN) approximation to inspiral dynamics Buonanno and Damour (1999). In particular, we provide a distinct late-time waveform model, alternative to the common “spectroscopic” model Dreyer et al. (2004); Berti et al. (2006) based on sums of quasinormal mode overtones.

## Ii Overview

We begin with some examples of gravitational strain waveforms as they might be observed by gravitational-wave instruments. In observational work, and PN analysis, it is customary to describe the radiation in terms of gravitational-wave strain, . In representing the strain, it is convenient to combine the two real waveform polarization components, and , into one complex strain waveform,

(1) |

We decompose the strain waveforms measured on a sphere of radius , into spin-weighted spherical harmonic components, . The details of the decomposition, and how the waveform information is extracted from the numerical simulations, are given in Appendix A.

The waveforms in this section are aligned in time and phase so that the maximum strain amplitude occurs at . The remaining figures of this paper will be aligned in a similar way, but with marking the time of peak (2,2) mode energy flux, (unless stated otherwise).

Fig. 1 shows waveforms from mergers of nonspinning black holes for various mass ratios, as observed at distance on the rotational/orbital axis of the system. The figure shows for each of the four mass ratios 1:1, 2:1, 4:1, and 6:1. For these observers the observed waveforms will be circularly polarized, so that is out of phase with . We use units in which and and express both time and spatial distances in terms of the total mass , where .

More typically, the observer will not be located on the system’s
orbital axis. The left panel of Fig. 2 shows
for the 4:1 case. The strain is measured at an azimuthal angle
of and various
inclinations^{1}^{1}1The inclination angle is
defined here as the angle between the line-of-sight with respect to
the detector and the orbital axis of the binary. This is the same
angle referred to as “inclination” in the PN/NR literature, and most
equations are constructed using that definition. However, the
astronomical literature has often defined inclination to be the angle
between the line of sight and the orbital plane of the binary,
resulting in a inconsistency.. The detailed shapes of
the waveforms change as the system is reoriented so that the observer
moves off the system’s rotational axis. For larger inclinations
(closer to being viewed edge-on) there are notable modulations at half
the base gravitational-wave frequency.

The right panel of Fig. 2 shows for different mass ratios oriented at an inclination of and an azimuthal angle of . For this orientation, constitutes the full strain waveform. For larger mass ratios, the lower frequency modulation increases.

For gravitational-wave observations of sufficiently strong binary black-hole sources, the types of differences shown in Fig. 2 could be exploited to estimate the inclination and mass-ratio of the source system.

For observational purposes, the combined waveform information encoding all possible source orientations can be conveniently represented in terms of spin-weighted spherical harmonic components [see Eqs.(27)-(29)], providing a neat description of the leading waveform features, The multipolar decomposition is even more valuable as a tool for exposing the hallmark simplicity of the merger radiation. The readily apparent simplicity in the waves viewed from the system’s orbital axis in Fig. 1 extends to each of the spherical harmonic components. Viewed off axis, these components linearly combine to yield the more complex appearance of the waveforms in Fig. 2.

This characterization of the gravitational radiation from a merging black-hole binary in terms of circular polarization was first recognized in the Lazarus project studies Baker et al. (2002a). In this picture, the radiation can be represented by a slowly varying amplitude and a polarization angle; see Eq. (2) below. This description relies on how the radiation appears to distant observers located on the rotational axis of the system. Other observers will typically see elliptically polarized waves, having a generally simple pattern that conforms to the rotational nature of the source. In the equatorial plane, the radiation reduces to the plus polarization, corresponding to the observer seeing no circulation in the source. Looking along the negative rotation axis, the observer sees circular polarization with the opposite helicity.

Each of the spherical harmonic waveform components exhibits circular polarization with steadily varying phase and amplitude, providing a natural framework for developing a practical and intuitive understanding of binary black-hole merger radiation. Our basic waveform description, in Sec. IV, and the more detailed analysis that follows, is based on this spherical harmonic decomposition.

As we describe the waveforms, we will also suggest a heuristic
interpretation of what the radiation tells us about the motion and
structure of the binary black-hole source. In the weak-field
description of radiation from a rotating object, the multipolar
waveform components of the gravitational radiation can be associated
with dynamics of multipolar moments of the radiating source
Thorne (1980). It is useful, in conceptualizing the full
coalescence radiation from inspiral through merger and ringdown, to
think of the multipolar radiation description as providing information
about the motion of a changing source object, described as a sum of
several multipolar mass moments. This source object is what we will
interpret as an effective rigid rotator radiation source, with a
slowly changing structure. We refer to this as the *implicit
rotating source* (IRS).

In the process of coalescence, the source begins as a separated black-hole binary system and ends as a single distorted black hole. For nonspinning binary mergers, numerical and PN results consistently indicate that the radiation is circularly polarized, in the sense first recognized in the Lazarus project studies Baker et al. (2002a), not only in the inspiral, but uniformly through the merger and ringdown. In our conceptual source description, this pattern of circular polarization is consistent with radiation generated by rotational motion of each source multipolar moment, where the polarization phase is tied to the instantaneous orientation of the source. Similarly, we think of the amplitude of the radiation multipole as related, through some generalization of the quadrupole formula, to the magnitude of the implicit source multipole.

We can write each multipolar component in a specific polar form natural for circularly polarized radiation:

(2) |

Each amplitude is expected to be a slowly varying
function of time and can be conceptually considered as a function only
of the magnitude of the source multipolar moment and its rotational
frequency. The additional sign for some cases allows a
consistent interpretation for the component phases and the component
amplitudes. The waveform phase is given here by , but in most of our analysis we refer to ,
which we call the *rotational phase*. In our terms of our
implicit source heuristic, the rotational phase for some particular
mode can be thought of as the azimuthal orientation of the
specified multipolar source component. In the inspiral, where the
source can be considered as a separated binary,
should coincide with the orbital phase, independent of . In
the post-Newtonian expansion Blanchet (2006), all the
defined here agree to at least 2PN order, while
the amplitudes remain real and non-negative. The
equatorial-plane symmetry of these mergers ensures that and that , so that we need consider only modes
for this analysis.

The expansion (2) is not appropriate for the modes. This points to an important caveat to our implicit rotating source interpretation of the radiation, that it applies to the degree that the radiation is circularly polarized. While not strictly vanishing, the waveform components, and other deviations from circular polarization are generally extremely small, and largely unmeasurable at the resolutions of the simulations we study. For the most part, we will not address deviations from circular polarization, and the consequent limitations of our implicit-rotating-source interpretation in this paper, focusing for now on the dominant features of the radiation.

## Iii Simulations

Our analysis is based on four simulations, representing mass ratios 1:1, 2:1, 4:1 and 6:1. Results from the 1:1 and 4:1 simulations have appeared in previous publications (Baker et al. (2007a) and Buonanno et al. (2007b), respectively). More recently, higher accuracy simulations have been presented by other groups for the 1:1 case Husa et al. (2007); Boyle et al. (2007). Our older waveform is sufficiently accurate for our present purpose, to examine the general features of the nonspinning merger waveforms.

Our numerical simulations are carried out with the hahndol evolution code Imbiriba et al. (2004), which uses finite-differencing methods to solve a 3+1 formulation of Einstein’s equations on a Cartesian grid. For initial data we solve the elliptic equation given by Brandt and Brügmann for conformally flat data in which the black holes are represented by punctures Brandt and Brügmann (1997). This is performed numerically using the multigrid solver amrmg Brown and Lowe (2005), which is second-order-accurate but tuned to give truncation errors typically much smaller than those produced by the evolution code. The momentum parameters are chosen according to the 2PN-accurate quasicircular approximation given by Kidder Kidder (1995), which has been found to result in low eccentricity. We evolve these data using the moving puncture method Campanelli et al. (2006); Baker et al. (2006c) with a modified version of the Baumgarte-Shapiro-Shibata-Nakamura equations Shibata and Nakamura (1995); Baumgarte and Shapiro (1998). Specifically, as suggested in van Meter (2006), we replaced the conformal factor variable with , which vanishes at the punctures. Further, we added the constraint-damping terms suggested in Duez et al. (2004), and the dissipation terms suggested in Kreiss and Oliger (1973); Hübner (1999). For the gauge we use the specific 1+log lapse and Gamma-freezing shift conditions recommended for moving punctures in van Meter et al. (2006).

Accurate simulations require adequate spatial resolution near the black holes (length scales ) as well as in the wave zone where the gravitational waves are extracted (length scales ). To this end, the grid has multiple refinement levels, determined adaptively near the black holes, but fixed in regions farther away (typically, ) where the waves are extracted; all grid refinement is handled within the framework of the software package paramesh MacNeice et al. (2000). The adaptive mesh refinement criterion near the black holes is designed to keep the scale of the square root of an invariantly defined curvature component, known as the Coulomb scalar Beetle et al. (2005); Burko et al. (2006), roughly constant with respect to the grid spacing. Interpolation in guard-cells between refinement regions is fifth-order-accurate, coupling with differencing stencils to yield at least fourth-order accuracy in the bulk.

Spatial derivatives are taken by sixth-order-accurate differencing
stencils, with the exception of advection derivatives, which are
handled by fifth-order-accurate mesh-adapted differencing for greater
stability Baker and van Meter (2005)^{2}^{2}2Sixth-order center-differenced
advection is unstable, and sixth-order lopsided advection is too
costly in terms of paramesh guardcells, which motivated our
particular modification.. Time integration is performed with a
fourth-order Runge-Kutta algorithm.

The initial configurations of the simulations we analyze are given in Table 1. In each case, the initial separation was chosen to be large enough to result in at least five orbits. The finest resolution, , ranged from to , as required to adequately resolve the black hole with the smaller mass in each case. The outer boundary was typically at , far enough away to prevent reflections from reaching the wave-extraction region during the simulation.

Mass ratio | ||||||||||

1:1 | 0.4872 | 0.4872 | 10.800 | 0.09118 | 0.9847 | 0.9907 | 1.0005 | 0.2500 | ||

2:1 | 0.3202 | 0.6504 | 8.865 | 0.09330 | 0.8271 | 0.9889 | 0.9989 | 0.9990 | 0.2228 | |

4:1 | 0.1890 | 0.7900 | 8.470 | 0.06957 | 0.5893 | 0.9929 | 1.0003 | 1.0004 | 0.1601 | |

0.1890 | 0.7900 | 8.470 | 0.06957 | 0.5893 | 0.9929 | 1.0003 | 1.0004 | 0.1601 | ||

0.1890 | 0.7900 | 8.470 | 0.06957 | 0.5893 | 0.9930 | 1.0003 | 1.0005 | 0.1601 | ||

6:1 | 0.1338 | 0.8490 | 8.003 | 0.05559 | 0.4449 | 0.9942 | 1.0000 | 1.0001 | 0.1226 |

We have measured the individual black-hole masses using the
*apparent-horizon mass* , a quantity calculated from the
area of each hole’s horizon, which we locate using the
AHFinderDirect code Thornburg (2004). From these
horizon masses, we calculate the symmetric mass ratio . This gives the most precise specification of the
actual mass ratio attained in our simulations. In the text we will
refer to the simulations by the mass ratio (*e.g.* 4:1).

We define the total, infinite-separation, mass of the system as an analogue for the total rest mass parameter used in PN studies. We measure in two ways:

(3) |

the sum of the individual BH horizon masses, and

(4) |

defined as the difference between , the total energy of the initial data, and the (negative) binding energy of the binary. The binding energy is estimated from an effective-one-body PN treatment Buonanno and Damour (1999), given the initial angular momentum . The result shows a very close correspondence between and , with differences at the level . For the rest of this paper, we use , except for the 1:1 simulation data, where was not available for technical reasons.

In interpreting the late-time radiation, it is valuable to know the mass and spin of the final Kerr black hole formed by the merger. We discuss the state of the final black hole, determined consistently by several means, in Appendix C.

For the 4:1 mass ratio case, we have carried out runs at three different resolutions in order to assess the quality of the simulations. The convergence of the constraints and waveforms is discussed in Appendix B.

The most important products of our simulations are the gravitational radiation waveforms, which we extract from the evolved simulation data as explained in Appendix A. Strain-rate waveforms for the 4:1 case at various resolutions are shown in Fig. 3, where the times and phases have been shifted to agree at the moment of peak energy flux, as is generally done in our analysis below. We can get some measure of the error in the waveforms by comparing the difference between the high and medium resolution simulations. Fig. 4 shows the relative differences in amplitudes, scaled by the high-resolution result. Ignoring the high frequency noise, the (2,2)-mode differences (upper panel) indicate a combination of a secular amplitude difference and a sinusoidal effect, which results from the combination of the eccentricity in the orbital dynamics and the difference in peak time due to limited resolution. These combine to give differences generally at the 3% level, somewhat smaller at late times. The eccentricity plays less of a role in the (4,4) differences (lower panel), as the relative secular error is much ( times) larger. Sinusoidal eccentricity effects are also visible in the phasing error [see Fig. 5]. Overall, we find waveform amplitude and phase errors to be consistent with between fourth- and fifth-order convergence [see Appendix B].

Assuming fourth-order convergence, and using Richardson extrapolation, our nominal expectation for these simulations leads to an error estimate for the high-resolution simulation applied in our analysis of the difference shown in Fig. 4. To be conservative, we could instead assume second-order convergence, which would lead to an error estimate of the difference shown in Fig. 4.

The errors for the 2:1 case should be comparable to the 4:1 case. The resolution for the 6:1, scaled by the smaller black hole’s mass is about 15% lower than lowest resolution 4:1 simulation, suggesting errors eight times larger, if we conservatively assume fourth-order convergence and that the errors around the smaller black hole dominate. The errors for the 1:1 mass ratio case are discussed in Ref. Baker et al. (2007a).

## Iv Descriptive analysis of waveforms

In this section we provide a descriptive analysis of the waveforms from our simulations. We try to serve two purposes in analyzing the radiation. In the first place, we are hoping to provide material for gravitational-wave observers, and others outside the field of numerical relativity, which makes clear some of the general characteristics of the radiation from these mergers. Beyond that, we also push the analysis in more detail, hoping to generate deeper insight into the physics which generates the radiation. Through this analysis we explore the similarities and differences for the various mass-ratio simulations, and among the different multipole components of the 4:1 case case. In this way, we examine the waveform amplitudes and energy, and the waveform phasing. As we proceed, we will interpret the results in terms of our implicit-rotating-source model, building up a heuristic description that applies through the inspiral, merger and ringdown of the binary.

Following a brief discussion of strain rate in IV.1, we study the waveform amplitudes and the associated energetics of the merger in IV.2. In IV.3 we address the polarization phase of several waveform modes, relating them to a common implicit source phase. Next, in Sec. V, we will examine the late-time frequency evolution and the relation to amplitude in more quantitative detail through the inspiral-merger-ringdown transition.

### iv.1 Strain-rate

In the Introduction, we motivated the spherical harmonic phase and
amplitude waveform decomposition with a discussion in terms of strain
. In analyzing our numerical simulation results, however, we can
work more directly with the *strain rate*
[see Appendix A for a more detailed
discussion]. As with the strain decomposition
(2), we will expand the strain rate as

(5) |

with real and non-negative.^{3}^{3}3A similar
expression, differing from (2) by an overall
sign, would be equally appropriate for direct interpretation of
numerically derived waveforms. Direct differentiation of
(2) reveals the relationships between the phases
and amplitudes defined in (2) and
(5), with while
differs from only at 2.5PN
order. Note that the differentiation produces a phase shift of
, so that strain-rate phases should be defined in the slightly
unconventional form (5) if we wish to
preserve the property that all are equal to orbital
phase in the limit of well-separated binaries. This representation
allows a more meaningful comparison of the phases of different
multipolar modes, ensuring that instantaneous phase corresponds to the
orientation of the binary system during the inspiral.

The strain-rate amplitude is directly related to the radiation power for the mode by Eq. (32), . Henceforth we shall use the strain-rate-derived rotational phase , rather than . We also limit our presentation to the modes, as equatorial symmetry implies and .

Unless otherwise indicated, in the remainder of the paper, the time axis of each plot will be shifted so that the peak of (and hence of the strain-rate amplitude ) occurs at . As the (2,2) mode is strongly dominant, this will closely approximate the peak time of the total .

### iv.2 Amplitude and energetics

We first study wave amplitudes across modes and mass ratios. Since we are examining strain rate waveforms, the modal energy flux is effectively equivalent to the square of the mode amplitude, as in Eq. (32). Preferring the most physical language we will express our modal amplitude comparisons as energy flux comparisons. In terms of the implicit rotating source model, we can think of the energy carried by the radiation as energy lost by the source.

In Fig. 6, we plot the actual peak values of the
dominant (2,2) energy flux contribution from
Eq. (32) as a function of symmetric mass ratio
. As this mode contribution is proportional to
, and we expect (the test-particle limit), we fit it to
a quadratic-quartic form, obtaining the fit^{4}^{4}4A quadratic-cubic
form is equally plausible, but fit the numerical data worse in this
case.:

(6) | |||||

We also plot the peak of the total energy flux for each mass ratio, scaled by one-half, since the (2,2) and (2,-2) modes contribute equally to . The difference between and increases as , reflecting the increased importance of other modes for unequal masses.

Aside from the value of the radiation maxima, it is interesting to see how the radiation power evolves in time near the peak. In Fig. 7, we show shapes of the dominant (2,2) contributions to the peak energy fluxes (32) for each mass ratio. These are scaled to the same peak height to allow shape comparison, and shifted in time so that the peaks of are aligned. We note the striking similarity of the peak shape and duration across all mass ratios. During the late-inspiral phase, the more extreme mass ratios appear to radiate more energy; however, since we have normalized each curve by peak height, this only means that the equal-mass binary experiences a steeper climb to its peak power rate. Nevertheless, the different mass ratios follow similarly shaped tracks approaching merger, and differences have been effaced by before the peak power. The post-peak portion of the curve is determined by the dominant quasinormal mode (QNM) damping time, which varies only slightly with underlying Kerr spin for moderate spins [see, for instance, tables in Berti et al. (2006)]. The lower spins on the black holes formed by smaller- mergers should cause power to fall off faster in these cases. The inset in Fig. 7 shows the difference in fall-off rate relative to the equal-mass case from after the peak.

In Fig. 8, we concentrate on the 4:1 mass ratio, plotting the strongest modal contributions to the flux, again scaled to the same peak height. Here we have applied only an overall time shift, so that the total flux peak is at . The energy profiles of the radiation burst in modes with are similar, all peaking at approximately the same time. The subdominant modes are relatively stronger in the burst than in the inspiral, so that they show up here as weaker in the approach to the peak even after scaling to the peak radiation. We note that the subdominant modes with , and , are particularly similar in this regard. Similarity in modes, and distinction in the other modes is a general feature of the bursts in several ways.

The shape of the peaks with are particularly distinct. The (2,1) mode peaks particularly late, and the burst is much stronger than the inspiral. We note that the (3,2) mode shows a double-bump in its contribution to the energy flux. From Fig. 9, this appears to be robust in its gross shape over resolution and extraction radius. From Fig. 10 we note, however, that the extent of this double-bump effect is very dependent on mass ratio; it is less in the 6:1 case, and not evident at all in the 2:1 case. Later we will also note irregularities in the late-time frequency evolution of this mode, apparently indicating a deviation from circular polarization in this case.

In both Fig. 7 and Fig. 8 we have normalized the mode-flux peaks for the purposes of shape comparison. It is also important to understand the relative strengths of each mode. We show in Fig. 11 the relative mode contributions to for several dominant modes over the final inspiral and merger of the 4:1 case. In the merger-ringdown peaks, as in the inspirals, the modes dominate the energy flux, followed by the modes. More discussion of the relative mode strengths for general nonspinning mergers is given in Ref.Berti et al. (2007).

In addition to the direct mode contributions, we plot the PN-derived energy “partitions” – the power emitted in each mode as a fraction of the total, according to the leading-order “restricted” PN expressions found in Eqs. (30)-(36) of Buonanno et al. (2007b), where the underlying orbital frequency was derived from the (2,2) mode. These are shown in Fig. 11 using dashed lines. Until near the peak time, where we discontinue the PN curves, we see that the partitioning tracks the numerics very well for all modes except (2,1), which grows visibly faster than the restricted amplitude prediction after before merger. A similar study Buonanno et al. (2007b) showed that the restricted (leading-order) approximation for the amplitudes consistently overestimates the strength of the radiation. This shows that the leading-order partitioning of energy can provide a simple, but more accurate, approximation of the mode amplitudes. We will take advantage of this in Sec. VI to provide more accurate amplitudes in a variation on the analytic EOB-based waveform model studied in Buonanno et al. (2007b).

### iv.3 Waveform phasing

In gravitational-wave observations, the waveform phase provides most of the time variation in the signals, and consequently is critically important in encoding observable information about the source. Here we consider the phasing of the leading spherical harmonic waveform components. Following the discussion above, we conceptually interpret each waveform phase as describing the orientation of a particular multipole of an implicit rotating source of the gravitational waves.

Direct comparative analysis of phases provides a stronger probe of the phase relations among the multipolar modes than the comparative analysis of frequencies conducted in a number of previous studies of numerical simulations. Here we will discuss waveform phasing in terms of as it appears in (5), which we will compute from each strain-rate mode . As noted in Sec. II, we expect all phases to agree in the large-separation limit.

We first compute the strain-rate waveform phase using the conventional decomposition:

(7) |

with real and non-negative. Then, setting this equal to (5), and solving for the rotational phase , we find

(8) |

The term results from the factor in (5), while the factor there produces the term for . The terms express the ambiguity in defining waveform phase. Considering any mode in isolation leads to an -fold degeneracy in the associated rotational phase. We resolve this degeneracy by choosing the pair that yields the closest consistency between and at early times (near ). We then determine the remaining for closest early consistency with . This gives us a phase for each mode that can be interpreted as the rotational phase of the implicit rotating source that produced that component of the radiation.

Fig. 12 shows the rotational phase from several modes of the strain-rate waveforms from the highest-resolution () 4:1 run, together with the rotational phase calculated from the tracks of the punctures. The left panel shows that for the inspiral portion of the evolution, all waveform phases, extracted at , agree extremely well, except for the and modes, which differ by a significant part of a radian. The relative difference of each mode from the mode is shown in the right panel. The differences between the modes are , much smaller than the ambiguity in defining rotational phase from the waveforms.

Note that the early part of the waveforms contain, by far, the longest wavelength radiation present, suggesting a greater potential for problems caused by extracting the waveforms too close to the source, not yet in the wave zone. For the relatively deviant and modes, we performed Richardson extrapolation with respect to extraction radius , using the values extracted at and and assuming an error for each mode. These Richardson-extrapolated phases are also shown in the right panel of Fig. 12, subtracted from the (2,2) rotational phase. Richardson extrapolation evidently reduces the early phase deviations in these modes considerably.

The rotational phases calculated from the modes agree to within 0.025 radians during the inspiral for several hundred before the peak , effectively identical within the uncertainties of the numerical approach. This phase agreement is consistent with expectations based on the PN analysis, for which all phases should agree up to 2PN order. These rotational phases also agree to within about radians with the coordinate-dependent rotational phases measured from the puncture tracks, after shifting the puncture track phase by an overall factor of . Heuristically, we can think of each radiation waveform mode as having been generated by the rotation of its own implicit source component. The phase agreement would then be interpreted as indicating that these implicit source components remain aligned through the inspiral. This is to be expected for a system which can be effectively described as an orbiting pair of point particles.

It is, perhaps, more remarkable that a very tight agreement among the mode persists throughout the merger and even into the ringdown, remaining within about 1 radian until after the merger, when the amplitude has already diminished significantly. According to our interpretation, this phase agreement suggests that a significant portion of the implicit rotating radiation source maintains some structural integrity throughout the coalescence. That is, the implicit source we have considered appears to exhibit considerable “rigidity” through merger. This is only possible because of the close relationship among the fundamental quasinormal ringdown frequencies [see Sec. V.1], mimicking the harmonic frequency relationship that holds during the inspiral. For the modes this quasinormal frequency relation does not hold and the phases must separate in the merger. In terms of our implicit source picture, these components of the source seem to shear away from the main source structure to rotate at a faster rate.

Heuristically, the puncture motion is strongly tied to the rotation of the implicit source for most of the evolution. During this period, it is natural to think of the implicit source as an inspiralling pair of pointlike objects moving on timelike world lines. At late times the orientation phase angle of the puncture track disassociates from the waveform rotational phase. The punctures veer away from the implicit source at a late times as they fall into the final black hole. At this point, though we can continue to consider an implicit rotating radiation source, it no longer makes sense to think of that source as a pair of pointlike objects.

Having compared the phases of different multipolar modes for the 4:1 case, we now consider how the phase evolution depends on mass-ratio. There are various reasonable approaches to comparing the phases among simulations of the different mass-ratio cases. Having established above the rotational phase consistency for the different modes, we will compare only the dominant phases. An obvious approach is to compare phases directly against time, scaled by the total PN mass . In the early-time well-separated limit, however, the leading-order PN analysis indicates that phases for different mass ratios should evolve at similar rates when time is scaled by the chirp mass

Fig. 13 shows the rotational phase computed from the (2,2) mode of strain rate for different mass ratios. In the left panel, we align the rotational phases at an early time and scale time by the chirp mass . For this plot, we shift the rotational phases in (chirp) time so that at , the chirp frequency, which is rotational frequency multiplied by chirp mass, is 0.033 and the rotational phase is 0. Following this approach, we would expect good phase agreement at times sufficiently early that only the leading-order PN effects are significant. However, for the late portion of coalescence that we have simulated, we find that the different mass ratios remain roughly in phase for several hundred before and after , peeling away in order at late times.

In the right panel we compare phases in a manner common for numerical relativity waveform comparisons, we shift the curves in time so that each peak energy flux occurs at , and we rotate the phases so that the phases are 0 at this time. We scale the time by . For the equal-mass case , and for the other mass ratios it is smaller, so all of the curves in the left panel are stretched by at least a factor of in time relative to the curves in the right panel. In the -scaled right panel, the different mass ratios again remain approximately in phase for several hundred before and after . At sufficiently late times, and particularly for small , we might expect this manner of consistent phasing as the evolution of the system eventually [after the innermost stable circular orbit] may become dominated by the course of unstable geodesic trajectories around the larger black hole [or an effective black hole in the effective-one-body (EOB) framework]. In that case the frequency evolves independently of the more strongly -dependent rate of energy or angular momentum loss.

## V Detailed late-time analysis

In Sec. IV we have presented general information about the phasing and amplitudes of the radiation components. Our analysis has stayed close to the standard numerical relativity waveform analysis, though we have emphasized an interpretation in terms of an implicit rotating source model. In this section we go beyond the standard waveform presentation, exploring the radiation with the hope of developing a deeper understanding of the simple characteristics of the radiation as described above. Those features and our heuristic interpretation suggest a new approach to examining the structure of the late-time phasing and apparent relationships between frequency and amplitude evolution.

In Sec. V.1, we examine the phasing again, seeking a quantitative understanding of the late-time evolution of the polarization frequency. We introduce a practical model that captures the merger-ringdown transition without the need for multiple quasinormal mode overtones. We investigate the implications of this model in relating the frequency and amplitude close to merger in Sec. V.2.

This section is more technical than Sec. IV, with some subtle discussion of late-time radiation characteristics. For readers who may wish to jump ahead to Sec. VI, we note two results that we will carry forward: (1) a simple quantification of the peak chirp rate for modes, and (2) the idea that becomes approximately constant at late times, which may serve as summary of relationships between frequency and amplitude near the radiation peak.

### v.1 Waveform frequency evolution

First we will study the phasing in the merger and ringdown in more quantitative detail by comparing the polarization frequency evolution for each mode with a simple empirical model. Based on our heuristic model, we interpret the polarization frequency as corresponding to the rotational frequency component of the implicit source.

Fig. 14 shows the evolution of the rotational frequency for several modes of the 4:1 mass-ratio case. Similar frequency evolution curves are a common feature in papers on numerical relativity waveforms [see, for example, Refs. Baker et al. (2006a); Buonanno et al. (2007a); Berti et al. (2007)]. These curves are time derivatives of the phase evolution curves shown in Fig. 12, zoomed in on the late-time behavior, near the elbows in the phase curves. The striking similarity in phasing for the various modes implies similar frequency evolution, which has been noted in previous studies Schnittman et al. (2008).

At late times, this similarity in frequency is made possible because of a special approximate relationship among the fundamental quasinormal ringing frequencies, that they are nearly equal after dividing by the azimuthal mode number to get what we call the rotational frequency. This has been considered in Mashhoon (1985), which pointed to a connection between the quasinormal ringing frequencies and the frequencies of stable null orbits of a black hole at the “light ring” with frequency . The association extends to charged Kerr-Newman black holes and has been compared with recent precise quasinormal ringing frequency calculations in Berti et al. (2006) and Berti and Kokkotas (2005). Conceptually, this allows us to think of the rotational frequency of the modes at late times as corresponding to the rotational rate of gravitational perturbations orbiting at the “light ring”. This suggests a heuristic description of our implicit rotating source at late times as a gravitational distortion of the forming final black hole which predominantly revolves around the black hole on null orbits at the light ring.

Returning to Fig. 14, we note that the mode is different from all the others, showing two spikes, near and . Comparisons of waveforms extracted at different radii, and from simulations of different resolutions, suggest some sensitivity to extraction radius, but do not suggest that the features will vanish in more accurate simulations or with more distant wave extraction. These anomalies may be related to the unusual shape in the amplitude peaks noted in Fig. 9 above. We will discuss this mode’s behavior further in Sec. VII.

For all other modes the frequency evolution follows a simple smoothly evolving curve, qualitatively similar in each and mass-ratio case [see Fig. 15 below]. In particular, we note that, except for small noise contributions, each curve shows that the frequency increases monotonically, ultimately saturating at a frequency set by the fundamental quasinormal ringdown mode. This monotonic frequency development is a universal characteristic of the radiation from inspiral, through merger, and up to ringdown. In the PN analysis of quasicircular inspiral, this characteristic makes it possible to describe the changing structure of the hardening binary as a function of frequency instead of the more coordinate-specific separation. This allows us, for instance, to write the waveform amplitude as a function of frequency.

Subsequently, we will assume monotonic frequency development throughout the coalescence process. This principle underlies our empirical curve fitting of the frequency evolution, allowing more quantitative analysis of the late-time phasing evolution. In Sec. V.2 we will further apply this idea as we study relationships between late-time frequency and amplitude evolution.

To produce an empirical curve for describing the late-time frequency evolution, we assume that each mode has a monotonically increasing polarization frequency, which approaches the fundamental ringdown frequency at late times. The general expectation that the frequency decays exponentially toward the ringdown frequency suggests that we model frequency evolution based on the hyperbolic tangent function.

Specifically, we will compare frequency evolution of the strain-rate waveforms with a model of the form , where

(9) |

This provides a curve that first grows exponentially, with e-folding time , from some initial frequency , then decays exponentially, with e-folding time , to the final frequency . The presence of the exponent allows the early exponential growth rate to differ from that at late times. The early part of this model is a coarse approximation to the growth in frequency near the end of the inspiral. This approximation must, therefore, fail to fit the data if we look back sufficiently long before the time of peak radiation, but, as we will show, it provides a fair approximation of the approach to the peak. The rate at which the frequency grows, the chirp-rate, increases to some maximum, then decreases to zero on approach to the final frequency . For a more meaningful parametrization, we set

(10) |

With this choice, will peak with value at time . The model then depends on the five parameters , , , and .

As shown in Figs. 15 and 16, we find that fits to the model provide an excellent approximation to the numerical data for the strain-rate rotational frequencies in most significant cases. For the modes of all mass ratios, and for all but one of the significant modes for the representative 4:1 mass-ratio case, we find agreement within a few percent after , with the primary differences coming from apparent noise in the numerical simulations, within the uncertainties in the numerical results.

Fig. 15 shows the comparison for the (2,2) mode rotational frequencies for several mass ratios. A glance at the curves shows that the unequal-mass cases are quite similar to the equal-mass frequency evolution, previously examined in Refs. Buonanno et al. (2007a); Hinder et al. (2007). The dominant difference for the unequal-mass cases is that the final frequency decreases with , consistent with the decrease in the spin of the final black hole produced. We expect to correspond to , where is the fundamental () quasinormal ringing frequency for the specific mode for a black hole with the appropriate spin. In the infinite-mass-ratio limit (), should correspond to half the Schwarzschild quasinormal mode frequency , indicated by the horizontal dashed line in the left panel of Fig. 15.

We show a few examples of polarization frequency curves, for subdominant modes in Fig. 16. The mode is very similar to the modes shown in Fig. 15, as are the other modes (not shown). The mode is of similar shape, also well approximated by our fit.

The quantitative fit results are summarized in Table 2. The error bars are based on statistical fit estimates, also incorporating the ranges of best fit results obtained by varying the fit range starting between and and ending between and . The final frequencies approached in the fit curves in Figs. 15 and 16 were robustly determined by the fits within a fraction of a percent. The frequencies from the fits in Fig. 15 were applied in Table 5 to find final black-hole parameters consistent with those determined by conservation of energy and angular momentum.

Mass ratio | |||||||
---|---|---|---|---|---|---|---|

1:1 | (2,2) | ||||||

2:1 | (2,2) | ||||||

4:1 | (2,2) | ||||||

(2,1) | |||||||

(3,3) | |||||||

(4,4) | |||||||

6:1 | (2,2) |

The peak in the chirp-rate is a particularly significant quantity in determining the shape of curves of our general form . As shown in Table 2, our fits determine to within a few percent in all cases. An interesting relationship is apparent among the fits for all the cases for all values of studied. In the last column of the table, we show the peak chirp-rate values scaled by the mass of the final black hole and the final frequency . In each case where we find , consistent within the fit uncertainties. This scaling makes some sense, since the height of the rise in frequency through the final radiation burst is largely determined by , while the time scale over which this rise occurs seems to be similar when time is scaled by . In Sec. VI we will use this result to predict the phase evolution in an analytic waveform model.

Our model for late-time frequency evolution (9)
describes exponential decay toward at an e-folding rate given by
half our fitting parameter . For all cases, the values of are
within about 30% of . In some cases, the fits for are
rather sensitive to the initial starting time, varying by up to 20%
or 30% in the (2,1) and (4,4) modes of the 4:1 case. At this coarse
level, we note that the values for are similar to the exponential
decay rates for quasinormal ringdown mode *amplitudes* listed in
Table 4. We will consider this relationship
further in Sec. V.2 below.

The other parameters in our fit are , relating to the shape of our fit curve at early times, and , giving the time at which the frequency peak occurs. The parameter is not very precisely determined; as we would expect, it depends sensitively on the starting time of the fit interval, since the early exponential frequency growth is only a coarse approximation of the expected behavior. The values for show that the peaks in generally occur roughly before the total energy peaks at . As was the case for the power peaks in Fig. 8, the chirp-rates of the different spherical harmonic modes peak at slightly different times.

We have supplemented our general implicit rotating source picture with the additional idea that the rotational frequency for each mode grows monotonically, not only in the inspiral, but also through the merger and ringdown. Based on this expectation we have identified an analytic fit model for the late-time frequency evolution that precisely matches the data for all cases but the (3,2) mode. These fits provide a quantitative understanding of the late-time phasing yielding, in particular, a robust result for the peak chirp rate for all modes. We will apply this information in Sec. VI.

### v.2 Late-time frequency and amplitude relationships.

The last step in our waveform analysis is to consider relationships between the frequency evolution and the waveform amplitude.

In the PN description of the quasicircular inspiral, the orbital frequency not only tells us the rotational rate, but can also serve as a label for describing the momentary state of the rotating object (in the inspiral case this means that we can reference the state of the system in terms of ). The PN generalization of the quadrupole formula, describing radiation from the rotating system, then leads to an expression for amplitude as a function of frequency. Our description of the gravitational radiation suggests an implicit source rotating with monotonically increasing frequency as it continues to “harden”, as the system evolves smoothly into merger and ringdown. In this section we seek to further unify this picture of the full coalescence process, considering an analogue of the PN description of amplitude as a function of frequency that can describe the radiation in the merger and ringdown.

In Sec. IV.2 we emphasized that the radiation power provides essentially the same information as the strain-rate amplitude (31). If the wave frequency is known, then the modal contribution to the total radiative angular momentum can similarly provide information about the gravitational wave amplitude. Eq. (34) gives an expression for angular momentum flux in terms of wave amplitude and phase. The relation simplifies to

(11) |

up to 5PN order, indicating that the radiation carries maximal angular momentum , as we generally expect for circularly polarized radiation. If we know how the rotational frequency evolves, we can derive the mode amplitudes from the mode-by-mode relationship of either energy or angular momentum with frequency, or .

To approach an understanding of the late-time relationships between amplitude and frequency, we compare how the system’s energy and angular momentum approach their final values, with how the system’s frequency approaches its final value. In Fig. 17 we examine the radiative loss of energy and angular momentum as the coalescing system approaches its final quiescent state, comparing these with the remaining difference of gravitational-wave frequency from its late-time limit , as determined in Table 2. As well as our standard final rotational frequency defined in (5), we also show the frequency: based on the strain, as defined in (2). We have also rescaled the energy and angular momentum by a constant, selected so that the value matches that of at the time of peak radiation power.

Fig. 17 indicates a general correspondence between how the angular momentum approaches its final state, and how the gravitational-wave frequency approaches its final state. Of the several cases of spherical harmonic modes and mass ratios that we have examined in this manner, we show two examples: the case for the equal-mass simulation, where the evolution of angular momentum and frequency correspond most closely (left panel), and the mode from the 4:1 mass ratio simulation, with the weakest correspondence (aside from the nonconforming case) (right panel). The correspondence is closest, holding to a better approximation over a longer period of time, in associating with . The association with energy is slightly weaker.

Though these plots suffer significantly from small modulations in frequency that we have not resolved numerically, the results suggest an approximate relationship between rotational frequency and angular momentum in particular. If the frequency evolution is otherwise known, then the late-time evolution of angular momentum for each mode could be approximately described by , where is a case-dependent constant that we will not attempt to specify generically. As we may refer to it as the dynamical moment of inertia of the implicit rotating source. This late-time expectation creates the possibility of extending our PN-based understanding of angular momentum flux into the late time waveforms, giving us the additional information we need for a full (approximate) description of .

This relation between frequency, angular momentum, and amplitude provides a connection between the peaks in modal radiation power, shown in Fig. 11, and the peaks in chirp rate given in Table 2. Assuming , a constant, yields

(12) | |||||

(13) |

Since the curve is steeply increasing, we would expect the peak in to be near the peak in but slightly delayed. Expanding (13) linearly about , when , we find that reaches its peak at a time

(14) |

where and are the frequency and its third time derivative at . These can be evaluated from (9) and the fit parameters , , , and in Table 2; we find that the values for all lie in the approximate range , in rough agreement with the direct fit for in Table 2.

Conceptually, the monotonic evolution of frequency might lead us to hypothesize a relationship between the structure of the implicit rotating source and its rotational frequency, such that the changes in the structure of this radiating rotator will be associated with finite changes in frequency. Knowing that the frequency growth must be limited by the quasinormal ringing frequency implies a peak in the chirp rate. This leads to the expectation that finite changes in energy and angular momentum are associated with finite changes in frequency, so that and approach constants once the evolution in frequency slows. The peak in radiation power might then be viewed as a consequence of the peak in chirp-rate. More quantitatively we find that seems to be roughly constant even before the evolution in frequency slows down (i.e. before ). This correspondence will be applied in the next section to provide a model for amplitude evolution through the peak, based on information about the frequency evolution.

If we could postulate the constancy of , we might also apply that assumption to “explain” some of what we have seen above. In Sec. V.1 we noted a rough agreement between the timescale in our frequency fitting curve (9) and the quasinormal ringdown amplitude decay rates for the corresponding quasinormal modes. Following the discussion above, this relationship could be derived, in the limit, from the constancy of either energy or angular momentum losses with respect to change in frequency. For constant , Eqs. (31), (5) and (12) imply that

(15) |

In the limit, the strain-rate amplitude decays at the rate predicted from black-hole perturbation theory, , where is the e-folding rate for the amplitude decay for the fundamental quasinormal mode. In this limit, our frequency evolution fit model reduces to

(16) |

Applying these limiting expressions for amplitude and frequency in (15) in the limit yields

(17) |

where the left-hand-side derives from the amplitude, and the right-hand-side from frequency. Ignoring the constant coefficients, this implies that .

If we adventurously assume the constancy of on approaching the ringdown, and expand the amplitude in powers of , the implied amplitude frequency relation might also provide more information about the amplitude evolution. Since , and consequently the right-hand-side of (17), contains only even powers of , the next term in the expansion for amplitude should be . This suggestion motivated our expansion (38) applied in fitting the late-time amplitudes in Sec. III.

## Vi Variations on the EOB model

An approach to modeling black-hole binary radiation known as the effective-one-body (EOB) model has been presented in the literature Buonanno and Damour (1999); Damour et al. (2000); Buonanno and Damour (2000, 2002); Damour (2001); Buonanno et al. (2006); Damour and Nagar (2007). The late-time waveforms in these models are based on a now-common description of the merger process as an epoch of radiation from spiraling particlelike trajectories, followed, in a sudden transition, by black-hole ringdown dynamics with waveforms described by a superposition of quasinormal frequencies. Our waveform analysis provides a complementary description of black-hole binary merger radiation that can be applied in an alternative late-time waveform model.

A recent promising approach along these lines is the so-called
*pseudo-4PN* (p4PN) EOB model Buonanno
et al. (2007b). This model extends
the 3PN-accurate EOB metric with a term of 4PN order with a tunable
multiplier . The phasing obtained from this expansion is
combined with leading-order PN strain amplitudes (the restricted
approximation) to obtain mode-by-mode waveforms valid for
inspiral. For the merger and ringdown, the phasing and amplitude are
derived simultaneously from a superposition of quasinormal modes. For
each angular mode, the fundamental ringdown modes and a few
overtones are summed in proportions as required for continuity with
the late end of the radiation from the inspiral phase. The value of
the p4PN multiplier is then chosen to match the pre- and
post-merger waveform portions, optimizing the agreement with
full-numeric waveforms.

In this section, we show that it is possible to develop variations on the p4PN EOB model that usefully encode some of the waveform phase and amplitude relationships we have described above. A key difference with the new variant is our prescription for the transition from inspiral to merger-ringdown radiation. In contrast with Buonanno et al. (2007b), for each mode we consider the entire wave train as that of a slowly varying instantaneously rigid rotator, consistent with the dominant “circular-polarization” waveform pattern encoded in the radiation. The phase evolution might be thought of as arising from the rotation rate of the corresponding source structure, which is continued through ringdown by continuous matching to a function of the form (9). The wave amplitude will be derived directly from expectations for the energy or angular momentum content of the radiation.

We present two specific models of this nature. Model 1 is based on exactly the same EOB-based prescription for inspiral-plunge trajectories in Buonanno et al. (2007b), while Model 2 shows the effect of a slight variation in the underlying EOB model. In both variants, as in Buonanno et al. (2007b), we derive the waveform phasing directly from the EOB trajectories (with for the strength of the p4PN term) up until some matching time, which we take simply as the time at which the (2,2) wave frequency is half the frequency of the fundamental (2,2) ringdown mode. After this point we use our fit model (9) to describe the subsequent phase evolution.

Recall that this model depends on several parameters: , , , , and . The results of our analysis in Sec. V.2 guide us in producing a fully specified model for these parameters. We take from the fundamental ringdown frequency of the radiation. For the time-constant for frequency decay we use the fundamental quasinormal mode amplitude decay time constant. While this is not clearly implied by our fits in the last section, it will lead to the correct amplitude fall-off, as specified below. For the strongest modes, our fits indicate . Lacking any better model we simply increase this by a factor when , roughly consistent with the higher value of found for the (2,1) mode. For these models we derive the quasinormal modes using the fit for the final black hole mass and spin described in Buonanno et al. (2007b). The remaining parameters and are chosen to provide continuity, up to the second time-derivative of phase, with the direct EOB-based phasing at the matching time.

Our prescription for the wave amplitudes differs from the restricted
amplitude description applied in Buonanno
et al. (2007b). As we showed
in Fig. 11, we can improve on the restricted
amplitude approximation by using the leading-order PN expressions for
waveform mode amplitudes only to fix the
*partitioning* of radiation power into angular modes. We set the
total power independently, from the full-order EOB model description
of the radiation power; the resulting waveforms are then energetically
consistent with the EOB description of the dynamics, and also show
better agreement with the numerical results (note that in this model
there is no radiation in the nonrotational modes). After
reaching the matching frequency, we continue the amplitude evolution
based on the assumption, suggested in Sec. V.2, that
the amplitude is roughly constant through the radiation peak.

In our new model, we set the late-time amplitude by asserting the approximate relationship (15), written this time in terms of the polarization frequency :

(18) |

setting the value of for amplitude continuity at the match-frequency. With this model for the amplitude, the peak in the gravitational-wave amplitude is a direct consequence of the peak in the time-derivative of gravitational frequency, fixed by in (9). The exponential decay in amplitude also follows directly from the exponential approach of the wave frequency to .

We compare the frequency and amplitude of the modeled waveform component with the corresponding numerical result for the 1:1 and 4:1 mass ratio cases in Fig. 18. A similar comparison is shown for some of the subdominant modes in the 4:1 case, in Fig. 19. The model we have so far described, based on the p4PN EOB trajectories is labeled Model 1 in the figures. The matching frequency is indicated as a vertical line in each plot. The frequency curves indicate very good phasing agreement for all cases except the mode of the 4:1 mass ratio case. In that case the sharp rise in frequency occurs a few too late, several times worse than the agreement shown for the other modes. The difference is a consequence of the slightly higher frequency early on, as the numerical frequency begins to grow already before the matching point as compared with the rotational frequency consistent with the numerical modes.

The generally good amplitude agreement shown before merger represents an improvement over the simple restricted amplitude model employed in Buonanno et al. (2007b). For the modes, Model 1 overestimates the amplitudes just before and consequently after matching. This indicates that the PN-based flux expression applied to the p4PN EOB model, based on Buonanno et al. (2003) and Damour et al. (1998), overestimates the flux at high frequencies. This is perhaps not surprising, since the flux formula is formulated with a pole at finite frequency. Physically we expect the flux to decrease at late times, when the frequency approaches the quasinormal ringdown frequency.

To correct for this, we show also results for a second variation, Model 2, in which we have introduced a zero in the flux function at ringdown frequency. Following the notation of Buonanno et al. (2003), we modify their Eq. (45) for the flux to read

(19) |

and then respecify the coefficients in their Eq. (50) to again provide consistency to 3.5PN order with the Taylor series expansion for the flux. Note that the flux now depends on and via the quasinormal mode frequency and . The modified flux function anticipates that the radiation will cut off at the ringdown frequency. We find that the new EOB model provides a very good approximation to the original p4PN EOB phasing (with ) if we choose for the new version with modified flux. The figures show that, for the waveform modes, the amplitudes based on Model 2 with the modified flux show better agreement with the numerical results leading up to the matching point, and also at the peak, while the frequency evolution is nearly identical. For the (2,1) case, Model 2 suffers the same problems as Model 1.

We have applied the observations made in previous sections to successfully predict late-time waveforms. This provides an alternative description of the transition to ringdown, distinct from the widely-applied approach of summing quasinormal modes Buonanno et al. (2007a); Damour et al. (2007); Berti et al. (2007). In Ref. Buonanno et al. (2007b), we applied it to nonspinning mergers.

This raises the question: how is it that each of these quite distinct approaches to approximately describing the late-time radiation can be simultaneously effective? We can explore this question by considering the polarization frequency evolution in a waveform constructed from a sum of quasinormal modes (suppressing and labels)

(20) |

with and where and correspond to the quasinormal overtone mode. Next we restrict to the first two terms and to linear order in , which vanishes at late times. This yields

(21) | |||||

(22) |

where and . Taking the derivative of the expression in square brackets gives the polarization frequency of

(23) |

where . Note that the expression in brackets is periodic with period , where is the difference between the fundamental quasinormal ringing frequency and its first overtone [see Berti et al. (2006) for a table of these overtones].

For the waveform modes we consider here is quite small, generally 1 or 2%, which means that the period of the expression in brackets is . We rewrite , only keeping terms linear in , as

(24) |

Considering the second term, we note that is generally just under for the cases we’ve studied so that . The other factor, , is not predictable without knowledge of the initial conditions, but most values of leave the second term somewhat less than one.

If these conditions hold near the onset of ringdown, we might neglect the time dependence in the second term. The late-time frequency evolution is approximated by a exponential decay with a time constant of . Working under these assumptions, and comparing with the late-time frequency evolution model in this paper, yields the association , where is the fitting parameter in our late time frequency evolution model. Again, looking at the quasinormal mode values, we comment that so that , the approximate relationship noted in Sec.V.1, which was also approximately derived from the assumption that is constant at late times. Consistently, applying the relation in Eq. (22) also give an expression for late-time amplitude consistent with Eq. (38).

## Vii Discussion

With hope of reaching out to a wide range of researchers interested in gravitational radiation from black-hole binary mergers, we have provided a descriptive walk-through of many of the general features of late-time waveforms from generic mergers of nonspinning binary black hole systems, based on a series of numerical simulations covering systems from equal-mass up to mass ratio 6:1. In this basic waveform description we have examined waveform phase and amplitude, comparing results among different mass ratios, as well as among the different spin-weighted spherical harmonic component modes.

In our presentation, we have attempted to describe the radiation in the simplest physical terms, pointing out traits in the waveforms that are similar through the inspiral, merger and ringdown stages. Throughout the coalescence, we find simple waveforms in each mode, each exhibiting strong circular polarization and monotonically increasing polarization frequency.

In our amplitude comparisons, we find that the leading-order PN-prediction for energy-partitioning provides a good estimate of the amplitude until late in the merger for modes. In astrophysical units, our fit for the peak-power in the dominant mode is

(25) |

Scaling the amplitudes by the peak values yields a very similar shape through the peak for the mode amplitudes for all mass ratios we have studied. For a particular mass-ratio case, the peak-widths remain similar, though there is some variation in peak-time.

For each mass ratio, the phase (and frequency) of the different components are strongly related. While this should be expected for the early inspiral, where the waveform phase is directly connected to the orbital phase, we show that, for modes, the same relationship holds through the merger and into the inspiral. We compare the phasing among simulations with different mass ratios in two ways, with time scaled by chirp-mass, as is appropriate in the early-time limit, and then with time scaled by total system mass . With the latter scaling, the dominant waveforms are similar in phase (and -scaled amplitude) for the last .

In the near-peak waveform comparisons the mode does not exhibit the same simple behavior as the other modes. This mode is easily subject to coupling with the much stronger mode. It has been seen in Buonanno et al. (2007a); Schnittman et al. (2008) that the (3,2) mode demonstrates significant mode mixing with the dominant (2,2) mode – the QNM ringdown part of the (3,2) waveform contains the fundamental frequencies of both the (3,2) and the (2,2) modes. We speculate that this mixing might be partly due to the use of coordinate extraction spheres that are systematically warped from the areal-radius spheres appropriate for correct radiation extraction. A more refined choice of extraction spheres and perhaps a better tuned decomposition basis (e.g. spheroidal harmonics) would make it possible to represent the waveform content in a manner which is like that seen with the other modes. Thereby we expect that a similar simple physical representation of the waveform content can be extended to all modes.

We suggest a simple conceptual interpretation that applies through the
full coalescence. We think of the radiation as being generated by an
*implicit rotating source*, with each mode generated
separately by the moment of some implicit source (which we
understand here only in the context of the radiation). The nearly
fixed relationship among the phase moments is interpreted
to indicate that the implicit source maintains some structural
integrity throughout the coalescence, without shearing among the
various modal components. For the modes, this rigidity is
maintained through the merger and into the ringdown, a relationship
made possible by the approximate equality for each of the
quasinormal modes , where
is the orbital frequency of unstable circular prograde graviton (or
photon) orbits.

The following physical picture may underlie these relationships. For well-separated black holes, the fields that embody the implicit source object evidenced in the radiation may be tied directly to the pointlike centers of the orbiting black holes. The source rotation frequency is the orbital frequency of the timelike trajectories traced out by the black holes. As the binaries spiral together, the pair can continuously be viewed as a shrinking, distributed dumbbell-like rotator. Eventually, most of this dumbbell shrinks inside the light-ring, which roughly coincides with the potential barrier in the wave mechanics of gravitational perturbation theory. From inside, little radiation can escape to a distant observer, and the timelike motion of the black hole centers disconnects from the radiation. At late times, the effective radiation source becomes a gravitational disturbance orbiting the forming black hole at the light ring. This is a seamless transition, with nearly consistent rotational phasing among all modes throughout the process. For the modes the associated quasinormal-ringing dynamics are somewhat distinct, and the phasing and amplitude relationships begin to peel away from the main trend through the merger process.

For the late-time portions of the waveforms, including the approach to the peak and the ringdown, we have introduced a quantitative fitting model based on a monotonically increasing polarization frequency for each mode, which decays exponentially toward the expected fundamental quasinormal ringdown frequency at late times. These fits provide an excellent match for the frequency evolution beginning before the peak, and allow precise estimates of , as well as the peak rate of change in frequency . Scaling the latter quantity by the final black hole mass and , we find for all waveforms we have looked at, including the modes of each mass-ratio, and modes up to for the 4:1 mass-ratio case.

Conceptually, the monotonicity of the frequency evolution suggests that, as is the case for inspiralling systems, the frequency can be taken to label the state of the adiabatically changing implicit rotating source that we interpret as the source of the radiation. If we suppose that changes in the source structure are tied to loss of angular momentum, then we would expect that finite changes in frequency would be associated with finite angular momentum loss, so that has a finite, nonzero value even at late times. Since approaches a nonzero constant at late times, we would likewise expect to approach a constant value. In Sec. V.2, we show that the late-time evolution of and are approximately consistent, mode-by-mode, with constant beginning about before merger.

Such a relation between frequency and angular momentum also implies a connection between frequency and amplitude. The moment of peak amplitude is expected to be near the peak in . At very late times, constant implies a connection between the rates at which the frequency and amplitude approach their quasinormal late time state, namely that our fitting parameter , where is the damping time of the quasinormal amplitude decay [see Eq. 9].

The simple relationships between the waveform modes and simple dependence on mass-ratio make it possible to specify much of the late-time waveform information developed in our numerical simulations in terms of just a few quantities. This information can then be combined with information from the PN approximation about the inspiral trajectories to provide analytical models for full-coalescence waveforms. In Sec. VI we have applied the PN-consistent EOB-based p4PN trajectory model presented in Buonanno et al. (2007b) together with assumptions asserting several of the approximate waveform features observed in Sec. IV. As in Buonanno et al. (2007b), early waveform phasing is derived from the p4PN EOB trajectories, with waveform amplitudes based on our PN-based power partitioning, together with the PN flux model. As the waveform approaches the anticipated peak, we match to a waveform phasing model based on our fit model, with parameters specified according to the approximate relationships identified in Sec. IV, and with amplitudes derived from =constant. The fits show excellent phase agreement with the numerical simulation results for the most significant modes. A variation in the PN flux model that enforces that the flux for each mode vanish as gives better late-time amplitude agreement.

Lastly, we observe that our description of the late-time phasing and amplitudes provides a picture complementary to another approach applied in several previous studies Buonanno et al. (2007b); Damour et al. (2007); Berti et al. (2007); Buonanno et al. (2007a), which successfully treat the late-time waveforms as a sum of quasinormal fundamental and overtone modes for each waveform component. This is motivated by the expectation that waveforms from generic initially compact distortions of the forming black hole will quickly reduce to a sum of these quasinormal harmonics Teukolsky and Press (1974). In Buonanno et al. (2007b), we have shown that in comparisons with some of the runs presented here, for mass ratios up to 2:1, this assumption can lead to a predictive waveform model with similar accuracy to our alternative model presented in Sec. VI. As a link between the two approaches, we have shown that at late times the combination of the fundamental and first QNM overtones for a particular mode may, under reasonable circumstances, mimic the amplitude decay properties of our model.

This work suggests several directions for further study. In the immediate future, we plan to assess the fidelity of the available nonspinning waveforms and models, and the impact of mass ratio on the overall detectability of the merger signal. We also plan to apply our implicit-rotating-source description as a baseline in analyzing future higher-precision numerical simulations. This might provide insight into understanding finer features of the merger physics, some of which could violate our simplified description. Further understanding the anomalous mode waveforms will be a first step in this direction. We must also investigate whether our description of the merger radiation applies also to spinning black holes. It is plausible that even precessing systems might be analyzed in this way using a spherical harmonic basis that tracks the orbital axis Gualtieri et al. (2008). This may make it possible to extend the analytic EOB-based waveform model presented here to include spin effects. Including spins in such analytic models will be necessary for observational data analysis applications.

###### Acknowledgements.

We thank Emmanuele Berti and Alessandra Buonanno for interesting discussions, and Luciano Rezzolla for useful comments. This work was supported in part by NASA grant 05-BEFS-05-0044 and 06-BEFS06-19. The simulations were carried out using Project Columbia at the NASA Advanced Supercomputing Division (Ames Research Center) and at the NASA Center for Computational Sciences (Goddard Space Flight Center). B.J.K. was supported by the NASA Postdoctoral Program at the Oak Ridge Associated Universities.## Appendix A Radiation Extraction

We extracted radiation using the outgoing Weyl scalar , defined as in Baker et al. (2002b), calculated with a symmetric tetrad. This is related to the complex gravitational-wave-strain via

(26) |

, the strain , and its time-derivative (which we
call the *strain rate*)
are all functions of time , extraction radius ,
and polar angles , . As is customary in
numerical relativity, we decompose the radiation into
spin-(-2)-weighted spherical harmonic components:

(27) | |||||

(28) | |||||

(29) |

With this, two time integrations of our measured quantity, , yields the more familiar gravitational-wave strain .

To extract the radiation from a simulation, we define a series of coordinate spheres of different radii ; here we use extraction spheres having radii between and . We extracted the radiation in modes by integrating against different over these coordinate spheres, using fourth-order interpolation onto each sphere followed by Newton-Cotes angular integration.

The gravitational waves produced by the binary carry both energy and angular momentum. Overall, the rate of energy emission is given by an angular integral of the squared strain rate over a coordinate sphere [see Eq. (5.1) of Baker et al. (2002a)]:

(30) |

Then using Eq. (28), we can express the total energy flux (30) as a sum over modes:

(31) | |||||

(32) |

where we have used the strain-rate decomposition (5) and taken the limit to go from (31) to (32).

Similarly, the rate of radiation of the -component of angular momentum can be expressed as a sum over modes Lousto and Zlochower (2007):

(33) |

## Appendix B Convergence

We carried out three runs of the 4:1 mass ratio model at different resolutions to study the convergence properties of our simulations. For these cases, the mesh spacing of the finest grids (the ones including the smaller puncture) was taken to be (low resolution), (medium resolution), and (high resolution). To facilitate comparisons among these cases, the overall grid structure of the runs was kept the same. In this Appendix, we discuss the convergence properties of the constraints and gravitational waveforms.

In comparing our medium and high resolutions for the 4:1 mass ratio case, the Hamiltonian constraint was found to be manifestly fourth-order-convergent in the dynamical strong-field, where the black holes move, evidently dominated by the expected fourth-order error from refinement interfaces. The convergence falls off to an apparent rate closer to first order in the coarsest regions. This seems to result from stronger dissipation of high frequency noise. Though less noise is generated in the higher resolution simulations, a greater portion of it survives propagation into the distant coarser regions. As shown below, this does not appear to affect the waveforms, which are well-resolved in the wave zone. The momentum constraint appeared to be at least second-order-convergent in the dynamical strong-field region, but also fell off to an apparent rate closer to first order in the coarsest regions – from the wave-extraction region outwards.

In our waveform mode analysis, we have typically time- and phase-shifted the data so that the amplitude peak fell at time Baker et al. (2006a). In Fig. 3, the real part of the (2,2) and (4,4) strain-rate harmonics is shown for all three resolutions. The agreement between the resolutions is then seen to be excellent. In Figs. 4 and 5, we showed the errors to be expected from the maquillaged strain-rate waveform data.

However, the presence of considerable eccentricity in the binary makes
it difficult to compare such time-shifted waveforms between
resolutions and establish an unambiguous order of convergence. In
In Figs. 20 and 21 we show the
amplitude and phase errors, respectively, for our three resolutions
*without* time-shifting – that is, we plot the data
starting from the initial time in each case, scaling the medium - high
differences for fourth and fifth-order convergence. It is clear that
we observe convergence between fourth and fifth order throughout the
evolution until around before merger, when the difference in
merger times among the runs becomes important. Additionally, it is
possible that at this higher-frequency stage of merger, the
lowest-resolution data is no longer in the convergence regime. We do
note, however, that the rate of growth of amplitude difference between
the medium and high-resolution runs is comparable here to that
observed closer to the merger for time-shifted data [see
Fig. 4]. From this we deduce that the medium and
high-resolution runs are still in the convergence regime, with errors
consistent with fourth-order convergence.

## Appendix C End states

In this Appendix we discuss the state of the final black hole formed in our simulations. This can be measured by several independent means, which we compare, finding agreement to within for and for . Table 3, we present “coarse” results for the simulations. The total radiated energy is obtained by integrating Eq. (30), and the total radiated angular momentum is obtained by summing over mode contributions (33). The final mass and angular momentum of the post-merger black hole are calculated using

(35) | |||||

(36) |

where, by the symmetries of the current simulations, we only deal with the component of angular momentum. For the 1:1 simulations, we only had the leading-order radiation modes available, so our estimate is significantly truncated, by as much as 11% (a conservative error estimate based on the effect of similarly truncating mode contributions past for the 4:1 case); we have marked this and derived values. We also quote the measured value of , the time at which reaches its peak. All waveform and derived plots in this paper have been time-shifted by subtracting this time, as an approximate marker of the time of merger, unless otherwise indicated.

Table 3 also contains data about the common apparent horizon (CAH) of the merged binary: the time at which the CAH was first detected, the CAH’s irreducible mass , and its full (horizon) mass, obtained from the Christodoulou Christodoulou (1970) formula:

(37) |

where we use for the final hole’s angular momentum.

Mass |
---|