High-accuracy numerical simulation of black-hole binaries: Computation of the gravitational-wave energy flux and comparisons with post-Newtonian approximants

Expressions for the gravitational-wave (GW) energy flux and center-of-mass energy of a compact binary are integral building blocks of post-Newtonian (PN) waveforms. In this paper, we compute the GW energy flux and GW frequency derivative from a highly accurate numerical simulation of an equal-mass, nonspinning black-hole binary. We also estimate the (derivative of the) center-of-mass energy from the simulation by assuming energy balance. We compare these quantities with the predictions of various PN approximants [adiabatic Taylor and Padé models; nonadiabatic effective-one-body (EOB) models]. We find that Padé summation of the energy flux does not accelerate the convergence of the flux series; nevertheless, the Padé flux is markedly closer to the numerical result for the whole range of the simulation (about 30 GW cycles). Taylor and Padé models overestimate the increase in flux and frequency derivative close to merger, whereas EOB models reproduce more faithfully the shape of and are closer to the numerical flux, frequency derivative, and derivative of energy. We also compare the GW phase of the numerical simulation with Padé and EOB models. Matching numerical and untuned 3.5 PN order waveforms, we find that the phase difference accumulated until $M\omega =0.1$ is $-0.12$ radians for Padé approximants, and 0.50 (0.45) radians for an EOB approximant with Keplerian (non-Keplerian) flux. We fit free parameters within the EOB models to minimize the phase difference, and confirm the presence of degeneracies among these parameters. By tuning the pseudo 4PN order coefficients in the radial potential or in the flux, or, if present, the location of the pole in the flux, we find that the accumulated phase difference at $M\omega =0.1$ can be reduced—if desired—to much less than the estimated numerical phase error (0.02 radians).

Figure 1

Some aspects of the numerical simulation. From top panel to bottom: the leading mode ${\dot{h}}_{22}$; the two next largest modes, ${\dot{h}}_{44}$ and ${\dot{h}}_{32}$ (smallest); the frequency of ${\dot{h}}_{22}$ [see Eq. (5)].

Figure 2

Lower panel: relative difference between flux $F(\varpi )$ computed with 99 different intervals $[t_1,t_2]$ and the average of these. Upper panel: relative change in the flux $F(\varpi )$ under various changes to the numerical simulation. The gray area in the upper panel indicates the uncertainty due to the choice of integration constants, which is always dominated by numerical error. The dashed line in the upper panel is our final error estimate, which we plot in later figures.

Figure 3

Contributions of various $(l,m)$ modes to the total numerical gravitational-wave flux. Upper panel: plotted as a function of time. Lower panel: plotted as a function of frequency $M\varpi$. The lower panel also contains the error estimate derived in Fig. 2.

Figure 4

Lower panel: difference between frequency derivative $\dot{\varpi }$ computed with 99 different intervals $[t_1,t_2]$ and the average of these. Upper panel: change in the frequency derivative $\dot{\varpi }$ under various changes to the numerical simulation. The gray area in the upper panel indicates the uncertainty due to choice of integration constants, which dominates the overall uncertainty for low frequencies. The dashed line in the upper panel is our final error estimate, which we plot in later figures.

Figure 5

Ratio of GW frequencies $\omega$ and $\varpi$ to orbital frequency, as a function of (twice) the orbital frequency, for different PN models. The GW frequencies $\omega$ and $\varpi$ are defined in Eqs. (4, 5). Solid lines correspond to 3.5PN; dashed and dotted lines correspond to 3PN and 2.5PN, respectively.

Figure 6

Effect of choice of frequency. Shown are the PN fluxes for two representative PN approximants, plotted (correctly) as function of $\varpi$ and (incorrectly) as function of $2\Omega$. Plotting as a function of $2\Omega$ changes the PN fluxes significantly relative to the numerical flux $F_{\mathrm{NR}}$.

Figure 7

Comparison between the numerical energy flux and several PN approximants at 3.5PN order versus GW frequency $\varpi$ extracted from ${\dot{h}}_{22}$ in the equal-mass case.

Figure 8

Comparison between the numerical energy flux and several PN approximants versus GW frequency $\varpi$ extracted from ${\dot{h}}_{22}$ in the equal-mass case. We show the relative difference between numerical flux and PN flux, as well as the estimated error of the numerical flux (blue bars; see Fig. 2). Solid lines represent 3.5PN models and NR; dashed and dotted lines correspond to 3PN and 2.5PN models, respectively. For notation see Table  and caption therein.

Figure 9

Comparison of normalized energy flux $F/F_{\mathrm{Newt}}$ [see Eq. (73)] for the equal-mass case. Solid lines represent 3.5PN models and NR; dashed and dotted lines correspond to 3PN and 2.5PN models, respectively. For notation see Table  and caption therein.

Figure 10

