Fast post-adiabatic waveforms in the time domain: Applications to compact binary coalescences in LIGO and Virgo

We present a computationally efficient (time-domain) multipolar waveform model for quasi-circular spin-aligned compact binary coalescences. The model combines the advantages of the numerical-relativity informed, effective-one-body (EOB) family of models with a post-adiabatic solution of the equations of motion for the inspiral part of the two-body dynamics. We benchmark this model against other state-of-the-art waveforms in terms of efficiency and accuracy. We find a speed-up of one to two orders of magnitude compared to the underlying time-domain EOB model for the total mass range $2 - 100 M_{\odot}$. More specifically, for a low total-mass system, such as a binary neutron star with equal masses of $1.4 M_{\odot}$, like GW170817, the computational speedup is around 100 times; for an event with total mass $\sim 40 M_\odot$ and mass ratio $\sim 3$, like GW190412, the speedup is by a factor of $\sim 20$, while for a binary system of comparable masses and total mass of $\sim 70 M_{\odot}$, like GW150914, it is by a factor of $\sim 10$. We demonstrate that the new model is extremely faithful to the underlying EOB model with unfaithfulness less than $0.01\%$ across the entire applicable region of parameter space. Finally, we present successful applications of this new waveform model to parameter estimation studies and tests of general relativity.


I. INTRODUCTION
Since 2015, the detections of gravitational waves (gws) have yielded a wealth of remarkable discoveries [1][2][3][4][5][6][7][8]. In the three observing runs [7,8] of the Advanced ligo [9] and Advanced Virgo [10] detectors, a total of 50 events have been observed and confirmed; among these are both binary black hole (bbh) mergers and binary neutron star (bns) mergers [11]. Particularly interesting are the discoveries of binaries GW190412 with mass ratio 3 [12] and GW190814 (which could be the first ever detected merger of a black hole and a neutron star) [13] with mass ratio 10. GW190521 is the most massive binary detected so far with a total mass of 150M [14].
Detections of compact binary mergers are expected to increase in the coming years [15,16]: during the upcoming ligo and Virgo observing runs [17], and with future ground-based detectors like the Einstein Telescope [18] and Cosmic Explorer [19], and the space-based mission lisa [20]. Extracting information from such gw detections relies on accurate and computationally efficient models of the gravitational waveforms, which are emitted during coalescence [21,22]. Firstly, the estimation of the binary parameters of a typical event (using Bayesian inference, Markov chains, or similar methods) requires on the order of several million evaluations of the waveform models [23,24]. On the other hand, the gw phase needs to be accurate to less than a cycle of the binary in order to avoid ambiguity in the estimations [25,26]. Upcoming runs, as well as future detectors, will require * Electronic address: deyan@aei.mpg.de even better waveform accuracy in order to reliably identify and analyse gw events [27]. Accurately identifying the properties of a large population of binaries will allow us to make inferences on scenarios of compact-object binary formation [28], and also carry out more stringent tests of General Relativity (gr) in the highly dynamical, strong-field regime [29]. For these reasons, work on more advanced and innovative waveform models continues for ligo, Virgo, and future gw missions.
The time-domain SEOBNR models [43,45,46] are routinely employed in data analysis for sufficiently high-mass binaries with the ligo Algorithm Library (lal) Inference codes [24], while, for generic-mass binaries, fast parameter-estimation codes are required [69]. Nevertheless, the time for generating a waveform in the lowmass regime ( 5M ) can be on the order of ∼ 100 s or even longer starting at 20 Hz. Thus, there are parts of the binary's parameter space, for which the time-domain SEOBNR models are not suitable for direct use in parameter estimations without further optimizations, like the ones discussed in the current publication, which was originally introduced in Refs. [70,71]. Alternative methods have been developed in order to afford speedy data analysis for gws. Reduced order modeling and surrogate techniques [72][73][74][75][76][77][78][79][80][81] have been successfully applied to eob waveform models [43,72,73,75,79] and to pure nr-waveforms [74,76,78,82,83]. These methods work by decomposing and interpolating the waveforms on a sparse grid in time or frequency domain, and then using interpolation or more sophisticated statistical methods to obtain the fitting parameters across the binary's parameter space under study. The resulting waveform model is then verified for accuracy against an independent testing set. However, although very successful, such models suffer from certain limitations. By construction, they are restricted to confined regions of the parameter space, and have to be developed from scratch if the underlying time-domain model is updated -for example when more physical effects are included or higher-order post-Newtonian parameters are added to make these waveforms more accurate.
Here, we develop the multipolar SEOBNRv4HM_PA waveform model for spin-aligned compact binaries moving on quasi-circular orbits. This model is a computationally cheaper version of the time-domain SEOBNRv4HM waveform model [45] and, as such, it includes higher-order harmonics (or higher modes, hm), which are important for asymmetric mass-ratio binaries, high-mass and highinclination systems [84][85][86][87][88][89]. In the SEOBNRv4HM_PA waveform model, the binary dynamics is solved using a postadiabatic (pa) approach. The latter was proposed and applied to the TEOBResumS model in Refs. [70,71,90] and used in all subsequent publications (see, e.g., Ref. [91] and references therein). It was also implemented in LALSimulation. In the pa method, the inspiral evolution (until the last few orbits before merger) is approximated by an adiabatic solution of the (ordinary differential) equations of motion, with post-adiabatic corrections added iteratively up to the order needed to achieve the desired accuracy. In this work we apply this technique to construct a fast and accurate aligned-spin dynamics based on the SEOBNRv4 model [43] and implement it in LALSimulation. The speed-up and accuracy benchmarks are supported by applications of the pa waveform model for parameter estimation studies and tests of gr.
The paper is organized as follows. Section II of this article reviews the eob dynamics in the pa approximation for arbitrary Hamiltonians. Section III presents the implementation of this method in the LALSimulation waveform model library (as approximant SEOBNRv4HM_PA), and benchmarks the model against other established eob models. Section IV presents two parameter estimation (pe) studies using the SEOBNRv4HM_PA model. In Section V the model is applied to a ringdown test of gr. Section VI concludes the article with a discussion of the significance of this work and its possible future directions.
We shall work in natural units G = 1 = c.

