An improved numerical fit to the peak harmonic gravitational wave frequency emitted by an eccentric binary

I present a numerical fit to the peak harmonic gravitational wave frequency emitted by an eccentric binary system in the post-Newtonian approximation. This fit significantly improves upon a previous commonly-used fit in population synthesis studies, in particular for eccentricities $\lesssim 0.8$.

Here, M ≡ m 1 + m 2 is the total binary mass, and G is the gravitational constant. As shown in the seminal work of Peters & Mathews (1963) that assumed the (lowest-order) post-Newtonian approximation, an eccentric binary (eccentricity e) emits GWs at additional harmonics. Specifically, the power of the n th harmonic with frequency f GW,n = nf orb (with integer n ≥ 1) is given by P n = 32 5 G 4 c 5 m 2 1 m 2 2 M a 5 g(n, e), where the function g(n, e) quantifies the factor to which more power is emitted in the n th harmonic compared to a circular orbit. The latter function is given by g(n, e) = n 4 32 J n−2 (ne) − 2eJ n−1 (ne) + 2 n J n (ne) where J i (x) is the i th Bessel function of the first kind. Fig. 1(a) plots g(n, e) as a function of n (interpreted as a real number) for several values of e. More eccentric binaries typically emit more power at higher harmonics. It is of interest to consider the peak harmonic, n peak , i.e., n peak (e) is the value of n for which g(n, e) is at a maximum. The corresponding peak GW frequency immediately follows from the relation f GW,peak (e) = n peak (e)f orb .
Equation (5) is plotted as dashed vertical lines in Fig. 1(a). Although a reasonable fit for high eccentricities, it does not capture the true maximum of g(n, e) very accurately for lower eccentricities. Nevertheless, the fit of W03 is commonly used to estimate the peak GW frequency of eccentric binaries, in particular in the context of population synthesis studies (e.g. Here, I present a new numerical fit of n peak (e) which significantly improves upon that of W03 for low eccentricities (e 0.8), whereas also accurate for high eccentricities: where c 1 = −1.01678, c 2 = 5.57372, c 3 = −4.9271, and c 4 = 1.68506. The fit retains the factor (1−e 2 ) −3/2 from W03 that dominates at high eccentricities, but the behavior at smaller eccentricities is modified. Equation (6) correctly states that n peak (0) = 2 for circular orbits.
The new fits are indicated in Fig. 1(a) with solid vertical lines, and show significantly better match with the true n corresponding to a maximum in g(n, e) (when n is interpreted as a real number). Fig. 1(b) plots the peak harmonic n peak as a function of e. Grey open circles show the 'exact' calculation of n peak by numerically calculating the maximum of g(n, e) for given eccentricity. The latter is carried out in practice for real n ≥ 1, although it should be understood that n is actually an integer. The exact (real) calculation yields that n peak decreases below 2 as e → 0. However, it is clear that, in the circular limit, g(n, 0) = 1 for n = 2, arXiv:2111.08033v2 [gr-qc]  and g(n, 0) = 0 for all other integer n. Therefore, when determining n peak , one should take the integer value and limit to n peak ≥ 2 in order to retain the correct behavior in the fitting function as e → 0. By enforcing that n peak (0) = 2 and given the limitations of the assumed functional form, the fit of W03 significantly overpredicts n peak for e 0.8. The new fit better captures the low-eccentricity regime, while also still satisfying n peak (0) = 2 and giving a good description at high eccentricities. This is particularly clear in Fig. 1(c), in which rounded integer numbers are plotted.
Lastly, Fig. 1(d) shows the fractional residuals in the integer n peak computed as the difference between the 'exact' integer-rounded calculation with n peak ≥ 2 and the integer-rounded fits, i.e., For e 0.8, the fit of W03 systematically overpredicts n peak (∆n peak < 0). This is especially the case for the range 0.17 e 0.30, where Equation (5) predicts n peak = 3, whereas it should be n peak = 2. For larger e, the discrepancies become less severe, although at e = 0.8 the peak harmonic is still overpredicted by ∼ 10%. For even higher e, e 0.97, Equation (5) starts to slightly underpredict n peak (∆n peak > 0).
In contrast, the new fit Equation (6) has typically zero fractional residuals, with only a few spikes occurring due to rounding effects at transitions where n peak advances by unity. At e = 0.999 (the highest eccentricity considered when determining Equation 6), the new fit has zero fractional residuals, i.e., n peak (0.999) = 51, 755 according to Equation (6) is consistent with the exact calculation. In contrast, Equation (5) has a fractional error of 1% with n (W03) peak (0.999) = 51, 216. When applying the fit presented here, it should be remembered that it relies on the results of Peters & Mathews (1963) based on the lowest-order post-Newtonian terms that describe dissipation due to GW emission (2.5PN). Higher-order post-Newtonian corrections (e.g., Tucker & Will 2021) are not included.