Comprehensive numerical relativity -- effective-one-body comparison to inform improvements in waveform models for binary neutron star systems

TL;DR

This paper compares various tidal effective-one-body models with high-precision numerical relativity simulations for binary neutron star systems, identifying strengths and limitations to improve gravitational waveform modeling.

Contribution

It provides the most extensive comparison to date, evaluating multiple EOB models against numerical relativity across a broad parameter space.

Findings

Dynamical tides EOB model best matches simulations

All models underestimate tidal effects for stiff equations or small masses

Discrepancies increase near merger, highlighting model limitations

Abstract

We present a detailed comparison between tidal effective-one-body (EOB) models and new state-of-the-art numerical relativity simulations for non-spinning binary neutron star systems. This comparison is the most extensive one to date, covering a wide range in the parameter space and encompassing the energetics of the binary, the periastron advance, the time and frequency evolution of the gravitational wave phase for the dominant mode, and several subdominant modes. We consider different EOB models with tidal effects that have been proposed, including the model with dynamical tides of [Phys.Rev.Lett. 116 (2016) no.18, 181101] and the gravitational self-force (GSF) inspired tidal EOB model of [Phys.Rev.Lett. 114 (2015) no.16, 161103]. The EOB model with dynamical tides leads to the best representation of the systems considered here, however, the differences to the GSF-inspired model are small. A common feature is that for systems where matter effects are large, i.e. stiff equations of state or small total masses, all EOB models underestimate the tidal effects and differences to the results from numerical relativity simulations become noticeable near the merger. We analyze this regime to diagnose the shortcomings of the models in the late inspiral, where the two neutron stars are no longer isolated bodies moving in vacuum. Our work will serve to guide further advances in modeling these systems.

Tables

Table 1. Table 1: Configurations. The first column defines the name of the configuration. The subsequent 9 columns describe the corresponding physical properties: EOS, gravitational mass of the individual stars M A , B superscript 𝑀 𝐴 𝐵 M^{A,B} , baryonic mass of the individual stars M b A , B superscript subscript 𝑀 𝑏 𝐴 𝐵 M_{b}^{A,B} , stars’ compactnesses 𝒞 A , B superscript 𝒞 𝐴 𝐵 \mathcal{C}^{A,B} , the total dimensionless quadrupolar tidal coupling constant κ 2 T subscript superscript 𝜅 𝑇 2 \kappa^{T}_{2} , the initial dimensionless GW frequency M ​ ω 22 0 𝑀 superscript subscript 𝜔 22 0 M\omega_{22}^{0} , the ADM mass M ADM subscript 𝑀 ADM M_{\rm ADM} and angular momentum J ADM subscript 𝐽 ADM J_{\rm ADM} from the initial data. The last two columns describe the resolutions employed and whether the outermost refinement level consists of Cartesian boxes or spherical shells, cf.  Hilditch et al. ( 2013 ); Dietrich et al. ( 2016a ) .

Name	EOS	$M^A$	${𝒞}^A$	$M^B$	${𝒞}^B$	${\kappa }_2^T$	$M{\omega }_{22}^0$	$M_{\mathrm{ADM}}$	$J_{\mathrm{ADM}}$	Resolutions	Grid Config.
MS1b-138138	MS1b	1.375008	0.1445	1.375008	0.1445	262	0.035	2.728593	8.15792	(R1,R2,R3)	shell
SLy-138138	SLy	1.375004	0.1771	1.375004	0.1771	65	0.036	2.727632	8.00631	(R1,R2,R3)	shell
MS1b-153122	MS1b	1.527797	0.1591	1.222228	0.1299	267	0.035	2.728700	8.03219	(R1,R2,R3)	shell
SLy-153122	SLy	1.527388	0.1976	1.222225	0.1573	69	0.036	2.728071	7.93386	(R1,R2,R3)	box
MS1b-165110	MS1b	1.650025	0.1708	1.100000	0.1182	282	0.035	2.729230	7.79917	(R1,R2,R3)	shell
SLy-165110	SLy	1.650019	0.2149	1.097885	0.1416	78	0.036	2.726620	7.70012	(R1,R2,R3)	box
MS1b-150100	MS1b	1.500016	0.1565	1.000001	0.1085	460	0.030	2.482511	6.65711	(R1,R2,R3)	shell
MS1b-150100-ecc	MS1b	1.500016	0.1565	1.000001	0.1085	460	0.030	2.482230	6.61058	(R1)	box
SLy-150100	SLy	1.500009	0.1938	1.000007	0.1293	138	0.031	2.482086	6.58692	(R1,R2,R3)	box

Equations

$$H_{\mathrm{EOB}}=M\sqrt{1+2\nu \left(H_{\mathrm{eff}}-1\right)}-M,$$

$$H_{\mathrm{eff}}^2=p_{r\ast }^2+A\left(1+\frac{p_{\varphi }^2}{r^2}+2\left(4-3\nu \right)\frac{p_{r\ast }^4}{\nu r^2}\right),$$

$$\frac{dr}{dt}$$

$$\frac{d{p}_{r\ast }}{dt}$$

$$\frac{d\varphi }{dt}=\frac{\partial H_{\mathrm{EOB}}}{\partial p_{\varphi }},\frac{d{p}_{\varphi }}{dt}={\mathcal{F}}_{\varphi }.$$

$${\cal F}_{\phi}=-\frac{1}{\Omega_{\phi}}\mathcal{\dot{E}}_{\rm rad},\ \ \ \ \ {\cal F}_{r}=\frac{p_{r*}}{p_{\phi}}{\cal F}_{\phi},\ \ \ \ \ \$$

$${\dot{\mathcal{E}}}_{\mathrm{rad}}=\frac{{\Omega }_{\varphi }^2}{8\pi }\sum _{\ell =2}^8\sum _{m=0}^{\ell }{m}^2|h_{\ell m}^{\mathrm{F}}{|}^2.$$

$${h_{\ell m}^{\mathrm{F}}}_{\mathrm{pp}}=h_{\ell m}^{\mathrm{N}}{S}_{\ell m}{T}_{\ell m}{\rho }_{\ell m}^{\ell }{N}_{\ell m}.$$

$$A^{\text{tidal}}$$

$${\alpha }_{\ell }^{(A)}=\frac{2X_B}{X_A}{k}_{\ell }^{(A)}{R}_A^{2\ell +1}{u}^{2\ell +2}$$

$${\hat{A}}_{\mathrm{adPN}}^{(A)}=1+\frac{5}{2}{X}_{A}u+\left(\frac{337}{28}{X}_A^2+\frac{1}{8}{X}_A+3\right)u^2$$

$$h_{\ell m}^{\mathrm{F}}={h_{\ell m}^{\mathrm{F}}}_{\mathrm{pp}}+h_{\ell m}^{\mathrm{tidal}}.$$

$${\hat{k}}_{\ell \mathrm{dyn}}$$

$${\mathcal{Q}}_{\ell m}$$

$$\hat{t}=\frac{8}{5\sqrt{{ϵ}_m}}\left(1-\frac{r^{5/2}{\omega }_{0\ell }^{5/3}}{m^{5/3}{M}^{5/6}}\right),$$

$${ϵ}_m=\frac{256\mu M^{2/3}{\omega }_{0\ell }^{5/3}}{5m^{5/3}}.$$

$${\hat{k}}_{2\mathrm{dyn}}^{(A)\mathrm{diss}}=\frac{{\hat{k}}_{2\mathrm{dyn}}^{(A)}\left({\omega }_{0\ell }^2+6X_B{\Omega }^2\right)}{3{\Omega }^2\left(1+2X_B\right)},$$

$${\hat{A}}_{\mathrm{adGS}{\mathrm{F}}^{+}}^{(A)}$$

$$\tilde{A}(r_{\mathrm{LR}})-\frac{r_{\mathrm{LR}}}{2}{\tilde{A}}^{\prime }(r_{\mathrm{LR}})=0,$$

$${\kappa }_2^T:=2\left(\frac{q^4}{(1+q{)}^5}\frac{k_2^A}{C_A^5}+\frac{q}{(1+q{)}^5}\frac{k_2^B}{C_B^5}\right),$$

$$M=M^A+M^B$$

$$q=M^A/M^B.$$

$$E=(H_{\mathrm{EOB}}-M)/\mu$$

$$\ell =p_{\varphi }$$

$$E$$

$$\ell$$

$${\mathcal{E}}_{\mathrm{rad}}$$

$${\mathcal{J}}_{z\mathrm{rad}}$$

$$\Delta E_{\mathrm{NR}}(\ell )=\sqrt{\Delta E_{\Delta M}^2+\Delta E_{\Delta J}^2(\ell )+\Delta E_r^2(\ell )+\Delta E_{\Delta x}^2(\ell )}.$$

$$K^{-2}=D^{-1}\bigg{[}A\left(3+\frac{rA^{\prime\prime}}{A^{\prime}}\right)-2rA^{\prime}\bigg{]},$$

