Gravitational waveforms capturing binarys evolution through the earlyinspiral phase play a critical role in extracting orbital features that nearly disappear during the lateinspiral and subsequent merger phase due to radiation reaction forces for instance the effect of orbital eccentricity. Phenomenological approaches that model compact binary mergers rely heavily on combining inputs from both analytical and numerical approaches to reduce the computational cost of generating templates for data analysis purposes. In a recent work Chattaraj et al. Phys. Rev. D 106 124008 2022 arXiv2204.02377grqc we demonstrated construction of a dominant quadrupole mode inspiralmergerringdown IMR model for binary black holes BBHs on elliptical orbits. The model was constructed in timedomain and is fully analytical. The current work is an attempt to improve this model by making a few important changes in our approach. The most significant of those involves identifying initial values of orbital parameters with which the inspiral part of the model is evolved. While the ingredients remain the same as in arXiv2204.02377grqc resulting waveforms at each stage seem to have improved as a consequence of new considerations proposed here. The updated model is validated also against an independent waveform family resulting overlaps better than sim 96.5 within the calibrated range of binary parameters. Further we use the prescription of the dominant mode model presented here to provide an alternate but equivalent model for the dominant quadrupole mode and extend the same to a model including the effect of selected nonquadrupole modes. Finally while this model assumes nonspinning components we show that this could also be used for mildly spinning systems with component spins anti aligned w.r.t the orbital angular momentum.
['General Relativity and Quantum Cosmology', 'Astrophysics - High Energy Astrophysical Phenomena']
An improved IMR model for BBHs on elliptical orbits
{'An improved IMR model for BBHs on elliptical orbits': "Pratul Manna, 1, ∗ Tamal RoyChowdhury, 2, † and Chandra Kant Mishra 1, 3, ‡ 1 Department of Physics, Indian Institute of Technology Madras, Chennai 600036, India 2 University of Wisconsin Milwaukee, Milwaukee, WI 53201, USA 3 Centre for Strings, Gravitation and Cosmology, Department of Physics, Indian Institute of Technology Madras, Chennai 600036, India (Dated: September 18, 2024)  \nGravitational waveforms capturing binary's evolution through the early-inspiral phase play a critical role in extracting orbital features that nearly disappear during the late-inspiral and subsequent merger phase due to radiation reaction forces; for instance, the effect of orbital eccentricity. Phenomenological approaches that model compact binary mergers rely heavily on combining inputs from both analytical and numerical approaches to reduce the computational cost of generating templates for data analysis purposes. In a recent work, Chattaraj et al., Phys. Rev. D 106, 124008 (2022) [1], we demonstrated construction of a dominant (quadrupole) mode inspiral-merger-ringdown (IMR) model for binary black holes (BBHs) on elliptical orbits. The model was constructed in time-domain and is fully analytical. The current work is an attempt to improve this model by making a few important changes in our approach. The most significant of those involves identifying initial values of orbital parameters with which the inspiral part of the model is evolved. While the ingredients remain the same as in Ref. [1], resulting waveforms at each stage seem to have improved as a consequence of new considerations proposed here. The updated model is validated also against an independent waveform family resulting overlaps better than ∼ 96.5% within the calibrated range of binary parameters. Further, we use the prescription of the dominant mode model presented here to provide an alternate (but equivalent) model for the (dominant) quadrupole mode and extend the same to a model including the effect of selected non-quadrupole modes. Finally, while this model assumes non-spinning components, we show that this could also be used for mildly spinning systems with component spins (anti-) aligned w.r.t the orbital angular momentum.", 'I. INTRODUCTION': "Since the discovery event, GW150914 [2], the LIGOVirgo-KAGRA (LVK) collaboration has reported nearly hundred compact binary mergers observed during the first three observing runs (O1-O3) [3-6]. These numbers have doubled since and the list continues to be dominated by signals identified as mergers of black holes in a binary; see for instance, the GWOSC 1 page [7] that lists all reported events. While these (now also routine!) observations continue to help improve our understanding of compact binary physics and astrophysics, their origins remain unknown [8-10]. Astrophysical environments where a binary is formed and processes through which it is formed may leave imprints on binary's mass and spin parameters. However, the available statistics needs to grow in order to make inferences concerning binary's origin based on mass and spin measurements alone [11]. Eccentricity, on the other hand, can be a unique tool to identify binary's origins, as dynamically formed binaries may still retain residual orbital eccentricities [12-14] when observed in ground based detectors currently operating.  \nCurrent template-based search pipelines make use of circular templates due to the expected circularisation of most compact binary orbits caused by radiation reaction forces [15] as they enter the sensitivity bands of ground based detectors such as LIGO [16] and Virgo [17]. However, binaries formed through the dynamical interactions in dense stellar environments are likely to be observed with residual eccentricities e 20Hz ∼ 0 . 1 [10, 12, 18]. In fact, the first binary merger event involving an intermediate mass black hole, GW190521 [19], is likely an eccentric merger [9, 20] (see also Refs. [21-26] discussing events with eccentric signatures). While quasi-circular templates should be able to detect systems with initial eccentricities e 10Hz ≲ 0 . 1, binaries with larger eccentricities would require constructing templates including the effect of eccentricity [27, 28]. Moreover, the presence of even smaller eccentricities ( e 10Hz ∼ 0.01 - 0.05) can induce significant systematic biases in extracting source properties [29](see also, Refs. [30, 31]). Furthermore, next generation ground-based detectors, Cosmic Explorer [32-34] and Einstein Telescope [35, 36], due to their low frequency sensitivities, should frequently observe systems with detectable eccentricities [37, 38].  \nEven though inspiral (I) waveforms from eccentric binary mergers involving non-spinning compact components are sufficiently accurate [39-45], waveform models including contributions from merger and ringdown (MR) stages compared to quasi-circular counterparts  \nare significantly less developed. 2 Numerous efforts toward constructing eccentric inspiral-merger-ringdown (IMR) waveforms, useful for data analysis purposes, are underway [1, 48-52]. However, these efforts do not include important physical effects such as spins (pointing along or away from binary's orbital angular momentum) or higher order modes. Dominant mode (or quadrupole mode) models for eccentric BBHs with component spins (anti-)aligned w.r.t the binary's orbital angular momentum were recently developed in Refs. [53-58]. Since most mergers observed so far are consistent with a zero-effective spin [59-61], models neglecting spin effects can still be useful [49]. 3 In addition, modeling of higher order modes also seems necessary as Refs. [1, 62] argue. Eccentric versions of the effective-one-body (EOB) waveforms including higher modes [63, 64] and an eccentric numerical relativity (NR) surrogate model [65, 66] also became available in the past couple of years. Alternatively, sub-optimal methods (with little or no dependence on signal model being searched) may be used for detecting an eccentric merger. Although, these methods are sensitive to high mass searches (typically ≳ 70 M /circledot ) [10, 18], while most observed events have a mass smaller than this limit [5, 60]; see, for instance, Fig. 3 of [67].  \nThe present work follows our first paper [1] and can be viewed as an update to the same; referred to as Paper I here onwards. In Paper I, construction of hybrid waveforms by combining post-Newtonian (PN) waveforms with NR simulations through a least-squares minimization was demonstrated. These hybrids were then used as target models to produce a fully analytical dominant (or quadrupole) mode model by matching an eccentric PN inspiral with a quasi-circular mergerringdown waveform. Subsequently, the performance of the model was checked both against the target hybrids used in training the model as well as against an independent family of waveforms [50]. It was shown that overlaps between target hybrids and the model significantly improved compared to those against a circular template, at least at the low end of the binary masses and for small eccentricities considered there (See for instance, Figure 9 of Paper I). The current work aims to improve the model presented in Paper I in the view of efforts such as those of Refs. [68, 69]. References [68, 69] develop a standard, gauge-independent prescription for defining orbital eccentricity using a suitable combination of dominant mode GW frequency. It should be noted that, in GR, eccentricity is not  \nuniquely defined (see for instance Ref. [70]) although at the leading (Newtonian) order there is consensus. The definition of Refs. [68, 69] reduces to the Newtonian value both in small and large eccentricity limit and is also independent of gauge-ambiguities. This motivates us to employ this definition of eccentricity in our model and investigate into the improvements in its performance.", 'A. Gauge invariant definition of eccentricity': 'Eccentricity is not uniquely defined in GR and thus templates computed within the framework of GR may have forms different from the observed data. Both perturbative and numerical solutions describing the compact binary dynamics use gauge-dependent constructs including the very definitions of orbital parameters that are evolved. These choices are almost never identical in any two approaches which naturally leads to inconsistencies between different models. To get rid of the ambiguity associated with the definition of eccentricity, Refs. [68, 69] proposed a new definition of eccentricity based on the GW frequency data. At 0PN order, this definition exactly reduces to the Newtonian definition of eccentricity ( e t ). We reproduce the necessary relations below. This new eccentricity for an observed GW ( e gw ) signal is defined as  \ne gw = cos( ψ/ 3) -√ 3 sin( ψ/ 3) , (1)  \nwith,  \nwhere,  \ne ω 22 = √ ω p 22 -√ ω a 22 √ ω p 22 + √ ω a 22 , (3)  \nwhere, ω p 22 and ω a 22 refer to the ( /lscript = 2 , | m | = 2) mode periastron and apastron frequencies, respectively and are functions of time. Since this new eccentricity is written in terms of frequencies that can be measured, it is also free from any gauge-ambiguities. This definition is employed in the gw eccentricity package provided by Ref. [68].', 'B. Summary of the current work': "Considerations related to a change in the definition of eccentricity such as the ones proposed by Refs. [68, 69] demand revisiting each element of model construction presented in Paper I which we intend to closely follow here. We start by comparing the PN prescription and NR simulations used in constructing hybrids in Paper I. The new definition of eccentricity introduced in  \nψ = arctan ( 1 -e 2 ω 22 2 e ω 22 ) , (2)  \nRefs. [68, 69] is used to identify a set of reference values for eccentricity ( e ref ), mean anomaly ( l ref ) and GW frequency ( f ref ) with which the overlap between the PN and NR data is maximum in a time-window where the two are expected to give accurate predictions. Figure 1 plots the two, PN and NR model, together with the window returning maximum overlap between them.  \nThe PN model is evolved using reference values, ( e ref , l ref , f ref ), as an initial set (i.e at the start of the PN model) and matched with NR simulations in a hybridization window following the prescription of Refs. [1, 71]. Figure 3 displays one of the hybrids and compares it with corresponding NR simulation. Table I lists all hybrids constructed here along with starting values of eccentricity ( e 0 ), mean anomaly ( l 0 ) and a frequency dependent PN parameter ( x 0 ) which is related to the GW frequency of the dominant mode ( f 0 ) via x 0 = ( πMf 0 ) 2 / 3 . Finally, a fully analytical dominant mode model is obtained by matching an eccentric PN inspiral [41] with a quasi-circular merger-ringdown model [72] in Sec. III. Note however, as in Paper I, a new set of hybrids (including only /lscript = 2 , | m | = 2 mode), with PN model purely based on Ref. [41] (termed EccentricTD within LIGO Algorithmic Library [73]) are used in training the model to minimize the difference between the target and the model which in turn uses the PN prescription of Ref. [41]. (See Fig. 4.)  \nFigure 7 evolves the model using the initial set of parameters for three different hybrids listed in Table I and simply plots it against the amplitude (top-left), frequency (top-right) and plus polarisation data from the hybrids. While these provide a visual proof of closeness of the model with target hybrids, right panel of Fig. 8 plots the overlap (maximized over a simple time and phase shift) between the model and the set of hybrids used in calibrating the model (labelled as training set in Table I) (thin lines). Overlaps with two of the hybrids not used in building the model (labelled as testing set in Table I) are also plotted (thick lines). For comparison, we also plot (left panel) the overlap of (all training and two testing) hybrids with quasi-circular templates of Ref. [72] (termed SEOBNRv5 ). Clearly, the model outperforms its circular counterpart and recovers the target eccentric hybrids with an accuracy better than ∼ 96 . 5% for almost the entire range of parameter space spanned by the training set hybrids. Note that both the target and template waveforms used in Fig. 8 involve only the dominant ( /lscript = 2 , | m | = 2) mode. Additionally, Figure 8 only shows overlap between the hybrids and the model for fixed values of parameters, ( e 0 , l 0 , x 0 ), associated with each hybrid. . Like not in a seperate folder. But we can de not in a folder.e nontfolfter a The expectation, that the model shall produce reliable waveforms for any value of the parameters in the range it is calibrated, should also be verified. To test this we compare our model against an independent family of waveforms  \ndiscussed in Ref. [58] (termed TEOBResumS-Dali ) for randomly sampled values for this set in the range 0 ≲ e 0 ≲ 0 . 3, 1 ≲ q ≲ 3, and -π ≤ l 0 ≤ π , respectively. The range of parameter values are chosen to match the parameter space spanned by the hybrids chosen for calibration, except for orbital eccentricity, which is conservatively chosen to have a maximum value of e 0 = 0 . 3 and also explores near circular cases to see if the model correctly (and gradually) reproduces the circular limit. The results are displayed in Fig. 9 and as can be seen there, model reproduces the waveforms of TEOBResumS-Dali [58] with overlaps better than 96.5% for almost the entire range of parameters considered there and confirms the suitability of the waveforms for data analysis purposes. Note however, the difference between the eccentricity scale in Fig. 8 and in Fig. 9.  \nSection IV A presents an alternate model to the one obtained in Sec. III. The primary motivation here is to use PN input waveforms of higher PN accuracy (compared to EccentricTD [41] used in constructing the model in Sec. III) to maximize on the overlaps with target models at a small cost of losing sensitivity to high eccentricity cases within the calibration range. Comparisons with TEOBResumS-Dali [58] waveform show that such a model may provide a suitable alternative to the model constructed in Sec. III. Figure 10 displays this comparison. It is interesting to note that the mismatches seem to have visibly improved compared to those with the model based on EccentricTD for small eccentricity cases. This is likely due to higher PN accuracy of the amplitude and phase used in constructing the alternate model and thus matches better with the waveform TEOBResumS-Dali [58]. For larger eccentricities, the performance of the two seem similar (see a discussion in Sec. IV A).  \nFinally, while a higher mode model can be constructed following the methods used in constructing the dominant mode model discussed in Sec. III, one may simply use the prescription for the ( /lscript = 2 , | m | = 2) mode to combine an (eccentric) inspiral and a (quasi-circular) merger-ringdown prescription for each mode to obtain an ad hoc higher mode (HM) model as was done in Paper I. The alternate model presented in Sec. IV A is extended to include selected /lscript = | m | and /lscript -1 = | m | modes and is validated against TEOBResumS-Dali [58]. These are the modes which are included in SEOBNRv5HM [72]; the quasi-circular model used for the merger-ringdown part. Note also, these are also the modes for which we construct hybrids; see Fig. 3. The model is discussed in Sec. IV B and its performance is displayed in Fig. 11 for mildly inclined systems (30 · ).  \nThe paper is structured in the following manner. In Sec. II we start by comparing waveforms from PN and NR approaches and discuss the construction of target hybrids. Next, in Sec. III we construct the waveform model  \n<!-- image -->  \nFIG. 1. Amplitude and frequency of selected modes from an eccentric NR simulation (SXS:BBH:1364) together with an eccentric PN model are plotted. The PN model is evolved assuming a fixed value for initial eccentricity ( e 0 = 0 . 172), mean anomaly ( l 0 = 2 . 681) and the frequency dependent PN parameter ( x 0 = 0 . 0391). The merger time of the NR waveform is set to zero. The initial set ( e 0 , l 0 , x 0 ) is obtained by maximizing the overlap between the PN and NR waveform in a 1000 M wide time-window, in a region where PN model is expected to agree with NR simulation. Additionally, a time shift is performed on the PN inspiral. The time-window of maximum overlap is shown as the shaded region. Binary's component mass ratio ( q ) is 2, while the total mass ( M ) and the luminosity distance ( D L ) of the binary are set to M =1 M /circledot and D L =1Mpc respectively, following the convention of SXS simulations.  \n<!-- image -->  \nby combining an eccentric PN inspiral model with a quasi-circular merger-ringdown model at a suitable point obtained by performing comparisons with target models. Subsequently, the model is validated against the target models not used in calibrating it as well as against an independent family of waveforms. In Sec. IV, we discuss an alternate model based on the prescriptions for attachment times obtained in Sec. III and subsequently extend this alternate model to include higher order modes. Finally, Sec. V presents summary of results and conclusions.", 'A. PN and NR comparisons': 'Paper I compares the PN inspiral waveforms with NR simulations. The inspiral mode amplitudes constituting 3PN inspiral waveforms assuming non-spinning binary systems on quasi-elliptical orbits were computed in Refs. [39, 42, 43]. The orbital phase was taken from Ref. [41]. GW frequency for each PN mode was obtained using the following scaling relation [69]  \nω /lscriptm ∼ m 2 × ω 22 . (4)  \nNR simulations used in comparison with PN waveforms were produced using the Spectral Einstein Code (SpEC) developed by SXS collaboration and are publicly available [48, 74]. Comparison of PN and NR prescriptions for a specific simulation was shown in Figure 1 of Paper I. Subsequently, a common region of validity was identified in which the two could be matched suitably  \nto obtain hybrids listed in Table I there. (See Sec. III A-III C of Paper I for technical details).  \nHere too we aim to compare the PN and NR prescriptions leading to the construction of hybrids. While the hybridization method is same as in Paper I, comparison (which leads to identification of a suitable window for hybridization) is performed following a slightly different approach. In paper I, PN models were simply evolved to match a set of reference orbital parameters ( e ref , l ref ) computed at a reference GW frequency ( x ref ) computed in Ref. [48] for each of the 20 simulations considered in Paper I. Here, we simply choose to compare the two waveforms in a 1000 M wide time-window, slide it over the overlapping data and compute overlaps (the match maximized over a reference time and phase shifts; see Eq. (5) below) by varying ( e ref , l ref , f ref ) trio for the PN model, where f ref is related to the PN parameter x ref via x ref =( πMf ref ) 2 / 3 . The match ( M ) between two waveforms is defined as an inner product given as  \nM ( θ 1 , θ 2 ) ≡ max φ c ,t c 〈 h ( θ 1 ) , h ( θ 2 ) e i (2 πft c -φ c ) 〉 , (5)  \nwith,  \n〈 h 1 , h 2 〉 ≡ 4 Re [ ∫ ∞ 0 df ˜ h 1 ∗ ( f ) ˜ h 2 ( f ) S h ( f ) ] , (6)  \nwhere 〈 h 1 , h 2 〉 represents the inner product between two waveforms h 1 and h 2 having unit norm and are functions of an intrinsic set of binary parameters ( θ 1 , θ 2 ). The phase φ c , time t c are measured at coalescence and S h ( f ) represents the noise in the detector (see Ref. [75]).  \nWe find the overlap is optimal for a time-window of ( -2000 M to -1000 M ) and for a given ( e ref , l ref , f ref )  \nFIG. 2. Same as Fig. 1, except that only dominant mode data is plotted and the comparison is shown for few other simulations with varying mass ratios and orbital parameters (see Table I for details). Initial values ( e 0 , l 0 , x 0 ) with which the PN model is evolved for comparisons with various simulations are listed in Table I.  \n<!-- image -->  \ntrio which becomes the starting reference ( e 0 , l 0 , f 0 ) and may be different for different NR simulations (see Table I). 4 Compared to the earlier approach, our current approach helps the construction of hybrids in at least two distinct ways. (1) The identification of hybridization window is done using quantitative measures such as overlaps and (2) the reference orbital parameters such as orbital eccentricity ( e 0 ) and mean-anomaly ( l 0 ) at a given frequency ( f 0 ) will be free from gauge-ambiguities due to the use of the definition of Refs. [68, 69]. With the suitable hybridization window identified, we can now proceed to reconstruct the hybrids.  \nFigure 1 compares the data corresponding to the NR simulation bearing simulation ID SXS:BBH:1364, for a selected set of modes chosen based on their relative significance compared to the dominant mode. 5 PN model is evolved using a set of initial parameters ( e 0 , l 0 , f 0 ) obtained using the procedure discussed above and then plotted together with the NR simulation after performing a time shift. The window giving maximum match is also displayed. For completeness we also show similar comparisons for few other simulations in Fig. 2 albeit for only the dominant mode. It is interesting to also note the time-window returning maximum match is common  \nfor all simulations.', 'B. Construction of hybrid waveforms': "Complete IMR waveforms are constructed by matching PN and NR prescriptions for set of modes included in Fig. 1, in a region where the PN prescription closely mimics the NR data following the method of Ref. [71]. These are traditionally referred to as 'hybrids'. As discussed in Ref. [71], construction of hybrids including higher modes (in the circular case) is possible by performing at least two rotations (and a time shift) so as to align the frames in which PN/NR waveforms are defined. 6 This argument was simply extended to the case of eccentric orbits in Paper I, assuming that the effect of marginalising over parameters such as eccentricity and mean anomaly will not significantly affect the hybridization. As discussed above, for the current work we simply adopt the hybridization procedure of Paper I. 7 The prescription for construction of hybrids is discussed in detail in Ref. [71] as well as in Paper I, nevertheless, we reproduce some of the steps here for completeness.  \nFIG. 3. A hybrid model constructed by matching a NR simulation (SXS:BBH:1364) with a PN model (evolved using parameters consistent with the simulation) in a time-window where the two are expected to correctly predict the binary dynamics. For comparison the NR data is also plotted. The black dashed line marks the beginning of the NR waveform and the shaded light-red region t ∈ ( -2000 M, -1000 M ) shows the matching window. Overlapping hybrid and NR waveforms on the left of the matching window hint at the quality of hybridization performed here. Table I lists details of all hybrids considered in the current work.  \n<!-- image -->  \nt/M  \nTABLE I. Set of time-domain hybrids constructed by matching NR simulations from the SXS Catalog and state-of-the-art PN prescriptions for BBHs on eccentric orbits are listed. SXS simulation IDs are retained for identification with NR simulations used in constructing the hybrids. Each hybrid assumes a fixed value for initial eccentricity ( e 0 ), mean anomaly ( l 0 ) and the frequency dependent PN parameter ( x 0 ) obtained using the gw eccentricity package based on Ref. [68]. The PN parameter ( x ) is related to the GW frequency ( f gw ) of the dominant mode as x =( πMf gw ) 2 / 3 with M representing binary's total mass. Mass ratio ( q ) and number of orbits prior to the merger are also listed. N orb is computed by taking the phase difference between the start of the waveform and the peak of the dominant mode ( /lscript =2, | m | =2) amplitude.  \nA least-squares minimization of the integrated difference between the GW modes from the PN and NR waveforms in a time interval ( t i , t f ), in which the two approaches give similar results, is performed and can be defined as  \nδ = min t 0 ,ϕ 0 ,ψ ∫ t f t i dt ∑ /lscript,m ∣ ∣ ∣ h NR /lscriptm ( t -t 0 ) e i ( mϕ 0 + ψ ) -h PN /lscriptm ( t ) ∣ ∣ ∣ , (7)  \nwhere the minimization is performed over a time shift ( t 0 ) and the two angles ( ϕ 0 , ψ ) as discussed above. The hybrid waveforms are then constructed by combining the NR data with the 'best matched' PN waveform in the following way:  \nh hyb /lscriptm ( t ) ≡ τ ( t ) h NR /lscriptm ( t -t ' 0 ) e i ( mϕ ' 0 + ψ ' ) +(1 -τ ( t )) h PN /lscriptm ( t ) , (8)  \nwhere ( t ' 0 , ϕ ' 0 , ψ ' ) are the values of ( t 0 , ϕ 0 , ψ ) that minimize the integral of Eq. (7). In the above equation, τ ( t )  \nis a weighting function defined by  \nτ ( t ) ≡    0 if t < t i t -t i t f -t i if t i ≤ t < t f 1 if t f ≤ t . (9)  \nThe hybrids corresponding to a representative NR simulation (SXS:BBH:1364) for all relevant modes are shown in Fig. 3. The two waveforms are aligned at merger and the shaded grey region t ∈ ( -2000 M, -1000 M ) highlights the matching window where hybridization was performed. Overlapping hybrid and NR waveforms outside (on the left of) the matching window hint at the quality of hybridization performed here.  \nWe reconstruct IMR hybrids corresponding to all 20 eccentric NR simulations listed in Ref. [48] as was done for Paper I. These are listed in Table I and the SXS simulation IDs have been retained to identify the hybrids with the corresponding NR simulation. Each simulation starts with a specific initial eccentricity ( e 0 ), mean anomaly ( l 0 ) and frequency ( x 0 ) obtained following the procedure discussed in Sec. II A.", 'III. THE WAVEFORM MODEL': 'Until now we focused on constructing a (PN-NR) hybrid model which could be used as a target for building a fully analytical IMR model for eccentric binary black hole mergers. 8 While these hybrids could be used to construct the model following the procedure adopted in Paper I, as was done there, we construct an independent set of hybrids using PN model EccentricTD [41] (see Fig. 4). This is primarily done to minimize the difference between the target hybrids and the model (being constructed) that also uses the waveform EccentricTD [41] for the inspiral part. 9 Note that, these new hybrids only include the dominant mode ( /lscript = 2 , | m | = 2) since it is the dominant mode which we wish to model first. Note that, these are also the hybrids used in validating the model constructed in this section.  \nAs in Paper I, here too we obtain a fully-analytical dominant ( /lscript = 2 , | m | = 2) mode model by matching an eccentric PN inspiral [41] with a quasi-circular prescription for the merger-ringdown phase [72]. Here too we stick to the procedures adopted in Paper I which involves identifying attachment times for both amplitude and frequency data together with an overall shift and the  \nFIG. 4. Same as Fig. 3 except, the PN part of the hybrid is purely based on EccentricTD model of Ref. [41]. These are also the hybrids that are used in constructing and validating the model in Sec. III.  \n<!-- image -->  \nFIG. 5. The (numerical) amplitude and frequency model produced by combining an eccentric inspiral with a quasi-circular merger-ringdown waveform for a total mass of 30 M /circledot is plotted against the hybrid used in calibrating the model. The eccentric inspiral ( EccentricTD [41]) and the (quasi-circular) merger-ringdown ( SEOBNRv5 [72]) models are also plotted for comparison. The amplitude (frequency) model transition smoothly from the inspiral to the merger-ringdown stage inside the shaded region(s), at t A match ( t ω match ) values, maximizing the overlap between the model and the target hybrid.  \n<!-- image -->  \noutput is a coherent IMR model suitable for generating desired signals. Exact details concerning this model are outlined in Secs. III A 1-III A 3. Note however, that one must map the data for attachment times and time shifts to a set of physical parameters of the binary such as mass ratio and other relevant parameters. For M ∼ 25 M /circledot , the highest x 0 value of Table I corresponds to a frequency of f ∼ 20Hz (low frequency cutoff for advanced LIGO design [16]). This motivates us to work with a conservative choice of 30 M /circledot system when constructing the model. We find only 13 out of 20 cases reasonably agree (with overlaps > 97%) for a 30 M /circledot system with the hybrids and thus the final fits are obtained only using these 13 performing cases. Section III B discusses the details of construction of the analytical model.', '1. Time-shift': 'The process of generating a numerical model is similar to the procedure adopted in Paper I and we reproduce it here for clarity and completeness. As described in Sec. II B, hybridization involves a minimization over a time shift, so when producing the amplitude model, we  \nhave to first perform a time shift of the inspiral waveform relative to the circular IMR waveform, because the time to merger is not known. This is done by first setting the merger time for the circular IMR waveform to zero and then time sliding the eccentric inspiral about the merger. We start by making a trial choice of t shift and then generate an amplitude and a phase model by the methods described in Secs. III A 2 and III A 3, respectively.', '2. Amplitude model': "As can be seen in Fig. 1, the waveforms tend to circularize near merger. 10 Hence, in order to model this effect, we can join the eccentric inspiral to the circular IMR by suitable choice of an appropriate time t A match . The amplitude model is obtained by joining the eccentric inspiral with the circular IMR using a transition function over a time window of 500 M which ends at t A match . Given a target hybrid, and a trial choice of t shift , we start with a trial choice of t A match roughly 500 M before the merger and produce the amplitude model as given below,  \nFIG. 6. Numerical fits for t A match and t ω match as well as for t shift are mapped into the physical parameter space for eccentric systems characterised by the binary's eccentricity, mean anomaly at a reference frequency and the mass ratio parameter q or η depending upon the model best-fitting the data. Circles represent the numerical data points while crosses represent the value returned by the analytical best-fit model. The best-fit model(s) predict the numerical estimates for attachment times and the time-shift within ± 1 M .  \n<!-- image -->  \nA model 22 ( t ) ≡ τ a ( t ) A IMR 22 ( t ) +(1 -τ a ( t )) A inspiral 22 ( t ) , (10)  \nwhere τ a ( t ) is defined as  \nτ a ( t ) ≡    0 if t < t i t -t i t f -t i if t i ≤ t < t f 1 if t f ≤ t . (11)  \nWe set t i = t A match -500 M and t f = t A match as the bounds of the time interval over which the two waveforms are joined. Figure 5 demonstrates the process. The grey region is the time interval ending at t A match where the inspiral and circular IMR is joined.  \nAfter the amplitude model is obtained for a particular choice of trial t shift and t A match , we combine it with the target hybrid phase to obtain the polarizations and then calculate the match with the target hybrid. 11 We then change the trial choice of t A match by 5 M , bringing it closer to the merger, and repeat the process of producing the amplitude model, and calculating the match. This variation of t A match is done until roughly 30 M before merger. We thus obtain a set of match values for varying t A match but for a single trial t shift and pick the one that has the highest value of match. We repeat the exercise for other choices of t shift (trial t shift varies between -400 M and 400 M in steps of 5 M ) and find out the corresponding t A match with the highest value of match. Thus, we obtain a set of t shift and t A match pairs with a match value for each pair. From this set, the pair with the highest value of match is chosen as the  \nnumerical estimate for t shift and t A match for a particular target hybrid. We obtain numerical estimates using the same process for all 20 target hybrids.", '3. Frequency model': 'For the frequency model, we follow a similar procedure as described in Sec. III A 2 with the only difference being the duration of the time interval where the inspiral frequency is joined with the circular IMR frequency. The value of t shift is fixed to the one that was obtained while producing the amplitude model. Similar to the amplitude model procedure, we determine an appropriate t ω match for joining the inspiral frequency with the circular IMR frequency. However, the time interval where the two are joined, starts at t ω match and ends at a time close to 30 M before merger. 12 Just like the amplitude model, we start with the choice of a trial value of frequency t ω match roughly 6000 M before merger and obtain the frequency model as given below,  \nω model 22 ( t ) ≡ τ a ( t ) ω IMR 22 ( t ) +(1 -τ a ( t )) ω inspiral 22 ( t ) , (12)  \nwhere τ a ( t ) is as defined in Eq. (11) with the difference being t i = t ω match and t f ≲ -30 M . Figure 5 demonstrates the process.  \nOnce the frequency model is obtained for the choice of trial t ω match , we calculate the phase by integrating the frequency model. This is then combined with the amplitude model obtained for the same target hybrid (generated using the numerical estimate of t shift and t A match already obtained) to produce the polarizations and a match with the target hybrid is calculated. We  \nFIG. 7. Top: Dominant mode ( /lscript =2, | m | =2) amplitude and frequency model, obtained by stitching an eccentric inspiral ( EccentricTD [41]) with a quasi-circular merger-ringdown model ( SEOBNRv5 [72]), at times predicted by the best-fit values of the (analytical) model are plotted against three representative target hybrids used in training the model. Bottom: one of the polarizations, obtained by combining the amplitude and the frequency model shown in the top panel for the q = 2 case, is shown as a visual proof of the quality of the model being presented. The transition from inspiral to merger-ringdown is shown (for visual clarity) at t A match given by Eq. (14) and has no impact on the plotted data.  \n<!-- image -->  \nthen change the trial choice of t ω match by 1 M , bringing it closer to the merger and repeat the process of producing the frequency model, and calculating the match. 13 Once again, we do this variation until roughly 30 M before merger to obtain a set of match values for varying t ω match and pick the one that has the highest value of match. The corresponding value of frequency t ω match is the numerical estimate for a particular target hybrid. We obtain numerical estimates for all 20 target hybrids using the same process.', 'B. Analytical model': "We have described the procedure of producing (numerical) time-domain model fits for the dominant mode model, where we used a set of 20 eccentric hybrids as targets to calibrate our model. For each hybrid, we obtained a numerical estimate for t shift , t A match , and t ω match . In order to be able to generate waveforms for an arbitrary configuration these numerical fits need to be mapped into the physical parameter space for eccentric systems characterised by binary's eccentricity, mean anomaly at a reference frequency and the mass ratio parameter. In this section, we determine a functional form by performing  \nanalytical fits to these numerical estimates. For analytical fits, we consider only those simulations (13 of the 20) for which the match between numerical model and the corresponding eccentric hybrid is greater than 97% and collectively refer to them as training set and the remaining (7) simulations are categorized as testing set although only two of these (HYB:SXS:BBH:1357, 1371) can really be used to test the model as other simulations have initial eccentricities significantly larger than any of the training set hybrids and thus outside the calibration range for the model (See for instance, Table I). For this reason, when validating the model against hybrids we only include these two hybrids from testing set. The fitted functions obtained are of the form as mentioned below.  \nt shift ( q, e, l ) = ∑ α,β,γ,δ A αβγδ e α q β cos( γ l + δ e l + a αβγδ ) , (13)  \nfor time shift, where A αβγδ = a αβγδ = 0 for α + β > 3 and/or α > 2 and/or γ + δ > 1, and A α 001 = A 0 β 10 = A 0 β 01 = A 00 γδ = a αβ 00 = 0,  \nt A match ( η, e, l ) = ∑ α,β,γ,δ B αβγδ e α η β cos( γ l + δ e l + b αβγδ ) ,  \n(14)  \nfor amplitude, where B αβγδ = b αβγδ = 0 for α + β > 3 and/or α > 2 and/or γ + δ > 1, and B 0 β 10 = B 0 β 01 = B 00 γδ = b αβ 00 = 0, and  \nFigure 6 shows comparisons between the numerically obtained values for t shift , t A match , and t ω match , with  \n<!-- image -->  \n/circledot  \n<!-- image -->  \n/circledot  \nFIG. 8. Mismatch as a function of total mass for all 13 training set hybrids (thin lines) and for the two testing set hybrids ( HYB:SXS:BBH:1357 and HYB:SXS:BBH:1371 ; thick lines) when compared with dominant mode ( /lscript =2, | m | =2) quasicircular templates ( SEOBNRv5 [72]; left panel) and eccentric templates (model constructed in Sec. III; right panel). The two horizontal lines indicate a match of 96.5% (or mismatch of 3.5%) and 99% (or 1% mismatch), respectively, and the color bar displays initial eccentricity ( e 0 ) value at an initial GW frequency ( x 0 ) for the respective target hybrid (see Table I).  \n<!-- image -->  \n/circledot  \n<!-- image -->  \n/circledot  \nFIG. 9. Same as Fig. 8 except mismatches are computed against an independent family of waveforms ( TEOBResumSDali [58]). Additionally, the parameter space explored by randomly sampling values for initial eccentricity ( e 0 ) and mean anomaly ( l 0 ) as well as of the mass ratio ( q ), in the range 0 ≲ e 0 ≲ 0 . 3, -π ≤ l 0 ≤ π , and 1 ≲ q ≲ 3, respectively.  \nt ω match ( η, e, l ) = ∑ α,β,γ,δ C αβγδ e α η β cos( γ l + δ e l + c αβγδ ) (15)  \nfor frequency, where C αβγδ = c αβγδ = 0 for α + β > 3 and/or α > 2 and/or γ + δ > 1, and C 0 β 10 = C 0 β 01 = C 00 γδ = c αβ 00 = 0. The values for the coefficients A αβγδ , B αβγδ , C αβγδ , a αβγδ , b αβγδ , and c αβγδ obtained by performing a fit to the numerical values are tabulated in Tables II-IV.  \n,  \nthe values predicted by our analytical fits. The predictions are within ± 1 M for numerical estimates for t shift , t A match and t ω match . We show amplitude and frequency comparison between the target hybrids and our models for three cases along with the full waveform for the q = 2 case in Fig. 7.", 'C. Validation': 'The performance of the model can be assessed from the plots against the new hybrids presented in Fig. 7, as well as from the mismatch (1 - M ) plots displayed in  \nFig. 8. Note again, the new hybrids are purely based on the PN prescription, EccentricTD [41] - the same model that constitutes the inspiral part of the model constructed here. Certainly, the dominant mode model outperforms the quasi-circular templates. On top of that, unlike Paper I, the model provides ≥ 96 . 5% match against all the training set hybrids for almost the entire range of parameters considered here. Furthermore, we also try to test our model against the testing set hybrids using Nelder-Mead down-hill simplex minimization algorithm of Scipy [77] over the three initial parameters set ( e 0 , l 0 , f 0 ) for each testing set hybrid. The mismatch against testing hybrids are plotted as thick lines in both panels of Fig. 8. As can be seen there, for q = 1 (testing) hybrid, match is ≥ 96 . 5% for the entire range of total mass, while it degrades a little for q = 3 (testing) hybrid for the low mass range ( M ≲ 60 M /circledot ) which is consistent with the trends observed for mismatch against q = 3 training set simulations.  \nNote that the target hybrids used in these mismatch computations include only the ( /lscript = 2 , | m | = 2) mode so as to assess the actual performance of the dominant mode model. Given the quality of analytical fits it is not surprising that the analytical model performs well against the set of hybrids used in training the model. Moreover, we retain the cases agreeing closely with the hybrids apart from two testing hybrids. (Only 13 out of 20 simulations were used in finding the analytical model.) Keeping this in mind we also try to test our model against an independent waveform family TEOBResumS-Dali [58]. For this comparison, we choose to sample a parameter space that is not identical to the training set hybrids. We choose to randomly sample the values of a reference eccentricity ( e 0 ), mass ratio ( q ) and reference mean anomaly ( l 0 ) in the range 0 ≲ e 0 ≲ 0 . 3, 1 ≲ q ≲ 3 and -π ≤ l 0 ≤ π , respectively. Note that the range for initial orbital eccentricity ( e 0 ) is slightly different from the ones spanned by the hybrids. While the upper value is chosen conservatively to take e 0 = 0 . 3 to reduce any systematic differences between the model and the waveform TEOBResumS-Dali [58] at large eccentricity values, near circular cases are also included to see if the model gradually produces the circular limit despite being trained on purely eccentric target models.  \nNext, the template (dominant eccentric model) is optimized against the target ( TEOBResumS-Dali [58]) using the same minimization algorithm [77] used in validating the model against testing set hybrids. The mismatch plot obtained is shown in Fig. 9. Additionally, for comparison, mismatches of TEOBResumS-Dali [58] with quasi-circular SEOBNRv5 [72] templates are displayed in the left panel. Clearly, our model seems to do better compared to the circular templates at the low mass end where the overlaps are ≳ 96.5% for nearly the entire range of parameters considered in the comparison.  \nThe mismatches are comparable for heavier systems as expected.', 'A. Eccentric model based on TaylorT2 phase and PN corrected amplitudes': 'In this section, we discuss the possibility of finding a suitable alternative to the model constructed in the previous section. As mentioned earlier in the Sec. I B, we propose to replace the PN model used in constructing the model in the previous section to include amplitude terms with higher PN accuracy in the model. We employ 3PN accurate expressions for the dominant mode amplitude of Refs. [39, 42, 43] and a 3PN accurate phasing (based on TaylorT2 approximant) [40] to construct this alternate model. Note that the inspiral part of the model presented in Sec. III was entirely based on the work of Ref. [41]. The EccentricTD approximant is 2PN accurate in phase and only Newtonian accurate in amplitude for the eccentricity related effects although is based on a superior (compared to TaylorT2 ) approximant namely TaylorT4 and should also work better for larger eccentricities as it includes corrections to 6th power in eccentricity while the TaylorT2 phase we use involves only leading order corrections of eccentricity although is 3PN accurate [40]. To test the performance of this model, we again compare this model with TEOBResumS-Dali [58] by computing overlaps on the same set of parameters used in validating the dominant mode model in Sec. III C. The results are shown in Fig. 10. Overlaps between the model and the target waveforms are better than ∼ 96 . 5% for almost the entire range considered here, making it a suitable alternative to the presented EccentricTD model.', 'B. Inclusion of higher order modes': 'In this section, we extend the dominant mode model of Sec. IV A to obtain a higher mode model by including ( /lscript , | m | )=(2, 2), (2, 1), (3, 3), (3, 2), (4, 4), (4, 3) and (5, 5) modes. These are precisely the modes which are included in SEOBNRv5HM (Ref. [72]) which we use for the merger-ringdown part and also included in the hybrids constructed in Sec. II B. The (eccentric) inspiral and (quasi-circular) merger-ringdown models are again attached using the analytical expressions for t shift , t A match and t ω match , obtained for the dominant mode model in Sec. III. (The higher mode model of Paper I was also obtained following the same strategy albeit, the HM model there only included the /lscript = | m | modes.) The inspiral part of the model for each nonquadrupole mode is obtained by combining the mode amplitudes obtained in Refs. [39, 42, 43] and orbital fre-  \nFigure 12 compares our nonspinning model with the spinning version of TEOBResumS-Dali [58]. Thick curves represent the mismatch when the inclination angle ( ι ) is set to zero for which only | m | = 2 modes survive (in our case (2, 2) and (3, 2)). Thin lines on the other hand, show mismatches for the inclination angle of 30 · and thus display mismatches with all /lscript = | m | and /lscript -1 = | m | modes included in our higher mode model. For the equal mass case ( q = 1), we find excellent agreements with target models (matches ≳ 99%). For q = 2 case, the matches are still ≳ 96 . 5% for χ eff ≤ 0 . 1 for all systems with M ≳ 60 M /circledot . This clearly shows that despite being nonspinning in nature, the model could be used to analyse mildly spinning events observed routinely by current generation detectors such as LIGO and Virgo [59-61]. In other words, systematics due to neglect of spin effects may be ignored.  \n<!-- image -->  \n/circledot  \nFIG. 10. Same as in Fig. 9 except target waveforms ( TEOBResumS-Dali [58]) are compared with the alternate dominant mode model presented in Sec. IV A.  \n<!-- image -->  \n/circledot  \nFIG. 11. Similar to Fig. 10 except, that both the target and the model now include all modes included in our hybrids and the inclination angle of the binary is chosen to be 30 · .  \nquency (multiplied by an appropriate factor involving mode number; see Eq. (4)). Figure 11 compares our higher mode model against the HM version of the waveform, TEOBResumS-Dali [58]. While all modes up to /lscript = 8 are included in TEOBResumS-Dali , we choose to include the same set of modes in the target and the model waveform to control the systematics. The orbital inclination angle is chosen to be 30 · . We find that the higher mode model recovers the target waveforms with accuracy better than 96 . 5% for nearly the entire range of parameter values considered here and thus may be used for analysing signals containing non-quadrupole modes and inclination angles ≤ 30 · .  \n<!-- image -->  \n/circledot  \nFIG. 12. Recovery of mildly spinning target waveforms ( TEOBResumS-Dali [58]) with the alternate nonspinning model presented in Sec. IV A is displayed for two different mass ratio cases ( q = 1, 2) and a set of effective spin values (with χ eff ≤ 0 . 1).', 'C. Recovery of aligned spin binaries': "Although our alternate model (including higher modes) is nonspinning, we also test its performance against a spinning model for a few mildly spinning cases. We assume spin-precession to be absent and thus binary's spin is described solely by the effective spin parameter, χ eff . When expressed in terms of dimensionless spin components, χ i = ( /vector S i · ˆ L ) /m 2 i , it reads [78-80]  \nχ eff = χ 1 m 1 + χ 2 m 2 m 1 + m 2 . (16)  \nHere, m i refers to the mass of the binary component having spin angular momentum /vector S i , and ˆ L denotes the unit vector along the direction of the orbital angular momentum of the binary.", 'V. DISCUSSION AND CONCLUSION': 'Our current work is a follow up of our earlier work of Ref. [1] (referred as Paper I through preceding sections) where we developed a fully analytical dominant mode model for nonspinning binary black holes on elliptical orbits. The model was obtained by stitching an eccentric inspiral EccentricTD [41] with a quasi-circular merger ringdown prescription SEOBNRv4 [81]. Here, we revisit the construction of the model presented in Paper I following a new definition of orbital eccentricity presented in Refs. [68, 69].  \nWe started by comparing 20 distinct NR simulations from SXS collaboration presented in Ref. [48] with a PN prescription obtained by combing the results of Refs. [41, 43, 82] for a selected set of modes chosen based on their relative significance compared to the dominant modes. Figure 1 shows a comparison of the amplitude and frequency data for one particular dataset. The PN model is evolved using an initial set of binary orbital parameters such as orbital eccentricity consistent with the definition of Refs. [68, 69]. For completeness, we also show dominant mode comparisons for few other simulations in Fig. 2. We find that the PN and NR prescriptions have maximum overlap in a time-window of ( -2000 M, -1000 M ) shown by shaded regions of Fig. 1-2. In fact, we observe a minimum match of ∼ 96% for all 20 simulations sets against PN models in this window, allowing hybridization of the two in this window. Table I lists all hybrids constructed here. Figure 3 displays one of these hybrids and as shown there the hybrid are constructed for all modes included in the comparison presented in Fig. 1. The hybridization procedure is same as in Paper I and is reproduced in Sec. II B. However, the model constructed in Sec. III is trained, following Paper I, with a new set of dominant mode hybrids, which uses the same PN prescription used in our dominant mode model. Details of model construction is discussed in Sec. III A-III B and validation is performed in Sec. III C.  \nWe validate the model against the hybrids used in training the model (training set in Table I) as well as two simulations (HYB:SXS:BBH:1357, 1371) from the testing set. 14 While consistency with training set was expected (but also confirms the quality of fits), better than 96 . 5% match with the two testing set simulations for most of the parameter space confirm the reliability of the model (see thick curves in the right panel of Fig. 8). Clearly the model also outperforms quasi-circular model shown in the left panel of Fig. 8. We also validate the model against an independent waveform family  \nTEOBResumS-Dali [58] for values randomly sampled for initial values of eccentricity ( e 0 ), mass ratio ( q ), and mean anomaly ( l 0 ) in the range 0 ≲ e 0 ≲ 0 . 3, 1 ≲ q ≲ 3, and -π ≤ l 0 ≤ π , respectively. Note the difference in eccentricity scale explored here compared to the one in Fig. 8 (See, Sec. III C for a discussion). Figure 9 shows this comparison and the model recovers the target family with overlaps better than 96 . 5% for mass ratios ≤ 1 . 25. For q > 1 . 25, mismatches are slightly poorer for cases with M ≤ 50 M /circledot .  \nWe also provide an alternate model by combining PN inputs from Ref. [40] and Refs. [39, 42, 43] for phase and amplitude part of the model respectively. The performance of this model against the waveform used in validating the model in Fig. 9 has been shown in Fig. 10. As can be seen, the performance of this alternate model is at the same level as of the model presented in Sec. III. Further, following Paper I, we also extended this alternate model to include few leading /lscript = | m | and /lscript -1= | m | modes (up to /lscript = 5) and compare it against waveforms of TEOBResumS-Dali [58] in Fig. 11 keeping the same set of modes in both target and template. As we can see the model at least reliably reproduces the target model for inclination angles ≤ 30 · .  \nFinally, while our model(s) assumes spinless binary constituents, we also tried testing its suitability for analysing signals with small spin magnitudes in the absence of spin precession. Figure 12 compares our nonspinning model (including higher modes) presented in Sec. IV A-IV B against the HM model of TEOBResumS-Dali [58] (with spins switched on). As shown there, our model seems to be able to extract the equal mass spinning target waveforms with accuracy better than 99% as shown by the thick-solid lines. On the other hand, our model recover target model with matches larger than than 96.5% for the q = 2 case (thin-dashed lines). For higher spin magnitudes as well as higher mass ratios, mismatches are larger than 3.5% in the low mass range, at least for unequal mass cases.  \nNote that, both the alternate model as well as the HM model were not obtained by calibrating against hybrids but rather we simply used the prescriptions presented in context of dominant mode model in Sec. III B and thus could be improved, however we restrict ourselves here to a proof of principle demonstration that alternate prescriptions for dominant mode model as well as simple extensions like the one proposed here could be easily achieved and perform reliably and leave such updates for a future work. Apart from being able to construct an improved model compared to the one presented in Paper I, in the current work we also address a few concerns with the model there. First, it was found that for nearly 10-15% cases, the analytical fits for times for merger-ringdown attachment produced nonphysical values (beyond merger at t > 0). Since the NR simulations  \nused in this work essentially circularize by 30 M before the merger [48], we impose this condition for those nonphysical scenarios. This removes the discrepancy regarding the attachment times and ensures the validity of the model in the entire parameter space explored. Secondly, the model presented in Sec. IV A is significantly faster than the dominant mode model of Paper I. This is likely due to a speedup with merger-ringdown model ( SEOBNRv5 [72] instead of SEOBNRv4 [81]). The waveform generation rate of TaylorT2 model is ∼ 1 . 5 times higher than the model based on EccentricTD . Moreover, when compared with the waveform TEOBResumS-Dali [58], these waveforms nearly have a ∼ 2 times speedup and thus likely to be useful for parameter estimation studies.', 'VI. ACKNOWLEDGMENTS': "We thank Kaushik Paul and Prayush Kumar for sharing invaluable insights in validating our model. We thank the authors of Ref. [58] for making the implementation of the TEOBResumS-Dali available for public use [83], and Divyajyoti for helping us with technical details of the implementation and waveform generation. We are thankful to the SXS Collaboration for making a public catalog of numerical relativity waveforms. P.M. thanks the members of the gravitational wave group at the Department of Physics, IIT Madras for organizing the weekly journal club sessions and the insightful discussions. T.R.C. acknowledges the support of the National Science Foundation award PHY-2207728. C.K.M. acknowledges the support of SERB's Core Research Grant No. CRG/2022/007959. This document has LIGO preprint number LIGO-P2400355. We thank Md Arif Shaikh for useful comments and suggestions on our manuscript.", 'Appendix A: Coefficients of the analytical fit': "Here we tabulate the coefficients of the analytical fits for the parameters t shift , t A match and t ω match obtained in Sec III B. The expressions are given in Eq. (13), (14) and  \n6  \n(15).  \nTABLE II. Table of coefficients for the analytical expression of t shift in Eq. (13). All other coefficients not included in the table are zero.  \nTABLE III. Table of coefficients for the analytical expression of amplitude t match in Eq. (14). All other coefficients not included in the table are zero.  \nC  \nα  \n= 0  \n-  \n1  \n.  \n87711  \n×  \n10  \n1  \n1  \n.  \n30607  \n×  \n10  \n2  \n-  \n2  \n.  \n51977  \n×  \n10  \n7  \n7  \n8  \n.  \n92879  \n×  \n10  \n-  \n5  \n.  \n21053  \n×  \n10  \n1  \n.  \n04979  \n×  \n10  \n6  \n7  \n8  \n-  \n5  \n.  \n75326  \n×  \n10  \n6  \n0  \n0  \nTABLE IV. Table of coefficients for the analytical expression of frequency t match in Eq. (15). All other coefficients not included in the table are zero.  \n- [4] R. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. X 11 , 021053 (2021), arXiv:2010.14527 [gr-qc].\n- [5] R. Abbott et al. (LIGO Scientific, VIRGO), Phys. Rev. D 109 , 022001 (2021), arXiv:2108.01045 [gr-qc].  \nαβ  \n00  \nβ  \n= 0  \n1  \n2  \n- [8] B. P. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J. Lett. 818 , L22 (2016), arXiv:1602.03846 [astro-ph.HE].\n- [9] R. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J. Lett. 900 , L13 (2020), arXiv:2009.01190 [astro-ph.HE].\n- [10] A. G. Abac et al. (LIGO Scientific, VIRGO, KAGRA), (2023), arXiv:2308.03822 [astro-ph.HE].\n- [11] B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 116 , 221101 (2016), [Erratum: Phys.Rev.Lett. 121, 129902 (2018)], arXiv:1602.03841 [gr-qc].\n- [12] M. Zevin, J. Samsing, C. Rodriguez, C.-J. Haster, and E. Ramirez-Ruiz, Astrophys. J. 871 , 91 (2019), arXiv:1810.00901 [astro-ph.HE].\n- [13] A. Ramos-Buades, S. Tiwari, M. Haney, and S. Husa, Phys. Rev. D 102 , 043005 (2020), arXiv:2005.14016 [grqc].\n- [14] G. Fumagalli, I. Romero-Shaw, D. Gerosa, V. De Renzis, K. Kritos, and A. Olejak, (2024), arXiv:2405.14945 [astro-ph.HE].\n- [15] P. C. Peters, Phys. Rev. 136 , B1224 (1964).\n- [16] J. Aasi et al. (LIGO Scientific), Class. Quant. Grav. 32 , 074001 (2015), arXiv:1411.4547 [gr-qc].\n- [17] F. Acernese (Virgo), J. Phys. Conf. Ser. 610 , 012014 (2015).\n- [18] B. P. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J. 883 , 149 (2019), arXiv:1907.09384 [astro-ph.HE].\n- [19] R. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 125 , 101102 (2020), arXiv:2009.01075 [gr-qc].\n- [20] R. Gamba, M. Breschi, G. Carullo, S. Albanesi, P. Rettegno, S. Bernuzzi, and A. Nagar, Nature Astron. 7 , 11 (2023), arXiv:2106.05575 [gr-qc].\n- [21] C. Kimball et al. , Astrophys. J. Lett. 915 , L35 (2021), arXiv:2011.05332 [astro-ph.HE].\n- [22] I. M. Romero-Shaw, P. D. Lasky, and E. Thrane, Astrophys. J. Lett. 921 , L31 (2021), arXiv:2108.01284 [astroph.HE].\n- [23] E. O'Shea and P. Kumar, Phys. Rev. D 108 , 104018 (2023), arXiv:2107.07981 [astro-ph.HE].\n- [24] I. M. Romero-Shaw, P. D. Lasky, and E. Thrane, Astrophys. J. 940 , 171 (2022), arXiv:2206.14695 [astroph.HE].\n- [25] H. L. Iglesias et al. , (2022), arXiv:2208.01766 [gr-qc].\n- [26] N. Gupte et al. , (2024), arXiv:2404.14286 [gr-qc].\n- [27] D. A. Brown and P. J. Zimmerman, Phys. Rev. D 81 , 024007 (2010), arXiv:0909.0066 [gr-qc].\n- [28] E. A. Huerta and D. A. Brown, Phys. Rev. D 87 , 127501 (2013), arXiv:1301.1895 [gr-qc].\n- [29] M. Favata, C. Kim, K. G. Arun, J. Kim, and H. W. Lee, Phys. Rev. D 105 , 023003 (2022), arXiv:2108.05861 [gr-qc].\n- [30] B. P. Abbott et al. (LIGO Scientific, Virgo), Class. Quant. Grav. 34 , 104002 (2017), arXiv:1611.07531 [grqc].\n- [31] Divyajyoti, S. Kumar, S. Tibrewal, I. M. Romero-Shaw, and C. K. Mishra, Phys. Rev. D 109 , 043037 (2024), arXiv:2309.16638 [gr-qc].\n- [32] D. McClelland, M. Evans, R. Schnabel, B. Lantz, I. Martin, and V. Quetschke, (2016).\n- [33] S. Dwyer, D. Sigg, S. W. Ballmer, L. Barsotti, N. Mavalvala, and M. Evans, Phys. Rev. D 91 , 082001 (2015), arXiv:1410.0612 [astro-ph.IM].\n- [34] B. P. Abbott et al. (LIGO Scientific), Class. Quant. Grav. 34 , 044001 (2017), arXiv:1607.08697 [astro-ph.IM].\n- [35] M. Punturo et al. , Class. Quant. Grav. 27 , 194002 (2010).  \n- [36] S. Hild et al. , Class. Quant. Grav. 28 , 094013 (2011), arXiv:1012.0908 [gr-qc].\n- [37] M. E. Lower, E. Thrane, P. D. Lasky, and R. Smith, Phys. Rev. D 98 , 083028 (2018), arXiv:1806.05350 [astroph.HE].\n- [38] S. Tibrewal et al. , In preparation (2022).\n- [39] C. K. Mishra, K. G. Arun, and B. R. Iyer, Phys. Rev. D 91 , 084040 (2015), arXiv:1501.07096 [gr-qc].\n- [40] B. Moore, M. Favata, K. Arun, and C. K. Mishra, Phys. Rev. D 93 , 124061 (2016), arXiv:1605.00304 [gr-qc].\n- [41] S. Tanay, M. Haney, and A. Gopakumar, Phys. Rev. D 93 , 064031 (2016), arXiv:1602.03081 [gr-qc].\n- [42] Y. Boetzel, C. K. Mishra, G. Faye, A. Gopakumar, and B. R. Iyer, Phys. Rev. D 100 , 044018 (2019), arXiv:1904.11814 [gr-qc].\n- [43] M. Ebersold, Y. Boetzel, G. Faye, C. K. Mishra, B. R. Iyer, and P. Jetzer, Phys. Rev. D 100 , 084043 (2019), arXiv:1906.06263 [gr-qc].\n- [44] C. Konigsdorffer and A. Gopakumar, Phys. Rev. D 73 , 124012 (2006), arXiv:gr-qc/0603056 [gr-qc].\n- [45] B. Moore and N. Yunes, Class. Quant. Grav. 36 , 185003 (2019), arXiv:1903.05203 [gr-qc].\n- [46] G. Carullo, S. Albanesi, A. Nagar, R. Gamba, S. Bernuzzi, T. Andrade, and J. Trenado, Phys. Rev. Lett. 132 , 101401 (2024), arXiv:2309.07228 [gr-qc].\n- [47] G. Carullo, (2024), arXiv:2406.19442 [gr-qc].\n- [48] I. Hinder, L. E. Kidder, and H. P. Pfeiffer, Phys. Rev. D 98 , 044015 (2018), arXiv:1709.02007 [gr-qc].\n- [49] E. A. Huerta et al. , Phys. Rev. D 95 , 024038 (2017), arXiv:1609.05933 [gr-qc].\n- [50] Z. Chen, E. A. Huerta, J. Adamo, R. Haas, E. O'Shea, P. Kumar, and C. Moore, Phys. Rev. D 103 , 084018 (2021), arXiv:2008.03313 [gr-qc].\n- [51] Y. Setyawati and F. Ohme, Phys. Rev. D 103 , 124011 (2021), arXiv:2101.11033 [gr-qc].\n- [52] R. Gamba, D. Chiaramello, and S. Neogi, (2024), arXiv:2404.15408 [gr-qc].\n- [53] A. Ramos-Buades, S. Husa, G. Pratten, H. Estell'es, C. Garc'ıa-Quir'os, M. Mateu-Lucena, M. Colleoni, and R. Jaume, Phys. Rev. D 101 , 083015 (2020), arXiv:1909.11011 [gr-qc].\n- [54] D. Chiaramello and A. Nagar, Phys. Rev. D 101 , 101501 (2020), arXiv:2001.11736 [gr-qc].\n- [55] A. Albertini, R. Gamba, A. Nagar, and S. Bernuzzi, Phys. Rev. D 109 , 044022 (2024), arXiv:2310.13578 [grqc].\n- [56] A. Nagar and S. Albanesi, Phys. Rev. D 106 , 064049 (2022), arXiv:2207.14002 [gr-qc].\n- [57] A. Nagar, R. Gamba, P. Rettegno, V. Fantini, and S. Bernuzzi, (2024), arXiv:2404.05288 [gr-qc].\n- [58] A. Nagar, S. Bernuzzi, D. Chiaramello, V. Fantini, R. Gamba, M. Panzeri, and P. Rettegno, (2024), arXiv:2407.04762 [gr-qc].\n- [59] B. P. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J. Lett. 882 , L24 (2019), arXiv:1811.12940 [astro-ph.HE].\n- [60] R. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J. Lett. 913 , L7 (2021), arXiv:2010.14533 [astro-ph.HE].\n- [61] R. Abbott et al. (KAGRA, VIRGO, LIGO Scientific), Phys. Rev. X 13 , 011048 (2023), arXiv:2111.03634 [astroph.HE].\n- [62] A. Rebei, E. A. Huerta, S. Wang, S. Habib, R. Haas, D. Johnson, and D. George, Phys. Rev. D 100 , 044025 (2019), arXiv:1807.09787 [gr-qc].\n- [63] A. Ramos-Buades, A. Buonanno, M. Khalil, and S. Ossokine, Phys. Rev. D 105 , 044035 (2022), arXiv:2112.06952 [gr-qc].\n- [64] A. Nagar, A. Bonino, and P. Rettegno, Phys. Rev. D 103 , 104021 (2021), arXiv:2101.08624 [gr-qc].\n- [65] T. Islam, V. Varma, J. Lodman, S. E. Field, G. Khanna, M. A. Scheel, H. P. Pfeiffer, D. Gerosa, and L. E. Kidder, Phys. Rev. D 103 , 064022 (2021), arXiv:2101.11798 [grqc].\n- [66] T. Islam, (2024), arXiv:2403.15506 [astro-ph.HE].\n- [67] Divyajyoti, P. Baxi, C. K. Mishra, and K. G. Arun, Phys. Rev. D 104 , 084080 (2021), arXiv:2103.03241 [gr-qc].\n- [68] M. A. Shaikh, V. Varma, H. P. Pfeiffer, A. RamosBuades, and M. van de Meent, Phys. Rev. D 108 , 104007 (2023), arXiv:2302.11257 [gr-qc].\n- [69] A. Ramos-Buades, M. van de Meent, H. P. Pfeiffer, H. R. Ruter, M. A. Scheel, M. Boyle, and L. E. Kidder, Phys. Rev. D 106 , 124040 (2022), arXiv:2209.03390 [gr-qc].\n- [70] K. Arun, L. Blanchet, B. R. Iyer, and M. S. Qusailah, Class. Quant. Grav. 21 , 3771 (2004), [Erratum: Class. Quant. Grav. 22, 3115 (2005)], arXiv:gr-qc/0404085 [grqc].\n- [71] V. Varma and P. Ajith, Phys. Rev. D 96 , 124024 (2017), arXiv:1612.05608 [gr-qc].  \n- [72] L. Pompili et al. , Phys. Rev. D 108 , 124035 (2023), arXiv:2303.18039 [gr-qc].\n- [73] LIGO Scientific Collaboration, 'LIGO Algorithm Library - LALSuite,' free software (GPL) (2018).\n- [74] M. Boyle et al. , Class. Quant. Grav. 36 , 195006 (2019), arXiv:1904.04831 [gr-qc].\n- [75] B. J. Owen and B. S. Sathyaprakash, Phys. Rev. D 60 , 022002 (1999), arXiv:gr-qc/9808076.\n- [76] V. Varma, S. E. Field, M. A. Scheel, J. Blackman, L. E. Kidder, and H. P. Pfeiffer, Phys. Rev. D 99 , 064045 (2019), arXiv:1812.07865 [gr-qc].\n- [77] T. E. Oliphant and S. . Contributors, Nature Methods 17 , 261 (2020).\n- [78] Divyajyoti, N. V. Krishnendu, M. Saleem, M. Colleoni, A. Vijaykumar, K. G. Arun, and C. K. Mishra, Phys. Rev. D 109 , 023016 (2024), arXiv:2311.05506 [gr-qc].\n- [79] P. Ajith et al. , Phys. Rev. Lett. 106 , 241101 (2011), arXiv:0909.2867 [gr-qc].\n- [80] L. Santamaria et al. , Phys. Rev. D 82 , 064016 (2010), arXiv:1005.3306 [gr-qc].\n- [81] A. Boh'e et al. , Phys. Rev. D 95 , 044028 (2017), arXiv:1611.03703 [gr-qc].\n- [82] Y. Boetzel, A. Susobhanan, A. Gopakumar, A. Klein, and P. Jetzer, Phys. Rev. D 96 , 044011 (2017), arXiv:1707.02088 [gr-qc].\n- [83] https://bitbucket.org/teobresums/teobresums/ branches ."}