Frequency domain reduced order model of aligned-spin effective-one-body waveforms with higher-order modes

We present a frequency domain reduced order model (ROM) for the aligned-spin effective-one-body (EOB) model for binary black holes (BBHs) SEOBNRv4HM that includes the spherical harmonics modes $(\ell, |m|) = (2,1),(3,3),(4,4),(5,5)$ beyond the dominant $(\ell, |m|) = (2,2)$ mode. These higher modes are crucial to accurately represent the waveform emitted from asymmetric BBHs. We discuss a decomposition of the waveform, extending other methods in the literature, that allows us to accurately and efficiently capture the morphology of higher mode waveforms. We show that the ROM is very accurate with median (maximum) values of the unfaithfulness against SEOBNRv4HM lower than $0.001\% (0.03\%)$ for total masses in $[2.8,100] M_\odot$. For a total mass of $M = 300 M_\odot$ the median (maximum) value of the unfaithfulness increases up to $0.004\% (0.17\%)$. This is still at least an order of magnitude lower than the estimated accuracy of SEOBNRv4HM compared to numerical relativity simulations. The ROM is two orders of magnitude faster in generating a waveform compared to SEOBNRv4HM. Data analysis applications typically require $\mathcal{O}(10^6-10^8)$ waveform evaluations for which SEOBNRv4HM is in general too slow. The ROM is therefore crucial to allow the SEOBNRv4HM waveform to be used in searches and Bayesian parameter inference. We present a targeted parameter estimation study that shows the improvements in measuring binary parameters when using waveforms that includes higher modes and compare against three other waveform models.

We present a frequency domain reduced order model (ROM) for the aligned-spin effective-onebody (EOB) model for binary black holes (BBHs) SEOBNRv4HM that includes the spherical harmonics modes ( , |m|) = (2, 1), (3,3), (4,4), (5,5) beyond the dominant ( , |m|) = (2, 2) mode. These higher modes are crucial to accurately represent the waveform emitted from asymmetric BBHs. We discuss a decomposition of the waveform, extending other methods in the literature, that allows us to accurately and efficiently capture the morphology of higher mode waveforms. We show that the ROM is very accurate with median (maximum) values of the unfaithfulness against SEOBNRv4HM lower than 0.001%(0.03%) for total masses in [2.8, 100]M . For a total mass of M = 300M the median (maximum) value of the unfaithfulness increases up to 0.004%(0.17%). This is still at least an order of magnitude lower than the estimated accuracy of SEOBNRv4HM compared to numerical relativity simulations. The ROM is two orders of magnitude faster in generating a waveform compared to SEOBNRv4HM. Data analysis applications typically require O(10 6 − 10 8 ) waveform evaluations for which SEOBNRv4HM is in general too slow. The ROM is therefore crucial to allow the SEOBNRv4HM waveform to be used in searches and Bayesian parameter inference. We present a targeted parameter estimation study that shows the improvements in measuring binary parameters when using waveforms that includes higher modes and compare against three other waveform models.

I. INTRODUCTION
In the past five years gravitational-wave (GW) observations [1][2][3] have opened up a new window to the Universe. In the first two observing runs of the advanced LIGO [4] and Virgo [5] detectors ten confident detections of binary black-holes (BBHs) and one binary neutronstar (BNS) were made [3] and tens of detection candidates [6] have been identified so far in the third observing run of this network, among them another confident detection of a binary neutron star system [7]. Both the detection and inference of binary parameters for these compact binaries rely heavily on our knowledge of the gravitational waveform emitted in these coalescences as encapsulated in parametrized models of GWs. The construction of stochastic template banks and the Bayesian inference of binary parameters routinely require tens to hundreds of millions of waveform evaluations [8][9][10][11]. At the same time the phasing of the GWs needs to be tracked to an accuracy better than a fraction of a wave cycle to avoid missing events or mis-measuring binary parameters. Therefore, waveform models need to be fast and accurate to extract the binary properties imprinted in the emitted GWs.
Surrogate or reduced order modeling techniques [46][47][48][49][50][51][52][53][54] provide established methods for accelerating slow waveforms while retaining very high accuracy. These techniques have been successfully applied to EOB [23,[46][47][48]52] and NR [49-51, 55, 56] waveforms. They work by decomposing waveforms from a training set in orthonormal bases on sparse grids in time or frequency, and interpolating or fitting the resulting waveform data pieces over the binary parameter space. The result is a smooth, accurate (as tested against an independent validation set), and fast to evaluate (compared to the original waveform data) GW model. These surrogate models have proven invaluable for GW data analysis.
This paper is organized as follows. In Sec. II we give a brief description of the time-domain SEOBNRv4HM model. In Sec. III we discuss various techniques we use to build the ROM, notably waveform conditioning in Sec. III A. We continue with a summary of the basis construction and decomposition in Sec. III B, tensor product interpolation in Sec. III C. Domain decomposition in frequency and in parameter space is discussed in Sec. III D and Sec. III F, respectively. We summarize how we connect the ROM with PN solutions at low frequency in Sec. III E. We present results in Sec.
IV where we demonstrate the accuracy of the ROM in Sec IV A, and its computational efficiency in Sec. IV B. We showcase a parameter estimation application in Sec. IV C. Finally we conclude in Sec. V.