$$K=\frac{{\Omega }_{\varphi }}{{\Omega }_r}(t_k)=\frac{{\varphi }_{\Omega }(t_{k+1})-{\varphi }_{\Omega }(t_k)}{2\pi }=\frac{{\varphi }_{22}(t_{k+1})-{\varphi }_{22}(t_k)}{4\pi },$$

$$K_{\mathrm{NR}}(M\Omega )=\frac{c_0+c_1(M\Omega )+c_2(M\Omega {)}^2}{\sqrt{1-6(M\Omega {)}^{\frac{2}{3}}}},$$

Full text

I Introduction

The recent first detection of gravitational waves (GWs) of merging binary black holes (BBHs) Abbott et al. (2016a, b) has initiated a new observational era in astronomy. Within the next few years GW detections of neutron star (NS) binaries are also anticipated Aasi et al. (2016); Abbott et al. (2016c). These observations will have a unique potential to probe the equation of state (EOS) of the ultradense matter in NS interiors Read et al. (2009a); Damour et al. (2012a); Read et al. (2013); Del Pozzo et al. (2013); Lackey and Wade (2015); Flanagan and Hinderer (2008); Hinderer (2008); Damour and Nagar (2009); Binnington and Poisson (2009); Hinderer et al. (2010), which remains one of the major unsolved problems in modern astrophysics.

To successfully extract the EOS information from the GW data requires highly accurate waveform models (templates) in order to perform matched-filtered searches, where the data are cross-correlated with a template bank covering all physical values of the parameters. Constructing sufficiently accurate models requires a detailed understanding of late stages of the binary NS (BNS) coalescence, where matter effects on the GW signal become most prominent. In this regime, one must solve Einstein’s equations together with the equations of general relativistic hydrodynamics, a challenging task that can only be accomplished through numerical relativity (NR) simulations. Over the last few years the physical realism of the NR simulations of BNSs has been significantly advanced by including a refined treatment of magnetic fields and microphysical processes, see e.g. Giacomazzo et al. (2011); Sekiguchi (2010); Sekiguchi et al. (2011); Galeazzi et al. (2013); Kiuchi et al. (2014); Dionysopoulou et al. (2015); Palenzuela et al. (2015); Lehner et al. (2016); Sekiguchi et al. (2015); Foucart et al. (2016); Sekiguchi et al. (2016), and by improving the numerical accuracy of the predicted GW signals Bernuzzi et al. (2012a); Radice et al. (2014a, b); Hotokezaka et al. (2013); Radice et al. (2015); Hotokezaka et al. (2015); Bernuzzi and Dietrich (2016).

Unfortunately, NR simulations are too expensive to create template banks for BNS GW signals. Furthermore, analytical descriptions such as the Post-Newtonian (PN) approach become inaccurate close to the merger and are therefore inadequate for analyzing the regime where tidal effects are largest Favata (2014); Yagi and Yunes (2014); Lackey and Wade (2015). A successful way to recast PN results and combine them with information from NR simulations is the effective-one-body (EOB) formalism Buonanno and Damour (1999, 2000). This framework has been refined to devise an accurate, semi-analytical description of the dynamics and GW signals of coalescing BHs with arbitrary spins and mass ratios Taracchini et al. (2014); Pan et al. (2014); Damour and Nagar (2014). Tidal effects have also been included in the EOB model, starting with the first analysis Damour and Nagar (2010) that used Newtonian and partial 1PN tidal information Vines and Flanagan (2010) and was compared against NR simulations in Baiotti et al. (2010). Further theoretical developments including 2PN tidal effects in the conservative dynamics Bini et al. (2012) and 1PN tidal corrections to the waveform amplitudes Damour et al. (2012a) have led to an improved tidal EOB model that was tested against results from improved NR simulations with robust error estimates in Bernuzzi et al. (2012b); Hotokezaka et al. (2013). More recent work has focused on further advancing the EOB description of tidal effects in different ways. The work of Ref. Bini and Damour (2014a) suggested an EOB model based on the Schwarzschild tidal field together with gravitational self-force (GSF) corrections. The pure GSF-based model underestimates matter effects in the late inspiral, but was augmented by an additional term to develop a model that remains close to NR predictions Bernuzzi et al. (2015a); Hotokezaka et al. (2015). Another line of recent work has focused on the enhancement of tidal effects due to the NS’s finite fundamental ($f-$) mode frequency that also leads to an improved performance of the model in comparisons with NR results Hinderer et al. (2016); Steinhoff et al. (2016).

To date, EOB – NR comparisons for BNSs have been limited to equal mass systems. Although nearly equal mass systems seem to be common Swiggum et al. (2015); Lattimer (2012), recent observations have also discovered systems with mass ratios of $q=M^{A}/M^{B}\sim 1.3$ Martinez et al. (2015); Lazarus et al. (2016). The range of pulsar masses that have been accurately determined ($1.17-2.01M_{\odot}$) and ‘in situ’ models, see e.g. Dominik et al. (2012); Dietrich et al. (2015a), also suggest a larger range of the mass ratio. Therefore, in this paper20170207, we go beyond previous work in covering the BNS parameter space by considering mass ratios in the range $q=1.0-1.5$. Since for unequal-mass systems spherical harmonic modes in the GWs besides the dominant $(\ell,m)=(2,2)$ mode are larger than for equal masses, we also compare the EOB predictions for the (2,1), (3,3), and (4,4) modes with NR simulations. In addition, we perform a comprehensive study of other important diagnostic quantities to test the performance of the EOB models and gain additional insights. Specifically, we investigate the phase evolution of the (2,2) mode as a function of the GW frequency, analyze the periastron advance (PA), and study the dynamics of the binary by tracing the binding energy vs. specific orbital angular momentum curve over the inspiral.

In general, as also found in previous EOB-NR comparisons, for most scenarios tidal effects in NR simulations are larger than those predicted by EOB models. The origin of this trend is twofold. On the one hand, tidal EOB models generally underestimate tidal effects present in the last stages of the BNS coalescence because they lack higher PN order tidal terms and physical effects that become relevant once the NSs come into contact shortly before the merger. On the other hand NR simulations tend to overestimate tidal effects since numerical dissipation accelerates the inspiral similar to matter effects. This shows that a reliable error estimate is crucial for a proper numerical relativity- analytical relativity comparison of BNSs.

The paper20170207 is organized as follows. In Sec. II we briefly review the tidal EOB models used in this article. Sec. III describes the configurations we investigate (see in particular Tab. 1) and briefly summarizes the NR methods employed in the simulations. We then discuss quantities characterizing the dynamics by considering binding energy vs. angular momentum curves in Sec. IV and the PA in Sec. V. In Sec. VI we compare the EOB and NR waveforms. Finally, we discuss current problems of the EOB models during the last stages of the inspiral for two particular cases and outline aspects requiring further improvements to the EOB model in Sec. VII, where we also show results from combining the models of Bernuzzi et al. (2015a) and Hinderer et al. (2016) as a possibility to further enhance tidal effects in the late inspiral regime. We conclude in Sec. VIII.

Throughout this article, we use geometrical units $c=G=M_{\odot}=1$ unless otherwise stated.

II Tidal Effective-One-Body models

The EOB model Buonanno and Damour (1999, 2000) meshes results from the PN approximation, valid for any mass ratio but limited to weak gravitational fields, with knowledge of strong-field effects from the test-particle limit. The analytical results are further improved by calibrating parameterized higher-order PN terms in the model to NR data to obtain an accurate description of GW signals from BBH systems Taracchini et al. (2014); Pan et al. (2014); Damour and Nagar (2014). The conservative EOB dynamics are described by the Hamiltonian

[TABLE]

with $M=M^{A}+M^{B}$ being the total mass and $\nu=M^{A}M^{B}/M^{2}$ the symmetric mass ratio. We use $A$ and $B$ to label the two bodies. The effective Hamiltonian $H_{\rm eff}$ describes a particle of mass $\mu\equiv\nu M$ moving in an effective spacetime. For nonspinning binaries the motion is in a plane and the effective Hamiltonian is given by

[TABLE]

where $p_{\phi}$ is the canonical azimuthal angular momentum per unit reduced mass, $p_{r*}=p_{r}/\sqrt{D}$ is a rescaled radial momentum, and $A$ and $D$ are potentials describing the effective spacetime. The EOB potentials considered in this paper20170207 have the form $A=A_{\rm pp}+A_{\rm tidal}$ and $D=D_{\rm pp}$, where ”${\rm pp}$” denotes the point-mass part that is used for BBH evolutions (EOBBBH). For $A_{\rm pp}$ we use the result given in Eqs. (A1)-(A2h) of Ref. Steinhoff et al. (2016) with the calibration parameter $K=1.712-1.804\nu-39.77\nu^{2}+103.166\nu^{3}$ determined in Taracchini et al. (2014). The potential $D$ is taken from Eq. (A4) of Ref. Steinhoff et al. (2016). These choices correspond to the model SEOBNRv2111Note that $A_{\rm pp}$ used for the EOB${}_{\rm adGSF^{+}}$ model introduced in Sec. II.3 is slightly different than described above since the public available code of EOB is used, see also Barausse et al. (2012); Bini and Damour (2014b, 2013); Bernuzzi et al. (2015a). The effect of the underlying point particle model on the results will be described in more detail in Appendix C..