Cauchy convergence test of $F/F_{\mathrm{Newt}}$ for T and P approximants. We plot $\Delta F_{n+m}\equiv F_{n+m+1}-F_{n+m}$, and $\Delta F_n^m\equiv F_{n+1}^m-F_n^m$ for different values of $v_{\Omega }$. The T and P approximants are given by Eqs. (19, 39), respectively. Note that the P approximant has an extraneous pole at 1PN order at $v_{\Omega }=0.326$. We use $v_{\mathrm{lso}}=v_{\mathrm{lso}}^{2\mathrm{PN}}$, and $v_{\mathrm{pole}}=v_{\mathrm{pole}}^{2\mathrm{PN}}$.

Figure 11

Fitting several PN approximants to the numerical flux. The $x$ axis denote the orbital frequency $\Omega$. Because the numerical flux is computed as function of the GW frequency, we use for the numerical flux $\Omega \equiv \varpi /2$. The blue bars indicate estimated errors on the numerical flux; see Fig. 2. For notation see Table  and caption therein.

Figure 12

GW frequency derivative $\dot{\varpi }$ for the numerical relativity simulation and various PN approximants at 3.5PN order. For notation see Table  and caption therein.

Figure 13

Comparison of $\dot{\varpi }$ for the numerical results and various PN approximants. Dotted, dashed, and solid lines correspond to 1.5PN, 2.5PN, and 3.5PN models, respectively. For notation see Table  and caption therein.

Figure 14

Comparison of PN $\dot{\varpi }$ with a heavily smoothed version of the numerical $\dot{\varpi }$. Solid lines represent 3.5PN models and NR; dashed and dotted lines correspond to 3PN and 2.5PN models, respectively. For notation see Table  and caption therein.

Figure 15

We compare $dE/d\varpi$ versus GW frequency $\varpi$ for numerical relativity [see Eq. (76)] and PN approximants. Solid lines represent 3.5PN models and NR; dashed and dotted lines correspond to 3PN and 2.5PN models, respectively. For notation see Table  and caption therein.

Figure 16

Phase differences between the numerical waveform, and untuned, original EOB, untuned Padé, and Taylor waveforms, at two selected times close to merger. The E approximants are $F_n^m/H_p$, while the P approximants are $F_n^m/E_p^q$ (see Table  and caption therein). Waveforms are matched with the procedure described in Sec. 6a and phase differences are computed at the time when the numerical simulation reaches $M\omega =0.063$ (left panel) and $M\omega =0.1$ (right panel). Differences are plotted versus PN order. Note that at 1PN order the Padé flux has an extraneous pole at $v=0.326$ causing a very large phase difference. The thick black line indicates the uncertainty of the comparison as discussed in Sec. 6, $|{\Phi }_{\mathrm{PN}}-{\Phi }_{\mathrm{NR}}|\le 0.02$ radians.

Figure 17

Phase differences between untuned and tuned P approximants and NR. The untuned P approximant is $F_4^3/E_2^4$ ($v_{\mathrm{lso}}=v_{\mathrm{lso}}^{2\mathrm{PN}}$, $v_{\mathrm{pole}}=v_{\mathrm{pole}}^{2\mathrm{PN}}$). The tuned P approximants are $F_4^4/E_2^4$ and tunable $v_{\mathrm{pole}}$ ($v_{\mathrm{lso}}=v_{\mathrm{lso}}^{2\mathrm{PN}}$) and $p{F}_4^4/E_2^4$ ($v_{\mathrm{lso}}=v_{\mathrm{lso}}^{2\mathrm{PN}}$, $v_{\mathrm{pole}}=v_{\mathrm{pole}}^{2\mathrm{PN}}$) with tunable $A_8$ and $B_8$. In all cases, waveforms are matched over $t-r^{*}\in [1100,1900]M$.

Figure 18

The upper panel shows phase differences versus time (lower $x$ axis) and versus GW frequency (upper $x$ axis) for several tuned and untuned E approximants. For the tuned models, the optimal $a_5$ and $v_{\mathrm{pole}}$ values displayed in Table . In the lower panel, we show phase differences between numerical and E approximants computed at $t_{0.063}$, $t_{0.1}$, and the end of the numerical simulation $t_{0.16}$, as functions of $a_5$. For the same color and style, the curve with the steepest slope corresponds to $t_{0.16}$ and the curve with the smallest slope corresponds to $t_{0.063}$. (For notation see Table  and caption therein.)

Figure 19

Normalized energy flux $F/F_{\mathrm{Newt}}$ versus GW frequency $2\Omega$ in the test-mass limit. For notation see Table  and caption therein. For comparison, both panels also include the result of the numerical calculation of Poisson [45], labeled with NR.

Figure 20

Convergence of the PN approximants in the test-mass limit. Plotted are differences of $F/F_{\mathrm{Newt}}$ from the numerical result. The P approximants do not converge faster than the Taylor series.

Figure 21

Cauchy convergence test of $F/F_{\mathrm{Newt}}$ in the test-mass limit for the T and P approximants. We plot $\Delta F_{n+m}\equiv F_{n+m+1}-F_{n+m}$, and $\Delta F_n^m\equiv F_{n+1}^m-F_n^m$ at three different frequencies. At high frequencies, the 4.5PN and 5PN Padé approximants are contaminated by the extraneous pole of the 5PN Padé series; for low frequencies ($v_{\Omega }=0.1$), the pole is apparently irrelevant.