Consistent eccentricities for gravitational wave astronomy: Resolving discrepancies between astrophysical simulations and waveform models

Detecting imprints of orbital eccentricity in GW signals promises to shed light on the formation mechanisms of binary black holes. To constrain formation mechanisms, distributions of eccentricity derived from numerical simulations of astrophysical formation channels are compared to the estimates of eccentricity inferred from GW signals. We report that the definition of eccentricity typically used in astrophysical simulations is inconsistent with the one used while modeling GW signals, with the differences mainly arising due to the choice of reference frequency used in both cases. We also posit a prescription for calculating eccentricity from astrophysical simulations by evolving ordinary differential equations obtained from post-Newtonian theory, and using the dominant ($\ell = m =2$) mode's frequency as the reference frequency; this ensures consistency in the definitions. On comparing the existing eccentricities of binaries present in the Cluster Monte Carlo catalog of globular cluster simulations with the eccentricities calculated using the prescription presented here, we find a significant discrepancy at $e \gtrsim 0.2$; this discrepancy becomes worse with increasing eccentricity. We note the implications this discrepancy has for existing studies, and recommend that care be taken when comparing data-driven constraints on eccentricity to expectations from astrophysical formation channels.


INTRODUCTION
Definitive measurement of orbital eccentricity in a compact binary coalescence (CBC) signal is one of the most exciting prospects of gravitational wave (GW) astrophysics.Despite O(100) binary black hole (BBH) signals being detected by the LIGO-Virgo-KAGRA (LVK) network (LIGO Scientific Collaboration et al. 2015;Acernese et al. 2015;Aso et al. 2013;Akutsu et al. 2021;Abbott et al. 2018) to date (Abbott et al. 2023a), the distinct pathways that led to the formation and merger of these BBH signals remains a mystery.Measurements of masses, spins, and redshift evolution have helped to identify broad constraints on the underlying physical processes that lead to BBH formation (Abbott et al. 2023b), although measurement uncertainties paired with the inherent uncertainties embedded within formation aditya@utoronto.cachannel modeling make it difficult to use these measurables to robustly pin down the formation pathway for a single BBH or the relative fraction of BBHs that results from one formation channel compared to another (see e.g.Zevin et al. 2021b;Cheng et al. 2023).
The former of these possibilities is something that can only be attained through strong dynamical encounters within dense stellar environments, such as globular clusters, nuclear clusters, and the disks of active galactic nuclei.The synthesis of eccentric BBHs mergers in globular clusters is a robust and unavoidable result of simple physical processes such as two-body relaxation and (chaotic) gravitational encounters between three or more black holes within the cluster core (see e.g.Rodriguez et al. 2018a).Since the relative fraction of systems that proceed through this eccentric subchannel is relatively robust to uncertainties in cluster environments, the detection and measurement of even a small number of eccentric BBH mergers can lead to stringent constraints on the relative fraction of merging BBH systems that result from dynamical environments entirely (Zevin et al. 2021a).
Although eccentricity is a powerful constraint on formation pathways in theory, in practice the measurement and interpretation of an eccentric GW signal is much more muddled.Current matched-filter CBC searches applied to GW data (Abbott et al. 2023b) only use aligned-spin and quasi-circular templates (Sachdev et al. 2019;Davies et al. 2020;Aubin et al. 2021), and thus a selection effect exists against detecting eccentric GW signals.The selection function for eccentric systems due to this method for searching the data has yet to be adequately characterized (see discussion in Zevin et al. 2021a).Waveform-independent searches for eccentric mergers have also been carried out (Abbott et al. 2019;Abac et al. 2023) with no significant triggers; this is unsurprising, since these searches are typically only sensitive to short duration (i.e.high mass) systems.A number of eccentric waveform approximants have been built (Tanay et al. 2016;Huerta et al. 2018;Nagar et al. 2021;Ramos-Buades et al. 2022;Liu et al. 2022), although much more work needs to be done to make them accurate over a large range of eccentricities.In particular, simultaneously modeling eccentricity and spin precession has proved to be challenging (see Liu et al. 2023 for a recent attempt1 ), and these properties have degenerate effects on the GW waveform making it difficult to distinguish between precessing and eccentric hypotheses especially for massive signals (Romero-Shaw et al. 2023b;Xu & Hamilton 2023;Romero-Shaw et al. 2020).Lastly, the definition of eccentricity differs between GW waveform models.Eccentricity evolves during the inspiral of a binary, and hence can only be defined with reference to an epoch in the orbit of that binary.However, the definition of this reference epoch and of eccentricity itself is inconsistent across different state-of-the-art eccentric waveform approximants, in part because eccentricity is a gauge-dependent quantity within general relativity.There have been attempts towards bringing all these definitions on an equal footing (Knee et al. 2022;Ramos-Buades et al. 2022;Bonino et al. 2023;Shaikh et al. 2023), allowing consistent comparison of eccentricities derived from diverse waveform approximants.
A robust and consistent definition for eccentricity that is shared across waveform models and astrophysical models is imperative to yield any astrophysical conclusions in the event of an eccentric BBH detection.In this paper, we show that the definition of eccentricity currently employed in astrophysical models is different from the definition used in waveform models.This discrepancy between these definitions is relevant for eccentricity e ≳ 0.2 with a factor of ∼ 2 discrepancy at high e.We also paint a clearer picture as to where these discrepancies in the definition of eccentricity lie, and present a framework for determining a robust definition of eccentricity that can be utilized for self-consistent comparisons between predictions from the astrophysical modeling community and the GW waveform community.
Our paper is structured as follows.In Section 2, we review current methods that are used for defining eccentricity, both in astrophysical simulations and waveform modeling.Section 3 presents a standardized definition of eccentricity and a scheme for extracting this definition of eccentricity from astrophysical simulations for a direct comparison to GW waveforms.We also show how quoted eccentricity values significantly differ when a selfconsistent definition of eccentricity is not used, investigate the robustness of our improved definition to higherorder post-Newtonian (PN) effects and spin, and discuss how our self-consistent picture of eccentricity affects the inferred fraction of measurably eccentric sources from astrophysical models.Finally, we provide a summary and concluding remarks in Section 4. Example python scripts (Vijaykumar et al. 2024a) and data products (Vijaykumar et al. 2024b) are available on Zenodo.