II. POST-ADIABATIC APPROXIMATION TO THE INSPIRAL DYNAMICS
The eob formalism provides an analytical description of the gw emission from the process of binary coalescence, including inspiral, merger, and ringdown [30,31]. The accuracy of this description can be further improved by calibrating against nr simulations.
A binary system composed of two bhs moving on a quasi-circular orbit with spins aligned or anti-aligned (henceforth, spin-aligned for short) with the orbital angular momentum is described by four parameters: the component masses m 1 and m 2 , and the (dimensionless) spins χ 1 = S 1 /m 2 1 and χ 2 = S 2 /m 2 2 . In the eob approach the (center-of-mass) two-body dynamics is mapped onto the dynamics of an effective body of mass µ = m 1 m 2 /(m 1 + m 2 ), which moves in a deformed Kerr spacetime of mass M = m 1 +m 2 , the deformation parameter being the symmetric mass ratio ν = µ/M . The conservative two-body dynamics is obtained from the eob Hamiltonian [30,39]: where H eff is the Hamiltonian that describe the motion of the effective body of mass µ and spin S * = [(m 2 /m 1 )S 1 + (m 1 /m 2 )S 2 ]/M 2 in the (deformed) Kerr spacetime of mass M with spin S = S 1 + S 2 . For aligned-spin binaries, the motion is constrained to a fixed plane. Thus, we use polar coordinates and introduce the phase-space dimensionless variables (r, ϕ, p r * , p ϕ ) related to the physical ones through the following expressions The radial momentum p r * is conjugate to the tortoise coordinate of the deformed spacetime r * [38,92]. The dissipative effects in the eob formalism are described by the radiation-reaction force [31,[35][36][37] where Ω is the angular orbital frequency, L is the orbital angular momentum, d L is the luminosity distance and h m are the gravitational modes far from the source. In this setup, the equations of motion read [38]: dr dt = dp r * dp r ∂H ∂p r * , 4b) dp r * dt = − dp r * dp r ∂H ∂r Here, t = T /M is a dimensionless time variable. The usual procedure employed in eob waveform models involves solving the eqs. (2.4) numerically, using an ordinary differential equation (ode) integrator with a suitable time step and initial conditions. This is often computationally expensive (especially for longer waveforms), and is one of the bottlenecks for efficiently generating the eob waveform. The post-adiabatic (pa) approximation [31,70,90] converts the ode equations of motion into a set of non-linear algebraic equations which need to be solved numerically, but have a lower computational cost associated with them.
The adiabatic approximation assumes that the dynamics is comprised of a sequence of circular orbits. As such, there is no radiation reaction, hence F ϕ = 0 and p r * vanishes. Hence, the leading-order orbital angular momentum p ϕ can be calculated at a given radius from Eq. (2.4c): The post-adiabatic approximation assumes that the radiation reaction F ϕ which can be used to furnish p r * through a combination of Eqs. (2.4a) and (2.4c): At the post-post-adiabatic level, one can use the newly obtained approximation for p r * to additionally correct the orbital angular momentum p ϕ , this time utilising Eqs. (2.4b) and (2.4d): This approximation procedure can be iterated further, and the procedure for obtaining the corrections to the leading-order solution (p r * , p ϕ ) = (0, j 0 (r)) can be formalised in the following way, as described in [70]. For each value of the radial coordinate r, the radiation reaction can be written as an expansion in a formal parameter Therefore, the solutions of the eob equations of motion can also be written as an expansion in powers of this fictitious parameter: The pa procedure allows for the two momenta to be calculated with arbitrary precision by adding more terms in the expansions above. The corrections at (n, n + 1)th pa order can be found by iteratively solving Eq. (2.6) for p r * and Eq. (2.7) for p ϕ . In solving these two equations, one must remember that all other variables (apart from the unknown one) must be kept at their most recent pa order. This procedure can be repeated as many times as necessary, until the desired accuracy in terms of powers of is achieved [70].
In practice, we proceed as follows. As a start, a radial grid is constructed for the part of the two-body dynamics where the pa approximation is to be applied, between two radii r max and r min . At each node in this grid, the adiabatic solution j 0 (r) is obtained through Eq. (2.5) -it serves as the leading-order uncorrected values for the orbital angular momentum (the uncorrected value for p r * is chosen as 0 everywhere on the grid). To obtain the N -th order pa approximation, the momenta p r * (r) and p ϕ (r) are computed through Eqs. (2.6) and (2.7), respectively, at each point in the grid and this part is repeated up to the chosen pa order N . Whenever radial derivatives of the corrected quantities need to be computed (for instance, dp ϕ /dr in Eq. (2.6)), this is performed numerically on the grid. Finally, the time t and the orbital phase ϕ are obtained through numerical integration The waveform is built from the pa dynamics using the same prescription as the standard eob waveform model. The waveform strain h(t) = h + (t) − ih × (t) can be decomposed into multipoles according to where d L is the distance from the detector to the source, and −2 Y m (θ, φ) are the spin-weighted spherical harmonics for s = −2. max is the highest-order multipole which is calculated. More detailed accounts of the procedure for generating eob waveforms can be found in Ref. [42,47,48]. A robust implementation of the pa dynamics for an arbitrary spin-aligned eob Hamiltonian is presented in the Sec. III.