II. THE SEOBNRV4HM MODEL
The gravitational wave signal emitted by a coalescing binary black hole is usually divided into three different regimes: inspiral, merger and ringdown. During the inspiral the two black holes move at a relative speed v that is small compared to the speed of light c, therefore the solution to the two body problem can be found using a perturbative expansion in the small parameter v/c, the so-called PN expansion [64]. At some point, during the evolution of the binary system, the parameter v/c ceases to be small and the PN expansion is not valid anymore. This happens roughly at the innermost stable circular orbit (ISCO) and demarcates the end of the inspiral and the beginning of the merger regime. The signal in this regime can only be computed using NR simulations that solve Einstein's equations for a BBH system, fully numerically. Finally, in the ringdown regime, the perturbed black hole formed after the merger of the binary emits gravitational waves at frequencies that can be computed within the black hole perturbation theory formalism [65].
The EOB formalism, first introduced by Buonanno and Damour in Refs. [12,13], provides a natural framework to combine these three regimes and produce a complete waveform with inspiral, merger and ringdown. Within the EOB formalism the PN conservative dynamics of a BBH system during the inspiral is resummed in terms of the dynamics of a test particle with an effective mass and spin around a deformed Kerr metric. This improved conservative dynamics is combined with a resummed energy flux [66][67][68] to produce an inspiral waveform that is close to NR solutions. To improve the agreement with NR waveforms the EOB conservative dynamics is also calibrated using information from NR simulations [18,19]. In the EOB waveform the merger and ringdown part is built using a phenomenological fit produced using informations from NR waveforms and black hole perturbation theory [23,69]. NR-tuned versions of EOB models are usually referred to as EOBNR.

III. TECHNIQUES FOR BUILDING THE ROM
In this Section we describe the construction of our ROM, from the preparation of the waveforms to the reduced basis and interpolation techniques. We use techniques developed for previous ROM models [23,47,48], which we generalized to the higher-harmonics case.
We start with a general discussion of how to prepare and decompose waveform data for higher mode waveforms in Sec. III A. In particular, we discuss time domain conditioning in Sec. III A 1, stationary phase approximation in Sec. III A 2, the orbital carrier phase in Sec. III A 3, the introduction of coorbital modes in Sec. III A 4, scaling of frequencies in Sec. III A 5. We summarize basis construction in Sec. III B and tensor product interpolation in Sec. III C. We also explain how we decompose the model in both frequency range patches (see Sec. III D) and parameter space patches (see Sec. III F). Hybridization with inspiral waveforms is discussed in Sec. III E.

A. Preparation and decomposition of waveform data
The waveform polarizations h + , h × are decomposed in spin-weighted spherical harmonics as The h m are called the harmonics or simply the modes of the gravitational wave, with h 22 and h 2,−2 the dominant harmonics corresponding to quadrupolar radiation. These modes h m are affected by convention choices: first, by the choice of polarization vectors defining h + , h × , and secondly by the definition chosen for the source frame in which the waveform is described. For nonprecessing systems, the z-axis of the source frame is taken to be the normal to the orbital plane, with a residual freedom in choosing the origin of phase. One can take two points of view for the definition of phase: either fix the definition of the source frame (for instance, imposing that the initial separation vector is along x) and call "phase" the azimuthal angle of the observer in the source frame, or fix the direction to the observer (for instance in the plane (x, z)) and call "phase" the binary's orbital phase at a given time. We can also consider the definition of the origin of time as part of the source frame definition. During the inspiral, the individual harmonics obey a simple overall scaling with the orbital phase as but this scaling does not apply post-merger where the modes are driven by their respective ringdown frequencies.
There are several challenges regarding the conditioning of higher-harmonics waveforms for the purpose of reduced order modelling. We recall that one relies on two kinds of interpolation here: one is the interpolation of waveform pieces along the tracking parameter, either time or frequency, used to compress data; the other is the interpolation across the parameter space (masses and spins) used in the internals of the ROM, either of waveform quantities directly (as in [46, 49-52, 55, 56] in the empirical interpolation formalism), or as in our case, of reduced basis projection coefficients [23,47,48]. Both these interpolations require smoothness, and discrete jumps can cause significant (and non-local) errors.
As a result, zero-crossings in the subdominant harmonics h m (as noticed in Refs. [57,70]) cause difficulties for the usual amplitude/phase representation: if the envelope of a mode crosses zero with a positively defined amplitude, the phase jumps by π, a discontinuity that will harm the interpolation performed when reconstructing the waveform. Among other advantages, this is alleviated by the procedure used in [50,56] of modelling the waveform in a coorbital frame where the dominant phasing of Eq. (3.2) has been scaled out, so that a more robust real/imaginary representation can be chosen instead; here we will use the same kind of coorbital quantities, but built in the Fourier domain.
The natural 2π degeneracy in phase also requires care when interpolating across parameter space. Discrete 2π phase jumps leave the waveforms themselves invariant, but can break the interpolation in-between waveforms. This issue is particularly relevant when dealing with numerical Fourier transforms of time-domain waveforms: when phase-unwrapping the output of the Discrete Fourier Transform starting from f = 0, numerical noise causes the 2π interval to be essentially random. In [23,47,48] a linear fit of the Fourier-domain phase was removed. Here we will keep time and phase alignment information throughout the conditioning procedure, so instead we will impose a given 2π range for the phase at a reference point, corresponding to the time of alignment.
Other difficulties are caused by the relative alignment of the different harmonics. Dividing the phase of the dominant h 22 mode by 2, whether in time or frequency domain, comes with an ambiguity of π then propagated as mπ to the other modes. Such an ambiguity is not necessarily a problem if the phase alignment is done as a last step when generating the waveform (as is the case in the IMRPhenomHM model [32]); giving up the geometric interpretation of the source-frame definition, it is sufficient that a [0, 2π] range in the "phase" input by the user corresponds to a [0, 2π] range in geometric phase. It becomes a problem, however, when we need to interpolate across parameter space to build a ROM. In particular, when working from the Fourier domain waveform alone, we do not have access to the orbital phase (as read from trajectories) to lift these kind of degeneracies. Here we will make sure that the alignment is performed in the time domain before taking Fourier transforms, and we will further introduce an artificial carrier signal to have access to a proxy of the orbital phase in the Fourier domain.
We detail below our conditioning procedure, chosen to circumvent these issues.