Definitions of eccentricity used in this work
This work refers to eccentricity in many different contexts.For ease of reading, we provide a list of these below.
1. e W03 : Eccentricities calculated from the prescription outlined in Wen (2003).This prescription has been commonly used to extract eccentricities from astrophysical simulations of binary formation (e.g.Samsing & Ramirez-Ruiz 2017;Samsing 2018;Rodriguez et al. 2018a;Zevin et al. 2019Zevin et al. , 2021a;;Kritos et al. 2022;Chattopadhyay et al. 2023).As we shall see in Section 2.2, this prescription uses the GW peak frequency as the reference frequency.
2. e kPN t : The eccentricity computed by evolving PN equations accurate to the k-th order (see e.g.Tanay et al. 2016;Klein et al. 2018).This uses the dominant, quadrupolar mode (ℓ = m = 2) GW frequency as the reference frequency.

e f22
W03 : Eccentricities calculated using the W03 prescription, but extracted at fixed 22-mode frequency.
5. e 0 : Eccentricity at a given reference point from astrophysical simulations, which is used as an initial condition (along with initial separation a 0 ) to evolve the system to a reference point suitable for GW detectors.

Peak harmonic of a GW signal in presence of eccentricity
At least in the inspiral phase of the binary, i.e. when the dynamics is well-described by PN expansion, a binary emits GWs at harmonics of its orbital frequency Note that the harmonic index n is different from the (ℓ, m) basis used to parametrize GW polarizations.At leading order, power emitted by the n-th harmonic is given by (Peters & Mathews 1963),  1) for the peak harmonic plotted along with the analytical fits from Wen (2003) and Hamers (2021).We see that the peak harmonic increases with eccentricity, and reaches values > 100 for e > 0.9.We also see that the fit in Hamers (2021) fits the numerical solution better, especially at lower values of eccentricity.
where the function g(n, e) quantifies the enhancement in the power emitted in the n th harmonic compared to a circular orbit (where there is no emission at n ̸ = 2), and is given by Here, J k (x) refers to Bessel functions of the first kind.For a given eccentricity, the peak power is emitted at a harmonic where g(n, e) has the maximum value.For ease of calculation, Wen (2003) provided a fit to the peak harmonic as a function of eccentricity, (1 + e) 1.1954  (1 − e 2 ) 3/2 . (3) Hamers (2021) noted that the above fit does not work well at low eccentricities, and provided an alternative fit to the peak harmonic,  Wen (2003) Since the value of eccentricity changes during the course of the inspiral, the value of eccentricity is often quoted at a specific reference frequency (or reference time).Wen (2003) (henceforth referred to as W03) uses the peak harmonic's frequency as a reference for defining eccentricity, and this definition has been employed by many other studies (e.g.Samsing & Ramirez-Ruiz 2017;Samsing 2018;Rodriguez et al. 2018a;Zevin et al. 2019Zevin et al. , 2021a;;Kritos et al. 2022;Chattopadhyay et al. 2023).Given semimajor axis a and total mass M , the average orbital frequency f orb is