III. IMPLEMENTATION IN LIGO ALGORITHM LIBRARY
We have implemented the post-adiabatic (pa) inspiral dynamics model described in Sec. II in LALSuite [93] and it is available through the SEOBNRv4HM_PA waveform model approximant.
When this model is used, the dynamics of the binary system, starting from the initial separation r 0 = r max until some final separation r min is approximated with the pa procedure described in Sec. II. The radius at which the pa procedure is terminated, r min , as well as the size of the grid dr, are empirically chosen to ensure that the faithfulness of the waveform is maximised while keeping the computational cost minimal. Around 10 4 waveforms were generated, covering the space of binary parameters and exploring the effects of varying these two parameters. In each case we compute the unfaithfulness and choose values for these parameters that ensure the fastest waveform generation while still being sufficiently accurate. We find that the following prescription satisfies these requirements: r min = 1.6 r isco and dr = 0.3. (3.1) Furthermore, the pa order is a free parameter of our model, with the default being 8 th order. Our studies show that lower orders cannot always achieve the desired accuracy, while higher orders incur computational cost without further improving the solution, or can be prone to numerical noise (e.g. above 12 th order). The pa approach is independent of any particular form of the Hamiltonian, and here we focus on the SEOBNRv4HM Hamiltonian [39,43] 1 . Our procedure is set up to use either analytical or numerical derivatives of the eob Hamiltonian (e.g. ∂H/∂r), this giving additional flexibility to the user. The numerical derivatives in Eqs. (2.5), (2.6), and (2.7) are computed using an 8 th -order finite difference method [94][95][96], while numerical integration is performed using a standard cube-spline quadrature algorithm [97,98].
In order to provide an implementation of the waveform model that is maximally efficient while preserving the faithfulness at each point in parameter space, a number of further changes are introduced to the algorithms for calculating the binary dynamics and for computing the waveform modes. These changes are summarized below: 1. Analytic derivatives are used both during the pa routine and the final ode integration (where the approximant SEOBNRv4_opt is used [43,99,100]). Analytic derivatives are computationally more efficient than finite-difference methods, and therefore provide valuable speedup for computing the binary dynamics.