Time-domain conditioning
In building this ROM, we will carry along time and phase alignment information all the way to the final Fourier-domain waveforms. This is in contrast to previous Fourier-domain waveform models (SEOBNRv4 ROM, IMRPhenomD) where the time and phase are adjusted after generating the waveform, as will be described below.
Individual harmonics are decomposed as an ampli-tude 2 and phase, following with the scaling In the early inspiral regime, for low frequencies, the phases ∆φ m are approximately constant. We choose the same polarization convention as in [71], for which we have in the low-frequency limit. When getting closer to merger, deviations from Eq. (3.5) become more important. In the notations of the EOB factorized waveforms [16,17], these phase deviations come from the phases e iδ m and tail factors T m (see Eqs. (14) and (21) in [17]), and from non-quasicircular corrections (NQCs) close to merger (see Eq. (22) in [17]). We choose the source frame convention for our model by imposing that its direction x is along the separation vector between the two black holes n(t align ), with an arbitrary time of alignment in the late inspiral t align = −1000M (with t = 0 being defined as the amplitude peak of h 22 ). In practice, rather than using n(t align ) we simply impose and we use the orbital phase φ orb as read from the EOB dynamics to resolve the π-ambiguity and impose φ orb 0 rather than φ orb π. These alignment properties will be reproduced, up to small numerical errors, by the reconstructed ROM waveforms.

Stationary phase approximation
As we will use it to guide our conditioning procedure, we recall here the Stationary phase approximation (SPA) for waveforms with higher harmonics. First, we introduce the Fourier transform for a time-domain signal h as (3.7) Note the sign difference in the exponential with respect to the more usual definition (used in particular in LAL [72]). This choice is made for convenience, as it will ensure that Fourier-domain modesh m with m > 0 and m < 0 have support at positive and negative frequencies respectively. One can come back to the LAL Fourier convention by the simple map f ↔ −f , which for real signals h(t) ∈ R translates ash LAL (f ) =h * (f ). Let us first consider a generic signal with an amplitude and phase as h(t) = a(t)e −iφ(t) and define ω ≡φ. The SPA is applicable under the assumptions |ȧ/(aω)| 1, ω/ω 2 1 and (ȧ/a) 2 /ω 1, that are well verified in the inspiral. Defining a time-to-frequency correspondence t(f ) implicitly by Applying this to the individual h m modes (3.3), treating the ∆φ m as constants, defining ω m =φ m mω orb , each mode will have a separate time-to-frequency correspondence and the phase (3.11) This is the Fourier-domain equivalent to the time-domain relation (3.4). We note useful relations between different mode numbers. The various t m (f ) functions are related by while the phases satisfy (3.13) This last relation holds regardless of the time and phase alignment of the waveform: as a phase change δφ m = mδφ, or a time change δΨ m (f ) = −2πf δt would both leave 2Ψ m (mf /2) − mΨ 22 (f ) invariant. It is sensitive however to the quantities ∆φ m (that we treat here as constants in the early inspiral), that depend on the choice of polarization convention.
Finally, we recall that we can build a time-to-frequency correspondence directly from the Fourier-domain waveform as (3.14) Note that this definition of time is strictly speaking only accurate in the inspiral phase, where the SPA applies and the two definitions (3.14) and (3.8) coincide. However, it can be used as a proxy for time everywhere, as we can evaluate (3.14) for any frequency 3 f .

Orbital carrier
In order to carry over information about the alignement of the respective mode from the time domain to the Fourier domain, we find it convenient to introduce a fictitious carrier signal k(t), that evolves with the orbital phase instead of twice the orbital phase as the h 22 mode does.
The choice made here of keeping the same amplitude as the h 22 mode is quite arbitrary, but will ensure that it decays in the ringdown, giving us a smooth Fourier transform for this carrier. Note that this construction is artificial, as the carrier does not correspond to any physical signal.
As mentioned before, the carrier half-phase φ 22 /2 comes with a π-degeneracy. We can alleviate this by forcing the carrier phase to be within π of the orbital phase, as read from the SEOB dynamics, at the time of alignment. This is, in fact, our main motivation for building this carrier in the time domain: it allows us to avoid the issues listed above, with all the conditioning being ultimately tied to the orbital phase, a quantity that is smooth across parameter space.
The Fourier transform of the carrier signal introduced in (3.15) is decomposed as an amplitude and phase as where A k = |k| will be discarded in the rest of the conditioning. When the SPA applies, we have approximately Sincek(f ) is computed via an FFT, nothing forbids arbitrary jumps of 2π of the phases between different points in parameter space. We use the relation above to remove this 2π-ambiguity in Ψ k . At the frequency f align such that t(f align ) = t align , we impose (3.19) In this way, Ψ k is directly tied to φ orb that is smooth in parameter space in our time-domain conditioning procedure.
We will factor out the Fourier domain phase of the carrier, build a ROM for the carrier separately, and then factor in the modelled carrier phase when reconstructing the waveform.

