Improving the NRTidal model for binary neutron star systems

TL;DR

This paper enhances the NRTidal gravitational waveform model for binary neutron star systems by incorporating improved analytical expressions, tidal amplitude corrections, and extended spin effects, resulting in more accurate and computationally efficient waveforms.

Contribution

The paper introduces a new closed-form tidal phase expression, tidal amplitude correction, and extended spin effects into the NRTidal model, and integrates these into existing waveform models for better BNS GW analysis.

Findings

The improved model matches numerical relativity waveforms more accurately.

The new phase and amplitude corrections enhance waveform fidelity.

Model validation shows consistency across a wide parameter space.

Abstract

Accurate and fast gravitational waveform (GW) models are essential to extract information about the properties of compact binary systems that generate GWs. Building on previous work, we present an extension of the NRTidal model for binary neutron star (BNS) waveforms. The upgrades are: (i) a new closed-form expression for the tidal contribution to the GW phase which includes further analytical knowledge and is calibrated to more accurate numerical relativity data than previously available; (ii) a tidal correction to the GW amplitude; (iii) an extension of the spin-sector incorporating equation-of-state-dependent finite size effects at quadrupolar and octupolar order; these appear in the spin-spin tail terms and cubic-in-spin terms, both at 3.5PN. We add the new description to the precessing binary black hole waveform model IMRPhenomPv2 to obtain a frequency-domain precessing binary neutron star model. In addition, we extend the SEOBNRv4_ROM and IMRPhenomD aligned-spin binary black hole waveform models with the improved tidal phase corrections. Focusing on the new IMRPhenomPv2_NRTidalv2 approximant, we test the model by comparing with numerical relativity waveforms as well as hybrid waveforms combining tidal effective-one-body and numerical relativity data. We also check consistency against a tidal effective-one-body model across large regions of the BNS parameter space.

Tables

Table 1. Table 1: The non-spinning BNS and BBH hybrids employed in the construction of the NRTidalv2 model. The columns refer to the name, the employed EOS, the individual masses of the stars M A superscript 𝑀 𝐴 M^{A} , M B superscript 𝑀 𝐵 M^{B} , the tidal deformabilities Λ A , Λ B superscript Λ 𝐴 superscript Λ 𝐵 \Lambda^{A},\Lambda^{B} , the tidal coupling constant κ eff T subscript superscript 𝜅 𝑇 eff \kappa^{T}_{\rm eff} [Eq. ( 8 )], and the ID in the CoRe and SXS databases.

Name	EOS	$M_A$	$M_B$	${\Lambda }_A$	${\Lambda }_B$	${\kappa }_{\mathrm{eff}}^T$	ID
SLy	SLy	1.350	1.350	392.1	392.1	73.5	CoRe:BAM:0095333Our work employs a higher resolution than currently available for this setup in the CoRe catalog.
H4	H4	1.372	1.372	1013.4	1013.4	190.0	CoRe:BAM:0037
MS1b	MS1b	1.350	1.350	1389.4	1389.4	288.1	CoRe:BAM:0064
BBH	–	1.350	1.350	0	0	0	SXS:BBH:0066

Table 2. Table 2: Overview of the existing NRTidal approximants. The individual columns refer to: the name of the approximant, the BBH baseline, the employed tidal phase, the employed spin-spin and cubic-in-spin contribution, employed tidal amplitude corrections, and the incorporation of precession, as well as the computational time Δ ​ T f min Δ subscript 𝑇 subscript f min \Delta T_{\rm f_{min}} of the model to produce a single waveform for a non-spinning, equal mass binary with individual masses M A , B = 1.35 subscript 𝑀 𝐴 𝐵 1.35 M_{A,B}=1.35 and Λ A , B = 400 subscript Λ 𝐴 𝐵 400 \Lambda_{A,B}=400 on an Intel Xeon E5-2630v3 processor for various starting frequencies.

LAL approximant name	BBH baseline	${\psi }_T$	spin-spin	cubic-in-spin	tidal amp.	precession	$\Delta T_{f_{\min }}\left[s\right]$
							$10\mathrm{Hz}$	$20\mathrm{Hz}$	$30\mathrm{Hz}$	$40\mathrm{Hz}$
IMRPhenomD_NRTidal	IMRPhenomD	NRTidal	up to 3PN (BBH)	✗	✗	✗	2.55	0.29	0.14	0.07
IMRPhenomD_NRTidalv2	IMRPhenomD	NRTidalv2	up to 3PN	✗	✗	✗	2.54	0.29	0.14	0.07
SEOBNRv4_ROM_NRTidal	SEOBNRv4_ROM	NRTidal	up to 3PN	✗	✗	✗	3.39	0.40	0.18	0.09
SEOBNRv4_ROM_NRTidalv2	SEOBNRv4_ROM	NRTidalv2	up to 3PN	✗	✗	✗	3.34	0.40	0.18	0.09
IMRPhenomPv2_NRTidal	IMRPhenomPv2	NRTidal	up to 3PN	✗	✗	✓	7.30	0.90	0.43	0.21
IMRPhenomPv2_NRTidalv2	IMRPhenomPv2	NRTidalv2	up to 3.5PN	up to 3.5PN	✓	✓	8.56	1.06	0.51	0.28

Table 3. Table 3: NR BNS configurations for validation of the time domain phasing. The columns refer to the CoRe -ID of the setup, the EOS (see Ref.  Read  et al. ( 2009 ) ), the NSs’ individual masses M A , B subscript 𝑀 𝐴 𝐵 M_{A,B} , the stars’ dimensionless spins χ A , B subscript 𝜒 𝐴 𝐵 \chi_{A,B} , the tidal deformabilities of the stars Λ A , B subscript Λ 𝐴 𝐵 \Lambda_{A,B} , the tidal deformability of the binary Λ ~ ~ Λ \tilde{\Lambda} , the grid resolution covering the NS, and the residual eccentricity of the configuration. In the last column we state whether we employ Richardson extrapolation for a better estimate of the phase. For those setups the errors shown in Fig.  8 present a conservative error measure and are shown as green shaded regions. We note that setups CoRe:BAM:0037 and CoRe:BAM:0064 have also been employed for the calibration of the model; see Eq. ( 5 ).

Name	EOS	$M_A$ $\left[M_{\odot }\right]$	$M_B$ $\left[M_{\odot }\right]$	${\chi }_A$	${\chi }_B$	${\Lambda }_A$	${\Lambda }_B$	$\tilde{\Lambda }$	$h_{\mathrm{fine}}$ $\left[M_{\odot }\right]$	$e$ $\left[{10}^{-3}\right]$	Richardson
CoRe:BAM:0001	2B	1.371733	1.371733	0.000	0.000	126.73	126.73	126.73	0.0930	7.1	✗
CoRe:BAM:0011	ALF2	1.500006	1.500006	0.000	0.000	382.77	382.77	382.77	0.1250	3.1	✗
CoRe:BAM:0037	H4	1.371733	1.371733	0.000	0.000	1006.2	1006.2	1006.2	0.0833	0.9	✓
CoRe:BAM:0039	H4	1.372588	1.372588	0.141	0.141	1001.8	1001.8	1001.8	0.0833	0.5	✓
CoRe:BAM:0062	MS1b	1.350398	1.350398	-0.099	-0.099	1531.5	1531.5	1531.5	0.0970	1.8	✓
CoRe:BAM:0064	MS1b	1.350032	1.350032	0.000	0.000	1531.5	1531.5	1531.5	0.0970	1.8	✓
CoRe:BAM:0068	MS1b	1.350868	1.350868	0.149	0.149	1525.2	1525.2	1525.2	0.0970	2.3	✓
CoRe:BAM:0081	MS1b	1.500016	1.000001	0.000	0.000	863.8	7022.3	2425.5	0.1250	15.	✗
CoRe:BAM:0094	MS1b	1.944006	0.944024	0.000	0.000	182.9	9279.9	1308.2	0.1250	3.4	✗
CoRe:BAM:0105	SLy	1.350608	1.350608	0.106	0.106	388.2	388.24	388.2	0.0783	0.7	✓

Table 4. Table 4: BNS hybrid configurations. The columns describe: the name of the hybrid (CoRe database ID), the EOS, cf.  Read  et al. ( 2009 ) , the NSs’ individual masses M A , B subscript 𝑀 𝐴 𝐵 M_{A,B} , the stars’ dimensionless spins χ A , B subscript 𝜒 𝐴 𝐵 \chi_{A,B} , the stars’ compactnesses C A , B subscript 𝐶 𝐴 𝐵 C_{A,B} , the tidal deformabilities of the stars Λ A , B subscript Λ 𝐴 𝐵 \Lambda_{A,B} , the tidal deformability of the binary Λ ~ ~ Λ \tilde{\Lambda} , the effective dimensionless coupling constant κ eff T subscript superscript 𝜅 𝑇 eff \kappa^{T}_{\rm eff} , and the merger frequency f mrg subscript 𝑓 mrg f_{\rm mrg} .