Eccentricity extraction prescription of
This means that the peak GW frequency is, In astrophysical simulations, we are generally interested in obtaining eccentricity at a specified reference frequency, given an initial value of semi-major axis a 0 and eccentricity e 0 .These initial values can be thought of as arising either from some fiducial distribution of eccentricities and separations, or from stopping conditions in dynamical simulations.W03 uses the leading (Newtonian) order evolution (Peters 1964) For a given value of the reference frequency f peak,GW , the eccentricity can be obtained by solving Eq. ( 6) and Eq. ( 7) simultaneously.

Possible issues with the W03 prescription
While the W03 procedure has been used in deriving eccentricities at a particular GW frequency in many simulations, it has a couple of drawbacks: 1.The reference frequency is taken to be the peak frequency of the orbit.This is inconsistent with recommendations in the literature attempting to standardize the definition of eccentricity in waveforms (Ramos-Buades et al. 2022;Shaikh et al. 2023).Specfically, these works motivate using the orbit-averaged 22 (i.e.ℓ = 2, m = 2) mode frequency as a consistent reference point, on account of its smooth and monotonic evolution over the course of the inspiral.
2. Since all equations are based on Peters (1964), they are restricted to Newtonian (i.e.leading) order.This is not a good approximation of the evolution of the orbit, especially for systems close to merger.Specifically, the ordinary differential equations describing the evolution of eccentricities and the 22 mode frequencies have been calculated to several PN orders beyond the Newtonian order, also with the inclusion of spins.Furthermore, this equation breaks down at very high eccentricity.
In the rest of our paper, we shall formulate a procedure that can be applied to any set of scattering experiments to extract eccentricities, and conduct sanity checks of our procedure.In order to accomplish this, we shall solve the problems listed above2 .

A CONSISTENT PRESCRIPTION FOR CALCULATING ECCENTRICITY
As shown in Figure 1, the peak harmonic of an eccentric inspiral takes higher values for higher eccentricities.For low values of eccentricity (e ≲ 0.3), the peak harmonic is n peak = 2, and hence the 22-mode frequency is the same as frequency of the peak harmonic.However, the difference between the peak harmonic frequency and the 22-mode frequency increases very quickly, with f peak = 10f orb = 5f 22 at e ≈ 0.7.This discrepancy means that the eccentricities derived from astrophysical simulations of different binary formation channels is defined at a different reference frequency as compared to the definition of reference frequency built into waveforms as per recommendations of Shaikh et al. (2023).As is also evident, this discrepancy is worse for higher eccentricities, which are arguably more interesting from an astrophysical standpoint.
To remedy this discrepancy, the definitions of eccentricity from the simulations and waveforms need to be placed on the same footing.Hence, instead of defining eccentricity at the frequency of the peak harmonic, we define eccentricity at the orbit-averaged 22 mode frequency.Let M = m 1 + m 2 be the total mass of the binary and η = (m 1 m 2 )/M 2 be the symmetric mass ratio, where m 1 and m 2 are the component binary masses.Given an initial value of semimajor axis a 0 and eccentricity e 0 for a particular binary from a simulation, we evolve differential equations in e and v given by PN formulae, where ϵ = √ 1 − e 2 , and α i , β i are PN coefficients that depend on e and η (Peters 1964;Junker & Schaefer 1992;Rieth & Schäfer 1997;Gopakumar & Iyer 1997;Arun et al. 2009) 3 .For the purposes of this work, we will use equations accurate up to 2PN order (i.e.terms up to O v4 ), and will neglect the effect of spin.We note that the full expressions are known up to 3PN order (Arun et al. 2009) including the effect of spin (Henry & Khalil 2023), but we find that restricting to 2PN gives results accurate enough for our purposes 4 .Evolving the above equations in time starting from the initial conditions, we extract eccentricities at an appropriate reference frequency.This way of defining eccentricity makes it consistent with the one defined in waveforms, modulo systematics in modeling the GW emission itself.In the next subsection, we compare estimates derived from the prescription mentioned above (which we denote by e 2PN t ) to those derived from the W03 prescription (denoted as e W03 ) as used in the astrophysical simulations literature, and also investigate potential biases due to incomplete physics.