Fourier-domain coorbital modes and waveform building blocks
Next, we build coorbital modes by scaling out the Fourier-domain phase of the carrier following These modes are built so as to factor out the main contribution to the phase of the Fourier-domain modes, to leave the coorbital modesh c m with an approximately constant phase in the inspiral regime. Namely, for the inspiral regime, where the SPA is valid, t k (f ) = t m (mf ) and applying (3.11) and (3.17) gives Note that our Fourier-domain construction is approximate, and these "coorbital" quantitiesh c m do not correspond exactly to a coorbital frame defined in the time domain as in [50,56].
We stress that these modes are not strictly coorbital, in the sense that they are not built from a time-domain coorbital frame built from the orbital phase. Indeed, the definition (3.20) is rooted in the Fourier domain, and its physical meaning is unclear in the high-frequency range where the SPA does not apply anymore.
Thus, the basic building blocks for the ROM will be • Ψ k = −Arg k , the Fourier-domain carrier phase; • Re hc m , the real part of the coorbital modes; • Im hc m , the imaginary part of the coorbital modes.
Conversely, to rebuild the modesh m from these waveform pieces, it is enough to factor in the carrier phase as in (3.20).

Scaling of frequencies using ringdown frequency
One of the prerequisites of our ROM procedure is to represent the waveform on a common frequency grid. However, the frequency range covered varies with physical parameters, notably spin. In SEOBNRv4 ROM, this was alleviated by extending waveforms to higher frequencies.
Here, we choose to apply a scaling to the frequencies of the waveform building blocks, depending on the ringdown frequency. For every mode ( , m) and the carrier we define where ω QNM m is the quasi-normal mode frequency, and varies for different waveforms as it depends on the spin of the remnant black hole. We will then use for all waveforms a common grid of this rescaled parameter y. Given this scaling, we have to carefully adjust the starting frequency of the waveforms of our training set so that the frequency range of the carrier phase Ψ k covers all modes after undoing the scaling. The maximal values of y m , y k where we cut the data are (y max 22 , y max 21 , y max 33 , y max 44 , y max 55 ) = (1.7, 1.7, 1.55, 1.35, 1.25) and y max k = 2.5. This technique is only used for building the high-frequency ROM (see Sec. III F); for the lowfrequency ROM, the ringdown frequency is irrelevant and the scaling would induce an additional cost in generating the waveforms of the training set.

B. SVD decomposition
We decompose all waveform data pieces defined in Sec. III A 4 into respective singular value decomposition (SVD) bases and subsequently interpolate the projection coefficients in each SVD basis over the parameter space, as discussed in Sec. III C. This method follows earlier work in [23,47,48].
We start with a waveform data piece X(f i ; θ), given on a discrete grid in frequency {f i } m i=1 , and on a regular grid of points θ in the three-dimensional binary parameter space in mass-ratio q and aligned BH spins χ i , (q, χ 1 , χ 2 ). We flatten the parameter grid and arrange the data in matrix form X ij = X(f i ; θ j ) ∈ R m×n , where n is the total number of input waveforms.
We then resample the data in a log-spaced frequency grid of 300 points. The number of points used for resampling are based on previous studies (see Ref. [23,48]). We compute the SVD [73, 74] X = V ΣU T and obtain an orthonormal basis for the column space of the matrix X in the first r = rank(X) columns of V . The SVD provides us with a decomposition of the range space of X, range(X) = span{v 1 , . . . , v r }, where the v j are the left singular vectors of X.
Given the basis B X = V , we expand the waveform data pieces x j (f i ) that make up the columns of X in this basis and can write the expansion To construct a waveform model we need to predict the coefficients c X at a desired parameter space point θ * . To do that we need to fit or interpolate c X over the parameter space. This is discussed in Sec. III C.

C. Tensor-product spline interpolation
In the following we describe how we obtain projection coefficients at arbitrary parameter space points. In low dimensional spaces we can afford to use dense grids built from the Cartesian product of one-dimensional sets of points. We choose cubic splines as the univariate interpolants and obtain a tensor-product interpolant (TPI) [47,54,75] for a three-dimensional coefficient tensor c ijk which can be written as (3.23) Here, the Ψ are B-spline basis functions [76] of order 3 for the chosen one dimensional sets of parameter space points in each dimension. We use "not-a-knot" boundary conditions to avoid having to specify derivatives at the domain boundaries. We built the model using TPI, a Cython/C package [75] to provide tensor product spline interpolation in Python, and later implemented the model in LAL [72].