Name	EOS	$M_A$ $\left[M_{\odot }\right]$	$M_B$ $\left[M_{\odot }\right]$	${\chi }_A$	${\chi }_B$	$C_A$	$C_B$	${\Lambda }_A$	${\Lambda }_B$	$\tilde{\Lambda }$	${\kappa }_{\mathrm{eff}}^T$	$f_{\mathrm{mrg}}$ [Hz]
equal mass, non-spinning
CoRe:Hyb:0001	2B	1.3500	1.3500	0.000	0.000	0.205	0.205	127.5	127.5	127.5	23.9	2567
CoRe:Hyb:0002	SLy	1.3500	1.3500	0.000	0.000	0.174	0.174	392.1	392.1	392.1	73.5	2010
CoRe:Hyb:0003	H4	1.3717	1.3717	0.000	0.000	0.149	0.149	1013.4	1013.4	1013.4	190.0	1535
CoRe:Hyb:0004	MS1b	1.3500	1.3500	0.000	0.000	0.142	0.142	1536.7	1536.7	1536.7	288.1	1405
CoRe:Hyb:0005	MS1b	1.3750	1.3750	0.000	0.000	0.144	0.144	1389.4	1389.4	1389.4	260.5	1416
CoRe:Hyb:0006	SLy	1.3750	1.3750	0.000	0.000	0.178	0.178	347.3	347.3	347.3	65.1	1978
unequal mass, non-spinning
CoRe:Hyb:0007	MS1b	1.5000	1.0000	0.000	0.000	0.157	0.109	866.5	7041.6	2433.5	456.3	1113
CoRe:Hyb:0008	MS1b	1.6500	1.1000	0.000	0.000	0.171	0.118	505.2	4405.9	1490.1	279.4	1170
CoRe:Hyb:0009	MS1b	1.5278	1.2222	0.000	0.000	0.159	0.130	779.6	2583.2	1420.4	266.3	1301
CoRe:Hyb:0010	SLy	1.5000	1.0000	0.000	0.000	0.194	0.129	192.3	2315.0	720.0	135.0	1504
CoRe:Hyb:0011	SLy	1.5274	1.2222	0.000	0.000	0.198	0.157	167.5	732.2	365.6	68.6	1770
CoRe:Hyb:0012	SLy	1.6500	1.0979	0.000	0.000	0.215	0.142	93.6	1372.3	408.1	76.5	1592
equal mass, spinning
CoRe:Hyb:0013	H4	1.3726	1.3726	+0.141	+0.141	0.149	0.149	1009.1	1009.1	1009.1	189.2	1605
CoRe:Hyb:0014	MS1b	1.3504	1.3504	-0.099	-0.099	0.142	0.142	1534.5	1534.5	1534.5	287.7	1323
CoRe:Hyb:0015	MS1b	1.3504	1.3504	+0.099	+0.099	0.142	0.142	1534.5	1534.5	1534.5	287.7	1442
CoRe:Hyb:0016	MS1b	1.3509	1.3509	+0.149	+0.149	0.142	0.142	1531.8	1531.8	1531.8	287.2	1456
CoRe:Hyb:0017	SLy	1.3502	1.3502	+0.052	+0.052	0.174	0.174	392.0	392.0	392.0	73.5	2025
CoRe:Hyb:0018	SLy	1.3506	1.3506	+0.106	+0.106	0.174	0.174	391.0	391.0	391.0	73.5	2048

Equations

$$h(t)=A(t)e^{-\mathrm{i}\varphi (t)},$$

$$\varphi (\hat{\omega })={\varphi }_{\mathrm{pp}}(\hat{\omega })+{\varphi }_{\mathrm{SO}}(\hat{\omega })+{\varphi }_{\mathrm{SS}}(\hat{\omega })+{\varphi }_{\mathrm{T}}(\hat{\omega })+\cdots ,$$

$$\tilde{h}(f)=\tilde{A}(f)e^{-\mathrm{i}\psi (f)},,$$

$$\psi (\hat{\omega })={\psi }_{\mathrm{pp}}(\hat{\omega })+{\psi }_{\mathrm{SO}}(\hat{\omega })+{\psi }_{\mathrm{SS}}(\hat{\omega })+{\psi }_{\mathrm{T}}(\hat{\omega })+\cdots .$$

$$\begin{split}\phi_{T}^{\rm NR}&=\frac{1}{320+192+192}\bigg{[}320(\phi_{\rm SLy}-\phi_{\rm BBH})+192(\phi_{\rm H4}-\phi_{\rm BBH})\\ &\quad+192(\phi_{\rm MS1b}-\phi_{\rm BBH})\bigg{]}\\ &=\frac{1}{11}\bigg{(}5\phi_{\rm SLy}+3\phi_{\rm H4}+3\phi_{\rm MS1b}-11\phi_{\rm BBH}\bigg{)}.\end{split}$$

$${\varphi }_T$$

$$+c_2^A{x}^2+c_{5/2}^A{x}^{5/2})+[A\leftrightarrow B],$$

$$c_{\mathrm{Newt}}^A$$

$$c_1^A$$

$$c_{3/2}^A$$

$$c_2^A$$

$$c_{5/2}^A$$

$${\kappa }_{\mathrm{eff}}^T=\frac{2}{13}\left[\left(1+12\frac{X_B}{X_A}\right){\left(\frac{X_A}{C_A}\right)}^5{k}_2^A+(A\leftrightarrow B)\right],$$

$$\tilde{\Lambda }=\frac{16}{3}{\kappa }_{\mathrm{eff}}^T,$$

$${\Lambda }_{A,B}=\frac{2}{3}\frac{k_2^{A,B}}{C_{A,B}^5}.$$

$${\varphi }_T(x)=-{\kappa }_{\mathrm{eff}}^T\frac{13}{8\nu }{x}^{5/2}{P}_{NRTidalv2}(x),$$

$$P_{NRTidalv2}(x)=\frac{1+n_{1}x+n_{3/2}{x}^{3/2}+n_2{x}^2+n_{5/2}{x}^{5/2}+n_3{x}^3}{1+d_{1}x+d_{3/2}{x}^{3/2}+d_2{x}^2}.$$

$$n_1$$

$$n_{3/2}$$

$$n_2$$

$$d_{3/2}$$

$$c_1$$

$$c_2$$

$$n_{5/2}$$

$$d_1$$

$$\frac{{\text{d}}^2{\psi }_T(\omega )}{\text{d}{\omega }^2}=\frac{1}{\omega }\frac{\text{d}{\varphi }_T(\omega )}{\text{d}\omega }$$

$${\psi }_T(x)=-{\kappa }_{\mathrm{eff}}^T\frac{39}{16\nu }{x}^{5/2}{\tilde{P}}_{NRTidalv2}(x),$$

$${\tilde{P}}_{NRTidalv2}(x)=\frac{1+{\tilde{n}}_{1}x+{\tilde{n}}_{3/2}{x}^{3/2}+{\tilde{n}}_2{x}^2+{\tilde{n}}_{5/2}{x}^{5/2}+{\tilde{n}}_3{x}^3}{1+{\tilde{d}}_{1}x+{\tilde{d}}_{3/2}{x}^{3/2}+{\tilde{d}}_2{x}^2},$$

$${\tilde{n}}_1$$

$${\tilde{n}}_{3/2}$$

$${\tilde{n}}_2$$

$${\tilde{d}}_{3/2}$$

$${\tilde{c}}_1$$

$${\tilde{c}}_2$$

$${\tilde{n}}_{5/2}$$

$${\tilde{d}}_1$$

$${\tilde{A}}_T^{6PN}=\sqrt{\frac{160\pi \nu }{27}}\frac{M^2}{D_L}{\kappa }_{\mathrm{eff}}^T{x}^{-7/4}\left(-\frac{27}{16}{x}^5-\frac{449}{64}{x}^6\right),$$

$$\begin{array}{rl}{\tilde{A}}_T^{\mathrm{Kawaguchi}}& =\sqrt{\frac{160\pi \nu }{27}}\frac{M^2}{D_L}{\kappa }_{\mathrm{eff}}^T{x}^{-7/4}\\ & \times \left(-\frac{27}{16}{x}^5-\frac{449}{64}{x}^6-4251x^{7.890}\right).\end{array}$$

Full text

The upgrades are: (i) a new closed-form expression for the tidal contribution to the GW phase which includes further analytical knowledge and is calibrated to more accurate numerical relativity data than previously available; (ii) a tidal correction to the GW amplitude; (iii) an extension of the spin-sector incorporating equation-of-state-dependent finite size effects at quadrupolar and octupolar order; these appear in the spin-spin tail terms and cubic-in-spin terms, both at 3.5PN.

We add the new description to the precessing binary black hole waveform model IMRPhenomPv2 to obtain a frequency-domain precessing binary neutron star model. In addition, we extend the SEOBNRv4_ROM and IMRPhenomD aligned-spin binary black hole waveform models with the improved tidal phase corrections. Focusing on the new IMRPhenomPv2_NRTidalv2 approximant, we test the model by comparing with numerical relativity waveforms as well as hybrid waveforms combining tidal effective-one-body and numerical relativity data. We also check consistency against a tidal effective-one-body model across large regions of the BNS parameter space.