2.
A larger integration step is used for the final part of the dynamics calculation (using the SEOBNRv4_opt model), which speeds up the ode integration significantly.
3. Following Refs. [43,99], quantities which do not vary with the mode numbers ( , m) are precomputed instead of being repeatedly generated during each iteration. This helps to remove a large portion of the computational overhead in building the waveform modes.
4. Finally, the waveform is computed over a nonuniformly-spaced time grid, which is comprised of the sparse grid for pa approximation, and the denser grid for the final part of the dynamics (plunge and merger). This speeds up the waveform generation considerably as waveform generation on an equally spaced grid is expensive. To obtain the final modes on an equally spaced grid we follow the interpolation approach described in [81].
A. Computational performance of the model The SEOBNRv4HM_PA model was benchmarked against other relevant waveform models which are available in LALSimulation: against SEOBNRv4HM [45], SEOBNRv4HM_ROM [81] and SEOBNRv4T_surrogate [79] in computational efficiency, and against SEOBNRv4HM in accuracy. Figure 1 shows the time for generating a waveform for total masses between 2 and 100M , with starting frequency of 10 Hz, and for three values of the mass ratio, q = {1, 3, 10}. The SEOBNRv4HM_PA model performs significantly faster than the SEOBNRv4HM model across all values of the total mass M . The speedup is most significant for lower total mass (∼ 60×), and drops to its lowest for high total mass (∼ 10×). Most importantly, the speedup is substantial (∼ 30×) for M ∼ 40 − 60M , where many of the events of current interest lie. For M 10M , the time to generate a waveform using SEOBNRv4HM_PA is less than 1 s.
Comparing to the frequency-domain SEOBNRv4HM_ROM model, we see that, as expected, the reduced order model is faster in almost all cases, except for very low total mass (M 10M ) where the two models are comparable in terms of computational cost (see Fig. 1).
As a further test, the SEOBNRv4HM_PA model was benchmarked against the SEOBNRv4T_surrogate model, which only includes the ( , |m|) = (2, 2) mode. Figure 2 compares the time for generation of a waveform with the SEOBNRv4HM_PA and SEOBNRv4T_surrogate models where only the quadrupole mode is computed, at different sampling rates, and with different starting frequencies, for a 1.4M + 1.4M bns with no spins.
Finally, we benchmark the time necessary to generate each higher-oder multipole of the model. Figure 3 shows the times to generate waveforms with the SEOBNRv4HM_PA model for which all ( , |m|) modes with ≤ max (see the legend) are resolved at initial frequency of f 0 = 20 Hz. As Compared to the SEOBNRv4HM model, the post-adiabatic model is between 10 and 10 2 times faster depending on the total mass (for a starting frequency of 10 Hz). The frequency-domain SEOBNRv4HM_ROM model is faster for high total mass M , but the two models have near-equal performance in the low-totalmass regime. In these tests, the compact objects have spins χ1 = 0.8 and χ2 = 0.3, and the sampling rate has been chosen so that it is large enough to resolve the (5, 5) mode for large total mass, but also to never exceed 8192 Hz. putational cost to the waveform generation.

