Surrogate model for an aligned-spin effective one body waveform model of binary neutron star inspirals using Gaussian process regression

TL;DR

This paper develops a fast, accurate frequency-domain surrogate model for aligned-spin binary neutron star waveforms using Gaussian process regression, enabling efficient parameter estimation for gravitational wave data.

Contribution

It introduces a Gaussian process regression-based surrogate for the SEOBNRv4T waveform model, significantly reducing evaluation time while maintaining high fidelity.

Findings

Maximum mismatch of 4.5e-4 for LIGO

Evaluation time reduced to 0.13 seconds

Accurate parameter recovery in tests

Abstract

Fast and accurate waveform models are necessary for measuring the properties of inspiraling binary neutron star systems such as GW170817. We present a frequency-domain surrogate version of the aligned-spin binary neutron star waveform model using the effective one body formalism known as SEOBNRv4T. This model includes the quadrupolar and octopolar adiabatic and dynamical tides. The version presented here is improved by the inclusion of the spin-induced quadrupole moment effect, and completed by a prescription for tapering the end of the waveform to qualitatively reproduce numerical relativity simulations. The resulting model has 14 intrinsic parameters. We reduce its dimensionality by using universal relations that approximate all matter effects in terms of the leading quadrupolar tidal parameters. The implementation of the time-domain model can take up to an hour to evaluate using a starting frequency of 20Hz, and this is too slow for many parameter estimation codes that require $O(10^7)$ sequential waveform evaluations. We therefore construct a fast and faithful frequency-domain surrogate of this model using Gaussian process regression. The resulting surrogate has a maximum mismatch of $4.5\times 10^{-4}$ for the Advanced LIGO detector, and requires 0.13s to evaluate for a waveform with a starting frequency of 20Hz. Finally, we perform an end-to-end test of the surrogate with a set of parameter estimation runs, and find that the surrogate accurately recovers the parameters of injected waveforms.

Equations

$${\Delta }_u={\Delta }_u^{\text{pm}}+{\Delta }_u^{\text{DT}},$$

$${\Delta }_u^{\text{DT}}=$$

$$-15{\Lambda }_{3,\text{dyn}}^A(u)X_A^6{X}_B{u}^8\left[1+\left(-2+\frac{15}{2}{X}_A\right)u+\left(\frac{8}{3}-\frac{311}{24}{X}_A+\frac{110}{3}{X}_A^2\right)u^2\right]+(A\leftrightarrow B).$$

$${\Lambda }_{\ell ,\text{dyn}}^{A,B}(u)={\Lambda }_{\ell }^{A,B}{\hat{k}}_{\ell ,\text{dyn}}(u;{\hat{\omega }}_{0\ell }^{A,B}),$$

$${\Lambda }_{\ell }^{A,B}=\frac{2}{(2\ell -1)!!}\frac{k_{\ell }^{A,B}}{C_{A,B}^{2\ell +1}},$$

$$x=\frac{{\hat{\omega }}_{0\ell }}{m\Omega },ϵ\equiv \frac{256\eta }{5}{\left(\frac{{\hat{\omega }}_{0\ell }}{m}\right)}^{5/3},\hat{t}\equiv \frac{8}{5\sqrt{ϵ}}\left[1-x^{5/3}\right],$$

$${\hat{k}}_{\ell ,\mathrm{dyn}}(u;{\omega }_{0\ell })=a_{\ell }+b_{\ell }\left[f(x)+\sqrt{\frac{\pi }{3ϵ}}{x}^2\mathcal{Q}\right],$$

$$f(x)=x^2\left[\frac{1}{x^2-1}+\frac{5}{6}\frac{1}{1-x^{5/3}}\right].$$

$$\mathcal{Q}$$

$$-\sin \left(\frac{3}{8}{\hat{t}}^2\right)\left[1+2F_C\left(\frac{\sqrt{3}}{2\sqrt{\pi }}\hat{t}\right)\right],$$

$$h_{\ell m}^{\text{insp-plunge}}=h_{\ell m}^{\text{pm}}+h_{\ell m}^{\text{tidal}},$$

$$H_{\mathrm{eff}}=H_S+{\beta }^i{p}_i+\alpha \sqrt{{\mu }^2+{\gamma }^{ij}{p}_i{p}_j+{\mathcal{Q}}_4(p)}+H_{\mathrm{extra}}^{\mathrm{BBH}},$$

$$H_{\mathrm{extra}}^{\mathrm{BBH}}=\frac{1}{2r^3}(3n_i{n}_j-{\delta }_{ij})\frac{\mu }{M}{S}_{\ast }^i{S}_{\ast }^j,$$

$$S_{\ast }^i=\frac{m_2}{m_1}{S}_1^i+\frac{m_1}{m_2}{S}_2^i.$$

$$H_{\mathrm{extra}}^{\mathrm{BNS}}$$

$$\displaystyle\;\cdot\bigg{[}\frac{\mu}{M}S_{*}^{i}S_{*}^{j}+(\kappa_{1}-1)\frac{m_{2}}{m_{1}}S_{1}^{i}S_{1}^{j}+(\kappa_{2}-1)\frac{m_{1}}{m_{2}}S_{2}^{i}S_{2}^{j}\bigg{]}$$

$$h_{22}^{\text{LO-SS}}=-\frac{8\eta M^3{\omega }^2}{R}\sqrt{\frac{\pi }{5}}{e}^{-2i\Phi }\mathcal{A},$$

$${\mathcal{A}}^{\mathrm{BBH}}=\frac{1}{M^2}\left(\frac{S_1^2}{m_1^2}+2\frac{S_1{S}_2}{m_1{m}_2}+\frac{S_2^2}{m_2^2}\right),$$

$${\mathcal{A}}^{\mathrm{BNS}}=\frac{1}{M^2}\left({\kappa }_1\frac{S_1^2}{m_1^2}+2\frac{S_1{S}_2}{m_1{m}_2}+{\kappa }_2\frac{S_2^2}{m_2^2}\right).$$

$${\mathcal{F}}^{\text{LO-SS}}=\frac{32{\eta }^2}{5}{v}^{14}\left[2\mathcal{A}+\frac{1}{16M^2}{\left(\frac{S_1}{m_1}-\frac{S_2}{m_2}\right)}^2\right],$$

$${\rho }_{22}^{\text{LO-SS}}=\frac{1}{2}\mathcal{A}{v}^4,$$

$$t_{\text{match}}=\min \left(t_{\text{peak}}^{\text{amp}},t_{\text{peak}}^{\text{freq}}\right).$$