I Introduction

The first gravitational wave (GW) signal associated with electromagnetic (EM) counterparts, detected on the 17th of August 2017, marks a breakthrough in the field of multi-messenger astronomy Abbott et al. (2017, 2017, 2017). Analyses of the GW and EM signatures favor a binary neutron star (BNS) coalescence, e.g., Tanaka et al. (2017); Nicholl et al. (2017); Tanvir et al. (2017); Perego et al. (2017); Waxman et al. (2018); Metzger et al. (2018); Abbott et al. (2019); Kawaguchi et al. (2018); Hinderer et al. (2018); Coughlin et al. (2018); Coughlin and Dietrich (2019); Siegel (2019). Due to the increasing sensitivity of advanced GW detectors, multiple detections of merging BNSs are expected in the near future Abbott et al. (2018).

A prerequisite to extract information from the data are theoretical predictions about the emitted GW signal. The properties of the system are typically inferred via a coherent Bayesian analysis based on cross-correlation of the measured strain with predicted waveform approximants, e.g., Veitch et al. (2015). These cross-correlations are done for a large number of target waveforms and require large computational resources. Thus, the computation of each individual waveform needs to be efficient and fast to ensure that the Bayesian parameter estimation of signals, containing several thousand GW cycles (as typical for BNS systems), is at all manageable. On the other hand, waveform models need to be accurate enough to allow a correct estimate of the source properties, such as the masses, the spins, and internal structure of the NSs.

Over the last years, there has been significant progress modeling the GW signal associated with the BNS coalescence, including the computation of higher-order tidal corrections or spin-tidal coupling, e.g., Refs. Damour et al. (2012); Pani et al. (2018); Banihashemi and Vines (2018); Abdelsalhin et al. (2018); Landry (2018); Jiménez Forteza et al. (2018), and improved accuracy of BNS numerical relativity (NR) simulations Hotokezaka et al. (2015); Bernuzzi and Dietrich (2016); Kiuchi et al. (2017); Dietrich et al. (2018, 2018); Foucart et al. (2019). However, although the analytical progress has improved the performance of post-Newtonian (PN) waveform approximants, PN models still become increasingly inaccurate towards the merger, e.g. Bernuzzi et al. (2012); Favata (2014); Wade et al. (2014); Hotokezaka et al. (2016); Dietrich et al. (2019); Dudi et al. (2018); Samajdar and Dietrich (2018).

Most of the current time-domain tidal waveform models Bernuzzi et al. (2015); Hotokezaka et al. (2015); Hinderer et al. (2016); Steinhoff et al. (2016); Nagar et al. (2018); Nagar and Rettegno (2019); Akcay et al. (2019); Nagar et al. (2019) are based on the effective-one-body (EOB) description of the general relativistic two-body problem Buonanno and Damour (1999); Damour and Nagar (2010). This approach has proven to be able to predict the BNS merger dynamics in large regions of the BNS parameter space, but recent numerical relativity (NR) data revealed configurations for which further improvements of the tidal EOB models are required Hotokezaka et al. (2015); Dietrich and Hinderer (2017); Akcay et al. (2019). While one can expect that over the next years, these issues will be overcome due to further progress in the fields of NR, gravitational self-force, and PN theory, the high computational cost for a single EOB waveform is yet another disadvantage. One possibility to speed up the EOB computation is the use of high-order post-adiabatic approximations of the EOB description to allow an accurate and efficient evaluation of the waveform up to a few orbits before merger Nagar and Rettegno (2019). The other possibility, and most common approach, is constructing reduced-order-models Lackey et al. (2017, 2019). Those models allow the fast computation of waveforms in the frequency domain and are well suited for a direct use in parameter estimation pipelines.

In addition to PN and EOB approximants, there have been proposals for alternative ways to describe tidal GW signals. Refs. Lackey et al. (2014); Pannarale et al. (2013) develop phenomenological black hole-neutron star (BHNS) approximants based on NR data. Ref. Barkett et al. (2016) transforms NR simulations of binary black hole (BBH) systems by adding PN tidal effects, and Refs. Lange et al. (2017, 2018) develop a method to employ NR waveforms or computationally expensive waveform approximants (such as tidal EOB waveforms) directly for parameter estimation.

Another approach to describe BNS systems was presented in Ref. Dietrich et al. (2017), in which BBH models have been augmented by an analytical closed-form expression correcting the GW phase to include tidal effects. This waveform model Dietrich et al. (2017, 2019), referred to as NRTidal, was implemented in the LSC Algorithm Library (LAL) LIGO Scientific Collaboration and Virgo Collaboration (2018) to support the analysis of GW170817 by the LIGO and Virgo Collaborations (LVC) Abbott et al. (2017, 2019, 2018, 2018, 2019) and has also been used outside the LVC, e.g. Dai et al. (2018); Radice and Dai (2019). In addition, Ref. Kawaguchi et al. (2018) developed an alternative tidal approximant in the frequency domain combining EOB and NR information following a similar idea as in Ref. Dietrich et al. (2017).

Studies showed that for GW170817, with its signal-to-noise ratio (SNR) of $\sim 30$, waveform model systematics are within the statistical uncertainties, i.e., that different employed tidal GW models give slightly different, but consistent constraints on the binary properties, e.g., Abbott et al. (2019). However, systematic effects will grow for an increasing number of detections or GW observations with larger SNRs Dudi et al. (2018); Samajdar and Dietrich (2018). Ref. Samajdar and Dietrich (2018) stated that for a GW170817-like event measured with the anticipated design sensitivity of the Advanced LIGO and Advanced Virgo detectors, systematic effects will dominate and the extracted equation of state (EOS) constraints between existing waveform approximants will become inconsistent. Furthermore, the analysis presented in Ref. Dietrich et al. (2019); Dudi et al. (2018) showed that the original NRTidal model could potentially underestimate tidal deformabilities, leading to possible biases for future detections with larger SNRs.

Therefore, to further push for the availability of a fast and accurate waveform model employable for the upcoming observing runs in the advanced detector era, after recalling the basics of NRTidal and discussing the NR simulations and hybrid waveform construction in Sec. II, we improve the NRTidal description by:

1. (i)

Recalibrating the closed-form phenomenological tidal description including additional analytical knowledge and using improved NR data (Sec. III.1); 2. (ii)

Adding a tidal GW amplitude correction to the model (Sec. III.2); 3. (iii)

Incorporating EOS-dependent 3.5PN spin-spin and cubic-in-spin effects proportional to the quadrupole and octupole moments of the NSs Marsat (2015); Bohé et al. (2015); Krishnendu et al. (2017); Nagar et al. (2019) (Sec. III.3).

We validate the new NRTidalv2 approximant with a set of $10$ high-resolution NR waveforms (Sec. IV.1) and $18$ hybrids of NR waveforms and the TEOBResumS tidal EOB model Nagar et al. (2018) (Sec. IV.2). Furthermore, we compare the model in a larger region of the parameter space than currently covered with NR simulations by computing the mismatch with respect to the SEOBNRv4T tidal EOB model Hinderer et al. (2016); Steinhoff et al. (2016) (Sec. IV.3). We note that for this waveform model Ref. Lackey et al. (2019) recently developed a reduced order model which can also be used directly for GW data analysis. We conclude in Sec. V. In the Appendices, we discuss possible extensions to the model, considering the tidal amplitude correction (Appendix A) and the mass ratio dependence of the tidal phase (Appendix B).

In this article geometric units are used by setting $G=c=M_{\odot}=1$. At some places units are given explicitly to allow a better interpretation. Further notations are $M=M_{A}+M_{B}$ for the total mass of the system, $\chi_{A},\chi_{B},\Lambda_{A},\Lambda_{B}$ for the individual dimensionless spins and tidal deformabilities of the stars. The mass ratio of the system is $q=M_{A}/M_{B}$ and the symmetric mass ratio is $\nu=M_{A}M_{B}/(M_{A}+M_{B})^{2}$. We define the labeling of the individual stars so their masses satisfy $M_{A}\geq M_{B}$.

II Basic ideas and improved numerical relativity data

II.1 The basic idea of NRTidal

During the BNS coalescence, each star gets deformed due to the gravitational field of the companion. These tidal deformations accelerate the inspiral and leave a clear imprint in the GW signal, e.g., Flanagan and Hinderer (2008). Consequently, the theoretical modeling of BNSs and the extraction of tidal effects from measured GW signals is an important way of determining the internal structure of NSs and thus the EOS of supranuclear dense matter.

The complex time-domain GW signal is given by

[TABLE]

with amplitude $A(t)$ and time-domain phase $\phi(t)$. Here we only consider the dominant $2,2$ [spin $(-2)$ weighted spherical harmonic] mode. We assume in the following that the phase can be decomposed into

[TABLE]