Comparison with the W03 prescription
As an illustration of our prescription, we consider binary mergers simulated with the Cluster Monte Carlo (CMC) code (Rodriguez et al. 2022).CMC performs long-term evolution of a globular cluster using a Hénon type Monte Carlo algorithm (Hénon 1971a,b;Joshi et al. 2000Joshi et al. , 2001;;Fregeau et al. 2003;Fregeau & Rasio 2007;Chatterjee et al. 2010Chatterjee et al. , 2013;;Pattabiraman et al. 2013;Rodriguez et al. 2015) given a set of initial conditions for the cluster.While we will only consider the CMC simulations as an example in this section, we again note that the W03 prescription has been in many other works (see e.g.Samsing & Ramirez-Ruiz 2017;Samsing 2018;Rodriguez et al. 2018a;Zevin et al. 2019;Kritos et al. 2022;Chattopadhyay et al. 2023).
A catalog of CMC simulations on a grid of initial cluster masses, metallicity, radii, and galactocentric distances is publicly available (Kremer et al. 2020;Kremer et al. 2020).This catalog also includes information about compact object mergers that take place in these clusters.Merging compact binaries are extracted from these simulations using the following procedure: 1.During strong GW encounters, the approximated N -body dynamics (including leading-order PN corrections that account for orbital energy dissipation from GW emission (Rodriguez et al. 2018b)) is evolved until two of the compact objects are in a bound state and reach a critical separation (typically 100M ), at which point the binary is assumed to definitively lead to a GW-driven merger and the properties of the binary are recorded.) at a reference frequency of 10 Hz set by each definition.The left panel shows a scatter plot of the two quantities.We see that there is considerable disagreement on the derived values of eccentricity, especially at e ≳ 0.2.This can be understood as the effect of the peak frequencies being shifted to higher n peak .The middle panel shows the histogram of (log) eccentricities for GW captures in globular clusters that occur during strong, few-body resonating interactions, derived using the W03 prescription and the prescription described in this work.The histograms do not deviate significantly from each other at low values of eccentricity, but are significantly different at e ≳ 0.2 (shaded region).The right panel shows the histogram of eccentricities in this e ≳ 0.2 region, where the discrepancies are more evident.Notably, the peak at e ≈ 1 is no longer present in the histogram; this peak is an artifact of the procedure used to extract eccentricity in astrophysical simulations.
2. With this pair of eccentricity and semimajor axis, the eccentricity at a given reference frequency is extracted using the prescription of W03.
Instead of relying on the W03 definition of eccentricity, we apply the procedure outlined in Section 3 to each binary extracted from the simulations.We then calculate eccentricities at a fixed 22-mode reference frequency of 10 Hz.In the left panel of Figure 2, we show a scatter plot of eccentricities derived from the two different definitions.Here, e W03 corresponds to the eccentricity as calculated by the W03 prescription, whereas e 2PN t corresponds to the one calculated using the prescription presented in this work.At low eccentricities, the two definitions agree with each other.However, as one increases the eccentricity, the disagreement between the two increases quite rapidly.This is a direct result of the peak frequency getting shifted to higher harmonics as one increases eccentricity; for higher eccentricity, a given peak frequency f peak,GW will correspond to a lower f 22 (since 1 2 n peak f 22 = f peak,GW ), hence describing the e 2PN t at a much lower reference frequency.This also means that the binary will lose more of its eccentricity by the time it reaches f 22 = 10 Hz.As a consequence, e 2PN t is always less than e W03 .
Notably, all the points that have e W03 ≈ 1 have e 2PN t that range from 0.4 − 0.8.We note that in many places in the literature, this e ≈ 1 peak has been interpreted to be binaries that form "in-band" and merge quickly after getting captured; these are classically hyperbolic sys-tems that are captured with a periapse frequency above the reference frequency considered (e.g., 10 Hz) via orbital energy loss through GW emission5 .While that interpretation still remains qualitatively correct, our analysis shows that the value of eccentricity that should be supplied to waveforms in comparison studies is not e ≈ 1 but rather is shifted to much lower values.
The middle panel of Figure 2 shows histograms of e W03 and e 2PN t calculated from the CMC simulations, with appropriate weights applied to each sample based on cluster mass, metallicity, and detection probability6 .As would be expected from the preceding discussions, the histograms match for low values of eccentricity, and start deviating at e ≳ 0.2.Most strikingly, the peak at e W03 ≈ 1 in the e W03 histogram is not present in the e 2PN t histogram.These observations are more evident in the right panel of Figure 2, where we show the histograms of eccentricities for e > 0.2.Since the waveform changes significantly as a function of eccentricity, one would expect these discrepancies to impact any analysis that assumes consistency between the definitions of ec-centricities in waveforms and astrophysical simulations.We will comment on this briefly in Section 3.4.
In the preceding discussion, we have assumed that the reference frequency is defined in the source-frame of the binary as opposed to the detector-frame.We make this choice since the original calculation of eccentricity using the W03 prescription with the CMC outputs made this choice.Strictly, this is incorrect and adds another significant level of inconsistency in the conventions for eccentricity.For completeness, we also show the eccentricity extracted at a reference frequency corresponding to a constant f 22 M = 1000 HzM ⊙ in Figure 3, as f 22 M does not depend on redshift.At a f 22 corresponding to constant f 22 M , all (black hole) binaries would be at the same epoch in their orbital evolution, i.e. they will all be at a fixed number of cycles before merger.This criterion is easy to implement while estimating parameters from BBH signals, and also allows for a more direct comparison of eccentricities for systems with different masses.Defining the eccentricity at a fixed M f22 has a few advantages: (i) M f22 is independent of redshift, (ii) it defines eccentricity consistently at a fixed number of cycles before merger, and (iii) it is easy to implement while estimating source parameters from BBH signals.We recommend that all astrophysical simulations provide e 2PN t estimates extracted at a fixed M f22 for easy comparison with estimates of eccentricity derived using GW waveform models.
It is natural to ask whether the discrepancy is primarily driven by the mismatch in conventions for the reference frequency, or due to only accounting for the leading order PN term in the W03 prescription.The results in the preceding paragraphs would lead us to believe that it is the former.We test this in another way: since we know n peak (e), we can also adjust the eccentricities calculated via W03 as a post-processing step by assuming that the e W03 is defined at a reference frequency f 22 = 2 × f peak /n peak (e), and evolve to a fixed reference (22-mode) frequency using PN equations (Eqs.( 8)-( 10)) but only up to the leading order.Let the eccentricity calculated using this modified version of the W03 procedure be e f22 W03 .We find that e f22 W03 ≈ e 2PN t for the binaries we consider7 , with a ∼ 5% scatter that could be attributed to non-inclusion of higher-order PN terms (see also Section 3.2).This further confirms that the discrepancy in the calculated eccentricities is primarily driven by the difference in conventions for the reference frequency.