We evolve the EOB dynamics using the equations of motion

[TABLE]

The factor $A/\sqrt{D}$ arises from the fact that $p_{r*}$ and $r$ are not canonically conjugate variables. The gravitational radiation reaction forces ${\cal F}_{\phi}$ and ${\cal F}_{r}$ are constructed from the energy flux through

[TABLE]

where the energy flux is computed from

[TABLE]

The sum here is only over positive $m$ since $|h_{\ell-m}^{\rm F}|=|h_{\ell m}^{\rm F}|$. The factorized EOB waveforms for point masses are computed from

[TABLE]

The various factors in Eq. (10), are described in detail in Taracchini et al. (2012) and given explicitly in Eqs. (16), (17), and (A1)-(A14h) therein.

II.1 Adiabatic tides with 2PN Taylor expanded potential: EOBadPN

Adiabatic tidal effects arise when the distorted NS remains in hydrostatic equilibrium. In this case the ratio of its induced multipolar deformation to the perturbing tidal field is a constant related to its Love number $k_{\ell}$. The deformation of the NSs changes the energy of the system, which can be captured in the EOB approach by including tidal contributions to the EOB potentials of the form

[TABLE]

Here $\alpha_{2}^{(A)}$ is the Newtonian tidal potential given by

[TABLE]

where $X_{A,B}=M^{A,B}/M$, $R_{A}$ is the radius, and $u=1/r$. The 2PN corrections to the tidal potential $\hat{A}^{(A)}$ from Eqs. (6.6) and (6.18) of Ref. Bini et al. (2012) are given by

[TABLE]

Tidal effects also influence the dissipative sector since the tidal bulges contribute to the gravitational radiation. In the EOB model this is accounted for by adding to the waveform modes of Eq. (6) a tidal contribution $h_{\ell m}^{\rm tidal}$ to the GW modes so that

[TABLE]

The explicit results for the adiabatic tidal terms $h_{\ell m}^{\rm tidal}$ were worked out in Eqs. (A14)-(A17) of Ref. Damour et al. (2012a) for $l\leq 3$. We use only those terms from Ref. Damour et al. (2012a) for which the PN knowledge of tidal effects is complete.

II.2 Dynamical tides with 2PN Taylor expanded potential: EOBdyPN

Dynamical tides become important when the tidal forcing frequency due to the companion’s orbital motion approaches an eigenfrequency of the NS’s normal modes of oscillation. The inclusion of the effects of dynamic tides from the NS’s $f-$modes of oscillation is discussed in detail in Ref. Steinhoff et al. (2016). Here, we consider only the most practical yet approximate of the models devised therein, where the potential from Eqs. (7) and (9) is used but with $k_{\ell}$ multiplied by a frequency-dependent enhancement factor such that $k_{\ell}\to k_{\ell}\hat{k}_{\ell\ {\rm dyn}}$ with

[TABLE]

Here, $\Omega^{\prime}=3/8$, $\Omega=M^{1/2}r^{-3/2}$, and $\omega_{0\ell}$ are the $f-$mode angular frequencies of the NS, where the subscript [math] indicates that the mode amplitude has no radial nodes. The quantity ${\cal Q}_{\ell m}$ is

[TABLE]

where the functions ${\rm F}_{\rm S}$ and ${\rm F}_{\rm C}$ are Fresnel sine and cosine integrals respectively using the conventions in Mathematica. We have defined

[TABLE]

and

[TABLE]

In Eq. (11) a body label $A,B$ on the quantities $\omega_{0\ell}$, $\hat{t}$, $\epsilon_{m}$, and ${\cal Q}_{\ell m}$ is implied. For each $\ell$-multipole only $m=\ell$ contributes in Eq. (11) because the effect of modes with $m<\ell$ has already been taken into account as adiabatic contributions. For the lowest order multipoles the coefficients are given by $(a_{2},a_{3},a_{4})=(\frac{1}{4},\frac{3}{8},\frac{29}{64})$ and $(b_{2},b_{3},b_{4})=(\frac{3}{4},\frac{5}{8},\frac{35}{64})$.

Similar to the treatment for the conservative dynamics, the effect of dynamic $f-$mode tides can be incorporated in the dissipative sector in an approximate way by multiplying the occurences of $k_{\ell}$ in $h_{\ell m}^{\rm tidal}$ in Eq. (10) by an effective function $k_{\ell}\to k_{\ell}\ \hat{k}_{\ell\ {\rm dyn}}^{\rm diss}$. For $\ell=2$ this function is given by

[TABLE]

where $\hat{k}_{2\ {\rm dyn}}^{(A)}$ is the enhancement function for body $A$ in the conservative dynamics from Eq. (11). We find that similar factors for other modes and from higher multipoles have only a small impact on the results and therefore neglect them in this study; they are further discussed in Ref. Taracchini, Andrea and others (2017).

II.3 Adiabatic tides with Gravitational Self-Force-inspired potential:

EOB${}_{\rm adGSF^{+}}$

Gravitational self-force calculations have recently computed tidal invariants that contain information about strong-field tidal effects in the limit of small mass ratios, to linear order in $X_{A}$ Dolan et al. (2014); Bini and Damour (2014a). These results have been augmented in Ref. Bini and Damour (2014a); Bernuzzi et al. (2015a) by a term $\propto X_{A}^{2}$ that would describe currently unknown second-order self-force effects. Specifically, in this model Eq. (7) is employed with the potential $\hat{A}^{(A)}$ given by

[TABLE]

The caveats with the choice of gauge for mapping GSF results to EOB potentials that contain explicit poles at the light ring $r_{\rm LR}$ are discussed in Steinhoff et al. (2016). In Eq. (16), the coefficient $a_{1}^{\rm GSF}(u)$ is obtained from Eqs. (7.24)–(7.27) of Ref. Bini and Damour (2014a). In the model of Ref. Bernuzzi et al. (2015a) the unknown parameter $a_{2}^{\rm 2GSF}$ is set to $a_{2}^{\rm 2GSF}(u)=337u^{2}/28$ and the unknown exponent $p$ is chosen to be $p=4$. The radius of the light ring $r_{\rm LR}$ is obtained from the conservative EOB dynamics by solving

[TABLE]

with the potential $\tilde{A}=A_{\rm pp}+A^{\rm tidal}_{\rm adPN}$ from the PN model in Eqs. (7) and (9).

For tests involving the EOB${}_{\rm adGSF^{+}}$ model we use the publicly available code of EOB that was also used to create the surrogate model in Lackey et al. (2016). Note that this code differs from the EOB model described above in the point-particle sector: (i) non-quasi-circular (NQC) corrections (the factor $N_{\ell m}$ in the waveform modes of Eq. (6) that is necessary to accurately describe the BBH merger) are omitted, (ii) a different resummation of the potentials $A_{\rm pp}$ and $D$ is used, (iii) the arguments of factors in $h_{\ell m}^{\rm F}$ such as $\rho_{\ell m}$ involve powers of $v_{\phi}=(\partial H_{\rm EOB}/\partial p_{\phi})^{-2/3}\Omega$ evaluated for circular orbits instead of $v=\Omega^{1/3}$, (iv) the higher-order logarithms in the tail terms are Pade-resummed, and (v) the calibration of the BBH model is different. The influence of these differences play only a marginal role compared to the tidal effects and how those are incorporated in the EOB model, therefore, we will present results obtained with EOB for EOB${}_{\rm adGSF^{+}}$. However, more details about the influence of the underlying point mass models will be given in Appendix C.

III Binary configurations and numerical methods

In this work we consider eight different configurations that have already been studied in Dietrich et al. (2016a) with lower resolution. The initial configurations are constructed with SGRID, see Tichy (2009, 2012); Dietrich et al. (2015a) for detailed descriptions. Dynamical evolutions are performed with the BAM code, see e.g. Brügmann et al. (2008); Thierfelder et al. (2011); Dietrich et al. (2015b).

We choose two zero-temperature EOSs modeled by piecewise polytropes Read et al. (2009b): SLy Douchin and Haensel (2001) and MS1b Müller and Serot (1996). The EOSs are chosen to cover a large range in the BNS parameter space with respect to the tidal coupling constant

[TABLE]

which characterizes Newtonian tidal effects in the EOB Hamiltonian. For the dynamical evolution we add an additional component to the pressure to describe thermal effects $p_{\rm th}=\epsilon\rho(\Gamma_{\rm th}-1)$ with $\Gamma_{\rm th}=1.75$ Bauswein et al. (2010). For both EOSs we focus on four different setups defined by their total mass

[TABLE]

and their mass ratio

[TABLE]

The parameters for our setups are $(M,q)=$ $(2.75,1.0)$, $(2.75,1.25)$, $(2.75,1.50)$, $(2.50,1.50)$. To illustrate the effect of changes in the parameters on the merger process we show in Fig. 1 the 2D density snapshots of these systems at merger, which corresponds to the peak in the GW amplitude. The density snapshots are gauge-dependent quantities, therefore, Fig. 1 enables only a qualitative interpretation.