where the dimensionless GW frequency is given by $\hat{\omega}=M\omega=M\partial_{t}\phi(t)$. Here $\phi_{\rm pp}$ denotes the nonspinning, point-particle, contribution to the overall phase, $\phi_{\rm SO}$ corresponds to contributions caused by spin-orbit coupling, $\phi_{\rm SS}$ corresponds to contributions caused by spin-spin effects (both self-spin and spin-interactions), and $\phi_{\rm T}$ denotes the tidal effects present in the GW phase.

Similar to Eq. (1), the waveform can be written in the frequency domain as

[TABLE]

with GW frequency $f$ and frequency domain amplitude $\tilde{A}(f)$ and phase $\psi(f)$. Here we assume again:

[TABLE]

Constraints on the supranuclear EOS governing the matter inside NSs rely on an accurate measurement of the tidal phase contribution. This contribution enters first at the 5th PN order.111 There is also the possibility of extracting EOS information from the spin-spin interaction first entering in the 2PN $\psi_{\text{SS}}$ contribution, where the individual terms of $\psi_{\text{SS}}$ are proportional to the square of the individual spins, i.e., $\chi_{A}^{2}$, $\chi_{B}^{2}$, or $\chi_{A}\chi_{B}$. Although the maximum NS spin in a BNS is not precisely known, the fastest spinning NS in a BNS system observed to date (PSR J1946+2052 Stovall et al. (2018)) will only have a dimensionless spin of $\sim 0.02$–$0.04$ at merger Kumar and Landry (2019). Thus, obtaining EOS information from the spin-spin phase contribution is extremely challenging.

The main idea of the NRTidal approach is to provide a closed-form approximation for the tidal phase $\phi_{\rm T}$ or $\psi_{\rm T}$. Because standard GW data analysis is carried out in the frequency domain, the frequency domain model is of particular importance, due to its efficiency. In addition to the tidal contribution, the final NRTidal approximant also incorporates EOS dependent effects in $\psi_{\rm SS}$, since the spin-spin contributions depend on the quadrupole and higher moments of the individual stars, and thus on the internal structure of the stars.

We note that there are higher-order spin-tidal coupling effects that have recently been computed Abdelsalhin et al. (2018); Landry (2018). However, as outlined in Jiménez Forteza et al. (2018), these terms will be unmeasurable in the advanced GW detector era. Therefore, we do not include them in the current description to avoid unnecessary computational costs.

II.2 High-precision NR simulations

The field of NR has made significant progress over the last years. Nevertheless, the production of highly accurate gravitational BNS waveforms remains challenging and there exist only a small number of simulations with low eccentricity and with phase errors small enough to allow GW modeling; cf. Refs. Hotokezaka et al. (2015); Bernuzzi and Dietrich (2016); Dietrich et al. (2018, 2018); Kiuchi et al. (2017); Foucart et al. (2019).

In addition to the dataset used for the original NRTidal calibration Dietrich et al. (2017, 2019), we performed one additional simulation for a non-spinning equal-mass BNS setup employing a piecewise-polytropic parametrization of the SLy Read et al. (2009) EOS. This EOS is in agreement with recent constraints extracted from GW170817 Abbott et al. (2018, 2018); Bauswein et al. (2017); Most et al. (2018); Radice and Dai (2019); Coughlin et al. (2018) and thus is a natural choice for our work.222While the maximum mass of $2.05M_{\odot}$ of the SLy EOS is slightly outside of the $68.3\%$ credible region of the recent heavy pulsar mass measurement in Cromartie et al. (2019) ($[2.07,2.28]M_{\odot}$), it is well inside the $95.4\%$ credible region of $[1.97,2.40]M_{\odot}$, which is why we still consider it here. The same physical configuration has already been used in the past for the construction of the NRTidal model Dietrich et al. (2017, 2018); cf. Table 1 for further details. In Dietrich et al. (2017, 2018), we have simulated this setup with the BAM code Brügmann et al. (2008); Thierfelder et al. (2011); Dietrich et al. (2015); Bernuzzi and Dietrich (2016) for $5$ different resolutions with $64$, $96$, $128$, $192$, and $256$ points in the finest refinement level covering the individual NSs. Here, we add one additional simulation with $320$ points in the finest refinement level. This corresponds to a spatial resolution of $0.047M_{\odot}\approx 70\ \rm m$ and computational costs of $\sim 5$ million CPU-hours for this single resolution.

The availability of six different resolutions and the presence of clean convergence across multiple resolutions allows us to employ Richardson extrapolation to obtain an improved GW signal and to provide an associated error budget; see Ref. Bernuzzi and Dietrich (2016) for more details. We present the GW signal for the different resolutions in Fig. 1 (top panel) and the phase difference and convergence properties in the middle and bottom panels.

Except for the lowest resolution, clean second order convergence is obtained throughout the inspiral. This becomes evident by comparison of the individual phase differences with the phase differences rescaled assuming second order convergence (dashed lines). For the lowest resolution setup ($n_{64}$), second order convergence is lost a few orbits before merger ($u\approx 1500M$). Merger times for each resolution are indicated by vertical solid lines in Fig. 1.

The phase difference between the highest ($n_{320}$) and second highest resolution ($n_{256}$) is $0.52\ \rm rad$ at the moment of merger. Performing a Richardson extrapolation Bernuzzi and Dietrich (2016), we obtain more accurate phase descriptions. We denote the Richardson extrapolated data obtained from the resolutions $n_{X}$ and $n_{Y}$ as $R(n_{X},n_{Y})$. We crosscheck the robustness of the procedure by presenting the phase differences $R(n_{320},n_{256})-n_{320}$ and $R(n_{256},n_{192})-n_{256}$ in the bottom panel of Fig. 1. Rescaling the phase difference of $R(n_{320},n_{256})-n_{320}$ assuming second order convergence shows excellent agreement with $R(n_{256},n_{192})-n_{256}$. This demonstrates that the leading error term scales quadratically with respect to the grid spacing/resolution.

Thus, we can estimate the uncertainty of the Richardson extrapolated waveform $R(n_{320},n_{256})$ to be the difference with the $n_{320}$ resolution. At the moment of merger, this gives an uncertainty of $0.37\ \rm rad$. At this time the estimated error due to the finite radius extraction is below $0.044\ \rm rad$, which leads to a conservatively estimated total error of $\Delta\phi_{\rm mrg}\lesssim 0.38\ \rm rad$ at merger.

An alternative, but not conservative, error measure is given by the difference between the two Richardson extrapolated waveforms (green line in the bottom panel). We find that throughout the inspiral the difference between the $R(n_{320},n_{256})$ and $R(n_{256},n_{192})$ is below $0.1\ \rm rad$ (at the moment of merger $\Delta\phi=0.087\ \rm rad$, which would lead to a total error of $\lesssim 0.1\ \rm rad$ once finite radius extraction is included).

In addition to this new setup, we also consider the additional two high resolution simulations available in the CoRe database Dietrich et al. (2018), cf. Table 1. These setups, CoRe:BAM:0037 and CoRe:BAM:0064, only employ 192 points across the star and have conservatively estimated phase uncertainties at merger of $1.20\ \rm rad$ and $2.27\ \rm rad$, respectively. We incorporated this accuracy difference by weighting the individual setups differently during the construction of the NRTidalv2 phase, as discussed in the next subsection.

II.3 Hybrid Construction

In the original NRTidal work, PN, EOB, and NR approximants have been separately used in different frequency intervals. Here, we start by constructing hybrid waveforms consisting of a time domain tidal EOB model (TEOBResumS) inspiral Nagar et al. (2018) connected to the high resolution NR simulation discussed above. The hybridization is performed as discussed in Refs. Dietrich et al. (2019); Dudi et al. (2018) to which we refer for further details. In addition to the BNS hybrid waveforms, we also create a hybrid between the non-tidal version of the TEOBResumS model and a binary black hole waveform computed with the SpEC code SpE , setup SXS:BBH:0066 of the public SXS catalog Mroué et al. (2013); Blackman et al. (2015). All the hybrids have an initial frequency of $20\ \rm Hz$.

We present the time domain phase evolution of the BBH and BNS hybrids in Fig. 2. For this plot we align the waveforms at $\sim 22\ \rm Hz$.

We emphasize that only the four hybrid waveforms listed in Table 1 are used for calibration of the NRTidalv2 model, where the dataset we are going to fit is

[TABLE]

The factors are obtained by linearly weighting the resolutions of the individual NR data, i.e., 320 points accross the star for the SLy setup and 192 for H4 and MS1b setups. We decided to use this minimal dataset since these are the available data with the highest accuracy. Note that a simple restriction to the highest resolution, i.e., the SLy data, leads to a phase description which does not accurately characterize binaries with large tidal deformabilities. Thus, it would be preferable to include in the future a larger number of NR simulations with varying masses, spins, mass ratios, and EOSs once these are available. However, while there are a small number of high quality waveforms Dietrich et al. (2018), these waveforms do not span a sufficiently large region of the parameter space to incorporate additional mass ratio, EOS, or mass dependencies in our phenomenological ansatz.