$$\displaystyle\hat{A}(t)=\left\{\begin{array}[]{ll}|h_{22}^{\rm insp-plunge}(t)|,&t\leq t_{\rm match}\,,\\ a+b(t-t_{\rm match}),&t>t_{\rm match}\,,\end{array}\right.$$

$$W(t)=\frac{1}{1+\exp [(t-t_{\text{match}}^{22}-15M)/\tau ]}.$$

$$\omega (t)={\omega }_{\mathrm{match}}-\Delta \omega \exp \left[-(t-t_{\mathrm{freq}})/(12M)\right],$$

$$\displaystyle\phi_{22}(t)=\left\{\begin{array}[]{ll}\phi_{22}^{\textrm{insp-plunge}}(t),&t\leq t_{\rm freq}\,,\\ \phi_{\rm freq}+\omega_{\rm match}(t-t_{\rm freq})+12M\Delta\omega\left\{\right.&\\ \left.\exp\left[-(t-t_{\rm freq})/(12M)\right]-1\right\}\,,&t>t_{\rm freq}\,.\end{array}\right.$$

$$\kappa=\left\{\begin{array}[]{ll}1+f\Lambda_{2}+g\Lambda_{2}^{2}+h\Lambda_{2}^{3}\,,&0\leq\Lambda_{2}\leq 1\\ e^{a+b\xi+c\xi^{2}+d\xi^{3}+e\xi^{4}}\,,&\Lambda_{2}>1\end{array}\right.,$$

$$\Lambda_{3}=\left\{\begin{array}[]{ll}\Lambda_{2}(f+g\Lambda_{2}+h\Lambda_{2}^{2})\,,&0\leq\Lambda_{2}\leq 10^{-2}\\ e^{a+b\xi+c\xi^{2}+d\xi^{3}+e\xi^{4}}\,,&\Lambda_{2}>10^{-2}\end{array}\right.,$$

$$\omega_{02}=\left\{\begin{array}[]{ll}f\,,&\xi<0\\ a+b\xi+c\xi^{2}+d\xi^{3}+e\xi^{4}\,,&1\leq\xi\leq 9\\ g\,,&\xi>9\end{array}\right.,$$

$$\omega_{03}=\left\{\begin{array}[]{ll}f\,,&\Upsilon<-1\\ a+b\Upsilon+c\Upsilon^{2}+d\Upsilon^{3}+e\Upsilon^{4}\,,&-1\leq\Upsilon\leq 10\\ g\,,&\Upsilon>10\end{array}\right.,$$

$$\Delta \ln (A)(Mf;\mathbf{x})$$

$$\Delta \Phi (Mf;\mathbf{x})$$

Full text

I Introduction

The detection of the merging binary neutron star (BNS) system GW170817 Abbott and et al. (2017) has demonstrated that the Advanced LIGO Harry and LIGO Scientific Collaboration (2010) and Advanced Virgo F. Acernese, et al. (2009) (Virgo Collaboration) gravitational-wave detectors can measure BNS properties such as the masses, spins, and tidal parameters The LIGO Scientific Collaboration et al. (2018a), and these results can even be extended to measure quantities such as the radii and equation of state (EOS) of NSs Raithel et al. (2018); De et al. (2018); The LIGO Scientific Collaboration et al. (2018b). Furthermore, an estimate of the BNS merger rate of 110–3840 events Gpc*-3* yr*-1* (90% confidence interval) O2C (2018) from this event indicates that, when second generation detectors reach design sensitivity, we may eventually observe a population of tens or hundreds of events after a couple of years of observation Abadie et al. (2010), and KAGRA Somiya (2012) and LIGO-India Iyer et al. (2011) could improve this further. These multiple events could then be “stacked” to significantly improve measurements of the EOS Del Pozzo et al. (2013); Lackey and Wade (2015), and also measure the population distribution of masses and spins.

The reliability of these measurements, however, depends on Bayesian analyses that require fast and accurate waveform models to match to the gravitational-wave data. Standard Bayesian parameter estimation tools such as Markov chain Monte Carlo (MCMC) and Nested Sampling often require $\mathcal{O}(10^{7})$ sequential waveform evaluations. Waveform evaluation times must therefore be significantly less than 1 s in order to run in less than a month for each event. To avoid large systematic errors, highly accurate waveform models are needed that include all relevant physical effects such as point-particle interactions, spins, and tides Favata (2014); Yagi and Yunes (2014); Wade et al. (2014). Most time-domain waveform implementations involve solving a set of ordinary differential equations, which depending on the model can take seconds to hours for long BNS waveforms. The results then need to be transformed into the frequency domain where the analysis takes place. Because of this, analytic approximations providing waveforms directly in the frequency-domain are often used instead.

For the analysis of GW170817 in Refs. The LIGO Scientific Collaboration et al. (2018a); O2C (2018), therefore, a set of four frequency-domain waveform models were used as templates for the main results. These models are modifications to spinning binary black hole (BBH) waveform models (see Table I of Ref. The LIGO Scientific Collaboration et al. (2018a) and references therein). They use either an aligned-spin post-Newtonian (PN) approximant in the Fourier domain (TaylorF2), the aligned-spin effective one body (EOB) formalism (SEOBNRT), or the phenomenological formalism for aligned-spin (PhenomDNRT) and precessing spin (PhenomPNRT). Tidal interactions were then included by simply adding a correction to the phase of each waveform. For the TaylorF2 waveform, the analytic 5PN and 6PN order terms were added. For the other three models, a fit that combined results from PN and numerical BNS simulations (referred to as the NRTidal fit in Ref. The LIGO Scientific Collaboration et al. (2018a)) was added Dietrich et al. (2017).

In this paper, we will use an implementation of the EOB formalism that treats matter effects consistently with the other parameters. It includes several dynamical effects beyond the standard adiabatic inspiral evolution and can also be tuned to numerical relativity (NR) simulations. The model, named SEOBNRv4T, includes the tidally induced $\ell=2$ and 3 multipole moments, the induced $\ell=2$ and 3 $f$-mode resonances Hinderer et al. (2016); Steinhoff et al. (2016), and the spin-induced quadrupole moment that can be important for large spins Poisson (1998); Harry and Hinderer (2018). This model agrees with NR simulations of BNS systems to a level consistent with the numerical error of the simulations (less than 1 rad during the last $\sim 10$ orbits before merger) Dietrich and Hinderer (2017); Kiuchi et al. (2017). It also includes a prescription for tapering the end of the waveform after merger in a manner consistent with numerical simulations. We also note there is an alternative EOB model, TEOBResumS Nagar et al. (2018a, b), that uses a different point-particle and spin prescription. For matter effects, this model includes the $\ell=2,3$, and 4 tidal terms and uses a method for re-summing these expressions. It also includes the spin induced-quadrupole term, but does not include the effect of dynamical tides from the induced $f$-mode resonances.

Although these EOB models are more expensive than frequency-domain models, parameter estimation results are still possible. In Refs. The LIGO Scientific Collaboration et al. (2018a); O2C (2018), an alternative parallelized code, RapidPE, was also used that approximately fits the posterior density function (PDF) for the intrinsic waveform parameters instead of directly sampling it Pankow et al. (2015); Lange et al. (2018). The measured tidal parameter using EOB waveforms and RapidPE were broadly consistent with the results using frequency-domain waveforms and traditional samplers (see Fig. 9 of Ref. O2C (2018)). However, a direct comparison using identical waveforms has not yet been published for BNS systems. As another possible solution for doing parameter estimation with EOB models, significant work has been done to optimize some EOB implementations, and speedups of a factor of several hundred Devine et al. (2016) have been achieved. Different formulations of the differential equations can also speed up integration Nagar and Rettegno (2018). However, many of these optimizations are specific to each waveform model, and it is not clear that sufficient speed ups can be achieved for the SEOBNRv4T model discussed here with dynamical tides. As an alternative to these methods, we will develop a surrogate model of SEOBNRv4T that bypasses the need for specialized parameter estimation tools and optimization techniques.

Surrogate modeling techniques have had significant success in gravitational-wave data analysis for rapidly evaluating waveform models. The essential idea is to construct a fitting function that approximates a waveform model as a function of frequency (or time) and waveform parameters $\mathbf{x}$. Previous works have focused on efficiently representing the space of all waveforms using a reduced basis of orthonormal functions. The two most common approaches for selecting a reduced basis have been singular value decomposition Pürrer (2014, 2016) and a greedy method Field et al. (2014); Lackey et al. (2017); Blackman et al. (2015, 2017a, 2017b). In this work, instead of using a reduced basis, we will instead use the fact that the Fourier-transformed EOB waveform is approximately known analytically via the fast, analytic TaylorF2 model. We then build a surrogate of the difference (or residual) between the EOB waveform and the TaylorF2 waveform. Because this residual is small, extremely high accuracy is not required, and we can efficiently represent the residual as a function of frequency using cubic splines between a small set of frequency nodes.

The final step in building a surrogate is to interpolate between waveform parameters $\mathbf{x}$. Many previous works focused on parameter spaces with three or fewer dimensions. This allows one to approximate the waveform as a function of $\mathbf{x}$ using standard interpolation techniques such as tensor spline Pürrer (2014, 2016) or Chebyshev interpolation Lackey et al. (2017). These interpolation techniques typically require waveforms to be evaluated on a rectangular grid, and thus suffer from the curse of dimensionality: for $N$ points per dimension $d$, $N^{d}$ waveforms are needed. For the EOB waveform here, we will only be able to reduce the parameter space to five dimensions, so grid based interpolation is unfeasible. Recent work on optimally choosing waveforms for analytic Doctor et al. (2017) and NR Blackman et al. (2017a, b) waveforms has shown that there are sufficiently accurate alternatives to grid-based interpolation that do not suffer from the curse of dimensionality.

In this paper we focus on a technique known as Gaussian process regression (GPR) Rasmussen and Williams (2006) that does not require a regular grid. Importantly, GPR also provides a convenient estimate of its own uncertainty. This allows us to iteratively add new points to the training set in such a way as to minimize the interpolation error over the parameter space for a given number of waveforms in the training set. Similar approaches have recently been used to optimally sample waveforms for a 2-dimensional aligned spin BBH surrogate using GPR Doctor et al. (2017) as well as for a surrogate of nonspinning, eccentric BBH mergers Huerta et al. (2017). GPR has also been used in GW data analysis to marginalize over waveform uncertainties in parameter estimation Moore and Gair (2014); Moore et al. (2016).

This work also takes a different approach from the surrogate developed in Ref. Lackey et al. (2017) for a nonspinning BNS EOB waveform. Whereas in Ref. Lackey et al. (2017) the authors constructed a time-domain surrogate with the goal of reproducing the original waveform model as accurately as possible, in this work, we construct a frequency-domain model that can be extended down to arbitrarily low frequencies and does not require an online Fourier transform for each waveform evaluation. It also enables additional techniques to accelerate likelihood evaluations such as reduced order quadrature Antil et al. (2013); Canizares et al. (2013, 2015); Smith et al. (2016), multi-banding Vinciguerra et al. (2017), and relative binning Zackay et al. (2018).

We organize the paper as follows. In Section II we provide an overview of the SEOBNRv4T waveform model and the approximations used to reduce the dimensionality of the parameter space. In Section III we describe the details of building a surrogate model, choosing training set waveforms, and evaluating the final model. We then compare the accuracy and speed of the surrogate to the original model in Section IV. In Section V we verify that the surrogate can correctly extract the parameters of injected waveforms. Finally, we discuss future improvements in Section VI and give the expressions for the TaylorF2 base model in the Appendix.

Conventions: Unless explicitly stated, we use units where $G=c=1$.

II Aligned-spin, dynamical tides EOB model

II.1 Inspiral-plunge waveform

The EOB approach to the general-relativistic 2-body problem, first described in Ref. Buonanno and Damour (1999), has proven successful in modeling the dynamics and GW emission of compact binaries. State-of-the-art aligned-spin EOB models Bohé et al. (2017); Nagar et al. (2017) can accurately match hundreds of NR simulations of aligned-spin BBH systems for mass ratios up to 8 and spin magnitudes up to 0.85 for unequal-mass (up to 0.98 for equal-mass) binaries. The EOB framework can also accommodate precessing-spin BBHs, showing good agreement to mildly precessing NR simulations Babak et al. (2017).

This research program has been extended to accommodate tidal effects Damour and Nagar (2010); Vines and Flanagan (2013); Damour et al. (2012); Bini et al. (2012); Bernuzzi et al. (2015); Hinderer et al. (2016); Steinhoff et al. (2016); Dietrich and Hinderer (2017) for binaries that contain NSs. The main effect of tides is to make the gravitational interaction more attractive with respect to the vacuum case. In particular, Refs. Hinderer et al. (2016); Steinhoff et al. (2016) built upon the aligned-spin BBH EOB model of Ref. Taracchini et al. (2014) and proposed a way to include the effect of dynamical tides. Neutron stars that are part of a compact-object binary will deform in the tidal field generated by the companion. The forcing tidal field varies at a multiple of the orbital frequency. Thus, in the late stages of the inspiral, the characteristic $f$-mode frequency of the neutron star can be dynamically approached, resulting in a resonant excitation of the $f$-mode. The net effect is an amplification of tidal effects as compared to the adiabatic limit, which assumes that the $f$-mode frequency is much larger than the frequency of the forcing tidal field.

Dynamical tidal effects are implemented in the EOB model through a modification of the potential $\Delta_{u}$, which is the $tt$-component of the metric of the effective spacetime. We adopt the tidally-augmented expression for $\Delta_{u}$ discussed in Appendix A of Ref. Steinhoff et al. (2016):

[TABLE]

where $\Delta_{u}^{\textrm{pm}}$ is the 4PN-accurate point-mass EOB term (Eq. (2.2) of Ref. Bohé et al. (2017)), and $\Delta_{u}^{\textrm{DT}}$ is the contribution due to dynamical tides. If either component of the binary is a black hole, then we set the tidal polarizabilities to zero. The tidal contribution, including quadrupolar and octupolar dynamical tides, is

[TABLE]

Here, $m_{A,B}$ are the masses of bodies $A$ and $B$, $M=m_{A}+m_{B}$ is the total mass, $X_{A,B}=m_{A,B}/M$, $u=1/r$ is the inverse of the ($M$-rescaled) EOB radial coordinate $r$, and $\Lambda_{\ell,\textrm{dyn}}^{A,B}(u)$ are the dimensionless $2^{\ell}$-polar dynamical tidal polarizabilities. Within the dynamical tides model, the tidal polarizabilities are not constant, but rather depend on the orbital separation and on the values of the $f$-mode frequencies $\hat{\omega}_{0\ell}^{A,B}$111Here, $\hat{\omega}_{0\ell}^{A,B}=m_{A,B}\omega_{0\ell}^{A,B}$, where $\omega_{0\ell}^{A,B}$ is the frequency in geometrized units.. In particular, the dimensionless dynamical tidal polarizability reads

[TABLE]

where $\Lambda_{\ell}^{A,B}$ is the dimensionless adiabatic tidal polarizability

[TABLE]

with $k^{A,B}_{\ell}$ the tidal Love number and $C_{A,B}=R_{A,B}/m_{A,B}$ the NS compactness, which depends on the NS radius $R_{A,B}$. Here, $\hat{k}_{\ell,\textrm{dyn}}(u;\hat{\omega}_{0\ell}^{A,B})$ is the separation-dependent, dimensionless enhancement factor Steinhoff et al. (2016); Dietrich and Hinderer (2017), which depends on the value of the $f$-mode $2^{\ell}$-pole dimensionless angular frequency, $\hat{\omega}_{0\ell}^{A,B}$.

Introducing a Keplerian orbital frequency defined from the radius variable as $M\Omega\equiv u^{3/2}$, using the notations (dropping the indices $A,B$)

[TABLE]

the dynamical enhancement factor of Steinhoff et al. (2016); Dietrich and Hinderer (2017) reads

[TABLE]

with $a_{\ell}$, $b_{\ell}$ constant coefficients with relevant values $\{a_{2},b_{2}\}=\{1/4,3/4\}$ and $\{a_{3},b_{3}\}=\{3/8,5/8\}$. These expressions are used in practice for $\ell=m$. Here, we have rewritten the first two resonant terms of Eq. (11) of Dietrich and Hinderer (2017) as

[TABLE]

Although the two terms are individually divergent222This requires to single out the near-resonance region in a numerical implementation. at the resonance $x=1$, the function $f(x)$ is actually regular, taking the value $f(1)=-1/12$. The last term in (6) above reads

[TABLE]

with $F_{S}$, $F_{C}$ the Fresnel sine and cosine functions. This term is regular at resonance $\hat{t}=0$. The enhancement factor (6) is 1 for low orbital frequencies, and increases to $\sim 2$ as $\ell$ times the orbital frequency approaches the $f$-mode frequency.

Tidal effects also enter the radiative part of the model. In particular, the point-mass $\ell=2,3$ inspiral-plunge waveform modes are corrected by tidal terms

[TABLE]

where the point-mass piece $h_{\ell m}^{\textrm{pm}}$ is discussed in Section II.B of Ref. Bohé et al. (2017) and the tidal piece $h_{\ell m}^{\textrm{tidal}}$ is given by Eqs. (A14)-(A17) of Ref. Damour et al. (2012). Following Ref. Dietrich and Hinderer (2017) in the computation of $h_{22}^{\textrm{tidal}}$, we include a dynamical enhancement factor that depends on the orbital separation and on the $f$-mode frequency, and is given in Eq. (15) of Ref. Dietrich and Hinderer (2017).

The waveform modes $h_{\ell m}^{\textrm{insp-plunge}}$ are then used to calculate the gravitational-wave flux (Eqs. (5) and (6) of Ref. Dietrich and Hinderer (2017)) from which the radiation-reaction force is derived. The Hamiltonian equations of motion with this radiation-reaction force are numerically integrated beginning with the quasicircular initial conditions described in Ref. Buonanno et al. (2006). The inspiral-plunge waveform is obtained by evaluating Eq. (9) using the solution to the orbital dynamics. We only include the dominant $(\ell,m)=(2,2)$ mode in constructing the final waveform, although we use all modes when evaluating the radiation reaction force.

II.2 Spin-induced quadrupole-monopole terms

For the model used in this paper, SEOBNRv4T, spin effects are based on the point-mass model SEOBNRv4 Bohé et al. (2017) that was calibrated to 141 NR simulations of BBH systems with aligned spins. It includes consistently spin-orbit interactions up to 3.5PN order and spin-spin interactions up to 2PN order. In addition, the SEOBNRv4T model accounts for the spin-induced quadrupole moment of neutron stars Poisson (1998), an extended-body effect distinguishing neutron stars from black holes. This effect is quadratic in the spins, appears at 2PN order as compared to the tidal effects that first appear at 5PN order, and is significant for systems with large spin Harry and Hinderer (2018). Although complete PN expressions for this contribution in the dynamics and waveform are known at 3PN Bohé et al. (2015)333And at 4PN for the dynamics only Levi and Steinhoff (2016)., at the moment the effect is only consistently included at the leading 2PN order. An extension to higher orders is left for future work444We note that TaylorF2 models do include all known spin-square terms at 3PN (see Appendix A). They are also present in the PhenomPNRT and SEOBNRT models. TEOBResumS was recently updated to include these next-to-leading terms Nagar et al. (2018b)..

In the compact binary system we consider, the two material bodies, e.g. two neutron stars, with masses and spins $m_{A}$ and $\bm{S}_{A}$ for $A=1,2$, have spin-induced mass quadrupole moments $Q_{A,\mathrm{SS}}^{ij}=-\kappa_{A}S_{A}^{\langle i}S_{A}^{j\rangle}/m_{A}$ Poisson (1998), with $\kappa_{A}$ the quadrupole-monopole parameter. The case of two black holes is recovered for $\kappa_{1}=\kappa_{2}=1$, while for a neutron star $\kappa$ can be larger, of order 10 for hard EOSs.

The SEOB model should be modified as follows in order to take this into account at leading order. The SEOB effective Hamiltonian is related to the real Hamiltonian by (7.2) of Barausse and Buonanno (2010). Its structure is given by (5.70) of Barausse and Buonanno (2010) as

[TABLE]

where

[TABLE]

in terms of the total mass $M$, the reduced mass $\mu=M\eta$, and the spin combination

[TABLE]

We refer to Barausse and Buonanno (2010) for the notations and meaning of the other terms. We find that, for systems including neutron stars, the latter contribution should be replaced with

[TABLE]

where we chose to use the spin variables $\bm{S}_{1}$ and $\bm{S}_{2}$ instead of $\bm{S}_{*}$ and $\bm{S}=\bm{S}_{1}+\bm{S}_{2}$ to avoid the occurence of the mass difference $m_{1}-m_{2}$ in denominators. This expression generalizes the black hole case and reduces to it for $\kappa_{1}=\kappa_{2}=1$.

In the waveform, at the leading 2PN order, the only contribution to consider will be in the mode $h_{22}$. The total leading-order spin-squared contribution, for aligned-spin circular orbits, is given by Kidder (1995); Will and Wiseman (1996); Buonanno et al. (2013)

[TABLE]

with $\Phi$ the orbital phase, $\omega=\dot{\Phi}$ the orbital frequency and $R$ the distance to the observer. Here, for a BBH,

[TABLE]

while for a BNS, this should be replaced with

[TABLE]

The total spin-squared contribution to the flux also includes squares of spin-orbit terms, and reads Kidder (1995); Will and Wiseman (1996); Poisson (1998)

[TABLE]

with $v=(M\omega)^{1/3}$.

When translating the above for the EOB factorized waveform (see e.g. (17) in Taracchini et al. (2012)), spin contributions in the effective source, tail factors, and phase contribution would enter only at higher order, so that the only modification to consider is in the $\rho_{\ell m}$ factor for $\ell=2$, $m=2$. We have simply

[TABLE]

with $\mathcal{A}$ the quantity given by (15) for a BBH or (16) for a BNS.

II.3 Waveform termination

We now look for a suitable time $t_{\rm match}$ to stop the numerical integration and match with an effective post-merger waveform. In the point-mass model, SEOBNRv4, nonquasicircular corrections to the inspiral-plunge signal guarantee that the waveform peaks, and $t_{\rm match}$ is chosen to be the time of the peak amplitude $t_{\rm peak}^{\rm amp}$. Here, although we include nonquasicircular corrections, they are computed as in the BBH case and are not tuned to BNS simulations, so we have less control over the behavior of Eq. (9) in the late inspiral. We therefore choose the following definition for $t_{\rm match}$. Let $t_{\rm peak}^{\rm amp}$ be the earliest time when the amplitude $|h_{22}^{\textrm{insp-plunge}}|$ of the (2, 2) mode peaks, and let $t_{\rm peak}^{\rm freq}$ be the time when the frequency $\omega_{22}^{\rm insp-plunge}$ of the (2, 2) mode peaks. Then,

[TABLE]

For some waveforms, the amplitude does not peak. In this case, we choose $t_{\rm peak}^{\rm amp}$ to be the earliest time when the slope of the amplitude $\partial_{t}|h_{22}^{\rm insp-plunge}|$ reaches a minimum after having reached a peak.

In BNS simulations, the waveform rapidly decreases in amplitude in 1–2 gravitational-wave cycles after reaching peak amplitude. The BNS system then either undergoes prompt collapse or post-merger oscillations of the remnant. We do not attempt to model post-merger oscillations, and instead model an approximate peak emission and a subsequenct tapering to zero. This signal will be represented in terms of analytic functions for the amplitude $A^{\textrm{post-mrg}}(t)$ and phase $\phi^{\textrm{post-mrg}}(t)$ such that $h_{22}^{\textrm{post-mrg}}(t)=A^{\textrm{post-mrg}}(t)\exp{[i\phi^{\textrm{post-mrg}}(t)]}$.

For the amplitude, we smoothly extend the waveform after $t_{\rm match}$ with a linear fit, and then taper the resulting amplitude. The linear extension is defined by

[TABLE]

where $a=|h_{22}^{\rm insp-plunge}(t_{\rm match})|$ and $b=\partial_{t}|h_{22}^{\rm insp-plunge}(t_{\rm match})|$. The tapering function is centered $15M$ after $t_{\rm match}$ and has a decay time $\tau=2\pi/\omega_{\rm match}$ (where $\omega_{\rm match}=\omega_{22}^{\rm insp-plunge}(t_{\rm match})$) of one gravitational wave period. It is given by

[TABLE]

The final amplitude after windowing is then given by $|h_{22}(t)|=\hat{A}(t)W(t)$.

For the phase, we smoothly extend the waveform frequency such that it agrees with the inspiral frequency at a time $t_{\rm freq}=t_{\rm match}-12M$ before the matching time $t_{\rm match}$, but then stretches out the frequency evolution such that it only approaches $\omega_{\rm match}$ asymptotically. We define this frequency evolution as

[TABLE]

where $\Delta\omega=\omega_{\rm match}-\omega_{\rm freq}$. Integrating the frequency, and requiring continuity at $t_{\rm freq}$, results in the final expression for the phase

[TABLE]

Although this effective post-merger model has not been fit to NR BNS simulations, it is in reasonable qualitative agreement with the post-merger behavior of the two equal-mass nonspinning NR BNS simulations that were analyzed in Ref. Hinderer et al. (2016). A more sophisticated, NR-informed model of the post-merger emission will be part of future investigations.

II.4 Reducing the number of matter parameters

After rescaling with the total mass $M$, the aligned-spin EOB waveform with dynamical tides depends on 13 intrinsic parameters: the mass ratio $q=m_{B}/m_{A}\leq 1$, the two dimensionless spin components along the orbital angular momentum $\chi_{A,B}$, the two spin-induced quadrupole-monopole parameters $\kappa_{A,B}$, the two adiabatic quadrupolar tidal polarizabilities $\Lambda_{2}^{A,B}$, the two adiabatic octupolar tidal polarizabilities $\Lambda_{3}^{A,B}$, and the four $\ell=2,3$ fundamental $f$-mode angular frequencies $\omega_{0\ell}^{A,B}$. To reduce the dimensionality of the intrinsic parameter space, we use nearly EOS-independent fits (universal relations) for these parameters in terms of $\Lambda_{2}$ as discussed below. This reduces the number of matter parameters from ten to two.

Eq. (15) of Ref. Yagi and Yunes (2017) provides a fit for $\kappa$ as a function of $\Lambda_{2}$. This relation has been fit for a sample of EOSs and for tidal parameters in the range $1\leq\Lambda_{2}\leq 10^{4}$. This formula, however, diverges at small values of $\Lambda_{2}$, so we replace it with a polynomial function that approaches the Kerr value of $\kappa=1$ at $\Lambda_{2}=0$. We require this extension to be continuous at $\Lambda_{2}=1$. The exact relation we use is

[TABLE]

where $\xi=\ln\Lambda_{2}$. The coefficients $\{a,b,c,d,e\}=\{0.194,0.09163,0.04812,-4.283\times 10^{-3},1.245\times 10^{-4}\}$ are from Table I of Ref. Yagi and Yunes (2017), and $\{f,g,h\}=\{0.42769,-0.32434,0.11074\}$.

Eq. (60) of Ref. Yagi (2014) provides a fit for $\Lambda_{3}$ as a function of $\Lambda_{2}$ with an accuracy of about 10% depending on the EOS. This relation has been fit for a sample of EOSs and for tidal parameters in the range $1\leq\Lambda_{2}\leq 10^{6}$. This formula, however, diverges at small values of $\Lambda_{2}$, so we replace it with a polynomial function that vanishes at $\Lambda_{2}=0$. We require this extension to be continuous at $\Lambda_{2}=10^{-2}$, and to fit the universal relation in the range $10^{-5}\leq\Lambda_{2}\leq 10^{-2}$. The exact relation we use is

[TABLE]

where $\xi=\ln\Lambda_{2}$, $\{a,b,c,d,e\}=\{-1.15,1.18,2.51\times 10^{-2},-1.31\times 10^{-3},2.52\times 10^{-5}\}$ from Ref. Yagi (2014), and $\{f,g,h\}=\{0.440649,-34.632322,1762.112913\}$.

Eq. (3.5) of Ref. Chan et al. (2014) gives relations for $\omega_{02}$ as a function of $\Lambda_{2}$ and $\omega_{03}$ as a function of $\Lambda_{3}$ to within a few percent error. For $\omega_{02}$, the fitting range used by Ref. Chan et al. (2014) was $0\leq\xi\leq 9$, and outside this range we require continuity. The relation we use is then

[TABLE]

where $\{a,b,c,d,e\}=\{0.182,-6.836\times 10^{-3},-4.196\times 10^{-3},5.215\times 10^{-4},-1.857\times 10^{-5}\}$ from Ref. Chan et al. (2014), and $\{f,g\}=\{0.182,0.161\}$. For $\omega_{03}$, the fit is given in terms of $\Upsilon=\ln\Lambda_{3}$ in the range $-1\leq\Upsilon\leq 10$, and outside this range we require continuity. The relation we use is then

[TABLE]

where $\{a,b,c,d,e\}=\{0.2245,-1.5\times 10^{-2},-1.412\times 10^{-3},1.832\times 10^{-4},-5.561\times 10^{-6}\}$ and $\{f,g\}=\{0.2379,0.1165\}$.

With these relations the waveform only depends on the five intrinsic parameters $\mathbf{x}=\{q,\chi_{A},\chi_{B},\Lambda_{2}^{A},\Lambda_{2}^{B}\}$. We note that this list does not include the total mass $M$ of the system. The point-mass part of the model is scale-invariant, and the tidal corrections only depend on $X_{A}=1/(1+q)$, $X_{B}=q/(1+q)$ and $\Lambda_{2}^{A,B}$. Finally, the simple model of the post-merger signal that we employ rescales with $M$ as well. It is not clear at what point in the transition from the inspiral to the post-merger this approximation breaks down. For example, whether the merging binary undergoes prompt collapse or forms a hypermassive remnant depends sensitively on the total mass. NR simulations will be needed to determine when after the merger this approximation is no longer valid.

III Surrogate model

In this section we describe how we decompose the EOB waveform into smooth, slowly-varying functions and train a surrogate model for these functions. We will work with the frequency-domain waveform $\tilde{h}(Mf;\mathbf{x})$. This is favorable for data analysis which is usually done in the frequency domain. It also allows us to extend the model down to arbitrarily low frequencies with an analytic, frequency-domain model. Unfortunately, Fourier transforming a finite length waveform leads to a surrogate with more noise. We will show below, however, that sufficient filtering can solve this problem. The SEOBNRv4T waveform model, as implemented, works for mass ratios in the range $q\in[1/3,1]$, spins $\chi_{1,2}\in[-0.5,0.5]$ and tidal parameters $\Lambda_{1,2}\in[0,5000]$. Our surrogate of this model, SEOBNRv4T_surrogate, will be valid for the same range of parameters.

III.1 Decomposition of the waveform

Because the waveform $\tilde{h}(Mf;\mathbf{x})$ is an oscillatory function of $Mf$ and $\mathbf{x}$, the waveform is usually decomposed into an amplitude $A(Mf;\mathbf{x})$ and phase $\Phi(Mf;\mathbf{x})$ as $\tilde{h}(Mf;\mathbf{x})=A(Mf;\mathbf{x})e^{i\Phi(Mf;\mathbf{x})}$. The amplitude and phase are smoother, mostly monotonic functions of frequency. This can be seen in Fig. 1 where we show the waveforms for the 32 corners of the 5-dimensional parameter space. Unfortunately, the amplitude and phase still span a wide range of values. The phase, for example, spans about $10^{4}$ rad between waveforms with different parameters $\mathbf{x}$ (see Fig. 1). To avoid systematic errors in the tidal parameters, for example, we need phase errors of $\lesssim 1$ rad over most of this frequency range, leading to a requirement on the fractional interpolation error of $\lesssim 10^{-4}$ rad. This is a difficult requirement to achieve for 5-dimensional interpolation.

For aligned-spin waveforms, we can solve this problem by using the fact that the waveform can be approximated with the analytic TaylorF2 waveform $\tilde{h}_{\rm F2}(Mf;\mathbf{x})=A_{\rm F2}(Mf;\mathbf{x})e^{i\Phi_{\rm F2}(Mf;\mathbf{x})}$. This allows us to write the EOB amplitude and phase in terms of small residuals, $\Delta\ln(A)(Mf;\mathbf{x})$ and $\Delta\Phi(Mf;\mathbf{x})$, relative to TaylorF2

[TABLE]

such that

[TABLE]

The exact functional form of the TaylorF2 waveform we use is given in Appendix A, and the residuals are shown in Fig. 1.

For the amplitude residual, we use a log-ratio instead of a ratio because it guarantees that interpolation errors will not lead to a negative amplitude for the reconstructed waveform (Eq. (35)). In addition, because the waveform amplitude spans several orders of magnitude at high frequencies, the log-ratio better captures this behavior. Comparing the phase $\Phi$ and phase residual $\Delta\Phi$ in Fig. 1, we find that the range in $\Delta\Phi$ is a few orders of magnitude smaller than the range in $\Phi$ except at very high frequencies.

Finally, we have found that the functions $\Delta\ln A(Mf;\mathbf{x})$ and $\Delta\Phi(Mf;\mathbf{x})$ sometimes vary rapidly as a function of the tidal parameters $\Lambda_{1}$ and $\Lambda_{2}$ for tidal parameters in the range [0, 1000]. We therefore perform a change of variables that stretches out the parameter space for small values of the tidal parameters:

[TABLE]

With $\mathbf{x}=\{q,\chi_{A},\chi_{B},\xi_{A},\xi_{B}\}$, the functions $\Delta\ln A(Mf;\mathbf{x})$ and $\Delta\Phi(Mf;\mathbf{x})$ are smoother functions of the parameters, making them easier to fit.

III.2 Conditioning the training set waveforms

The accuracy of the final frequency-domain surrogate depends on how well the finite-length, numerical waveform is Fourier transformed and filtered to remove numerical artifacts. We now describe the procedure to condition the training-set waveforms used to construct the surrogate.

We evaluate the EOB waveform with a starting frequency of $Mf_{{\rm win},i}=0.000197$, equivalent to a physical frequency of 20 Hz for a binary with total mass $M=2M_{\odot}$. In order to take a discrete Fourier transform, we window the start of the waveform with a Planck window McKechan et al. (2010) in the interval $[Mf_{{\rm win},i},Mf_{{\rm win},f}]=[0.000197,0.00021]$ to reduce Gibbs oscillations. The end of the waveform has zero amplitude, so the end does not need to be windowed. We then resample the waveform with a spacing $\Delta t/M=5$ and pad the end of the waveform with zeros such that all waveforms in the training set have the exact same time samples. After evaluating the discrete Fourier transform, we calculate the residuals $\Delta\ln(A)(Mf)$ and $\Delta\Phi(Mf)$ between the EOB and TaylorF2 waveforms using Eqs. (33) and (34).

Waveforms have free time and phase parameters $t_{c}$ and $\phi_{c}$, and in the frequency-domain this means that one can freely add a linear term $\phi_{c}+2\pi(Mf)(t_{c}/M)$ to the phase $\Phi(Mf)$. We use this freedom to match the EOB waveform to the analytic TaylorF2 waveform near the starting frequency. We do this by subtracting a linear fit to $\Delta\Phi(Mf)$ at the beginning of the waveform in the window $[Mf_{{\rm fit},i},Mf_{{\rm fit},f}]=0.00021[1,1.05]$. At $Mf=Mf_{{\rm fit},i}$, the resulting phase residual $\Delta\Phi(Mf)$ is zero and has zero slope, guaranteeing that the surrogate smoothly matches to TaylorF2 below this frequency.

The resulting waveforms still have some remaining Gibbs oscillations which can be seen as small-amplitude, high-frequency oscillations in the residuals $\Delta\ln(A)(Mf)$ and $\Delta\Phi(Mf)$. These come from two sources. The first is the fact that the Planck window at the beginning of the waveform was not sufficiently long. We could make this window longer, but that would require us to start the surrogate at a higher frequency. The second source comes from the end of the waveform where the amplitude rapidly drops to zero amplitude during the $\sim 1$ cycle after the peak amplitude. Some contribution to the Gibbs oscillations at high frequencies is therefore a genuine feature of the EOB model. We reduce these oscillations in the frequency domain using a moving average filter centered on $Mf$ with interval $[Mf(1-\alpha),Mf(1+\alpha)]$. We use a width of $\alpha=0.1$ for the amplitude residual and $\alpha=0.05$ for the phase residual. Smoothing these oscillations makes it significantly easier to interpolate the amplitude and phase residuals as functions of $Mf$ and $\mathbf{x}$.

Finally, we truncate the residuals outside the interval $[Mf_{{\rm trunc},i},Mf_{{\rm trunc},f}]=[0.00021,0.07]$. We note that the gravitational-wave frequency at the Schwarzschild innermost stable circular orbit (ISCO) is $Mf_{\rm ISCO}=1/(6^{3/2}\pi)\approx 0.022$, and higher frequencies represent only the last $\sim 1$ cycle of the waveform during the merger. Thus, the high-frequency cutoff of $Mf_{{\rm trunc},f}=0.07$ is a very conservative upper bound for data-analysis purposes. However, truncating the waveform at lower frequencies can visibly effect the shape of some waveforms if inverse Fourier transformed back into the time domain.

To validate that the waveforms are sufficiently conditioned, we plot cross-sections of the residuals in Fig. 2 for fixed frequencies as functions of one of the waveform parameters $\chi_{1}$. The fact that the residuals are smooth functions of the waveform parameters indicates that most numerical noise has been removed. The exception is the noisy amplitude residual at high frequencies, $\Delta\ln(A)(Mf=0.052;\mathbf{x})$, where the amplitude is very small and dominated by numerical noise.

III.3 Spline interpolation for frequency $f$

We now seek to interpolate the conditioned residuals $\Delta\ln(A)(Mf;\mathbf{x})$ and $\Delta\Phi(Mf;\mathbf{x})$. We begin by choosing a method for interpolating as a function of $Mf$ for fixed $\mathbf{x}$. The majority of previous papers on gravitational-wave surrogates have used an orthonormal basis of global functions $\hat{e}_{i}(Mf)$ for interpolating as a function of frequency (or time) Pürrer (2014, 2016); Field et al. (2014); Lackey et al. (2017); Blackman et al. (2015, 2017a, 2017b). This type of surrogate, built from a reduced basis, is usually referred to as a reduced order model. For a generic function $g(Mf;\mathbf{x})$, this decomposition can be written as $g(Mf;\mathbf{x})\approx\sum_{i=1}^{N}c_{i}(\mathbf{x})\hat{e}_{i}(Mf)$. The coefficients $c_{i}(\mathbf{x})$ are then interpolated as functions of $\mathbf{x}$ Pürrer (2014, 2016). Alternatively, using the empirical interpolation method Barrault et al. (2004); Chaturantabut and Sorensen (2010); Field et al. (2014), one can re-express these basis functions $\hat{e}_{i}(Mf)$ in terms of empirical interpolating functions $B_{j}(Mf)$ and the value of $g(Mf;\mathbf{x})$ at empirical nodes $MF_{j}$ as $g(Mf;\mathbf{x})\approx\sum_{j=1}^{N}B_{j}(Mf)g(MF_{j};\mathbf{x})$. The location of these nodes $MF_{j}$ is then optimized to minimize interpolation errors.

For the problem here, we have found that global basis functions are not necessarily the optimal solution. The residuals are small and smooth at low frequencies and large and noisy at high frequencies. If using the empirical interpolation method with global basis functions, for example, errors in evaluating $g(MF_{j};\mathbf{x})$ at high frequencies can propagate to large errors between the nodes $MF_{j}$ at low frequencies. Instead, we find that spline interpolation works significantly better. We use 40 frequency nodes $MF_{j}$ log-spaced in the interval $[Mf_{{\rm trunc},i},Mf_{{\rm trunc},f}]$ (see Fig. 1). We evaluate the residuals $\Delta\ln(A)(MF_{j};\mathbf{x})$ and $\Delta\Phi(MF_{j};\mathbf{x})$ at these nodes as discussed below. We then interpolate between these frequencies using third-order splines. The local, third-order polynomials, that are only connected by the requirement of smoothness, do not propagate high-frequency errors down to low-frequency errors as significantly as do global basis functions.

III.4 Gaussian process regression for parameters $\mathbf{x}$

Next we choose a method to interpolate the residuals $\Delta\ln(A)(MF_{j};\mathbf{x})$ and $\Delta\Phi(MF_{j};\mathbf{x})$ at each of the frequency nodes $MF_{j}$ as a function of the five waveform parameters $\mathbf{x}$. Most multivariate interpolation techniques (e.g. tensor spline or Chebyshev interpolation) require a function to be sampled on a rectangular grid, and thus suffer from the curse of dimensionality: the number of samples grows exponentially with the dimension $d$ ($N^{d}$ samples for $N$ points per dimension). For our 5-dimensional problem, $10^{5}$ waveform evaluations are needed for only 10 samples per parameter. At a starting frequency of 20 Hz, over an hour is needed to evaluate a single SEOBNRv4T waveform on a standard CPU, so this is at the limit of what is reasonable. However, there are other methods that do not require a rectangular grid, and this allows us to choose more efficient experimental designs as discussed in Sec. III.6 below. The method we choose is GPR Rasmussen and Williams (2006).

In GPR, the values of a function $g(\mathbf{x})$ at the points $\mathbf{x}$ are assumed to be a realization of a Gaussian process

[TABLE]

where $m(\mathbf{x})$ is a mean function at each point $\mathbf{x}$ and $k(\mathbf{x},\mathbf{x}^{\prime})$ is a covariance function between the points $\mathbf{x}$ and $\mathbf{x}^{\prime}$. When modeling data with GPR, one often attempts to subtract out the mean by first fitting the data with a parameterized function. Then, the residual is represented by a zero-mean Gaussian process. The covariance $k(\mathbf{x},\mathbf{x}^{\prime})$ is described by a kernel function with tunable hyperparameters. For the problem here, we have already subtracted the TaylorF2 waveform from the EOB waveform, so we model the residual in terms of a zero-mean Gaussian process.

Explicitly, we are interested in predicting the function value $y_{*}$ at the point $\mathbf{x}_{*}$ given the sampled function values $y_{i}$ at the points $\mathbf{x}_{i}$. In a zero-mean Gaussian process, this is represented by the following multivariate Gaussian distribution

[TABLE]

where $K_{ij}=k(\mathbf{x}_{i},\mathbf{x}_{j})$ is a matrix, $K_{*i}=k(\mathbf{x}_{i},\mathbf{x}_{*})$ is a vector, and $K_{**}=k(\mathbf{x}_{*},\mathbf{x}_{*})$ is a scalar.

From Eq. (38), the conditional probability for $y_{*}$ given the training set examples $y_{i}$ and kernel hyperparameters $\bm{\theta}$ is also a Gaussian

[TABLE]

where the mean and variance are

[TABLE]

Eq. (40) is the estimate of the function, and Eq. (41) is the estimate of the uncertainty.

We use radial kernels, which express the covariance in terms of a distance $r$ between points

[TABLE]

where we choose the matrix $M$ to be diagonal

[TABLE]

The tunable hyperparameters ${\bf\ell}$ represent the length scale over which the function $g(\mathbf{x})$ varies in each coordinate.

We try two classes of radial kernels. The first is the Matérn kernel

[TABLE]

where $K_{\nu}(x)$ is a modified Bessel function. The value of $\nu$ parameterizes the smoothness of the Gaussian process, and is $k$ times mean-square differentiable if $\nu>k$ Rasmussen and Williams (2006). For half-integer values of $\nu$, this kernel has a computationally cheap form without special functions, and we have had good results with $\nu=5/2$, resulting in a twice-differentiable function. The $\nu=5/2$ kernel is

[TABLE]

The second class we try is the more common squared exponential kernel that results in an infinitely differentiable function Rasmussen and Williams (2006)

[TABLE]

Our final kernel takes the form

[TABLE]

where $\sigma_{f}$ is a scale factor that describes the range of values that $g(\mathbf{x})$ takes over the domain, and $\sigma_{n}$ is a noise parameter. The white noise kernel $\sigma_{n}^{2}\delta_{ij}$ (also called a nugget) parameterizes the noise in the data $y_{i}$. In our case, the training set waveforms have numerical noise that we will estimate by optimizing the hyperparameters. The full set of hyperparameters is now $\bm{\theta}=\{\sigma_{f},\ell_{q},\ell_{\chi_{1}},\ell_{\chi_{2}},\ell_{\xi_{1}},\ell_{\xi_{2}},\sigma_{n}\}$ for our parameter space $\mathbf{x}=\{q,\chi_{1},\chi_{2},\xi_{1},\xi_{2}\}$.

In order to estimate the hyperparameters, we use the above assumption (Eq. (38)) that the joint distribution of the data $y_{i}$ is a multivariate Gaussian with the following distribution:

[TABLE]

This is the log-likelihood for $y_{i}$ given the hyperparameters ${\bm{\theta}}$, and we can find the posterior for ${\bm{\theta}}$ given $y_{i}$ using Bayes’ theorem

[TABLE]

The prior $p({\bm{\theta}})$ is typically uniform and used to set the bounds on ${\bm{\theta}}$. One can sample this posterior if interested in the distribution of hyperparameters $\bm{\theta}$. However, for the problem here, we simply want the maximum posterior. We do this using the gaussian_process module in the scikit-learn package Pedregosa et al. (2011). With these optimized hyperparameters, the final interpolating function is given by Eq. (40) and its uncertainty by Eq. (41).

III.5 Surrogate waveform evaluation

Given the GPR fits for the amplitude and phase residuals, we can now reconstruct the frequency-domain waveform. We label the set of GPR fits at the frequency nodes $MF_{j}$ by $\{\mathcal{I}_{\rm GPR}[\Delta\ln A_{j}](\mathbf{x})\}$ for the amplitude residuals and by $\{\mathcal{I}_{\rm GPR}[\Delta\Phi_{j}](\mathbf{x})\}$ for the phase residuals. The surrogates for the residuals $\Delta\ln A_{S}(Mf;\mathbf{x})$ and $\Delta\Phi_{S}(Mf;\mathbf{x})$ are constructed by interpolating between the nodes $MF_{j}$ with cubic splines:

[TABLE]

We set these functions to 0 below the first frequency node $MF_{0}=Mf_{{\rm trunc},i}$ so that the waveform transitions to TaylorF2 at lower frequencies. With the analytic expressions for $A_{\rm F2}(Mf;\mathbf{x})$ and $\Phi_{\rm F2}(Mf;\mathbf{x})$ and the interpolated expressions $\Delta\ln A_{S}(Mf;\mathbf{x})$ and $\Delta\Phi_{S}(Mf;\mathbf{x})$, the final surrogates for the amplitude and phase are

[TABLE]

In physical units, for an inclination angle $\iota$, the $+$ and $\times$ polarizations of the waveform are (for $f>0$)

[TABLE]

III.6 Iterative construction of training set and surrogate

One of the main aims of this paper is to build a surrogate with as few waveform evaluations as possible. Using the freedom provided by GPR to sample waveforms at arbitrary locations, we try a set of designs more efficient than uniform grids. One such design is a Latin Hypercube Design (LHD) McKay et al. (1979). An LHD with $N$ samples divides each of $d$ dimensions uniformly into $N$ grid points for a total of $N^{d}$ grid points. However, unlike a uniform grid, the $N$ values in each dimension are sampled exactly once instead of $N^{d-1}$ times. For an LHD there are $(N!)^{d}$ ways to choose these points, and we choose one randomly555There are additional ways to choose an LHD. A standard requirement is that the LHD be space-filling, meaning that the points are as far apart as possible from each other. (See Ref. Husslage et al. (2011) for a review of methods for optimizing the placement of samples.) One such definition of space filling is that the chosen locations maximize the minimum Euclidean distance between any two samples. We have not experimented with these alternatives here.. An LHD has the property that the samples are non-collapsing; a projection onto a subspace is still an LHD and no points are repeated in any dimension. This avoids wasting samples when one of the parameters has much less of an influence than the other parameters.

We build an initial training set with 128 waveforms sampled with an LHD for the five parameters $\mathbf{x}=\{q,\chi_{1},\chi_{2},\xi_{1},\xi_{2}\}$. In addition, we find empirically that the GPR uncertainty estimates (Eq. (41)) are largest at the corners of the parameter space, so we also sample the 32 corners. We construct our initial surrogate with these 160 waveforms shown in Fig. 3.

With this initial surrogate, we compare two methods for adding samples and improving the accuracy of the surrogate. The first is to simply sample the parameter space $\mathbf{x}$ with a uniform distribution. We choose 400 samples shown in Fig. 3. The second method is sometimes referred to as uncertainty sampling Brochu et al. (2010). In this method, we can choose new training-set samples by iteratively searching the parameter space for new points ${\bm{x}}$ that maximize some error criterion, then adding a new sample at that point. The quantity that we use is the root-mean-squared (RMS) phase error at the frequency nodes $MF_{j}$:

[TABLE]

where $\sigma_{\Delta\Phi_{j}}(\mathbf{x})$ is the GPR estimate (Eq. (41)) of the uncertainty in $\Delta\Phi(MF_{j},\mathbf{x})$ at node $MF_{j}$. We only include nodes below $Mf=0.03$ because there is very little signal at higher frequencies, and the RMS error would otherwise be dominated by higher frequencies. There are many alternative ways to construct a scalar that estimates overall waveform error. One alternative would be an approximation to the mismatch discussed in Sec. IV.1. However, the quantity $\epsilon_{\rm RMS}(\mathbf{x})$ describes the most important quantity, phase, and is fast to compute. For a training set with $N$ samples, $\epsilon_{\rm RMS}(\mathbf{x})$ has $\sim N$ local maxima located in the voids between samples $\mathbf{x}$. To find the global maximum, we use a basin hopping algorithm Wales and Doye (1998); sci to avoid getting stuck in local maxima.

For low starting frequencies, EOB waveforms are expensive enough that we would like a method to efficiently choose $N_{\rm new}$ waveforms and evaluate them in parallel. We note that if we hold the hyperparameters $\bm{\theta}$ fixed, the GPR error estimate (Eq. (41)), and therefore $\epsilon_{\rm RMS}(\mathbf{x})$, only depends on the samples $\mathbf{x}$ and not on the waveform data. The algorithm for choosing the $N_{\rm new}$ new points is as follows.

For $i=1,\dots,N_{\rm new}$:

Construct the GPR error estimate $\epsilon_{\rm RMS}(\mathbf{x})$ (Eq. (56)). In practice this can be done by specifying the samples $\mathbf{x}$, the hyperparameters $\bm{\theta}$ from the initial surrogate, and dummy data for $\mathbf{y}$ since Eq. (56) does not depend on $\mathbf{y}$. 2. 2.

Find the point $\mathbf{x}_{\rm max}$ that maximizes $\epsilon_{\rm RMS}(\mathbf{x})$ over the parameter space. 3. 3.

Add $\mathbf{x}_{\rm max}$ to the list of samples: ${\bm{x}}\to[{\bm{x}},{\bm{x}}_{\rm max}]$.

These $N_{\rm new}$ waveforms can now be evaluated in parallel. An updated surrogate can then be constructed with re-optimized hyperparameters $\bm{\theta}$ using the $N+N_{\rm new}$ waveforms.

Fig. 3 shows the parameters of the $N_{\rm new}=400$ new waveforms chosen by the uncertainty sampling method when using a Matérn kernel. We note that the edges of the parameter space are more often chosen than the inner region. We find this is true for both the Matérn and squared exponential kernels. Ref. Doctor et al. (2017) found similar behavior for 1- and 2-dimensional problems. As discussed there, this results because the edges have fewer nearby samples than the interior, resulting in a larger error estimate. The uncertainty sampling algorithm compensates by adding samples near the edges where the error estimate is largest.

Fig. 4 shows the estimated RMS error $\epsilon_{\rm RMS}(\mathbf{x})$ maximized over the parameters $\mathbf{x}$ for each new sample added to the training set. This is done for both the Matérn (black) and squared exponential (gray) kernels. For each kernel, we constructed an initial surrogate with optimized hyperparameters, then fixed the hyperparameters before choosing the 400 new samples using uncertainty sampling. The non-smoothness of the curves results because the basin-hopping algorithm does not always find the global maximum of $\epsilon_{\rm RMS}$ at each iteration. We find that increasing the number of waveforms from 160 to 560 decreases the global maximum of $\epsilon_{\rm RMS}$ by a factor of $\sim 4$ for the Matérn kernel and a factor of $\sim 10$ for the squared exponential kernel. In Section IV.1 we will examine how well this estimated error agrees with the true error.

Finally, after generating waveforms for the initial surrogate, the uncertainty sampling, and the uniform distribution, we combine all 960 waveforms. We then build the final surrogate with re-optimized hyperparameters.

IV Performance of the surrogate

IV.1 Accuracy

The accuracy of the surrogate can be assessed by comparing it to a test set of waveforms. We generate 1000 waveforms with the original parameters sampled uniformly in the ranges $q\in[1/3,1]$, $\chi_{1,2}\in[-0.5,0.5]$ and $\Lambda_{1,2}\in[0,5000]$. In Fig. 4, we show box plots of the RMS phase error of the test set waveforms for each of the 8 surrogate models. The first two surrogates are constructed using the initial design with either the Matérn or squared exponential kernel. The second two surrogates are constructed with the initial design and the 400 uniformly distributed samples and either the Matérn or squared exponential kernel. The next two surrogates use the initial design and the 400 samples chosen by the uncertainty sampling algorithm for each kernel. The final two surrogates use all 960 waveforms for each kernel. We find that for the Matérn kernel, the global maximum of $\epsilon_{\rm RMS}(\mathbf{x})$ is a reasonable estimate of the maximum test-set RMS phase errors. However, for the squared exponential kernel, the global maximum of $\epsilon_{\rm RMS}(\mathbf{x})$ significantly underestimates the true maximum test-set error. Additionally, we find that the surrogates constructed with the uniform distribution have smaller median test-set errors. However, the surrogates constructed with uncertainty sampling have smaller maximum test-set errors.

Another common measure of the surrogate model accuracy is the mismatch between the surrogate model and the test-set waveforms. The mismatch represents the loss in signal-to-noise ratio that results from using the surrogate model $h_{\rm Sur}$ instead of the original waveform $h$. It is defined as the deviation from a perfect overlap after aligning the two waveforms using the time and phase free parameters $t_{c}$ and $\phi_{c}$:

[TABLE]

The inner product between waveforms $h_{1}$ and $h_{2}$ is

[TABLE]

where the Fourier transformed waveforms are weighted by the noise power spectral density (PSD) $S_{n}(f)$ of the detector.

In Fig. 5, we show box plots of the mismatches $\mathcal{M}$ between the surrogate and the 1000 test-set waveforms for each of our 8 surrogates. We use the design sensitivity aLIGO PSD Aasi et al. (2013) and a total mass of $M=2.8M_{\odot}$. Our integration bounds are $Mf_{\rm low}=0.00021$ and $Mf_{\rm high}=0.07$, corresponding to physical units of 15.2 Hz and 5076 Hz respectively. As we increase the number of samples in our training set, the mismatches decrease. However, it is less clear which kernel and sampling strategy is optimal. For our final surrogate, we simply choose the surrogate that has the smallest maximum and median mismatches, and this is the surrogate with the Matérn kernel and all 960 waveforms. This surrogate has a maximum mismatch of $4.5\times 10^{-4}$.

We also show the surrogate amplitude and phase errors as a function of frequency for the parameters with the largest mismatch relative to the test-set in Fig. 6. As we add waveforms to the training set with the various sampling methods, the errors generally decrease. However, above a frequency of $Mf\sim 0.03$, the amplitude of the training set waveforms are noisy (see Figs. 1 and 2), and it is thus difficult to produce a good fit. This is not a problem, however, because aLIGO is insensitive to the waveform morphology at such high frequencies. With our final surrogate, we are able to achieve a maximum phase error of $\lesssim 1$ radian over the entire frequency range relevant to aLIGO.

Finally, to demonstrate that we can correctly recover the behavior of the original, time-domain EOB waveform, we inverse Fourier transform the surrogate and compare it to the time-domain waveform. In Fig. 7 we compare the final surrogate to the test set waveform that has the largest mismatch. We align the two waveforms in time and phase by maximizing the overlap at early times in the interval $(t-t_{\rm merger})/M\in[-10^{7},-10^{6}]$ using the method of Read et al. (2009a). The two waveforms remain nearly identical up to merger, far from the interval where they were aligned.

IV.2 Timing

The surrogate evaluation time is dominated by two parts. The first part is the time needed to evaluate the residuals $\{\mathcal{I}_{\rm GPR}[\Delta\ln A_{j}](\mathbf{x})\}$ and $\{\mathcal{I}_{\rm GPR}[\Delta\Phi_{j}](\mathbf{x})\}$ at each of the $N_{A}+N_{\Phi}$ interpolating nodes using GPR. Although optimizing the hyperparameters for GPR scales with the number of training set samples $N$ as $\mathcal{O}(N^{3})$, the evaluation time of a stored GPR scales as $\mathcal{O}(N)$ (see Eq. (40)). The evaluation time for all the residuals therefore has a cost of $\mathcal{O}(N(N_{A}+N_{\Phi}))$. The second part is the time needed to resample the final surrogate (Eqs. (52)– (III.5)) at uniformly spaced frequency samples in physical units beginning at a starting frequency $f_{\rm min}$.

In Fig. 8, we show the surrogate evaluation time as a function of the starting frequency $f_{\rm min}$ for an equal mass $1.4M_{\odot}$–$1.4M_{\odot}$ binary evaluated on a 3.5 GHz Intel Xeon processor. The waveform was sampled with a uniform frequency spacing corresponding to a sampling rate of 4096 Hz. At large starting frequencies the evaluation time is limited by the $\sim 0.01$ s needed to calculate the GPR fits at each frequency node. As $f_{\rm min}$ decreases, the waveform becomes longer and a smaller frequency spacing is needed. The evaluation time is then dominated by the spline interpolation needed to resample the waveform. We also compare the evaluation time to three BNS waveform models used recently in the analysis of GW170817 (TaylorF2, SEOBNRv4_ROM_NRTidal, IMRPhenomD_NRTidal) The LIGO Scientific Collaboration et al. (2018a) and the BBH model SEOBNRv4_ROM. Only the analytic TaylorF2 model is noticeably faster. In the bottom panel, we compare the evaluation time to the original time-domain SEOBNRv4T waveform model and find it is $\sim 1500$–6000 times faster.

V Parameter estimation

As a final test of the surrogate, we perform several parameter estimation runs where we inject the original, time-domain SEOBNRv4T waveform and recover the injected parameters with the SEOBNRv4T_surrogate template. We use the design sensitivity PSDs for the two advanced LIGO detectors and the advanced Virgo detector Aasi et al. (2013). We choose as our waveform parameters $m_{1}=m_{2}=1.4M_{\odot}$, $\chi_{1}=\chi_{2}=0.1$, and $\Lambda_{1}=\Lambda_{2}=1286$ corresponding to the EOS named MS1b in Read et al. (2009b). The sky location and inclination angle are fixed at an arbitrary value, and the distances are chosen such that the optimal network signal-to-noise ratios (SNRs) are either 30, 60, or 120. We inject the waveforms into zero noise data so that the recovered parameters do not depend on the specific noise realization Nissanke et al. (2010).

As discussed in more detail in The LIGO Scientific Collaboration et al. (2018a), the posterior $p(\vec{\vartheta}|d)$ for the parameters $\vec{\vartheta}$ given the detector data $\vec{d}$ is given by Bayes’ theorem

[TABLE]

where $p(\vec{\vartheta})$ is the prior and $\mathcal{L}(\vec{d}|\vec{\vartheta})$ is the likelihood. For the prior we choose the sky position and orientation to be uniform on the unit sphere; the distance to be uniform in co-moving volume; the masses, spins, and tidal parameters to be uniform in the ranges $m_{1,2}\in[0.7,2]M_{\odot}$, $\chi_{1,2}\in[-0.5,0.5]$ and $\Lambda_{1,2}\in[0,5000]$; and the mass ratio is restricted to the range $q\in[0.5,1]$. As was done in The LIGO Scientific Collaboration et al. (2018a), we sample the data and waveform at 4096 Hz such that the Nyquist frequency is 2048 Hz, and integrate the likelihood function in the interval $[20,2048]$ Hz. We sample this posterior with the parallel-tempered MCMC code in LALInference LIGO Scientific Collaboration, Virgo Collaboration (2018) which requires $\sim 10^{7}$ iterations to produce 10,000–20,000 statistically independent samples.

We note that the SEOBNRv4T waveform model still has some signal above 2048 Hz (see Fig. 1). We therefore also performed runs with a sampling frequency of 8192 Hz and integrated the likelihood function in the interval $[20,4096]$ Hz. However, this did not have a noticeable effect on the recovered posteriors.

We show marginalized 1-dimensional posteriors for four key parameters in Fig. 9. These are the chirp mass $\mathcal{M}=(m_{1}m_{2})^{3/5}/(m_{1}+m_{2})^{1/5}$, the mass ratio $q=m_{2}/m_{1}$, the effective spin $\chi_{\rm eff}=(m_{1}\chi_{1}+m_{2}\chi_{2})/(m_{1}+m_{2})$ and the effective tidal deformability parameter $\tilde{\Lambda}=(16/13)[(1+12q)\Lambda_{1}+(12+q)q^{4}\Lambda_{2}]/(1+q)^{5}$. As the SNR increases, the peaks of the marginalized 1-d posteriors become more closely aligned with the injected value. However, the distributions are noticeably asymmetric and are not centered on the injected values. This results mainly because prior boundaries on some parameters, such as the mass ratio $q$, shift the marginalized PDFs for correlated parameters. This effect is strongest for the interaction between $\chi_{\rm eff}$ and $q$ (Fig. 10), where the injected value $q=1$ is at the boundary of the prior. Although the injected parameter $(\chi_{\rm eff},q)=(0.1,1)$ lies on the ridge of maximum density in the marginalized 2-d PDF, the marginalized 1-d PDF for $\chi_{\rm eff}$ is noticeably offset. To counteract this effect in Fig. 9, we also take a thin cross section containing the injected value of $q=1$ by cutting samples outside the interval $q\in[0.9,1]$. This causes the peaks of the other 1-d distributions to be much closer to the injected values. We consider the difference between these peaks and the injected values to be a conservative bound on the size of the systematic errors due to errors in the surrogate model.

VI Discussion and future work

We have constructed a fast, frequency-domain surrogate of one of the most accurate BNS waveform models to date. This aligned-spin model, SEOBNRv4T, incorporates the tidally induced $\ell=2$ and 3 multipole moments as well as the effect of dynamical tides as the excitation approaches the $\ell=2$ and 3 $f$-mode frequencies. We have achieved mismatches of no more than $4.5\times 10^{-4}$ and phase errors of $\lesssim 1$ rad up to the merger frequency. These are sufficient to not bias results in any of the parameters.

The evaluation time has a flat cost of $\sim 0.01$ s to perform the GPR interpolation at each node. The rest of the time is spent resampling the waveform with spline interpolation. For a starting frequency of 20 Hz, this takes 0.13 s when the waveform is matched with data uniformly sampled at 4096Hz. For current parameter estimation codes, this is sufficient. However, one could further improve run times using reduced order quadrature Antil et al. (2013); Canizares et al. (2013, 2015); Smith et al. (2016) which requires frequency-domain waveforms, multi-band waveform interpolation Vinciguerra et al. (2017), or relative binning Zackay et al. (2018). Finally, we note that the production MCMC sampler used for the GW170817 analysis has an autocorrelation length of $\mathcal{O}(10^{3})$ for aligned-spin BNS models, so significant improvements to the parameter estimation runtime can be made through better samplers.

The hierarchical method presented here, where we begin with the analytic TaylorF2 reference model then make a surrogate of the residual, can be used to further improve the waveform model. For example, one could build a surrogate for numerical BNS simulations using SEOBNRv4T_surrogate as the base model, then constructing a surrogate of the residual. Such a model would have the accuracy of EOB below $\sim 400$ Hz and the accuracy of numerical simulations for the last several cycles before merger. Current state-of-the-art NR simulations have phase errors of several tenths of a radian over the last $\sim 20$ gravitational-wave cycles Dietrich and Hinderer (2017); Kiuchi et al. (2017). Using GPR and uncertainty sampling discussed above, one could optimally choose the waveform parameters for the numerical simulations, and run them in parallel to build the training set. If the difference between SEOBNRv4T and numerical BNS simulations is small, one would not need a high fractional accuracy for a surrogate of the difference, and 10–100 waveforms may be sufficient.

This model notably does not include precession for non-aligned spins. Although none of the EOB models currently available include both tidal effects and precession, there are ways to rotate this waveform model to a precessing frame and approximately incorporate precession. For example, Chatziioannou et al. have analytically solved the 2PN-accurate precession equations for generic spins Chatziioannou et al. (2017a, b), with the exception of transitional precession which is unlikely for low mass-ratio and spin BNS systems. Using shifted uniform asymptotics Klein et al. (2014), they can also analytically Fourier-transform the solution. Importantly, one is free to specify the frequency-domain amplitude and phase evolution, such as the SEOBNRv4T_surrogate here, for the waveform in the co-precessing frame. This approach would provide a fast, accurate model for BNS systems with tides and generic spins. We note, however, that including precession did not noticeably affect the parameter estimation results for GW170817 The LIGO Scientific Collaboration et al. (2018a).

Acknowledgements.

The authors sincerely thank Justin Vines and Tanja Hinderer for deriving the quadrupole-monopole corrections needed in the model. The authors are grateful to Zoheyr Doctor, Bhooshan Gadre, and Prayush Kumar for reviewing and testing the implementation of SEOBNRv4T_surrogate in LAL. BL also thanks the participants of the Bayesian Methods in Nuclear Physics workshop at the University of Washington for helpful advice at the start of this work.

Appendix A TaylorF2 reference waveform

The explicit expressions for the amplitude and phase of the TaylorF2 waveform are as follows. For the amplitude we use the 1PN correction to the leading order waveform,

[TABLE]

where $v=(\pi Mf)^{1/3}$ is the standard PN parameter. The sign is set following the conventions of LAL LIGO Scientific Collaboration, Virgo Collaboration (2018). The phase has the schematic form

[TABLE]

where $t_{0}$ and $\phi_{0}$ are constants related to the freedom of choosing the time and phase of coalescence. We separated point particle terms $\Phi_{\rm F2}^{\rm PP}(\eta)$, spin terms $\Phi_{\rm F2}^{\rm Spin}(\eta,\chi_{1},\chi_{2})$ and tidal terms $\Phi_{\rm F2}^{\rm Tidal}(\eta,\Lambda_{1},\Lambda_{2})$. We include point particle terms to the highest know order of 3.5PN (see e.g. Blanchet (2014)).

The spin terms can be decomposed into spin-orbit and spin-spin terms (ignoring cubic-in-spin contributions at 3.5PN) as $\Phi_{\rm F2}^{\rm Spin}=\Phi_{\rm F2}^{\rm SO}+\Phi_{\rm F2}^{\rm SS}$. We include spin-orbit terms up to 3.5PN Blanchet (2014). For spin-square terms, as the SEOBNRv4T waveform does not include yet 3PN spin-spin effects (see Sec. II.2), we stop at the leading 2PN in the TaylorF2 waveform and ignore the known 3PN terms. The main reason for doing this is because including only common physical effects in both SEOBNRv4T and TaylorF2 makes the residual $\Delta\Phi(Mf;\mathbf{x})$ smaller and have less variation, making it easier to fit at a given error requirement.

The tidal terms $\Phi_{\rm F2}^{\rm Tidal}(\eta,\Lambda_{1},\Lambda_{2})$ are known to 6PN order Vines et al. (2011). These tidal contributions to the phase are exactly as used in the LAL waveform TaylorF2.

For completeness, we give below explicit expressions of the PN coefficients entering the phasing, using the notations $\delta\equiv(m_{1}-m_{2})/(m_{1}+m_{2})$ for the mass difference, and $\chi_{s,a}=(\chi_{1}\pm\chi_{2})/2$, $\kappa_{s,a}=(\kappa_{1}\pm\kappa_{2})/2$ for the symmetrized and antisymmetrized combinations of $\chi_{A}$, $\kappa_{A}$.

The point particle terms coefficients up to 3.5PN can be found e.g. in Eq. (3.18) of Ref. Buonanno et al. (2009), and have the structure

[TABLE]

The non-zero coefficients read

[TABLE]

where $\gamma_{E}$ is Euler’s constant.

The spin-orbit and spin-spin corrections have the structure

[TABLE]

The non-zero spin-orbit coefficients are

[TABLE]

while the leading-order quadratic-in-spin term reads

[TABLE]

As explained above, to be consistent with the fact that this PN information has not yet been incorporated in SEOBNRv4T, we do not include the next-to-leading 3PN term, which for reference is given by Bohé et al. (2015)

[TABLE]

Finally, the tidal contributions to the phasing take the form Vines et al. (2011)

[TABLE]

Introducing convenient mass-weighted combinations of the tidal parameters, expressed using $X_{1}=m_{1}/(m_{1}+m_{2})$ and $X_{2}=m_{2}/(m_{1}+m_{2})$,

[TABLE]

the coefficients are given by

[TABLE]