Effects of ignoring higher order PN terms and spin effects
We now systematically investigate if there is a relative advantage in using the 2PN equations above just the leading order (0PN) evolution equations.For an equal mass binary of total mass M = 50 M ⊙ , we consider three different 22-mode starting frequencies f 0 = 5, 10, 15 Hz, and calculate e 2PN t and e 0PN t at a reference frequency of 20 Hz for a range of initial eccentricities e 0 .The fractional differences between e 2PN t and e 0PN t estimates are plotted in the top panel of Figure 4, and show an increasing trend for lower values of f 0 , with a discrepancy as high as 12% for f 0 = 5 Hz.For a given f 0 , the fractional differences are lower for higher e 0 .In typical applications, one might want to calculate eccentricities at a given reference frequency for much lower values of f 0 .While a ∼ 10% difference might not be relevant while comparing eccentricity estimates from data to simulations, we recommend using e 2PN t estimates since the costs of evaluating e 0PN t and e 2PN t are similar (see Appendix C).
Our prescription in Section 3 does not include the effects of spin.To see if this is a valid choice, we use the same mass, f 0 , f ref , and e 0 values as in the previous paragraph and calculate e 2PN t with and without including spin terms in the evolution for two representative dimensionless (aligned) spin values of χ 1,z = χ 2,z = χ z = [0.7,0.99].For the spin terms in the evolution, we use 2PN accurate terms from et (2018).
fractional between e 2PN t estimates with and without spin effects are shown in the lower panel of Figure 4. We see that the discrepancies are at worst ∼ 10 −5 , shifting to even lower values when considering smaller fiducial spins or higher f 0 .We hence conclude that the non-spinning approximation is sufficient.