III Improvements

The NRTidalv2 approach can be added to any BBH model: we focus our discussion here on the frequency-domain IMRPhenomPv2, IMRPhenomD, and SEOBNRv4_ROM models. We primarily concentrate on the extension of IMRPhenomPv2 Hannam et al. (2014); Khan et al. (2016) describing precessing systems. In addition, we have also added the improved tidal phase description to the SEOBNRv4_ROM Bohé et al. (2017) and the IMRPhenomD Khan et al. (2016) approximants.444See Pürrer (2014, 2016) for more details of the reduced order model technique used to construct SEOBNRv4_ROM from the time domain approximant SEOBNRv4. For SEOBNRv4_ROM and IMRPhenomD, we decided to include only the tidal phase description to reduce additional computational costs and allow a faster computation of waveforms than for IMRPhenomPv2_NRTidalv2.

We present an overview of all existing NRTidal models in Table 2.

III.1 Recalibrating the NRTidal phase

III.1.1 Ansatz for the NRTidal time-domain phase

Non-spinning tidal contributions start entering the GW phasing at the 5PN order and partially known analytical knowledge exists up to 7.5PN Damour et al. (2012):

[TABLE]

with the dimensionless EOB tidal parameter $\kappa_{A}$ (defined below) and $x(\hat{\omega})=(\hat{\omega}/2)^{2/3}$. The individual coefficients $c_{i}^{A}$ are

[TABLE]

and similarly with $A\leftrightarrow B$. Here $X_{A,B}=M_{A,B}/M$. We note that although analytic knowledge exists up to the 7.5PN order, some unknown terms are present at 7PN. As discussed in Ref. Damour et al. (2012), these terms are expected to be small and are set to zero in our definition of $c_{2}^{A}$, $c_{2}^{B}$.

As in the original NRTidal description Dietrich et al. (2017, 2018) we introduce the effective tidal coupling constant $\kappa_{\rm eff}^{T}$ which describes the dominant tidal and mass ratio effects:

[TABLE]

where $C_{A,B}\equiv M_{A,B}/R_{A,B}$ are the compactnesses of the stars at isolation, and $k^{A,B}_{2}$ the Love numbers describing the static quadrupolar deformation of one body in the gravitoelectric field of the companion Damour (1983); Hinderer (2008); Damour and Nagar (2009); Binnington and Poisson (2009). The parameter $\kappa^{T}_{\rm eff}$ is related to $\tilde{\Lambda}$ (the mass-weighted tidal deformability commonly used in GW analysis Wade et al. (2014)) by

[TABLE]

and the individual tidal deformability parameters are given by

[TABLE]

The EOB tidal parameter used in Eq. (6) is given by $\kappa_{A}=3X_{B}X_{A}^{4}\Lambda_{A}$.

In the following, we restrict the parameters $c_{1}^{A,B},c_{3/2}^{A,B},c_{2}^{A,B},c_{5/2}^{A,B}$ to their equal-mass values (due to the absence of a large set of high-quality unequal mass NR data), and therefore, discard the superscripts $A$ and $B$. For this case, an effective representation of tidal effects is obtained using

[TABLE]

with the Padé approximant

[TABLE]

To enforce consistency with the analytic PN knowledge [Eqs. (6)-(7e)], some of the individual terms are restricted

[TABLE]

with

[TABLE]

The remaining unknown $4$ parameters are fitted to the data:

[TABLE]

Figure 3 shows our findings. We show as a gray shaded region the parameter space in $\phi_{T}/\kappa_{\rm eff}^{T}$ covered by our simulations, where the gray dashed line refers explicitly to the SLy configuration. Here we do not include any error estimate in the BBH hybrid used to extract the tidal phase. In addition, we present the 6PN tidal contribution, which the old NRTidal approximant reduces to in the low frequency limit; the 7.5PN contribution, which the new NRTidalv2 reduces to in the low frequency limit; and the 7PN contribution, which is the PN approximant showing the best agreement to the NR data. We also show the tidal phase given in Kawaguchi et al. Kawaguchi et al. (2018),555We obtain the time domain tidal phase approximant from the frequency domain expression given in Ref. Kawaguchi et al. (2018) using the stationary phase approximation. which has been calibrated to NR simulations up to a frequency of $1000\ \rm Hz$ (thin dashed line). The model of Ref. Kawaguchi et al. (2018) loses validity outside its calibration region and overestimates tidal effects at the moment of merger, though this would not affect GW data analysis if a maximum frequency of $1000\ \rm Hz$ is employed, or the signal at frequencies $\gtrsim 1000\ \rm Hz$ is sufficiently suppressed by the detectors’ noise. In addition, we find good agreement between the Kawaguchi et al. fit and the new NRTidalv2 approximation below $1000\ \rm Hz$. We also show the estimated tidal phase extracted by comparing our BBH hybrid with a tidal EOB waveform computed for our SLy configuration using the TEOBResumS Nagar et al. (2018) model. The tidal phase estimate of the TEOBResumS model is slightly less attractive for frequencies around $1000\ \rm Hz$, but more attractive at higher frequencies. Finally, the original NRTidal model is shown as a green line. The tidal contribution is overestimated at about $f\sim 1000\ \rm Hz$, and later underestimated. This oscillatory behavior has been seen before, e.g., Dietrich et al. (2017, 2018, 2019), and could potentially lead to biases in the estimate of tidal effects from GW signals Dudi et al. (2018). For both NRTidal and NRTidalv2 the growth of the tidal phase around merger is much smaller than for any other approximant, which generally reduces possible pathologies in more extreme regions of the parameters, e.g., a cancellation of the point-particle and attractive tidal phase close to merger. As expected the NRTidalv2 model stays within the gray shaded region and, thus, close to the numerical relativity dataset used for the calibration.

III.1.2 Frequency domain phase

As in Ref. Dietrich et al. (2017) we will employ the stationary phase approximation (SPA), discussed in, e.g., Damour et al. (2012), to derive the tidal phase contribution $\psi_{T}$ in the frequency domain, i.e., we solve

[TABLE]

to obtain $\psi_{T}$. Although $\phi_{T}$ is given explicitly, we solve Eq. (16) numerically and approximate the result with a Padé approximant similar to Eqs. (11) and (12):

[TABLE]

with

[TABLE]

and

[TABLE]

where the known coefficients are:

[TABLE]

and the fitting coefficients are:

[TABLE]

We present the final tidal phase contribution in the frequency domain in Fig. 4 for a number of different GW approximants. Fig. 5 shows the corresponding phase differences with respect to the NRTidalv2 model on a double logarithmic scale.

III.2 Tidal amplitude corrections

The extraction of binary properties relies mostly on the GW phase, which makes an accurate description of $\psi$ the primary target of GW modeling. However, a realistic estimate of the GW amplitude is also of importance, e.g., for a precise distance measurement.

Therefore, we will discuss a possible extension of the NRTidal approach including a tidal amplitude correction in the frequency domain. An alternative time and frequency domain amplitude correction is presented in Appendix A.

Here, we will derive the frequency domain tidal correction from the frequency domain representation of the SLy and BBH TEOBResumS-NR hybrids, described in Table 1. We do not employ the H4 and MS1b setup for the amplitude correction since their lower merger frequencies add additional complications during the construction procedure. The top panel of Fig. 6 shows $\tilde{A}$ for our generic setups and also the BBH result augmented ($\tilde{A}=\tilde{A}_{\rm BBH}+\tilde{A}_{\rm T}$) by the 6PN expression:

[TABLE]

e.g., Ref. Hotokezaka et al. (2016), where $D_{L}$ is the luminosity distance of the source, which is the appropriate substitution for the effective distance used in Ref. Hotokezaka et al. (2016) for our case.666Note further that as before, we have restricted our analysis to the leading order mass ratio effect and do not incorporate further mass ratio dependence in the PN parameters. Furthermore, we restrict our consideration to gravitoelectric contributions and do not consider gravitomagnetic tidal effects recently computed in Banihashemi and Vines (2018). Kawaguchi et al. Kawaguchi et al. (2018) extended Eq. (22) to

[TABLE]

Based on the good agreement we have found between the results of Ref. Kawaguchi et al. (2018) and the new NRTidalv2 phase description below $1000\ \rm Hz$, we want to use Eq. (23) as baseline for a possible frequency amplitude extension of the NRTidalv2 approximant.777We note that the phase and amplitude extension presented in Kawaguchi et al. (2018) follow different approaches: While the tidal phase correction is based on an additional contribution due to non-linear tides, i.e., a higher order $\tilde{\Lambda}$ contribution, the amplitude correction only uses linear tidal effects, but adds an effectively higher order PN coefficient. Therefore, the proposed amplitude extension of Kawaguchi et al. (2018) can easily be incorporated in our approach.

For this purpose we employ the ansatz

[TABLE]