Building on results from our previous work Dietrich et al. (2016a), we have performed new simulations of the setups at higher resolution. Thus, three different resolutions denoted by (R1,R2,R3) are available for each configuration. The resolutions on the finest grid covering the NSs are $(0.23,0.15,0.12)$ for the SLy setups and $(0.25,0.16,0.13)$ for MS1b. With these resolutions the diameters of the stars are covered by approx. 64, 96, and 128 grid points, respectively. Details about the configurations are given in Tab. 1.

To compute results for the periastron advance we consider as a representative case the setup MS1b-150100. Reliably extracting this information from NR data, as described in Appendix A, requires an additional simulation with higher eccentricity that we perform following the approach of Moldenhauer et al. (2014); Dietrich et al. (2016a).

IV Dynamics

To study the dynamics and energetics of BNS mergers, we present results for the binding energy vs. orbital angular momentum of the binary. They were proposed as important diagnostic tools to assess the performance of EOB models in Damour et al. (2012b) and were later applied for BNS systems in Bernuzzi et al. (2014a).

IV.1 Binding energy curves in the EOB model

For the binding energy curves we evaluate

[TABLE]

and

[TABLE]

using the solutions to the evolution equations for $(r,p_{r*},p_{\phi})$ and obtain the results parametrically.

IV.2 NR error budget

We compute the dimensionless binding energy and angular momentum per reduced mass as

[TABLE]

respectively. Here, $M_{\rm ADM}$ and $J_{\rm ADM}$ describe the ADM-mass and ADM-angular momentum, and $\mathcal{E}_{\rm rad}$ and $\mathcal{J}_{\rm rad}$ are the energy and angular momentum emitted from the system via GWs. The radiated energy and angular momentum perpendicular to the orbital plane are calculated as

[TABLE]

with $l_{\rm max}=8$.

The binding energy vs. orbital angular momentum curves are directly accessible from our numerical simulations, i.e. no alignment or additional shift is needed. This makes $E(\ell)$ curves an unambiguous and effective tool for a comparisons with analytical models.

To assess the uncertainty of $E(\ell)$ for our simulations we consider four different sources of error:

1. (i)

Uncertainty of the binding energy at $t=0$. Because of the limited number of points in SGRID and the interpolation of the SGRID initial data on the BAM grid, the initial ADM-mass as well as the individual component masses are only accurate up to $\sim 10^{-5}$. We account for this uncertainty through $\Delta E_{\Delta M}=2\times 10^{-5}$. 2. (ii)

Uncertainty of the angular momentum at $t=0$. Similar to the ADM mass, the angular momentum also contains numerical uncertainties of the order of $\Delta J_{\rm ADM}(t)\sim 10^{-3}$. This leads to $\Delta\ell(t=0)=\Delta J_{\rm ADM}(t=0)/(M^{A}M^{B})\sim 10^{-3}$, where we estimate the uncertainty conservatively as $\Delta\ell=2\times 10^{-3}$. Since $\Delta\ell$ leads to a horizontal shift in the $E(\ell)$ curve we assign a time dependent error on the binding energy of $\Delta E_{\Delta J}(\ell)=E(\ell)-E(\ell-\Delta\ell)$. 3. (iii)

Finite radius extraction effects. Another source of error is caused by the finite size of computational domain and the finite extraction radii. As shown in Bernuzzi and Dietrich (2016) this error is small ($\lesssim 10^{-1}\rm rad$). To estimate the resulting uncertainty we use two extraction radii, $r_{1}=1000M$ and $r_{2}=750M$, and compute the difference in $E(\ell)$, i.e. $\Delta E_{r}(\ell)=E_{r_{1}}(\ell)-E_{r_{2}}(\ell)$. In future work with more accurate NR data the waveform has to be extrapolated to infinity, but for the scope of this paper20170207 the simple estimate is sufficient. 4. (iv)

Numerical discretization. We estimate this error by considering different resolutions. As shown in Bernuzzi et al. (2012a); Bernuzzi and Dietrich (2016) a conservative error estimate for the triplet of resolutions employed in this article is to compute the difference between resolution R3 and R2 $\Delta E_{\Delta x}(\ell)=E_{\rm R3}(\ell)-E_{\rm R2}(\ell)$. We find that $\Delta E_{\Delta x}(\ell)$ dominates the overall error in most cases.

The total uncertainty of the numerical data is approximated by

[TABLE]

IV.3 Results for $E(\ell)$

IV.3.1 General features

Our results are summarized in Fig. 2, where $E(\ell)$ is shown for the NR data (black), the EOBBBH model (green), and the EOBdyPN model (orange) in the top panels and the difference between the NR results and the EOB models in the bottom panels. The ordering of the panels is the same as in Fig. 1.

To aid in the interpretation of our results we first discuss important general aspects of the $E(\ell)$ curves. Because of the attractive nature of tidal effects BNS systems are more bound than BBH systems once systems with the same specific angular momentum are compared. Consequently, the binding energy for a given angular momentum is smaller for BNS setups than for the corresponding BBH case 222Note that for systems with the same frequency the opposite statement holds, i.e., the binding energy at a given frequency is larger for BBHs than for BNSs see the description in Dietrich et al. (2016b)..

The influence of the EOS on $E(\ell)$ is deduced by contrasting the MS1b and SLy curves in adjacent panels of Fig. 2. We observe that for stiffer EOS (MS1b) the merger occurs at larger orbital angular momentum and larger binding energy, so the merger remnant is less bound for MS1b than for SLy. This effect is also visible in Fig. 1, where MS1b setups merge at larger distances because of their lower compactness. The influence of the mass ratio is illustrated in Fig. 2 by comparing panels with the same EOS (but excluding the 150100 setups that have a lower total mass). We see that for increasing mass ratio the moment of merger also occurs at higher binding energy and angular momentum. Finally, the bottom panels with the same EOS in Fig. 2 show that an increase in the total mass of the system also leads to a more bound merger remnant. Most of these findings can be explained by the quasi-universal properties of BNS mergers reported in Bernuzzi et al. (2014b, 2015b); Takami et al. (2015) and the influence of $\kappa_{2}^{T}$ on the inspiral and merger dynamics.

IV.3.2 EOB-NR comparison

Comparing the EOB and NR curves in Fig. 2 shows that for most setups the tidal EOB models deviate from the NR predictions shortly before merger. This disagreement between NR and EOB is larger for stiffer EOSs, smaller total masses, and higher mass ratios, i.e. setups with larger tidal effects. Furthermore, it is evident that for all simulations EOBdyPN and EOB${}_{\rm adGSF^{+}}$ stay closer to the NR result than EOBadPN. This observation is independent of the EOS, mass, and mass ratio. The origin of the larger tidal effects for EOBdyPN and EOB${}_{\rm adGSF^{+}}$ compared to EOBadPN is different for each model: for EOBdyPN the dynamical tides enlarge $k_{\ell}$ when the orbital frequency approaches the resonance frequency of the individual stars, while enhanced tides for EOB${}_{\rm adGSF^{+}}$ are caused by approaching the light ring in the GSF+-potential, Eq. (16).

Another interesting aspect is the location of the merger point of the individual models: except for the EOB${}_{\rm adGSF^{+}}$ model in the cases SLy-138138 and SLy-153122 (where the EOB${}_{\rm adGSF^{+}}$ merger corresponds to the endpoint of the blue dashed curve) the merger for all EOB models occurs at smaller angular momentum and binding energy than in the NR simulations. This is interesting since it shows that in most cases the final merger remnant is less bound than predicted by the EOB models and that additional repulsive effects that are not included in the current EOB description become important near merger.

We emphasize that the NR error quoted here is chosen conservatively, i.e., deviations between the NR and the EOB results outside the NR error bars mark region where the EOB models must be further improved. But note that due to numerical dissipation, which acts as an artificially attractive force, the continuum solution will have smaller tidal effects than the NR results presented here and will therefore be closer to the EOB models. Thus, the correct solution will lie within the lower part of the shaded regions in the bottom panels of Fig. 2.

V Periastron Advance

V.1 Periastron Advance in the EOB models

The PA is characterized by the ratio of angular and radial frequencies $K=\Omega_{\phi}/\Omega_{r}$ and, in the case considered here, is influenced by a combination of relativistic and tidal effects. In the limit of circular orbits, as discussed e.g. in Le Tiec et al. (2011); Hinderer et al. (2013), the PA for nonspinning binaries and for the classes of tidal models considered here can be computed from

[TABLE]

V.2 Extracting the Periastron Advance from NR Simulations

The computation of the PA for NR simulations relies on a suitable indicator for the eccentricity. In this work, we test three different indicators $e_{\Omega},e_{\omega},e_{\phi}$ based on the orbital frequency $\Omega$, the GW frequency of the (2,2)-mode $\omega_{22}$, and the GW phase $\phi_{22}$ of the (2,2)-mode. A detailed description is given in Appendix A.

From the eccentricity indicator we calculate the PA following the work of Mroue et al. (2010). We determine the times $t_{k}$ of the extrema of $e$ corresponding to perihelia and aphelia and compute the orbital phase $\phi_{\Omega}$ and GW phase $\phi_{22}$ at these times $t_{k}$. Finally the ratio between the angular and radial frequency is given by