Comparison to estimates derived from gw_eccentricity package
As mentioned earlier, different GW waveform approximants have different conventions to define eccentricity.To solve this issue, Shaikh et al. (2023) proposed a standardized definition of eccentricity that can be applied to any waveform approximant, e gw , which is given by where and ω p,a 22 (t) are the 22-mode frequencies through the pericenters and apocenters, respectively.This standardized definition is proposed in order to extract eccentricity directly from observables, in this case the gravitational waveform.More details on the motivations for this definition and the various methods of calculating ω p,a 22 (t) can be found in Shaikh et al. (2023).We also briefly investigate the differences between e gw derived from different waveform approximants in Appendix B.
While e gw does asymptote to the PN definition of eccentricity in the small eccentricity limit as well as the large separation limit, these two definitions may differ in general.For instance, waveforms could include more accurate dynamics at high eccentricity and closer to merger as compared to our 2PN prescription.Instead of the 2PN accurate equations, we could hence imagine replacing the prescription in Section 3 by generating waveforms at the specified initial epochs, and let gw_eccentricity calculate the eccentricities at the required reference frequency.
However, since calculating e gw involves computing the waveform starting at a specified initial frequency/time, the calculation would be on the order of a hundred times slower (see Appendix C) than computing the PN eccentricity.The cost of generating waveforms also increases with decreasing starting frequency.Since astrophysical simulations require calculating eccentricites of a large number of binaries, calculating the Shaikh et al. (2023) definition of eccentricity would be computationally inefficient.Therefore, we now compare the differences be- tween the two eccentricity definitions, e gw and e 2PN t , and quantify the reliability of e 2PN t eccentricity estimates from astrophysical simulations.For the investigations that follow, we will use the SEOBNRv4EHM (Ramos-Buades et al. 2022) waveform approximant for the calculation of e gw .Using the same mass, f 0 , and e 0 values as in Section 3.2, we calculate e 2PN t and e SEOB gw at a reference frequency of 20 Hz. Figure 5 illustrates the differences between these.While small (a few tens of a percent at worst), the differences do accumulate as the binary evolves for a longer period of time with f 0 values closer to the reference frequency giving more agreement.Since the cost of computing e gw is a factor of ∼ 100 greater as compared to calculating e 2PN t (see Appendix C), we conclude that e 2PN t is sufficient, especially in regimes where PN theory is valid.

Implications for existing work
A number of works have ascertained the detectability of eccentric mergers within the LVK band as well as in other detector bands, using data from astrophysical simulations such as CMC.As we have shown in the previous sections, the waveform eccentricities are not in agreement with the eccentricities given by simulations.Hence, we would expect these detectability estimates to change when defining eccentricity consistently in the simulations.For instance, Zevin et al. (2021a) found that ≈60% of potentially detectable eccentric sources with e ≳ 0.2 would be missed due to non-inclusion of eccentric templates in search pipelines.Using the same dataset as Zevin et al. (2021a) and changing the eccentricities from the e W03 definition to the e 2PN t definition, we find that only 30% of potentially detectable eccentric sources with e ≳ 0.2 are missed-a discrepancy by a factor of ≈ 2. The decrease in the number of missed signals is due to the shifting of the histogram of eccentricities towards lower values, which are better approximated by non-eccentric templates and hence less likely to be selected against.As another illustration, we calculate limits on the cluster branching fraction β c i.e the fraction of events in the total population coming from clusters using the same method and dataset as used in Zevin et al. (2021b).We find that assuming N ecc = 0 eccentric detections in GWTC-3, we would have inferred β c < 86% (90% credible level) using the e W03 estimates, and infer β c < 63% using the e 2PN t estimates.On the other hand, assuming N ecc = 1 gives β c > 10% using the e W03 estimates and β c > 8% using the e 2PN t estimates.Similarly, any study that uses the detected eccentricities (or lack thereof) to place constraints on the overall merger rate in eccentric systems in astrophysical environments will also be affected by these discrepancies.
For instance, Abac et al. (2023) use the non-detection of orbital eccentricity from a waveform-independent search to place an upper limit of 0.34 Gpc −3 yr −1 on the merger rate of binaries in dense stellar clusters, assuming a population model from CMC simulations.This population model is taken from Zevin et al. (2021a) with eccentricities defined using the W03 prescription, with additional conditioning on the mass (M > 70) and eccentricity (e < 0.3) parameter ranges of their search.From Figure 2, we see that the relative number of mergers having e < 0.3 increases when accounting for the prescription in this work.This means that an upper limit derived using the e 2PN t estimates from this work would be lower than the ones quoted in Abac et al. (2023).
We also note that Romero-Shaw et al. ( 2022), which uses signatures of non-zero eccentricity in four GW events to constrain the fraction of BBH mergers formed dynamically, does account for the disparate conventions for eccentricity while making their comparisons.They do so in an approximate fashion, similar to the calculation of e f22 W03 described in Section 3.1.We do not expect conclusions reported there to change significantly.