B. Accuracy of the model
To assess the accuracy of our model with respect to SEOBNRv4HM we use the notion of unfaithfulness as outlined below.
In general, the gw signal from a non-precessing, quasicircular bbh is characterized by a total of 11 parameters: the binary companion masses m 1 and m 2 , the (dimensionless) component spins χ 1 and χ 2 , which can be aligned or anti-aligned with the orbital angular momentum, the orientation of the binary in the source frame (ι, ϕ 0 ), the sky location of the source in the detector frame (θ, φ), and the luminosity distance d L , the time of arrival t c , and the polarization angle ψ. The detector response can be written as where λ = (m 1 , m 2 , χ 1 , χ 2 , ι, ϕ 0 , d L ) are the binary parameters and the functions F {+,×} (θ, φ, ψ) are the antenna patterns (see [101,102]). This can be cast as Here κ ≡ κ(θ, φ, ψ) is the effective polarization [85], defined as where the overall amplitude function A(θ, φ) is For each mass ratio q and ( , m) mode, the initial frequency was modified so that this mode is included at f = 20 Hz. The compact objects have spins χ1 = 0.8 and χ2 = 0.3, and the sampling rate has been fixed so that it is large enough to resolve the (5, 5) mode for large total mass, but also to never exceed 8192 Hz.
We can now define the match between a gw signal h s (t) (which for us is SEOBNRv4HM) and a template waveform h t (t) (which is SEOBNRv4HM_PA) (see [103]): where the parameters denoted with s (t) refer to the signal (template) waveform. For the purposes of this calculation, we marginalize over the phase φ 0t , the effective polarization κ t , and the time of arrival t c , and use the familiar definition for the inner product between two waveforms [101,102]: where a ∼ denotes a Fourier transform, a * denotes a complex conjugate, and finally S n (f ) is the one-sided power spectral density (psd) of the detector noise.
Here, we use f low = 20 Hz, f high = 2048 Hz, and the Advanced ligo "zero-tuned high-power" design sensitivity curve [104]. The definition of the faithfulness in Eq. (3.6) depends on the signal parameters (ι s , ϕ 0s , κ s ) and therefore allows us to work with either the maximum or the averaged unfaithfulness (or mismatch) 1 − F(ι s , ϕ 0s , κ s ): Using these definitions of the unfaithfulness, we examine the accuracy of the post-adiabatic (pa) model in comparison with SEOBNRv4HM. Figure 4 shows the maximum unfaithfulness U max across the parameter space of the effective spin χ eff and the mass ratio q for 4 separate values of the total mass M . We find that in all cases that were presented, the unfaithfulness is below O(0.01%), which makes our model comparable in accuracy to SEOBNRv4HM_ROM [81]. The variations of the unfaithfulness for a fixed total mass can be explained in terms of the length of the waveform and the heuristic condition which we use to transition from pa inspiral to merger and ringdown.
The variations across the separate panels of Fig. 4 can be better expressed in the form of a histogram of the average unfaithfulness U, which is shown in Fig. 5. The distribution for each value of the total mass M peaks at a higher value of the unfaithfulness, which is consistent with the fact that lower-mass binaries undergo a longer coalescence inside the frequency band relevant to ligo. In all cases, however, the average unfaithfulness stays below O(0.01%).

IV. PARAMETER ESTIMATION STUDY
While the previous section deals with the absolute performance of the SEOBNRv4HM_PA model, it is important to demonstrate that these benchmarks translate into a faster and reliable parameter estimation (pe) analyses. For this purpose, two pe studies were performed with the SEOBNRv4HM_PA model, as well as some of the established waveform models (SEOBNRv4HM or SEOBNRv4HM_ROM). The first one is an injection study involving a synthetic signal. The data is then analysed using SEOBNRv4HM_PA and SEOBNRv4HM models. The second pe study was performed on data from the event GW190412. It was analysed using the SEOBNRv4HM_PA and the SEOBNRv4HM_ROM models.
In both cases, the parameter estimation was done using the Markov-chain Monte Carlo code LALinferenceMCMC [93]. The results of these studies are presented below in the following subsections. We also note that pe studies using the TEOBResumS waveform model with the pa approximation were done in Refs [105][106][107][108][109].