[TABLE]

where each equal sign corresponds to a different way of extracting $K$ from the NR data.

V.3 EOB-NR comparison

Extracting $K$ from NR simulations is only possible for sufficiently long inspirals. Furthermore, as shown in e.g. Mroue et al. (2010), it is difficult to determine $K$ for simulations with almost vanishing eccentricity. For these reasons we focus on the MS1b-150100 setup which has the lowest starting frequency. Further we perform a new simulation with higher eccentricity: MS1b-150100-ecc.

Figure 3 shows $K$ for all eccentricity measurements for MS1b-150100 in the upper panel, and the results for MS1b-150100-ecc in the middle panel333We have simulated MS1b-150100-ecc only for resolution R1. The resolution has only a small influence on the measurement of the PA, since resolution effects become dominant at later times of the simulation, where no eccentricity can be measured at all.. We find that the scattering of data points obtained from $K_{\rm NR}(e_{\phi})$ is smallest and that these results are in agreement with $K_{\rm NR}(e_{\omega})$. By contrast the results from $K_{\rm NR}(e_{\Omega})$ are more scattered and do not help to put constraints on $K$. Thus, we combine all data points obtained with $K_{\rm NR}(e_{\phi})$ and $K_{\rm NR}(e_{\omega})$ and fit the NR data for both configurations according to

[TABLE]

cf. Hinderer et al. (2013). We indicate the $1\sigma$ and $2\sigma$ region of the fit as dark and light shaded regions respectively in Fig. 3.

The bottom panel of Fig. 3 illustrates that the NR data contains larger tidal effects than those predicted by the EOB models. Because of the difficulties in reliably computing the NR error for the PA (as discussed in Appendix A, see also Hinderer et al. (2013) for uncertainty estimates of PA in BBH setups) the results presented here have to be considered with caution. However, it is clearly visible from Fig. 3 that for low frequencies $M\Omega<0.019$ the PA predicted from the EOB models and the results from the NR simulations agree very well, cf. dark shaded region. At higher frequencies the curves start to deviate and the point-mass result EOBBBH becomes distinguishable from the NR uncertainty region. The EOBBBH result leaves the $1\sigma$ region at $M\Omega\approx 0.019$ and the $2\sigma$ region at $M\Omega\approx 0.021$. This demonstrates for the first time that current NR simulations are able to resolve tidal effects on the PA. Up to frequencies of order $M\Omega\approx 0.023$ the results from the EOB models EOBdyPN and EOB${}_{\rm adGSF^{+}}$ remain within the $2\sigma$ region of the NR data. Beyond this frequency no additional NR data points (see upper and middle panel) are available and, thus, no further statement to discriminate between different tidal EOB models can be made.

VI Waveforms

VI.1 EOB waveforms

The EOB waveforms are computed by evaluating the functions $h_{\ell m}^{\rm F}$ in Eq. (10). They depend on the EOB trajectory and enter into the backreaction of GW losses on the dynamics through Eq. (4). To compute the time evolution of the modes we therefore solve the algebraic equations (10) jointly with the EOB equations of motion.

VI.2 NR error budget

Gravitational waves are extracted using the Newman-Penrose formalism Newman and R. (1962) from the curvature scalar $\Psi_{4}$. Individual modes can be obtained by a multipole decomposition on extraction spheres using spherical harmonics with spin weight $-2$. The algorithm and implementation is described in detail in Brügmann et al. (2008). We reconstruct the metric multipoles $rh_{lm}$ from the curvature multipoles $r\Psi_{4\ lm}$ using the frequency domain integration of Reisswig and Pollney (2011) with a cutting frequency of $m\omega_{22}^{0}/2$, see Tab. 1.

VI.2.1 $\phi(t)$-Evolution

The NR waveforms are presented as a function of the retarded time

[TABLE]

For the extraction radius we use the value $r_{\rm extr}=1000$. As shown in Bernuzzi and Dietrich (2016) for this radius the error due to using waveforms extracted at a finite radius is of order $\sim 0.1$rad and decreases during the simulation because the GW frequency increases. As a simple error estimate we compute $\Delta\phi_{r}=\phi_{r_{1}}-\phi_{r_{2}}$ with $r_{1}=1000$ and $r_{2}=750$, similar to the investigation of the binding energy discussed above.

In most simulations the dominant source of uncertainty are truncation errors during the last few orbits before merger. We estimate the uncertainty due to finite resolution by the difference between setups R3 and R2, $\Delta\phi_{\Delta x}=\phi_{\rm R3}-\phi_{\rm R2}$. As shown in Bernuzzi et al. (2012a); Bernuzzi and Dietrich (2016) for these resolutions this is a conservative measure of the error.

Finally, we also take into account that mass is not exactly conserved in the simulations. In particular, as discussed in Dietrich et al. (2015b), the treatment of the atmosphere can cause mass loss or gain and is a source of systematic errors. The dephasing caused by variations in the mass is taken into account as $\Delta\phi_{\Delta M}=\omega u\Delta M_{b}/M_{b}$, where $M_{b}$ denotes the baryonic mass of the system at the beginning of the simulation and $\Delta M_{b}$ is the change in the baryonic mass.

The total phase uncertainty is then given by

[TABLE]

We take into account that some errors are non-monotonic by using ${\rm{max}}_{\tilde{u}\in[0,u]}[...]$. This ensures that the total error is always increasing over time 444Note that we do not include the influence of the residual eccentricity in our error budget since this would require the use of eccentricity reduced waveforms. But as presented in Dietrich et al. (2015a) the remaining eccentricity can cause phase differences of the order of $\sim 0.2$rad..

VI.2.2 $\phi(\omega)$- Evolution

In addition to the time evolution of the GW phase, we also present the phase accumulated over a given frequency interval. For all setups we compute the accumulated phase with respect to the reference frequency $M\omega_{22}=0.04$, which corresponds to $470\rm Hz$ for $M=2.75$ and to $517\rm Hz$ for $M=2.50$. The details of the computation are discussed in Appendix B.

The error budget of $\phi(\omega)$ contains the following components, similar to the analysis of the time evolution $\phi(t)$:

1. (i)

We take into account the fact that the GWs are extracted at finite radii by considering $\Delta\phi_{r}(\omega)=\phi_{r1}(\omega)-\phi_{r2}(\omega)$ with $r_{1}=1000$ and $r_{2}=750$. 2. (ii)

The uncertainty caused by numerical discretization is estimated as $\Delta\phi_{\Delta x}(\omega)=\phi_{\rm R3}(\omega)-\phi_{\rm R2}(\omega)$. 3. (iii)

We include an estimate resulting from the non-conservation of the total mass $\Delta\phi_{\Delta M}(\omega)$. 4. (iv)

An additional source of error is the uncertainty caused by the low-pass filtering of the raw NR data (Appendix B). We account for this by varying the cutoff-frequency by a factor of $2$ and obtain the estimate $\Delta\phi_{\Delta\omega_{\rm cut}}(\omega)=\phi_{\omega_{\rm cut}}(\omega)-\phi_{2\cdot\omega_{\rm cut}}(\omega)$. 5. (v)

A further contribution to the error budget arises from the fact that the phase evolution is computed with respect to a phase at $M\tilde{\omega}_{22}=0.04$. The uncertainty in this reference value is estimated to be $\Delta\phi_{M\tilde{\omega}}(\omega)=\phi_{M\tilde{\omega}}(\omega)-\phi_{M\tilde{\omega}\pm 2.5\times 10^{-4}}(\omega)$.

The total error is given as:

[TABLE]

VI.3 EOB-NR comparison: $\phi(t)$

As shown in e.g. Dietrich et al. (2016a) more than $99\%$ of the total energy emitted in a BNS merger up to a mass ratio of $q=1.5$ is contained in the (2,2) mode. This motivated previous EOB and NR comparisons to focus solely on this mode. However, we want to go beyond previous work and include the (2,1), (3,3), and (4,4) mode in our comparison.

By contrast to the comparisons of results for $E(\ell)$, $K$, and $\phi(\omega)$ discussed below, we do not include the EOBBBH model in our analysis of $\phi(t)$. The reason is that a robust comparison of the dephasing over time requires aligning two waveforms over a time interval where they agree. However, our simulations cover only the last $\sim 10$ orbits before merger, where tidal effects cannot be neglected and a reliable alignment of BBH and BNS waveforms cannot be achieved.

VI.3.1 Dominant (2,2)-mode

Figure 4 presents the results for the time evolution of the dominant (2,2)-mode. In the upper panels the real part of $rh_{22}$ is shown for the NR data (black) and the EOBdyPN model (orange dashed), we also include the amplitude $|rh_{22}|$ as thin lines. The time interval used for the alignment is marked by vertical gray dash-dotted lines555For large $\kappa_{2}^{T}$ changing the alignment interval effects $\Delta\phi(t)$ and a more robust comparison would require longer NR simulations.. In the bottom panel we present the difference between the NR waveform and the tidal EOB models. The moment of merger in the NR simulation is shown as a vertical black dashed line in the bottom panels. The numerical uncertainty is shown as a shaded region.