Equation (24) ensures that for small frequencies Eq. (23) is recovered, but that the high frequency behavior ($f>1000\ {\rm Hz}$) can be adjusted. We obtain $d=13477.8$ by fitting the data presented in Fig. 6 (blue dashed line in the bottom panel).

As for the previous NRTidal implementation, we add a Planck taper McKechan et al. (2010) to end the inspiral waveform. The taper begins at the estimated merger frequency [Eq. (11) of Dietrich et al. (2019)] and ends at $1.2$ times the merger frequency. Thus, the final amplitude is given as:

[TABLE]

Because of the smooth frequency and amplitude evolution even after the moment of merger, this taper only introduces negligible errors and does not lead to biases in the parameter estimation of even SNR 100 signals, as shown using an injection of an SLy hybrid with the same parameters as those considered here in Dudi et al. (2018).

III.3 Incorporating higher-order spin-spin effects

While nonspinning NSs and black holes only have a nonzero monopole moment, spinning neutron stars and black holes have an infinite series of nonzero (Geroch-Hansen) multipole moments, e.g., Refs. Pappas and Apostolatos (2014); Yagi et al. (2014). The contributions from the stars’ (mass) monopole and (spin) dipole to the binary’s motion are explicitly accounted for in the BBH baseline. Additionally, contributions from higher spin-induced multipoles in the BBH baseline model are indirectly included due to the calibration to NR simulations. However, without further adjustment, all multipoles would be specialized to the black hole values, which (as shown in Dietrich et al. (2019)) noticeably reduces the accuracy for spinning BNS systems. Thus, to improve the performance of NRTidalv2 for spinning configurations, we include an EOS dependence in the quadrupole and octupole, as these are the moments that appear in current PN calculations. These two moments (in their scalar versions) can be written as $M^{(A,B)}_{2}=-C_{\rm Q}^{A}M_{A,B}^{3}\chi_{A,B}^{2}$, $S_{3}^{(A,B)}=-C_{\rm Oc}^{A,B}M_{A,B}^{4}\chi_{A,B}^{3}$, respectively, for star $A$ and $B$. Here $C_{\rm Q}^{A,B}$ and $C_{\rm Oc}^{A,B}$ are the quadrupolar and octupolar spin-induced deformabilities for the individual stars. Both $C_{\rm Q}$ and $C_{\rm Oc}$ are $1$ for a black hole.

In this paper, we extend the existing LALSuite implementation, which currently contains the EOS dependence of the quadrupole moment only up to 3PN Bohé et al. (2015); Mishra et al. (2016) to include the $3.5$PN spin-squared terms, completed using the recently computed 3.5PN tail terms Nagar et al. (2019). We also include leading order spin-cubed terms entering at 3.5PN order.

The contributions of the quadrupole and octupole deformations of the stars to the binary’s binding energy and energy flux have been computed through $3.5$PN, Refs. Bohé et al. (2015); Marsat (2015); Nagar et al. (2018), building on earlier work reviewed in Blanchet (2014). We compute the phase in the frequency domain using the SPA. These contributions to the phase were already presented in Krishnendu et al. (2017), except the $3.5$PN spin-spin terms, as Krishnendu et al. (2017) did not have the $3.5$PN spin-squared tail term from Ref. Nagar et al. (2018). Explicitly, the self-spin (i.e., $C_{\rm Q}$ and $C_{\rm Oc}$) terms in the phasing that we add to the BBH baseline are

[TABLE]

with

[TABLE]

Here we use $\hat{C}_{\rm Q}^{A}:=C_{\rm Q}^{A}-1$ and $\hat{C}_{\rm Oc}^{A}:=C_{\rm Oc}^{A}-1$ to remove the contribution from the black hole multipoles already present in the baseline BBH phase.

Finally, we relate $C_{\rm Q}^{A}$ to the tidal deformability $\Lambda_{A}$ and $C_{\rm Oc}^{A}$ to $C_{\rm Q}^{A}$ using the EOS-insensitive relations (Tables 1 and 2 from Yagi and Yunes (2017)):

[TABLE]

and

[TABLE]

To allow a better interpretation of the spin-spin terms discussed above, we present in Fig. 7 the individual contributions $\hat{\psi}^{(A)}_{\text{SS, 2PN}}$, $\hat{\psi}^{(A)}_{\text{SS, 3PN}}$, and $\hat{\psi}^{(A)}_{\text{SS, 3.5PN}}$. In addition, for better visibility, we also show explicitly the spin-cubed octupole term $\hat{\psi}^{(A)}_{C_{\rm Oc},3.5\text{PN}}=-440\hat{C}_{\rm Oc}^{A}X_{A}^{3}\chi_{A}^{3}$. For an equal mass setup with $\Lambda^{A,B}=350$ and $\chi^{A,B}=0.2$ the 2PN contribution dominates up to $\hat{\omega}\sim 0.06$, before the positive 3PN term becomes larger. Overall, we find that except the 3PN contribution all terms are negative for the chosen setup. We also see that throughout the inspiral the octupole term is about 1 order of magnitude smaller than other contributions. This observation remains valid even for spins close to break-up $\chi\sim 0.75$. Thus, we do not attempt to include additional higher order multipoles.

III.4 Precession dynamics

We conclude the discussion of the model by shortly describing the incorporation of precession. The precession dynamics in IMRPhenomPv2_NRTidalv2 is included as in the previous IMRPhenomPv2_NRTidal approach Dietrich et al. (2019). For this we assume that the spin-orbit coupling can be approximately separated into components parallel and perpendicular to the instantaneous orbital angular momentum, where the component perpendicular to the orbital angular momentum is driving the precessional motion Kidder (1995); Apostolatos et al. (1994); Buonanno et al. (2003); Schmidt et al. (2011, 2012, 2015).

Consequently, we construct a precessing tidal waveform approximant from the spin-aligned model after adding all tidal corrections to the underlying spin-aligned point particle model. We then rotate the waveform to account for precession, as discussed in Refs. Schmidt et al. (2012, 2015).

IV Validation

IV.1 Time domain comparison with NR simulations

As a first validation check, we compute the time domain phase difference between IMRPhenomPv2_NRTidalv2 and a selected set of NR data; see Table 3. All of the employed waveforms are publicly available in the CoRe database (www.computational-relativity.org Dietrich et al. (2018)). In addition to IMRPhenomPv2_NRTidalv2, we also present the phase difference with respect to SEOBNRv4T and IMRPhenomPv2_NRTidal in Fig. 8.

Waveform alignment: For comparison, we align all waveforms with respect to the NR data by minimizing the phase difference in the time interval $[t_{i},t_{f}]$

[TABLE]

where $\rm x$ denotes the individual waveform approximant. The alignment windows are marked by vertical dashed lines in Fig. 8.

NR data uncertainty: For a quantitative comparison with respect to the NR data, we assign each dataset with an uncertainty, where we generally distinguish between (i) setups employing the high-order flux scheme of Bernuzzi and Dietrich (2016) for which clean convergence is found throughout the inspiral, and (ii) setups whose behavior is monotonic, but no clean convergence is present. For the setups employing the high-order flux scheme, we obtain a better phase estimate and an error measure (green shaded region) due to Richardson extrapolation Bernuzzi and Dietrich (2016); Dietrich et al. (2018); cf. Sec. II.2. Other configurations are marked by blue shaded regions. For these cases, the uncertainty due to numerical discretization is estimated by the difference between the two highest resolutions, which is not necessarily a conservative error estimate. For both scenarios, we also include an error measuring the effect of the finite radius extraction of the GW from the numerical domain. This error measure is obtained by computing the difference in the waveform’s phase with respect to different extraction radii; see e.g. Refs. Bernuzzi and Dietrich (2016); Dietrich and Hinderer (2017) for a more detailed discussion.

NRTidalv2 dephasing: Considering the performance of the NRTidalv2 approximation, we find that for all cases with reliable error measure (green shaded regions), the dephasing between the model and the NR data is well within the error estimate and never exceeds $1\ \rm rad$. The performance is comparable with the SEOBNRv4T model which is shown as a blue dashed-dotted line.888We note that very recently Ref. Lackey et al. (2019) has constructed a reduced order model of SEOBNRv4T which can also be used directly for parameter estimation. Considering the difference with respect to IMRPhenomPv2_NRTidal, we find that as expected the new NRTidalv2 model is less attractive, which is caused by the slightly different behavior in the frequency range $\hat{\omega}\geq 0.05$.

For the NR setups which show no clear convergence throughout the inspiral, we find that for most cases the estimated uncertainty is larger than the phase difference between the NRTidalv2 and the NR data, the exceptions are BAM:0081 and BAM:0094. These setups are characterized by high mass ratios [BAM:0094 is to date the NR dataset with the largest simulated mass ratio ($q=2.1$)] and tidal deformabilities which are in tension with the observation of GW170817 Abbott et al. (2018).

Additional simulations with clean convergence for large mass ratios are needed to allow an overall improvement of BNS models in these regions of the parameter space (see Appendix B).

IV.2 Mismatch Computations with respect to EOB-NR hybrids