D. Patching in geometric frequency
Here we discuss dividing the waveform domain in geometric frequency into separate sub-domains, where we build a separate ROMs. In Sec. III F we will instead discuss how to tackle non-uniform resolution requirements over the binary parameter space.
In the early inspiral waveforms modes tend to be well approximated by the PN expansion and an accurate ROM can be built from a relatively small training set. In contrast, the high geometric frequency part of the waveform modes encapsulates the late inspiral, merger and ringdown part of the signal which is more complex and non-linear, and this is also the regime where EOB waveforms are tuned against NR where they are available. Building an accurate ROM for the high geometric frequency part of the waveform modes consequently requires a higher density of training set waveforms.
Therefore, it is natural to treat the low and high geometric frequency part of the waveform separately, following the construction of previous ROMs [48]. This allows us to make the training set for the low geometric frequency part of the waveforms significantly smaller and reduce the computational cost of the training. Waveforms for low mass binaries are the most costly waveforms to generate. The cost is exacerbated due to the presence of higher modes with |m| > 2, since they require a lower starting orbital frequency to cover the same gravitational wave frequency range as the dominant mode. In Fig. 1 we show the sub-division into low and high geometric frequency sub-domains. The low frequency sub-domain is connected with PN waveforms modes in the early inspiral, as discussed in Sec. III E. We generated the SEOBNRv4HM waveforms at a sufficiently low frequency (at 15 Hz and a total mass of 5M to allow for some tapering) to have the complete set of higher harmonics included in SEOBNRv4HM present at a frequency M f = 0.0005 * 5/2 ≈ 0.0012, where the low geometric frequency sub-domain starts. We choose the geometric transition frequency between the low and high frequency sub-domains to be M f = 0.003 * m, using the natural inspiral scaling of the frequency of the waveform modes with m.. For the high frequency sub-domain we generated waveforms choosing the starting frequency as described in Sec. III A 5, ensuring the generated waveforms after undoing the scaling (3.22a) will cover this transition frequency. The complete waveform modes are then generated by blending together the low and high frequency parts at the frequency using a variant of the Planck taper function described in Ref. [77].

E. Hybridization with TaylorF2
Here we describe how we carry out the hybridization of the ROM waveform with the TaylorF2 inspiral waveform.
After evaluating the ROM waveform for all modes, we generate the TaylorF2 amplitude and phase for the (2, 2) mode from the lowest frequency necessary to be able to start all inspiral modes at a user specified frequency. We blend the TaylorF2 and ROM amplitude and phase for the (2, 2) mode using the same Planck taper function used to connect high and low frequency ROM. We can obtain the higher mode PN inspiral waveforms by rescaling the (2, 2) amplitude and phase. For the phase we follow Eq. (3.21) and Eq. (3.5) to compute the carrier phase from the TaylorF2 (2,2) phase and rescale it to obtain the phase for each mode. We then align the inspiral phase with the ROM phase for each mode and blend them together on a common frequency grid. For the amplitude we rescale the TaylorF2 (2,2) amplitude according to the PN amplitudes given in Ref. [78] Eqs.(12a-12t).

F. Patching in parameter space
As already noted for the previous ROMs of EOBNR waveforms (see Refs. [23,47,48]) model features often require more resolution in particular parts of the parameter space. However, regular grids do not allow for local refinement in (q, χ 1 , χ 2 ). Therefore, we partition the binary parameter space into subdomains on which resolution requirements can be satisfied with a particular regular grid choice.
The low frequency ROM does not need any special treatment and was built using waveforms placed on a Cartesian grid with 64 points in q and 12 points in χ 1,2 as shown in Fig. 2. Here, the 1D grids in χ 1 and χ 2 were chosen to be identical. The grids for q and χ 1,2 are the same as the ones used for SEOBNRv4 ROM (see Sec.VII in Ref. [23]), except that we limit the grid to q = 50.
On the other hand, as already noted in Refs. [23,48], modeling the non-linear merger and ringdown part of the waveform in the high geometric frequency ROM requires a higher resolution when approaching large positive values of the primary spin. Therefore, we build two different high frequency ROMs based on the value of the primary's spin, with one ROM having a finer grid in the region of high χ 1 . The inclusion of higher modes in the SEOBNRv4HM ROM model also requires additional resolution near equal mass. The modes with odd m vanish by symmetry on the line q = 1 and χ 1 = χ 2 and their behavior in the vicinity is non-trivial to model (see Refs. [57,70]). Therefore, we build two different high geometric frequency ROMs in mass-ratio, one of which is   Fig. 3. The physical domain covered by each patch is defined by a Cartesian product of intervals in binary parameters (q, χ1, χ2). We also indicate the number of grid points in each parameter per patch.
covering the region q → 1 with a finer grid. In total we then have four high frequency ROMs to correctly model the merger and ringdown part of the signal.
The 2D projection of the grid in (q, χ 1 ) for these four ROMs is shown in Fig. 3. Since no special choice is made for the grid in χ 2 we have omitted plotting the grid in this dimension. We set domain boundaries at q = 3 and χ 1 = 0.8 for these four ROMs. In Table I we collect information on how the four patches are placed in parameter space and the number of gridpoints in each dimension.

IV. RESULTS
In this section we discuss the accuracy and the increase in efficiency of SEOBNRv4HM ROM compared to SEOBNRv4HM. Finally we also perform a parameter estimation study to demonstrate the potential of this model in data analysis applications.

A. Accuracy of the model
We start by defining the faithfulness function, used to assess the closeness between two waveforms when higherorder modes are included. We then use this faithfulness measure to determine how accurately the ROM reproduces SEOBNRv4HM waveforms.