We conclude from Fig. 4 that:

1. (i)

For setups with a small tidal coupling constant the dephasing between EOB and NR waveforms is small and remains within the numerical uncertainty almost up to the merger. 2. (ii)

Variations in the mass ratio have a negligible effect on the performance of the EOB models when considering $\phi(t)$. 3. (iii)

For setups with stiffer EOSs the EOB models start to deviate from the NR simulation a few orbits before the merger. 4. (iv)

In all cases the EOB models underpredict tidal effects, i.e. $\phi_{\rm EOB}(u)<\phi_{\rm NR}(u)$. 5. (v)

The EOB merger generally occurs after the NR merger. 6. (vi)

For all setups the performance of the tidal EOB models is similar, with the EOBdyPN and EOB${}_{\rm adGSF^{+}}$ models remaining closer to the NR results than the EOBadPN model.

VI.3.2 Higher modes

While for equal mass configurations the (2,1) and (3,3) modes are zero, they are non-zero for unequal masses. The energy emitted in the (2,1) and (3,3) mode increases almost linearly for an increasing mass ratio. By contrast the amplitude and consequently the emitted energy of the (4,4) mode is almost independent of the mass ratio, see Dietrich et al. (2016a) for further details. We present the (4,4) mode for all models employing the MS1b EOS in Fig. 5666Note that we do not study systems employing the SLy EOS since those have (as presented in Dietrich et al. (2016a)) a larger drift of the center of mass, which effects the computation of the subdominant modes.. The (2,1) and (3,3) mode for the unequal mass setups with MS1b EOS are shown in Fig. 6. We used the alignment based on the $(2,2)$ mode, where we obtained the time and phase shifts $(\delta t^{\rm align},\delta\phi_{22}^{\rm align})$. This determines the shifts for the other $(\ell,m)$ modes to be $\delta\phi_{\ell m}^{\rm align}=m\delta\phi_{22}^{\rm align})/2$, together with the same time shift. Another consideration is that the point-mass EOB model for higher modes has not been calibrated to NR results, e.g. there is no equivalent function to $N_{22}$. As shown in Appendix C, the influence of such calibrations of the EOB model is $\lesssim 0.03$rad for the MS1b setups and thus much smaller than the size of tidal effects.

Overall we find that the dephasing for the higher modes is similar to the results obtained for the dominant (2,2) mode. The performance of the EOB models depends only weakly on the mass ratio since $q$ has only a small influence on the tidal coupling constant. By contrast, the total mass impacts the accuracy of the EOB models since tidal effects strongly decrease with increasing total mass. For example, for MS1b-150100 where $\kappa_{2}^{T}$ is $\sim 180$ larger than for MS1b-165110 the dephasing from NR is larger since the larger tidal effects are inadequately captured by the models.

Considering the agreement of the amplitude between the EOB models and the NR results, we find that close to the moment of merger the amplitude predicted by the EOB model is smaller than those for the NR simulation. This is true for all setups and is a robust (qualitative) observation. However, since the computation of the subdominant modes is challenging and the numerical simulations are more inaccurate close to the merger a more quantitative investigation is postponed to future work, where simulations with higher resolution and waveforms corrected for the center of mass drift will be available.

VI.4 EOB-NR comparison: $\phi(\omega)$

Another important means to assess the performance of EOB models is the phase evolution as a function of the frequency. The results of this comparison are shown in Fig. 7. The top panels exhibit the accumulated phase $\phi(\omega)$ with respect to the reference frequency $M\omega_{22}=0.04$ computed for the EOBBBH model (green), the EOBdyPN model (orange), and the NR data (black). The difference between the NR and the EOB results is shown in the bottom panels. In contrast to the small tidal contributions in the time evolution of the phase $\phi(t)$, the accumulated phase as a function of frequency due to tidal effects is of the order of $10$ to $40$rad for our setups, which enables an easier assessment of tidal effects. Furthermore, because no alignment is needed for comparisons of $\phi(\omega)$, the point-mass results EOBBBH can readily be included in the comparison.

The qualitative behavior of $\phi(\omega)$-plots is the following. Because the frequency increases during the inspiral, the accumulated phase within a given frequency interval $\Delta(M\omega)$ decreases, i.e. $\partial_{\omega}^{2}\phi(\omega)<0$. The effect of the EOSs (seen in Fig. 7 by comparing SLy and MS1b configurations) is that MS1b setups merge at lower frequency due to their smaller compactness. This is in agreement with the conclusions drawn from the results of the binding energy in Sec. IV. Considering the difference $\phi_{\rm NR}-\phi_{\rm EOB_{\rm BBH}}$, shown as the green curves in the lower panels of Fig. 7, we see that smaller total masses lead to a larger accumulated phase due to the larger tidal effects. The oscillatory behavior of $\Delta\phi$ for low frequencies is caused by the residual eccentricity that is not entirely eliminated by the low-pass filter, cf. Appendix B.

In contrast to the results for $\phi(t)$ shown in Fig. 4 the results for $\phi(\omega)$ in Fig. 7 clearly exhibit a substantial improvement of the EOBdyPN and EOB${}_{\rm adGSF^{+}}$ models over the EOBadPN model that has a larger discrepancy to the NR results. As described before the agreement between EOB and NR is generally better for softer EOSs for which tidal effects are smaller. But when considering the evolution of $\phi(\omega)$, where differences are more pronounced, it is apparent that several models are outside the NR uncertainty. Only for the cases SLy-153122 and SLy-165110 are the EOBdyPN results in agreement with the NR waveforms to within the current uncertainty in NR data. Interestingly, we observe that for some setups the dephasing of the EOB${}_{\rm adGSF^{+}}$ model decreases close to the merger, and, in one case terminates before attaining the merger frequency. These features indicate that the EOB${}_{\rm adGSF^{+}}$ model likely overestimates tidal effects in the late inspiral for these particular cases.

VII Current problems and future prospects

VII.1 Test cases

Having presented the results of comparing EOB predictions to NR results for a variety of systems in the previous sections, we next perform a more detailed analysis to identify the source of the discrepancies. For simplicity, we focus on two complementary cases at opposite corners in the parameter space considered here: SLy-138138, where both stars have a soft EOS and equal masses, and MS1b-165110, where both stars have a stiff EOS and the mass ratio is $q=1.5$.

VII.1.1 SLy-138138

In Fig. 8 we summarize important quantities for SLy-138138. In the top panels we show three snapshots of the rest-mass density $\rho$ obtained from the NR simulation, where the color bar is identical to that used in Fig. 1. The snapshots correspond to times $t=1500M$ (blue), $2500M$ (green) and the moment of merger $2675M$ (red). The other panels show the real part and amplitude of the GW, the dephasing of the EOB models with respect to the NR simulation, the difference of the binding energy between the EOB models and the NR simulation, and the dephasing as a function of the frequency. The retarded time, angular momentum, and frequency corresponding the density snapshots are marked as columns of the same color in the plots.

For setup SLy-138138 both stars have a soft EOS and the same gravitational mass of $M=1.375$. The tidal coupling constant at $t=0$ is $\kappa_{2}^{T}=65$. Within our data set this configuration has the smallest tidal effects. Consequently, differences between the EOB models and the NR simulation are also small.

At time $t=1500M$ (left density snapshot, blue columns) the two stars are well-separated and no tidal deformation is visible. The low density regions ($\rho\lesssim 10^{-8}$, dark shades of blue) in the plot are affected by the artificial atmosphere in the NR treatment and no clear star surface exists. The imprint of the artificial atmosphere decreases for increasing resolution. At this time the tidal EOB models are in agreement with the NR simulation for all quantities considered.

As the system evolves and the separation decreases tidal effects become larger. At $t=2500M$ (middle density snapshot, green columns) the tidal deformation of the stars is visible. At this stage the GW phasing predicted by the EOBadPN model starts to deviate from the NR data as seen in the plots of $\Delta\phi(t)$ and $\Delta\phi(\omega)$. This indicates that the EOBadPN model underestimates tidal effects.

Approximately $1.5$ orbits later the GW amplitude reaches its maximum and the stars merge (red frame/columns). Because of the enhanced tidal effects in EOB${}_{\rm adGSF^{+}}$ and EOBdyPN these models are able to describe the GW phase as a function of time up to the last GW cycle before the merger. We note three important observations: (i) The model EOB${}_{\rm adGSF^{+}}$ predicts a merger at a higher specific orbital angular momentum (endpoint of the blue dashed curve in the fourth panel) than the model EOBdyPN. The merger of the NR simulation lies between both models. (ii) The dephasing $\Delta\phi_{\rm EOB_{\rm adGSF^{+}}}(\omega)$ decreases shortly before merger, which is caused by the fact that the binary separation approaches the light ring and tidal effects are overestimated at this stage. (iii) Since tidal effects are overall small, both models EOB${}_{\rm adGSF^{+}}$ and EOBdyPN give a reasonably accurate representation of the dynamics and waveform. Further improvements of tidal EOB models in this part of the parameter space would require more accurate NR simulations.

VII.1.2 MS1b-165110