To validate the new NRTidalv2 model, we compare our LALSuite implementation against a set of target waveforms combining TEOBResumS and NR data by computing the mismatch. Those waveforms have been constructed for Ref. Dietrich et al. (2019) and are publicly available under www.computational-relativity.org Dietrich et al. (2018). We refer to Dietrich et al. (2019) for further details. The main properties of these target waveforms are summarized in Table 4.

Mismatch computation: We compute the mismatch according to

[TABLE]

where $\phi_{c},t_{c}$ are an arbitrary phase and time shift. The noise-weighted overlap is given by

[TABLE]

where tildes denote the Fourier transform, $S_{n}(f)$ is the spectral density of the detector noise, and $f$ is the GW frequency (in the frequency domain). We used the Advanced LIGO zero-detuning, high-power (ZERO_DET_high_P) noise curve of adl for our analysis999We note that this noise curve has recently been updated slightly upd , but for consistency with Ref. Dietrich et al. (2019) we employ the old noise curve. with a fixed $f_{\rm min}=30\,{\rm Hz}$ and a variable $f_{\rm max}$ ranging from 500 Hz up to the merger frequency ($f_{\rm mrg}$) reported in Table 4.

Mismatch with respect to hybrid waveforms: We compute the mismatch for $18$ TEOBResumS-NR hybrid waveforms (Table 4) against a range of different phenomenological models: IMRPhenomD Husa et al. (2016); Khan et al. (2016) (no tidal effects), IMRPhenomD_NRTidal (incorporating tidal effects using the NRTidal model of Dietrich et al. (2019) but no quadrupole-monopole self-spin terms), IMRPhenomPv2_NRTidal incorporating tidal effects using the NRTidal model of Dietrich et al. (2019) including quadrupole-monopole self-spin terms up to 3PN, and the new model IMRPhenomPv2_NRTidalv2. In addition, we include the new SEOBNRv4_ROM_NRTidalv2 and IMRPhenomD_NRTidalv2 approximants. We evaluate the waveform models at the parameters of the hybrids reported in Table 4 with an initial frequency of $30\,{\rm Hz}$.

Generally, we find that IMRPhenomPv2_NRTidalv2 performs as well or better than IMRPhenomPv2_NRTidal, except for 2 cases. For all configurations the mismatch stays below $5\times 10^{-3}$ even for maximum frequencies at or above the merger frequency. In addition, our comparisons show again that the inclusion of the quadrupole-monopole terms is important even for astrophysically reasonable spins—see Harry and Hinderer (2018); Dietrich et al. (2019); Samajdar and Dietrich (2019) for previous studies. In most cases the mismatches between the hybrids and IMRPhenomD_NRTidalv2 are marginally smaller compared to SEOBNRv4_ROM_NRTidalv2. Even less notable are the differences between IMRPhenomD_NRTidalv2 and IMRPhenomPv2_NRTidalv2 which are dominantly driven by the additional 3.5PN spin-spin and cubic-in-spin contributions in IMRPhenomPv2_NRTidalv2. The additional tidal amplitude corrections have almost a negligible effect; cf. the non-spinning configurations in Fig. 9.

Our comparison shows that the TEOBResumS-NR hybrids are well described by the new approximant and that no additional pathologies (in the low frequency regime) are introduced during recalibration.

IV.3 Cross validation against SEOBNRv4T

As a final check of the approximant, we compute the mismatch between the IMRPhenomPv2_NRTidalv2 and the SEOBNRv4T model for a number of randomly sampled configurations. We compare these mismatches with mismatches between IMRPhenomPv2 and IMRPhenomPv2_NRTidalv2 to give an impression of the importance of tidal effects. The match computation is restricted to the frequency interval of $f\in[40,2048]\ \rm Hz$. We have tested starting frequencies of $25\ \rm Hz$ and $30\ \rm Hz$ for a smaller number of cases and obtained smaller mismatches than for the $40\ \rm Hz$ initial frequency. Therefore, to save computational costs and to provide a conservative estimate, we use a minimum frequency of $40\ \rm Hz$.

IV.3.1 Non-spinning Configurations

We start this analysis by considering non-spinning configurations. For this purpose, we select $1000$ samples with flat priors in $M_{A},M_{B}\in[1,3]\,M_{\odot}$ and $\Lambda_{A},\Lambda_{B}\in[0,5000]$. The final analysis is shown in Fig. 10, where we compare the IMRPhenomPv2_NRTidalv2 and SEOBNRv4T approximant. For non-spinning configurations the mismatches between IMRPhenomPv2_NRTidalv2 and SEOBNRv4T are below $0.034$ for our set of configurations. The largest difference is found for large mass ratios; cf. upper left and lower right corners of the top panel. For better visualization, we mark mass ratios of $q=1.25;1.5;2.0$ by diagonal gray, dark gray, and black lines, respectively. Restricting to mass ratios below $1.5$, we find a largest mismatch of $\bar{F}=0.024$

In addition, our analysis shows that for larger tidal deformabilities the mismatch between the two models tends to increase; cf. upper right corner of the right panel in Fig. 10. We mark in the plot $\Lambda_{A}+\Lambda_{B}=1250;2500;5000$ with gray, dark gray, and black lines. Restricting our analysis to $\Lambda_{A}+\Lambda_{B}<2500$ leads to a maximum mismatch of $\bar{F}=0.016$.

Overall, the average mismatch between IMRPhenomPv2_NRTidalv2 and SEOBNRv4T for our dataset is $0.009$. Interestingly, if we restrict our analysis to the more physical parameter space in which the more massive star has the smaller tidal deformability,101010For equations of state with no phase transition, the dimensionless tidal deformability is a monotonically decreasing function of the star’s mass—see, e.g., Fig. 1 in Chatziioannou et al. (2015) for an illustration. However, in cases with a phase transition that yield twin stars, the tidal deformability is no longer a monotonically decreasing (or even a single-valued) function of mass, as illustrated in, e.g., Refs. Sieniawska et al. (2019); Han and Steiner (2019). Of course, even in twin star cases, the deviations from monotonic decrease and single-valuedness are not large, while the parameters we generate by the aforementioned random sampling can have significant violations. the average mismatch decreases by roughly a factor of 2 to $0.0059$.

Consequently, we find for non-spinning configurations a good agreement between the tidal EOB model SEOBNRv4T and IMRPhenomPv2_NRTidalv2.111111We note that as cross-validation of the implementation of the other approximants, we also tested the mismatch between IMRPhenomPv2_NRTidalv2 and SEOBNv4_ROM_NRTidalv2 and find an average mismatch of $\sim 5\times 10^{-4}$.

IV.3.2 Spinning Configurations

We further consider spinning configurations using flat priors $M_{A,B}\in[1,3]\,M_{\odot}$, $\Lambda_{A,B}\in[0,5000]$, and $\chi_{A,B}\in[-0.30,0.30]$ as well as $\chi_{A,B}\in[-0.60,0.60]$. For both prior choices we select $3000$ randomly distributed samples.

If we consider spins within $\chi_{A,B}\in[-0.30,0.30]$ (upper panels of Fig. 11), we find a maximum mismatch of $\bar{F}=0.034$, which is comparable with the non-spinning result presented before. Overall, the average mismatch of our $3000$ samples for spins within $\chi_{A,B}\in[-0.30,0.30]$ is $\bar{F}=0.0072$. For the same set of configurations, the average mismatch with respect to IMRPhenomPv2_NRTidal is $\bar{F}=0.0092$, i.e., $25\%$ larger. Furthermore, we find that, as for the non-spinning cases, the largest mismatches are obtained for configurations which have large mass ratios and large tidal deformabilities. If only spin magnitudes up to $|\chi_{A,B}|\leq 0.3$ are considered, we do not find a noticeable spin effect.

However, spin effects become important for large spin magnitudes. For spin magnitudes up to $0.6$, the largest mismatches between SEOBNRv4T and IMRPhenomPv2_NRTidalv2 are found for large anti-aligned spins, i.e., the lower left corner of the right-most bottom panel of Fig. 11. The maximum mismatch is ${\bar{F}=0.167}$ for our randomly chosen set of configurations.

Comparing average values, we find that while the average mismatch is $0.043$ between the original IMRPhenomPv2_NRTidal model and SEOBNRv4T, the average mismatch decreases to $0.021$ between the IMRPhenomPv2_NRTidalv2 model and SEOBNRv4T, i.e., much better agreement is found within this large region of the parameter space.

The disagreement between SEOBNRv4T and IMRPhenomPv2_NRTidalv2 for large anti-aligned spins needs further investigation and requires additional NR simulations in regions of the parameter space which are currently not covered. Note that for the largest anti-aligned spin, high-quality NR setups (CoRe:BAM:0062) the NSs have only a spin of $\chi_{A,B}=-0.10$. For this physical configuration, both waveform approximants (SEOBNRv4T and IMRPhenomPv2_NRTidalv2) describe the data within the estimated uncertainty.

V Summary