A. Injection-based study
The choice of parameters for this pe study was made so that it could emphasise the speedup which the SEOBNRv4HM_PA model offers, but also to manifest higherorder modes in the signal. The injected signal has total mass M = 60M with mass ratio 1/q ≈ 0.4, aligned spins (χ 1 , χ 2 ) = (0.8, 0.3), at a distance of approximately 1364 Mpc (corresponding to signal-to-noise ratio (snr) of ∼ 22) and at an inclination of π/3. The injection was made using the SEOBNRv4HM model, while the subsequent analysis was performed using both the SEOBNRv4HM and the SEOBNRv4HM_PA models in order to judge the performance of the post-adiabatic (pa) model.
The results of the parameter estimation analysis can be summarized in the series of marginalized 2-dimensional posterior plots in Fig. 6. Figure 6a shows the plot for the component source-frame masses m s 1 and m s 2 ; Fig. 6b shows the plot for the mass ratio q and the effective spin χ eff . Finally, Fig. 6c shows the posterior for the luminosity distance d L (cf. Eq. (2.11)) and the inclination angle ι. The plots demonstrate that the distributions of the posterior samples are in very good agreement between the two models, which demonstrates the reliability of SEOBNRv4HM_PA in pe studies. The Jensen-Shannon (js) divergence between the 1-dimensional posteriors which are shown in Fig. 6 is always at O 10 −3 or below, which is consistent with the effects of stochastic sampling on the recovered quantities [110].