Figure 9 illustrates the results for the system MS1b-165110 in a similar way as Fig. 8 does for SLy-138138.

At time $t=1500M$ (left density snapshot) both stars are clearly separated, however, because of the small compactness of the lower-mass star and a stiff EOS, tidal effects are already visible at this early stage of the simulation. Furthermore, the lower density regions are more affected by the artificial atmosphere, showing that the simulation of MS1b-165110 is less accurate than SLy-138138 and leading to larger NR uncertainties. We find that at $t=1500M$ all EOB models are able to describe the dynamics and GW emission of this system.

Up to a time $t\approx 2300M$ all tidal EOB models remain within the NR error estimate and hence consistent with the simulation. At this time the NSs are still separated. Some of the low density material in the NR simulation is transferred from the secondary to the more massive NS, however, as shown in Dietrich et al. (2016a) this mass transfer is in part a numerical artifact that decreases with increasing resolution. Later, at $t=2500M$ or about two orbits before the merger (middle density snapshot) tidal deformations are significant and the EOB models deviate from the NR models.

Around merger the lower-mass star is tidally disrupted by the primary star (right density snapshot). Tidal effects are significantly larger than for SLy-138138. The dephasing caused by tidal effects ($\Delta\phi_{{\rm EOB}_{\rm BBH}}(\omega)$) is of the order of 30 radiants, i.e. more than a factor of two larger than for SLy-138138 (although the system merges at a smaller GW frequency). Overall, during the last four GW cycles before the merger none of the EOB models give a representation that is in agreement with the NR uncertainty. Neither the enhanced tidal effects due to the GSF+-potential nor the dynamical tides of EOBdyPN are strong enough for an accurate modeling, although again the EOBdyPN performs better than the EOB${}_{\rm adGSF^{+}}$ model. Our results suggest that an additional enhancement of tidal effects in the EOB description for such regions of the BNS parameter space is needed.

VII.2 Towards improved EOB models

The studies presented in this paper20170207 show that for the region of the BNS parameter space and length of NR simulations considered here both the EOBdyPN and EOB${}_{\rm adGSF^{+}}$ tidal EOB models provide a robust description of BNSs during most of the inspiral, however, our NR simulations reveal weaknesses of the models in the last few orbits before the merger. In particular, for cases where tidal effects are large i.e. for large $\kappa_{2}^{T}$, the EOB models tend to underestimate tidal effects irrespective of the mass ratio. To quantify the impact of these differences on measurements with advanced GW detectors requires parameter estimation studies that are the subject of ongoing work, see also Lackey et al. (2014); Wade et al. (2014); Hotokezaka et al. (2016). Nevertheless, understanding and improving waveform models in the regime close to merger will be essential for future GW detectors. Thus, we want to mention three possibilities for advancing existing models in this part of the BNS parameter space.

VII.2.1 Combining Dynamical tides with the GSF-inspired potential: EOB${}_{\rm dyGSF^{+}}$

As discussed above, the current description of tidal effects in the EOBdyPN model tends to underestimate tidal effects in the late inspiral. To illustrate the effect of a simple enhancement of tidal contributions near the merger, we exchange the PN $A$-Potential with the GSF inspired potential introduced in Sec. II.3. This leads to a combination of EOB${}_{\rm adGSF^{+}}$ and EOBdyPN: EOB${}_{\rm dyGSF^{+}}$. For this combination, the light ring $r_{\rm LR}$ is computed from Eq. (17) but with the potential $\tilde{A}=A_{\rm pp}+A^{\rm tidal}_{\rm dyPN}$. While EOB${}_{\rm dyGSF^{+}}$ allows a more accurate representation of systems with large $\kappa_{2}^{T}$ it also enhances tidal effects for systems with small $\kappa_{2}^{T}$ and leads to problems. The effect of using the combination EOB${}_{\rm dyGSF^{+}}$ is presented in Fig. 10, where the tidal part of the EOB $A$-potential is shown atop the phasing plots. The configurations are SLy-138138 (left panels) and MS1b-165110 (right panels), i.e. the two cases discussed in detail before.

As seen from the right panels, the EOB${}_{\rm dyGSF^{+}}$ model has a smaller dephasing from the NR results for systems with large $\kappa_{2}^{T}$ than the other models. In fact it is only for SLy-138138 that the dephasing of EOB${}_{\rm dyGSF^{+}}$ is comparable to EOB${}_{\rm adGSF^{+}}$ and EOBdyPN, while for all other setups EOB${}_{\rm dyGSF^{+}}$ represents the NR data better than the current models (not shown here).

However, overestimating tidal effects for SLy-138138 is a substantial weakness of this ad-hoc combination of dynamical tides and the GSF inspired potential. While several possibilities to overcome these issues exist, a preferable approach that we will pursue in the future is to first understand and incorporate additional physical effects.

VII.2.2 Including additional physical effects

As demonstrated by the detailed analyses in this paper20170207, improvements to tidal EOB models for the late inspiral will require including additional physical effects that are missing from the current descriptions. For example, existing tidal EOB models assume that the NSs are described as separated objects up to the merger. We showed in Fig. 1, Fig. 8, and Fig. 9 that the low density material of the two NSs comes into contact already well before the merger. Subsequently, yet still prior to the merger, the tidally deformed high-density cores of the NSs are moving inside the viscous material from the former NS outer layers. These physical effects must be taken into account in the tidal EOB models. However, a caveat with this conclusion is the fact that with increasing NR resolution and improved numerical methods the NS’s surface becomes less ‘smeared out’ and the interaction of the low density material decreases. Therefore, we postpone a detailed modeling of these near-merger effects to future work, where we will consider NR simulations with even higher resolution than presented here.

VII.2.3 Calibrating to NR waveforms

Finally, an improvement of the tidal EOB models could also be achieved by incorporating information from NR simulations as done for EOB models for BBH coalescences, see e.g. Taracchini et al. (2014); Damour and Nagar (2014); Bohé et al. (2016). Calibrating to NR data entails introducing additional terms into the models and determining their coefficients by maximizing the overlap between the EOB and NR waveforms. This approach will require a large set of highly accurate NR simulations and possibly eccentricity reduced data. After calibration the models have to be checked against a variety of waveforms covering a large region in the parameter space. This is a challenging task since the unknown EOS and the larger computational cost for general relativistic hydrodynamic simulations make a full coverage of the BNS parameter space substantially more difficult than for BBHs.

VIII Conclusion

In this article we have presented the largest study to date of the performance of tidal EOB models for BNS systems. We considered setups covering a large portion of the BNS parameter space and analyzed variations in the total mass, the mass ratio, and the EOS. Furthermore, we performed a variety of tests to assess the successes and shortcomings of the models by considering several diagnostic quantities, some of which had not been studied in the context of BNSs before.

Our setups encompassed eight different configurations with two different EOSs, three different mass ratios, and two different total masses. We presented comparisons of the energetics of the BNS coalescence characterized by the binding energy vs. specific angular momentum for all models, the periastron advance for one configuration, the time evolution of the GW phase $\phi(t)$ of the (2,2), (2,1), (3,3), and (4,4) modes, and the frequency evolution of the GW phase $\phi(\omega)$ of the (2,2)-mode.

The most effective tool to study the GW signal emitted during the late stage of the BNS coalescence was to compare the accumulated phase within a given frequency interval $\Delta\phi(\omega)$. Considering this quantity enabled studies of the GW phasing without having to align waveforms and emphasized details of the last few orbits before the merger, where tidal effects are strongest.

Overall we found that the EOBdyPN model currently gives the best representation of the dynamics and GWs in the late-inspiral phase, however, the performance of the EOB${}_{\rm adGSF^{+}}$ model is comparable and differences between the models are small. For all setups the EOBadPN model showed a significantly larger deviation from the NR results. We also found that the mass ratio has only a small influence on the performance of the tidal EOB models. Comparing results for higher modes led to similar results as for the dominant (2,2)-mode. We concluded that the main limitation of current tidal EOB models is that tidal effects are underestimated for systems with large tidal coupling constants, i.e. those for which tidal effects are largest. This is a characteristic feature of all tidal EOB models which we have tested.

In addition to comparisons between existing tidal EOB models and new state-of-the-art NR simulations, we also outlined possible avenues towards an improved waveform modeling. As a simple extension of the EOBdyPN and the EOB${}_{\rm adGSF^{+}}$ models we considered the combination of both, i.e. the use of dynamical tides and a GSF inspired potential as a simple way to mimic the larger tidal effects near the merger due to missing physics in the models. Although this method is ‘ad hoc’ and not motivated by the inclusion of additional physical effects, the EOB${}_{\rm dyGSF^{+}}$ model gives a phase evolution $\phi(t)$ closer to the NR results for almost all cases considered.

We emphasize that besides improving the non-spinning tidal EOB model, an important goal for future work on EOB tidal models is to fully include the effects of spins to achieve a better representation of the inspiral dynamics over a larger region of the BNS parameter space. While first NR simulations have characterized the effect of the NS spin on the GW phasing Dietrich et al. (2016b), higher accuracy is needed to place constraints on EOB models, where spin effects for point masses are already incorporated. Further, the extension of the EOB models beyond the moment of merger is of importance and would improve their usability. While first attempts for a post-merger waveform model have been made, e.g. Clark et al. (2016), it may also be possible to derive a post-merger model using the same parameters characterizing the inspiral Bernuzzi et al. (2015c).