Definition of faithfulness
A GW signal emitted by a spinning, non-precessing and non-eccentric BBH is characterized by 11 parameters in the detector frame. These parameters are the BH masses m 1 and m 2 , the (constant) projection of the spins in the direction perpendicular to the orbital plane χ 1 and χ 2 , the angular position of the line of sight measured in the source's frame (ι, ϕ 0 ), the sky location of the source in the detector frame (θ, φ), the luminosity distance D L , the time of arrival t c of the signal and finally the polarization angle ψ. The detector response can be written as: where masses and spins are combined in the vector ξ ≡ (m 1 , m 2 , χ 1 , χ 2 ), and the functions F + (θ, φ, ψ) and F × (θ, φ, ψ) are the antenna patterns [79,80]. This equation can be rewritten as: with κ(θ, φ, ψ) being the effective polarization [59] defined in the range [0, 2π) as: where the function A(θ, φ) is an overall amplitude and is defined as: Given a GW signal h s (SEOBNRv4HM in our case) and a template waveform h t (SEOBNRv4HM ROM in this context), we define the faithfulness (or match) as [57,59,81]: where parameters with the subscript "s" ("t") refer to the signal (template) waveform. The expression above does not depend on A(θ, φ), therefore the only dependance on (θ, φ, ψ) is encoded in κ(θ, φ, ψ). For the faithfulness calculation we optimize over the phases ϕ 0t and κ t and the time of arrival t c because they are not interesting from an astrophysical perspective. We use the standard definition of the inner product between two waveforms (see [79,80]): where˜denotes the Fourier transform, * indicates the complex conjugate and S n (f ) is the one-sided power spectral density (PSD) of the detector noise. For this computation we use the Advanced LIGO "zero-detuned high-power" design sensitivity curve [82]. The integral is computed between the frequencies f low = 20Hz and f high = 3kHz. The same definition of faithfulness has been used to determine the agreement between SEOBNRv4HM and numerical relativity waveforms (see [57]) 4 . Since the faithfulness given in Eq. (4.5) depends on the signal parameters (ι s , ϕ 0s , κ s ), we will summarize the results using the maximum and the average unfaithfulness (or mismatch) [1−F(ι s , ϕ 0s , κ s )] over these parameters, namely [59,81,84]:

Faithfulness against SEOBNRv4HM
In order to avoid biases in data analysis applications when using the ROM instead of SEOBNRv4HM, it is important to verify that the additional modeling error introduced in the construction of the ROM is negligible compared to the inaccuracy of the SEOBNRv4HM waveforms with respect to the NR simulations. Since the typical unfaithfulness between SEOBNRv4HM and NR waveforms is O(1%) (see Figs. (11) and (12) in Ref. [57]), it is therefore natural to require the unfaithfulness between SEOBNRv4HM and SEOBNRv4HM ROM to be O(0.1%) or less. To that end we have generated 10000 SEOBNRv4HM waveforms with random (uniformly distributed) values of (q, χ 1 , χ 2 ) and computed their match against the same waveforms produced with SEOBNRv4HM ROM.
We summarize these results in Fig. 4 where we show a histogram with the unfaithfulnessŪ computed between the ROM and SEOBNRv4HM waveforms for different values of the total mass. For each total mass we report in Tab. II the median and maximum values of these unfaithfulness The median of these mismatch distributions is weakly dependent on the total mass and it is always ≤ 0.002% while their maximum value is always ≤ 0.08%. In Fig. 5 we display the distribution of mismatches shown in Fig. 4 as a function of (q, χ 1 ) and for different values of the total mass. The largest mismatches between the ROM and SEOBNRv4HM are obtained for M = 300.0M and large negative χ 1 , as it is clear from Fig. 4 and Fig. 5 (bottom right panel). ROM GW modes are generated up to a maximum frequency (in geometric units) that scales with the inverse of the total mass of the system. For large total masses the lack of signal above this maximum frequency is the main source of inaccuracy of the ROM. This maximum frequency for each mode is proportional to its least damped quasi-normal mode frequency as described in Eq. 3.22a. The mismatch is larger for large negative spins because the least damped quasi-normal mode frequency decreases in this region of the parameter space. We highlight that in this region the ROM waveforms still have mismatches 0.1% against SEOBNRv4HM  waveforms as demanded at the beginning of this section. The results described above do not change substantially when considering the distribution of U max instead ofŪ. In Tab. III we report the median and maximum values of these distributions.
For total masses M ≤ 20M it is more convenient to summarize the results of the faithfulness calculations as an histogram with a fixed m 2 instead of the total mass. In Fig. 6 we show theŪ distribution when fixing m 2 = 1.4M and varying m 1 in the interval 1.4M ≤ m 1 ≤ 18.6M such that the total mass of the system is always M ≤ 20M . The median of this distribution is 0.0003% while its maximum is 0.01%. In Fig. 7 we report theŪ distribution in Fig. 6 as a function of (m 1 , χ 1 ). The accuracy of the ROM in this case degrades for large values of m 1 and large positive spins but it is still well within the requirements. Also in this case the results are not very different when looking at the U max distribution for which the median is still 0.0003% while the maximum increases to 0.02%.
These analyses demonstrate that the modeling error introduced by the ROM is negligible with respect to the difference between SEOBNRv4HM and NR waveforms. For this reason the mismatch of the ROM against the NR waveforms is essentially the same as SEOBNRv4HM (see Figs.11-12 in Ref [57] and Fig.6 in Ref. [55] 5 ).