In this article we have presented our most recent update of the NRTidal model. The model gives a closed analytical expression for tidal effects during the BNS coalescence and can be added to an arbitrary BBH baseline approximant. We added the new NRTidalv2 approximant to IMRPhenomPv2 Hannam et al. (2014); Khan et al. (2016) to obtain a frequency-domain precessing BNS approximant as well as to the (frequency-domain) SEOBNRv4_ROM Bohé et al. (2017) and the IMRPhenomD Khan et al. (2016) approximants to allow an improved and fast modeling of spin-aligned systems.

Our main improvements in comparison to the initial NRTidal model are:

1. (i)

a recalibration of the tidal phase to improved NR data incorporating additional analytical knowledge for the low frequency limit; 2. (ii)

the addition of a tidal amplitude correction to the model; 3. (iii)

incorporation of higher order (3.5PN) quadrupole and octupole information to the spin sector of the model.

We also hope to further improve the NRTidalv2 model for higher mass ratios to allow an accurate description of high mass ratio systems. Such an extension requires additional high-quality NR simulations for a variety of different mass ratios.

An additional improvement would be the incorporation of the effect of $f$-mode resonances as recently computed in Refs. Andersson and Pnigouras (2019); Schmidt and Hinderer (2019), the incorporation of an updated precession dynamics as used in Khan et al. (2018), or the incorporation of higher modes London et al. (2018); Cotesta et al. (2018).

We have compared the IMRPhenomPv2_NRTidalv2 model with high resolution numerical relativity data and found agreement within the estimated uncertainty for all NR data with clear convergence. Overall, the performance of IMRPhenomPv2_NRTidalv2 is comparable with state-of-the-art tidal EOB models.

This accuracy was verified by the mismatch computation between IMRPhenomPv2_NRTidalv2 and TEOBResumS-NR hybrid waveforms, for which mismatches are well below $5\times 10^{-3}$.

We concluded the performance test of the model with a mismatch computation with respect to the tidal EOB model SEOBNRv4T. For non-spinning cases (or cases with small spins as employed in the low spin prior of the LVC analysis) the mismatch computed from a starting frequency of $40\ \rm Hz$ never exceeds $\bar{F}\approx 0.034$ for $M_{A,B}\in[1,3]\,M_{\odot}$ and $\Lambda_{A,B}\in[0,5000]$. Considering spinning setups ($\chi_{A,B}\in[-0.6,0.6]$), the mismatch increases to a maximum of $\bar{F}=0.164$.

Acknowledgements.

We acknowledge fruitful discussions with Sebastiano Bernuzzi, Tanja Hinderer, Alessandro Nagar, Patricia Schmidt, Ka Wa Tsang, and Chris Van Den Broeck and thank Alessandra Buonanno, Michael Pürrer, and Ulrich Sperhake for comments on the manuscript. We thank Sarp Akcay and Nestor Ortiz for the BAM:0001 high resolution data. T. D. acknowledges support by the European Union’s Horizon 2020 research and innovation program under grant agreement No 749145, BNSmergers. A. S. and T. D. are supported by the research programme of the Netherlands Organisation for Scientific Research (NWO). S. K. acknowledges support from the Max Planck Society’s Independent Research Group Grant. N. K. J.-M. acknowledges support from STFC Consolidator Grant No. ST/L000636/1. Also, this work has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement No. 690904. R. D. acknowledges support from DFG grant BR 2176/5-1. W. T. was supported by the National Science Foundation under grant PHY-1707227. The computation of the numerical relativity waveform was performed on the Minerva cluster of the Max-Planck Institute for Gravitational Physics.

Appendix A Alternative formulation of a tidal amplitude correction

A.1 Tidal amplitude corrections in the time domain

As in the frequency domain, the BNS waveform time domain amplitude can be obtained by augmenting existing BBH-models with additional tidal corrections $A_{\rm T}$, i.e.,

[TABLE]

Refs. Damour et al. (2012); Banihashemi and Vines (2018) present the tidal amplitude corrections for the leading and next-to-leading order

[TABLE]

While Eq. (34) describes the tidal amplitude corrections for small frequencies, it loses validity close to the moment of merger. For extreme cases, i.e., stiff EOSs and low NS masses, the additional amplitude corrections can become larger than $A_{\rm BBH}$ causing the overall amplitude to be negative. Thus, a further calibration to NR or EOB data is required. We will employ the quasi-universal relations, which allow an EOS-independent description of important quantities at the moment of merger (merger frequency, merger amplitude, reduced binding energy, specific orbital angular momentum, and GW luminosity) Bernuzzi et al. (2014); Takami et al. (2015); Dietrich (2016); Zappa et al. (2018). As shown in Fig. 6.7 of Ref. Dietrich (2016), the GW amplitude at the merger follows a quasi-universal relation as a function of the tidal coupling constant

[TABLE]

namely,

[TABLE]

We note that a straightforward extension of Eq. (36) [Eq. (6.15d) of Ref. Dietrich (2016)], would be the incorporation of a larger number of NR simulations as publicly available under www.computational-relativity.org, Ref. Dietrich et al. (2018). However, we postpone this to future work Zappa et al. , in which a more general discussion about quasi-universal relations during the BNS coalescence will be given.

To incorporate Eq. (36) in Eq. (34), we extend the analytical knowledge with an additional, unknown higher order PN-term and define the NRTidal amplitude correction as

[TABLE]

where the individual terms $\hat{c}_{A},\hat{c}_{1}^{A},\hat{c}_{B},\hat{c}_{1}^{B}$ can be obtianed from Eq. (34) once we express $\Lambda_{A,B}$ in terms of $\kappa_{A,B}$.

Enforcing

[TABLE]

gives us the unknown parameter $d$ according to

[TABLE]

with $\Delta A=A^{\rm mrg}(x^{\rm mrg})-A_{\rm BBH}(x^{\rm mrg})$.

We note that although the outlined approach has been tested for a selected number of cases, we did not implemented it in LALSuite due to the large computational costs inherent to time domain waveform approximants.

A.2 Frequency domain amplitude corrections by SPA

In addition to the frequency domain amplitude correction presented in the main text, we also want to present a possible alternative way to augment the frequency domain binary black hole amplitude with tidal correction. For this purpose we use the SPA to obtain the frequency-domain amplitude from Eq. (37). Following the SPA approach,

[TABLE]

where $\ddot{\phi}$ refers to the second time derivative of $\phi$.

Using Eq. (16) [see also Eq. (14) in Damour et al. (2012); note that the published version is missing an equals sign], Eq. (40) can be rewritten as

[TABLE]

Inserting $A=A_{\rm BBH}+A_{\rm T}$ and $\psi=\psi_{\rm BBH}+\psi_{T}$ leads to

[TABLE]

Treating the tidal phase correction as a small change of the underlying BBH waveform, we rewrite the expression as

[TABLE]

Linearizing in $\epsilon=\mathcal{O}(\kappa)$ and neglecting terms proportional to $\kappa^{2}$ [noting that $A_{\rm T}=\mathcal{O}(\kappa)$] leads to

[TABLE]

Thus, the final expression is given as

[TABLE]

The approach outlined in this appendix, i.e., Eq. (45), leads to much larger computational cost than the Padé approximant [Eq. (24)], which is why we chose the easier and more straightforward implementation shown in the main body of the paper. However, this additional approach might become relevant for a potential improvement/extension in the future.

Appendix B Extension of the NRTidal phase incorporating additional mass ratio dependence

We now outline a possible extension of the NRTidalv2 model which incorporates additional analytically known mass-ratio dependence. Since we do not find such an extension to perform better in our tests and to reduce computational costs, we limited the mass ratio dependence in the tidal phase simply to the prefactor $\propto\kappa_{\rm eff}^{T}/\nu$ in the current implementation of the model.

However, it is possible to recast Eq. (17) as

[TABLE]

with $c_{0}^{A}=-\frac{3}{16X_{A}X_{B}}\left(12+\frac{X_{A}}{X_{B}}\right)$ (and similarly for $c_{0}^{B}$).

The individual Padé approximants $\tilde{P}_{\rm NRTidal}^{A}(x)$, $\tilde{P}_{\rm NRTidal}^{B}(x)$ are similar to Eq. (18) together with the constraints in Eqs. (19), but the known PN coefficients $\tilde{c}_{1}^{A,B},\tilde{c}_{3/2}^{A,B},\tilde{c}_{2}^{A,B},\tilde{c}_{5/2}^{A,B}$ have a mass ratio dependence as given in Ref. Damour et al. (2012).

Due to the limited set of high-quality NR data, the fitting coefficients in Eqs. (20) ($\tilde{d}_{1},\tilde{d}_{2},\tilde{n}_{5/2},\tilde{n}_{3}$) can only be determined for the equal mass case.

We find that while such a choice of the coefficients leads to a correct mass ratio dependence for the low frequency limit, the higher frequency phase is described worse compared to the NRTidalv2 approximation given in the main text. We suggest that this is caused by the inconsistency introduced by adding the mass-ratio dependence in only some of the Padé coefficients.