Finally, we note that determining the priorities for future work on BNS systems will require inputs from data analysis studies. While it is desirable to have an EOB model that exactly reproduces accurate NR results for the GW phasing and energetics up to the merger, it is important to assess the extent to which the discrepancies for some of the setups discussed in this paper20170207 would contaminate measurements with LIGO and other ground-based GW detectors. 6

Acknowledgements.

It is a pleasure to thank S. Bernuzzi, B. Brügmann, A. Buonanno, B. Lackey, A. Nagar, J. Steinhoff, A. Taracchini, W. Tichy, M. Ujevic for helpful discussion and useful comments on the manuscript. We also want to thank W. Tichy for developing and sharing his SGRID code, S. Bernuzzi and B. Brügmann for allowing us to use the BAM code, A. Taracchini for help with the EOB implementation, and J. Steinhoff for computing the $f-$mode frequencies. Computations where performed on SuperMUC at the LRZ (Munich) under the project number pr48pu, Jureca (Jülich) under the project number HPO21, Stampede (Texas, XSEDE allocation - TG-PHY140019), and the compute cluster Minerva of the Max Planck Institute for Gravitational Physics (Albert Einstein Institute).

Appendix A Eccentricity from NR data

A.1 Eccentricity from the orbital frequency

The orbital frequency can be approximated from the coordinate location of the two NSs by

[TABLE]

where $\vec{r}(t)$ is the coordinate distance between the stars’ centers. (The center of each star is defined as the minimum of the lapse function inside the star.) Because of the non-zero eccentricity and numerical noise $\Omega(t)$ oscillates during the inspiral. We remove high-frequency noise with a low-pass filter. To extract the eccentricity we fit the orbital frequency according to

[TABLE]

with

[TABLE]

where $a_{1}$, $a_{2}$, $b_{1}$, $b_{2}$, $c_{1}$, $c_{2}$, $t_{c}$, and $\tau_{0}$ are fitting parameters. Eq. (35) is based on post-Newtonian calculations, e.g. Hannam et al. (2010), but extended by $(a_{1}t+b_{1})^{2}+(a_{2}t+b_{2})^{4}$. We find that by adding those terms the final fit is more robust than the choice of Hannam et al. (2010), i.e. less dependent on the fitting interval and the initial guess of the parameters. In Figure 11 (upper panel) we show the time evolution of $\Omega(t)$ for MS1b-150100. The raw data are presented in black. The eccentricity induced oscillations and smaller high frequency oscillations are clearly visible. We apply a low-pass filter to eliminate the high frequency oscillations (red dashed line) and fit the filtered results with the model function $\Omega_{\rm fit}(t)$ (blue solid line). The fit depends on the fitting interval and needs to be chosen in a way that the initial junk radiation and gauge dynamics is cut out and that the signal does not extend too close to merger where $\Omega_{\rm fit}(t)$ loses validity. In Fig. 11 we restricted the fit and the analysis to times $t\in[500M,3100M]$. The final eccentricity is then given as

[TABLE]

and shown for MS1b-150100 as a blue dash-dotted line in Fig. 11. Due to the emission of GWs the eccentricity decreases over time which can be seen by an increasing amplitude of $e_{\rm\Omega}$.

A.2 Eccentricity from the GW frequency

In addition to the coordinate location of the two NSs we also use the GW signal to determine the eccentricity. For this purpose we compute the dimensionless GW frequency $M\omega_{22}$. [In the following part of this section, we drop the indices 22 for the (2,2)-mode for simplicity.] We proceed as for $\Omega(t)$ and apply a low-pass filter. Finally, we fit $M\omega_{\rm filtered}$ according to Eq. (35) (blue line), see middle panel of Fig. 11. The eccentricity is then given as

[TABLE]

We present $e_{\omega}$ for MS1b-150100 in the bottom panel of Fig. 11 as a solid red line.

A.3 Eccentricity from the GW phase

The last eccentricity indicator is based on the phase $\phi$ of the (2,2)-mode. To extract the eccentricity from the phase evolution, we fit $\phi$ according to

[TABLE]

after applying a low-pass filter. The eccentricity is then given as

[TABLE]

The eccentricity $e_{\phi}$ for MS1b-150100 is shown in the bottom panel of Fig. 11 as a green dashed line.

Appendix B Extracting $\phi(\omega)$ from NR data

In previous work we used a non-linear fit of the GW phase to compute the phase evolution as a function of time Dietrich et al. (2016b). This approach is similar to that used for the PA in Appendix A, but more delicate since the fit has to extend up to the merger. To overcome this problem, we use a different method in this work.

The biggest problem is that due to eccentricity and numerical noise $\omega(t)$ might be non-monotonic. We minimize the residual eccentricity with a low-pass filter applied to $\phi(t)$ and $\omega(t)$. Since the eccentricity decreases during the evolution while the phase and frequency evolution accelerates, it is not possible to apply a low pass filter with a fixed cutoff frequency. Therefore we evaluate two different filtered functions and perform a smooth transition between them centered around $M\omega_{22}=0.055$ (except for MS1b-150100, where we use $M\omega_{22}=0.053$ since the signal is shorter). As an exemplary case we show the setup MS1b-153122 in Fig. 12, where the raw data are shown as red crosses and the fit as a black line. Notice that the NR data cluster at low frequencies since they are equally spaced in time. Although the filtering reduces the eccentricity, some remaining eccentricity in particular at the transition between the two filtered function can be present. This is the reason for oscillations of $\Delta\phi(\omega)$ in the lower panels of Fig. 7.

Appendix C Influence of the underlying point mass EOB framework

Although this paper20170207 focuses on the performance of tidal EOB models using different ways of including tidal effects, we also want to briefly quantify the influence of the underlying point mass EOB model during the merger of BNS systems. For this reason we present again results for SLy-138138 and MS1b-165110 as in Sec. VII. The particular effects we study are:

1. (i)

the radial radiation reaction force $\mathcal{F}_{r}$, where we compare results with $\mathcal{F}_{r}=0$ as employed in Damour et al. (2013); Bernuzzi et al. (2015a) with those that use $\mathcal{F}_{r}=(p_{r}/p_{\phi}){\cal F}_{\phi}$ as in Taracchini et al. (2012); Hinderer et al. (2016), (shown red in Fig. 13). 2. (ii)

NQC corrections using the prescription of Taracchini et al. (2012) (shown orange in Fig. 13). 3. (iii)

calibration and resummation of the potentials, i.e. using simply the PN Taylor expanded EOB potentials (shown blue in Fig. 13). 4. (iv)

different resummation techniques for the EOB potentials, where we compare a Pade resummation with the calibration summarized in Bernuzzi et al. (2015a) with a logarithmic resummation of the EOB potentials as given explicitly in Steinhoff et al. (2016) with the calibration from Taracchini et al. (2014) (shown black in Fig. 13). 5. (v)

different arguments in the $h_{\ell m}^{F}$ modes and thus also the radiation reaction force, where we compare results from using, in some occurences, a different frequency-related variable that was advocated in Damour et al. (2009) and is obtained by replacing $\Omega\to\Omega^{3}/(\partial H_{\rm EOB}/\partial p_{\phi})^{2}$ evaluated for circular orbits (shown green in Fig. 13).

Figure 13 summarizes our results. The baseline is the EOBdyPN model as employed throughout the article. To quantify the effects (i)-(v), we compute the difference in the time evolution of the GW phase (top panels), the binding energy (middle panels), and the phase evolution as a function of the frequency (bottom panels).

Clearly visible is that for all quantities the underlying point mass EOB model effects the SLy-138138 configuration more than the MS1b-165110 setup. This is caused by the larger compactness and smaller radii of the NSs for SLy-138138, consequently the stars come closer and reach higher frequencies during the inspiral. Considering the differences $\Delta\phi(t)=\phi_{\rm dyPN_{\mathcal{F}_{r}=0}}-\phi_{\rm dyPN}>0$ and $\Delta\phi(t)=\phi_{\rm dyPN_{\rm\neg NQC}}-\phi_{\rm dyPN}>0$ we find that the omission of a radial reaction force or NQC corrections accelerates the inspiral, which mimics larger tidal effects. The opposite is true when using the PN-potentials. The particular resummation technique and argument of the $h_{\ell m}$ modes have only a minor effect on $\phi(t)$.

Considering the binding energy curves, we find that except when using a Pade resummation, all other changes of the underlying point mass EOB model lead to $\Delta E(\ell)>0$, i.e., the setups are less bound than the EOBdyPN model used in this work.

Finally, considering the influence of the underlying point mass model on $\phi(\omega)$ we find that all phase differences $\phi(\omega)$ are below $0.5\rm rad$. This is significantly smaller than the dephasing accumulated by tidal effects and also significantly smaller than the difference between the EOBdyPN model and the NR simulations (cf. e.g. bottom panel of Fig. 10).