B. GW190412 parameter estimation
It is more interesting to demonstrate the viability of the SEOBNRv4HM_PA model on a real gw-event. The GW190412 from the LIGO O3a catalogue is a good candidate for such a study due to its parameters -total mass of 40M (m 1 ≈ 30M and m 2 ≈ 8M , asymmetric mass ratio of about 3, and dimensionless spin of the massive companion between 0.22 and 0.6 [12,[111][112][113]. The total mass and mass ratio suggest that waveforms for this event would be computationally much more efficient using the SEOBNRv4HM_PA model -which means it is particularly suitable to study this event. The analyses were performed with both the SEOBNRv4HM_PA and SEOBNRv4HM_ROM models, since any analysis with the SEOBNRv4HM model would have been impractical in terms of the time required to complete it.
The 2-dimensional marginalized posterior plots are shown in Fig. 7. Figure 7a shows the posterior for the component source-frame masses m s 1 and m s 2 of the binary; Fig. 7b shows the plot for the mass ratio q and the effective spin χ eff ; and finally, Fig. 7c shows the plot for the luminosity distance d L and the inclination angle ι. The plots demonstrate the incredibly good agreement between the SEOBNRv4HM_ROM and the SEOBNRv4HM_PA models, which is further confirmed by the js divergence between the samples, which is below the O 10 −3 level.
In conclusion for this section, the SEOBNRv4HM_PA model may be reliably used for parameter estimation analyses in place of models like SEOBNRv4HM, with no evident caveats which could hinder the results of such  analyses. Furthermore, the average speedup in the generation of samples for the mcmc chains was found to be around 1 order of magnitude. The reason that this is less compared to the speedup demonstrated in Fig. 1 is the fact that for the pe studies, the waveforms are generated using starting frequency f 0 = 20 Hz.

V. TESTS OF GENERAL RELATIVITY
In previous sections we have seen that SEOBNRv4HM_PA is a robust, accurate and fast alternative to SEOBNRv4HM and can therefore be used as a drop-in replacement. An interesting further application of this is to use the SEOBNRv4HM_PA model for tests of General Relativity, where the SEOBNRv4HM model was previously used as a baseline gr model. In particular, we consider a parametrized black hole (bh) ringdown test that measures the deviations of quasi-normal mode emission from predictions of gr. We summarize the test briefly here, see [29,114] for more details. In gr, the no-hair conjecture predicts that the physical properties of a (uncharged) bh are completely determined by its mass and spin. Consequently, the quasi-normal modes (qnm) that describe the gravitational waves emitted by a perturbed bh are also uniquely determined by its mass and spin. Thus one can check the validity of gr by measuring or constraining any deviations in the complex qnm frequencies. The pSEOBNR ringdown analysis uses a parameterised version of a full inspiral-merger-ringdown gw signal model to measure and constrain the (complex) qnm frequencies. Consequently, unlike other ringdown stud-ies restricted to the post-merger signal [115][116][117][118], the pSEOBNR analysis makes use of the entire signal power and does not suffer from the ambiguity of a ringdown start-time definition. In pSEOBNR, one starts with the SEOBNRv4HM model, but then introduces deviations δf lm0 and δτ lm0 (which are treated as free parameters) to the qnm frequency and damping time, so that 2) and the ringdown signal is different from the gr prediction. Here, f gr lm0 and τ gr lm0 are computed from final mass and spin as predicted by nr fitting formulae. The goal of the test is to infer the values of the (δf lm0 , δτ lm0 ) by doing full parameter estimation.
Performing full parameter estimation taking into account deviations from gr is a challenging problem, since it involves sampling a higher-dimensional space (two new parameters for every ( , m) mode considered), which means any improvements to waveform speed can lead to even bigger impact on the overall runtime. To check the speed and accuracy of using SEOBNRv4HM_PA for this application, we compare it to the results obtained with SEOBNRv4HM for the first gravitational wave event ever detected, GW150914. We choose this event as it has significant signal-to-noise ratio (snr) in the mergerringdown regime. For a high (quasi) equal-mass binary like GW150914, contributions from higher multipoles of the gw signal are expected to be negligible. Hence, we restrict our analysis to just the least-damped dominant qnm, i.e., (f 220 , τ 220 ), keeping the other qnms fixed at their nominal gr values. In Fig. 8  the posterior distributions for the fractional deviations in the damping time and frequency as recovered by using SEOBNRv4HM and SEOBNRv4HM_PA. It is evident that the posteriors are in extremely good agreement with the results showing consistency with gr. More quantitatively, the Jensen-Shannon (js) divergence between the 1-dimensional posteriors is O(10 −3 ) for δf 220 and δτ 220 , which is again within the range of what is expected due to stochastic sampling. We find a similarly good level of agreement for the other parameters. Finally, we find a speed-up of ∼ 10 times when using SEOBNRv4HM_PA instead of SEOBNRv4HM as the base gr model, significantly accelerating inference. With focus shifting from analyses of individual events to population studies, demands on computational resources and person-power are everincreasing, as demonstrated by large-scale studies in lvk catalog papers [5,7]. Hence, such increases in computational efficiency is immensely important for the future of gw data analysis.

VI. DISCUSSION AND CONCLUSIONS
We developed the SEOBNRv4HM_PA model, which combines the multipolar eob nr-informed model SEOBNRv4HM with the post-adiabatic (pa) approach for solving the (spin-aligned) binary dynamics (developed and used in the TEOBResumS models [70,71,90,91]). The resulting model is computationally cheap (at most on the order of seconds and less than 1 second for most of the parameter space) and highly accurate (unfaithfulness less than O(10 −3 ) when benchmarked against the SEOBNRv4HM model). Therefore, it can be used as a drop-in replacement for parameter estimation studies and tests of gr for many gw events where the use of the SEOBNRv4HM model would have been impractical.
In Sec. II we presented the eob formalism and the derivation of the pa equations, together with a discussion on the specifics of transitioning from the pa regime to the merger part of the dynamics. Section III describes the practical implementation of this model in the LALSimulation library of waveform models, and provides the important benchmarks in terms of speed and accuracy of the SEOBNRv4HM_PA model. In particular, we demonstrated that the SEOBNRv4HM_PA model provides a speed-up of 10 to 100 times compared to the SEOBNRv4HM waveform model. Furthermore, we showed that the new model is accurate at a level which allows us to use it for the purposes of ligo data analysis.
Section IV makes use of the results of the previous section and illustrates the use of the SEOBNRv4HM_PA model in 2 separate parameter estimation studies: recovering an injected synthetic signal and applying the model to analyse the event GW190412. In both cases, we find an extremely good agreement between the results obtained with our model and other established models. Compared to pe with the SEOBNRv4HM model, the study with SEOBNRv4HM_PA completes around 10 times faster. Finally, Sec. V depicts the use of the waveform model to a test of gr and discusses the importance of the model for this type of analyses. We find that there is an excellent agreement between the results obtained using the SEOBNRv4HM_PA and the established SEOBNRv4HM and SEOBNRv4HM_ROM models.
There are several interesting future directions that can be pursued to apply the pa approximation to other SEOBNRv4HM models, notably models that include tidal deformabilities for binary neutron stars [119,120], for which, currently, only the surrogate method has been applied [79] (see Refs. [105,121] for pa models with tidal effects using TEOBResumS). A more important but challenging extension would involve the precessing SEOBNRv4PHM Hamiltonian [46], which would allow the model to be used more efficiently and for a broader range of gw events when doing parameter-estimation studies. Finally, the pa approximation could also be extended to binary systems on eccentric dynamics. The main challenges with spin-precession and eccentricity is the presence of time scales beyond the orbital and radiation-reaction ones.