CONCLUSION
Orbital eccentricity in a compact binary is believed to be a telltale sign of its formation history.Consequently, there has been progress on building waveform models for eccentric sources as well as on simulating expectations from astrophysical formation channels.In this work, we note that the conventions used to define eccentricity in the GW waveform and astrophysical simulation communities differ from each other.The major difference in these conventions lies in the choice of the reference frequency at which eccentricity is defined-GW waveforms typically define eccentricity at a fixed orbit averaged 22-mode GW frequency, while astrophysical simulations define eccentricity at a fixed peak GW frequency following the W03 prescription (described in Section 2.2).Further, we delineate an eccentricity calculation procedure that can be applied to data from astrophysical simulations to extract eccentricities consistent with conventions laid down by the GW waveform community.This procedure involves evolution of ordinary differential equations describing the time-evolution of eccentricity and frequency obtained from PN expansions, and extracting eccentricity consistently at the orbit-averaged 22 mode frequency.
Using compact binaries from the publicly-available CMC catalog as an illustration, we find that there is a discrepancy between eccentricity estimates derived using the W03 prescription and the one described in this work.The discrepancy is small for e ≲ 0.1, but is large for higher values of eccentricity.Notably, all binaries with e W03 ≈ 1 (thought to be GW captures happening in the LVK band) have e 2PN t between 0.4-0.8.We investigate the effect that adding 2PN corrections to the leading order evolution has on derived eccentricities, and also the effects of not including spins in our prescription.We also compare our e 2PN t estimates to e gw estimates derived from the gw_eccentricity package using the SEOBNRv4EHM waveform.We find that the e 2PN t estimates should be sufficiently accurate in extracting eccentricities, and also computationally efficient to implement as a post-processing step in simulations.However, as outlined in Shaikh et al. (2023), gw_eccentricity is the most accurate way of extracting eccentricity and should be used when high-precision eccentricity estimates are necessary.
The different conventions for defining eccentricity currently employed in the literature mean that one should be careful in comparing eccentricity estimates/upper limits derived from real data to astrophysical simulations, creating mock populations of eccentric sources based on astrophysical expectations, and other investigations that pair predictions from astrophysical models with GW data analysis.As an illustration, we show that the detectability estimates of eccentric sources and branching fractions from Zevin et al. (2021a) will change when ensuring a consistent definition of eccentricity.We recommend that catalogs created from astrophysical simulations include e 2PN t estimates (extracted at a constant M f 22 ) for clarity and ease of direct comparison to GW observations.APPENDIX

A. PN FORMULAE FOR ECCENTRICITY AND VELOCITY EVOLUTION
The PN formulae (accurate up to 2PN) for the evolution of eccentricity and velocity are given below where  The prescription of calculating eccentricity e gw , proposed by Shaikh et al. (2023) and implemented in the gw_eccentricity python package only uses the GW waveform as seen by the detector as input, as described in , should be consistent with the input eccentricity e0 (black line).EccentricTD (blue) does not agree for any value of e0, while SEOBNRv4EHM (yellow) does agree for all e0 values considered.We also see that the final eccentricity reported for the given reference frequency f ref between the two waveforms (blue and yellow dashed lines) do not agree.
Section 3.3.This definition is beneficial since eccentricity is now linked to a direct observable.However, different waveforms will, in general, yield different values of e gw for the same binary.This could be due to a variety of reasons: differences in the modeling of the dynamics, different internal definitions of eccentricity in the model, or even choices made while implementing the approximant in some codebase such as LALSimulation (LIGO Scientific Collaboration et al. 2018).Figure 6 illustrates this point by comparing e gw derived from two waveforms, EccentricTD (Tanay et al. 2016) and SEOBNRv4EHM (Ramos- Buades et al. 2022).For this illustration, we assume an equal mass binary with a total mass of M = 50 M ⊙ .We vary the initial eccentricity e 0 (defined at initial frequency f 0 = 5 Hz) supplied to the model between 10 −2 and 0.8, and extract e gw from both EccentricTD e EccTD gw and SEOBNRv4EHM e SEOB gw , at both the initial frequency f 0 and also at f ref = 20 Hz.We see that e SEOB gw (f 0 ) matches all input values e 0 , whereas e EccTD gw (f 0 ) differs from e 0 by a few percent.This is in agreement with investigations performed in Shaikh et al. (2023), where they found that the e SEOB gw (f 0 ) = e 0 trend is maintained for all e 0 > 10 −5 , whereas e EccTD gw deviates from the e EccTD gw (f 0 ) = e 0 trend.We also see that there is a large discrepancy (factor of a few) between e EccTD gw (f ref ) and e SEOB gw (f ref ), strengthening the evidence for differences in the ways each of these waveform approximants define eccentricity.We use the SEOBNRv4EHM approximant for our comparisons in Section 3.3.