B. Computational performance
In this section we discuss the computational performance of the ROM in terms of walltime for generating a waveform. We first compare the ROM to SEOBNRv4HM and then to other waveform models that include higherorder modes.
FIG. 5: UnfaithfulnessŪ between the ROM and SEOBNRv4HM as a function of (q, χ1) and for different values of the total mass. For M = 20M there are no data in the region q > 19 because for these system m2 would have an unphysical subsolar mass.

Speedup with respect to SEOBNRv4HM
The speedup of the ROM with respect to SEOBNRv4HM is computed by the ratio of the walltimes of the two models for generating a frequency domain waveform at the same parameters. Since SEOBNRv4HM is a time domain model, we first generate the waveform in the time domain at a sample rate of 16384 Hz and then compute its Fourier transform. The ROM waveform is already in the frequency domain and it is generated using the sampling interval set to 1/T where T is the duration in seconds of the associated time domain waveform. The maximum frequency of the SEOBNRv4HM ROM waveform is set to 8192 Hz.
In Fig. 8 we show this speedup as a function of the total mass and for different values of the mass ratio. The speedup is of order 100. It increases with mass ratio and decreases with total mass. The maximum speedup is found around a total mass of 50M . Since the spins have only a limited effect on the waveform duration, the speedup depends only weakly on them.

C. Parameter estimation study
In this section we use the SEOBNRv4HM ROM model in a parameter estimation application. For this purpose we create two mock signals (or injections) with the same binary parameters, using either SEOBNRv4HM or NRHybSur3dq8 to generate the waveform. We then use SEOBNRv4HM ROM, SEOBNRv4 ROM, and, as a comparison between waveform models for the second case, IMRPhenomHM [32] 6 and NRHybSur3dq8 to compute posterior distributions from the mock signals. The analysis of the first mock signal will demonstrate the improvements in measuring binary parameters when using a model with higher harmonics with respect to a model including only the dominant ( , |m|) = (2, 2) mode. The analysis of the second mock signal will give us a sense of possible biases due to modeling errors in the original SEOBNRv4HM model with respect to NR-surrogate waveforms, which are close to NR simulations. In creating the mock signals we do not add detector noise. This choice is made to avoid additional uncertainty and bias introduced by a random noise realization. It is the natural choice given that the goal of this parameter estimation analysis is to check for possible biases due to inaccuracies in waveform models.

Setup
We choose parameters for the mock signals in order to emphasize the effect of higher-modes in the waveform. Since the contribution of higher-order modes in the emitted GWs increases with the mass ratio, we use for the mock signals q = 8, the largest mass ratio available for the model NRHybSur3dq8. For the total mass we use M = 67.5M such that the values of the component masses m 1 = 60M and m 2 = 7.5M are consistent with the masses of BBH systems observed during O1 and O2 (see [3] and [86]). The waveform models are restricted to non-precessing spins and we pick χ 1z = 0.5 and χ 2z = 0.3. To maximize the effect of the higher-modes we inject the signal at edge-on inclination (ι = π/2) with respect to the observer. The coalescence phase φ c is chosen to be 1.2 rad, while the polarization phase ψ is set to 0.7 rad. The signal has been injected at gps-time 1249852257 s with a sky-position defined by its right ascension of 0.33 rad and its declination of -0.6 rad. Finally the distance of the mock signal is chosen by demanding a network-SNR of 21.8 in the three detectors (LIGO Hanford, LIGO Livingston and Virgo) when using the Advanced LIGO and Advanced Virgo PSD at design sensitivity [82]. The resulting distance is 627 Mpc. We used PyCBC's pycbc generate hwinj [87] to prepare the mock signal. To carry out Bayesian parameter estimation we used the Markov chain Monte Carlo code LALInferenceMCMC [10,88]. We choose a uniform prior in component masses in the range [3,100]M . Aligned component spins are assumed to be uniform in [−1, 1]. The priors on the other parameters are the standard ones described in Appendix C.1 of Ref. [3].

Results
Let us first focus on the case in which the mock signal is generated with SEOBNRv4HM. In Fig. 10 we summarize the results of the parameter estimation analysis for some relevant binary parameters. The top left panel shows the marginalized 2D posterior for the component sourceframe masses, and the top right panel the marginalized 2D posterior for the mass ratio q and the spin parameter χ eff = (m 1 χ 1 + m 2 χ 2 )/(m 1 + m 2 ). In the bottom left panel we present the marginalized 2D posterior with inclination ι and luminosity distance d L and, finally, in the bottom right panel, we report the matched filter SNR. The star in the plots corresponds to the true value used for the mock signal, while the 2D contours of the posterior distributions represent 90% credible regions. The waveform templates used to infer binary parameters are SEOBNRv4 ROM (blue curve) and SEOBNRv4HM ROM (red curve). It is clear from the plots that all the parameters reported in Fig. 10 are more precisely measured when using SEOBNRv4HM ROM instead of SEOBNRv4 ROM. The posterior volume represents the degeneracy of the gravitational wave signal, and, in the absence of detector noise, this degeneracy is intrinsic to the waveforms. The inclusion of higher harmonics in SEOBNRv4HM ROM breaks the degeneracy between the parameters q−χ eff and ι−D L and allows to measure them more precisely. These results are consistent with what was previously found in the literature [63,89,90]. As expected SEOBNRv4HM ROM also measures a larger matched filter SNR.
Let us now consider the case in which the mock signal is represented by NRHybSur3dq8. In Fig. 11 we show the marginalized 2D posterior for the mass ratio q and the spin parameter χ eff (left panel) and the marginalized 2D posterior with inclination ι and luminosity distance d L (right panel) as measured by the waveform models SEOBNRv4HM ROM (red curve), SEOBNRv4 ROM (blue curve), IMRPhenomHM (green curve), and NRHybSur3dq8 (orange curve). As before the star in the plots corresponds to the true value used for the mock signal, while the 2D contours of the posterior distributions represent 90% credible regions. From the plots in Fig. 11 it is clear that, as before, SEOBNRv4HM ROM recovers the binary parameters more precisely than SEOBNRv4 ROM. It is important to highlight that with SEOBNRv4HM ROM the binary parameters are recovered inside the 90% credible regions. This means that for this quite asymmetric system at a moderately high SNR of ∼ 20 the bias due to modeling errors in the original SEOBNRv4HM model compared to NR waveforms is negligible with respect to the statistical uncertainty. In contrast, the marginal posterior distributions recovered for IMRPhenomHM are in general broader compared to the ones recovered by SEOBNRv4HM ROM, are notably bimodal in mass-ratio and effective spin, and extend a lot further along the line of q -χ eff degeneracy. In distance and inclination the IMRPhenomHM posterior shows little improvement over SEOBNRv4 ROM which does not include higher harmonics. Finally, the marginal posteriors for NRHybSur3dq8 are quite similar in size to those for SEOBNRv4HM ROM, but better centered around the true parameter values. This is as expected since the likelihood should peak at the true parameter values when the signal and template use the same waveform. The mock signal is also sufficiently loud for the posteriors to be likelihoodrather than prior-dominated, resulting in an unbiased parameter recovery.
This study shows that using SEOBNRv4HM ROM for parameter estimation yields unbiased measurements of the binary parameters at moderately high SNR even in a configuration where the effect of higher harmonics in the waveform is large. We defer a more comprehensive analysis to future studies.
While the construction of this Fourier domain ROM broadly follows previous work [47,48] we have introduced the following new features to accurately represent the higher harmonics and make the model more flexible (see Sec. III). While previous models used an amplitude / phase decomposition of the Fourier domain waveform, we here define a carrier signal (see Eq. (3.15)) based on the time domain orbital phase. Subsequently we extract the carrier phasing from each Fourier domain waveform mode (see Eq. (3.20)). This essentially makes the phase of modes almost constant in the inspiral, defining what we here call "coorbital modes". This choice allows us to avoid zero-crossings in the subdominant harmonics which could spoil the smoothness of the training data and make accurate interpolation of the waveform data over parameter space very difficult. We perform alignment in the time domain to keep track of the time of coalescence of the training set waveforms and this information is preserved in the ROM. We use here an alternative approach to dealing with the fact that the ringdown frequency varies over the parameter space, but waveform data needs to be given on a common frequency grid to build a ROM. We rescale the geometric frequency parameter so that the ringdown is reached before a fixed termination frequency which demarcates the end of the frequency grid. We use the inverse rescaling during the evaluation of the ROM. We extend the ROM to arbitrarily low frequencies by splicing it together with multipolar PN waveforms. Therefore, it can in principle be used for arbitrarily light compact binary systems. We decompose waveform input data in orthonormal bases using the SVD, and build a model by constructing a tensor product spline over the 3-dimensional parameter space of mass-ratio and the two aligned spins of a binary. To increase model accuracy and efficiency we use domain decomposition in frequency and in parameter space.
In Sec. IV A we demonstrate that the ROM has a very high faithfulness (or match) with SEOBNRv4HM. Maximizing over inclination, reference phase and effective polarization of the source waveform (see Eq. (4.7)) the maximum mismatch over the remaining source parameters is below 0.03% for binaries with a total mass below 100M and below 0.2% for binaries at 300M (see Table III). Even for this very conservative choice the mismatch is at least an order of magnitude lower than the unfaithfulness of SEOBNRv4HM against NR simulations. Therefore the additional modeling error introduced in building the ROM is strongly subdominant and the ROM very accu-rately represents the SEOBNRv4HM waveform model. In Sec. IV B we show that our ROM accelerates waveform evaluation by a factor 100 -200 compared to SEOBNRv4HM and favorably compares against other higher mode waveform models for BBH systems, being about an order of magnitude faster. We showcase in Sec. IV C (see Figs. 10, and 11) that our ROM can recover component masses and spins, and especially distance and inclination angle for a quite asymmetric and spinning BBH with increased precision compared to the SEOBNRv4 ROM waveform which only models the dominant mode. In addition we show that the ROM accurately recovers binary parameters, irrespective of whether the source is represented by a SEOBNRv4HM or a NRHybSur3dq8 waveform. Our ROM gives a significantly more accurate parameter recovery compared to the phenomenological IMRPhenomHM waveform and is close to the NRHybSur3dq NR-surrogate model, while being more versatile and covering a significantly larger parameter space.
This ROM should prove a very useful tool for GW data analysis to describe systems where the contribution of higher harmonics is important in terms of additional signal-to-noise-ratio and discriminating power for detection and parameter inference. We stress that the ROM is very fast and reproduces the SEOBNRv4HM model with a great accuracy over the widest range in parameter space of all inspiral-merger-ringdown higher mode models available to date, from mass-ratio 1 to 1:50 where aligned spins can take values in the full range allowed for Kerr BHs, up to extremal spins.