C. COMPUTATIONAL COST OF DIFFERENT ECCENTRICITY PRESCRIPTIONS
Astrophysical simulations keep track of many evolving binary systems, where it is often computationally infeasible to compute waveforms for all systems.In Figure 7, we demonstrate the speed-up in eccentricity calculation at the reference frequency, f ref , by computing the PN eccentricities, e kPN t , as opposed to the standardized waveform definition of eccentricity, e gw .The top left of Figure 7 demonstrates the average time to compute the 2PN eccentricity at various initial eccentricities, e 0 , and frequencies, f 0 .The top right plot illustrates that computing the waveform eccentricity, e gw using the SEOBNRv4EHM waveform model takes on the order of a hundred times longer on average.The bottom plots in Figure 7 demonstrate that using the 0PN eccentricity has a comparable cost to using the 2PN eccentricity, whereas including spin effects in the PN equations is ∼1.5 times slower.

Figure 1 .
Figure 1.Numerical solution of Eq. (1) for the peak harmonic plotted along with the analytical fits fromWen (2003) andHamers (2021).We see that the peak harmonic increases with eccentricity, and reaches values > 100 for e > 0.9.We also see that the fit in Hamers (2021) fits the numerical solution better, especially at lower values of eccentricity.

Figure 2 .
Figure 2. Comparison of orbital eccentricities in CMC simulations derived using the W03 prescription (e W03 ), to ones derived using the prescription in this work (e 2PN t

PDFFigure 3 .
Figure 3. Histogram of eccentricities extracted at M f22 = 1000 M⊙Hz.Defining the eccentricity at a fixed M f22 has a few advantages: (i) M f22 is independent of redshift, (ii) it defines eccentricity consistently at a fixed number of cycles before merger, and (iii) it is easy to implement while estimating source parameters from BBH signals.We recommend that all astrophysical simulations provide e 2PN

Figure 4 .
Figure 4.The top panel compares e t

Figure 5 .
Figure 5.Comparison between the eccentricity definition given byShaikh et al. (2023), calculated using the waveform approximant SEOBNRv4EHM (e SEOB gw ; dotted lines), and e 2PN t (dashed lines), as a function of the initial eccentricity, e0, for three initial frequencies, f0 = 5, 10, 15 Hz.We assume an equal mass non-spinning binary with total mass M = 50 M⊙ for this purpose.

Figure 6 .
Figure6.Comparison of internal definitions of eccentricity for different waveform models as a function of the initial input eccentricity, e0 with a fiducial input frequency of 5 Hz.Initial eccentricities using the standardized definition, egw(f0) (dotted lines), should be consistent with the input eccentricity e0 (black line).EccentricTD (blue) does not agree for any value of e0, while SEOBNRv4EHM (yellow) does agree for all e0 values considered.We also see that the final eccentricity reported for the given reference frequency f ref between the two waveforms (blue and yellow dashed lines) do not agree.

Figure 7 .
Figure7.Average time to compute the 2PN eccentricity left) and the relative time to compute the eccentricities using the SEOBNRv4EHM waveform (top right), the eccentricity (bottom left), and the 2PN eccentricity for spinning binaries (bottom right).The grid covers various input eccentricity, e0, and frequency, f0, values.Calculating eccentricities using waveforms takes ∼ 100× longer than calculating eccentricities using our PN approximation.