Sub-radian-accuracy gravitational waves from coalescing binary neutron stars II: Systematic study on the equation of state, binary mass, and mass ratio

TL;DR

This study performs high-precision numerical relativity simulations of binary neutron star mergers, validating gravitational waveform templates and examining universal relations, with implications for gravitational wave data analysis.

Contribution

It provides systematic, high-resolution simulations of binary neutron star mergers, validating waveform templates and critically assessing universal relations for tidal deformability.

Findings

Phase error less than 0.5 rad over 200 rad phase evolution

Waveform accuracy within 0.1 rad in phase and 20% in amplitude up to 1000 Hz

Universal relations for post-merger signals show significant systematics

Abstract

We report results of numerical relativity simulations for {\it new} 26 non-spinning binary neutron star systems with 6 grid resolutions using an adaptive mesh refinement numerical re\ lativity code {\tt SACRA-MPI}. The finest grid spacing is $\approx 64$--$85$ m, depending on the systems. First, we derive long-term high-precision inspiral gravitational waveforms and show that the accumulated gravitational-wave phase error due to the finite grid resolution is less than $0.5$ rad during more than $200$ rad phase evolution irrespective of the systems. We also find that the gravitational-wave phase error for a binary system with a tabulated equation of state (EOS) is comparable to that for a piecewise polytropic EOS. Then we validate the SACRA inspiral gravitational waveform template, which will be used to extract tidal deformability from gravitational wave observation, and find that accuracy of \ our waveform modeling is $\lesssim 0.1$ rad in the gravitational-wave phase and $\lesssim 20 \%$ in the gravitational-wave amplitude up to the gravitational-wave frequency $1000$ Hz.\ Finally, we calibrate the proposed universal relations between a post-merger gravitational wave signal and tidal deformability/neutron star radius in the literature and show that th\ ey suffer from systematics and many relations proposed as universal are not very universal. Improved fitting formulae are also proposed.

Figures

Tables

Table 1. Table 1: List of the systems for which we performed new simulations. The names of the systems are given in the 1st column. The 2nd and 3rd columns show gravitational mass of less massive NS, m 1 subscript 𝑚 1 m_{1} , and massive NS, m 2 subscript 𝑚 2 m_{2} , respectively. The 4th column shows EOS. Dimensionless initial orbital angular velocity, m 0 ​ Ω 0 subscript 𝑚 0 subscript Ω 0 m_{0}\Omega_{0} , with the total gravitational mass of the binary systems, m 0 = m 1 + m 2 subscript 𝑚 0 subscript 𝑚 1 subscript 𝑚 2 m_{0}=m_{1}+m_{2} , is given in the 5th column. The 6th, 7th, and 8th columns show chirp mass, ℳ c = ( m 1 ​ m 2 ) 3 / 5 ​ ( m 1 + m 2 ) − 1 / 5 subscript ℳ 𝑐 superscript subscript 𝑚 1 subscript 𝑚 2 3 5 superscript subscript 𝑚 1 subscript 𝑚 2 1 5 {\cal M}_{c}=(m_{1}m_{2})^{3/5}(m_{1}+m_{2})^{-1/5} , symmetric mass ratio, η = m 1 ​ m 2 ​ ( m 1 + m 2 ) − 2 𝜂 subscript 𝑚 1 subscript 𝑚 2 superscript subscript 𝑚 1 subscript 𝑚 2 2 \eta=m_{1}m_{2}(m_{1}+m_{2})^{-2} , and binary tidal deformability, Λ ~ ~ Λ \tilde{\Lambda} , respectively. Location of outer boundary in a computational domain, L 𝐿 L , and grid spacing of a finest AMR level, Δ ​ x 9 Δ subscript 𝑥 9 \Delta x_{9} , are given in the 9th and 10th columns, respectively. The grid spacing with N = 182 , 150 , 130 , 110 , 102 , 𝑁 182 150 130 110 102 N=182,150,130,110,102, and 90 90 90 is shown in the parenthesis in the 10th column. The final column shows the extraction radii of gravitational waves.

System	$m_1\left[M_{\odot }\right]$	$m_2\left[M_{\odot }\right]$	EOS	$m_0{\Omega }_0$	${ℳ}_c$	$\eta$	$\tilde{\Lambda }$	$L\left[\mathrm{km}\right]$	$\Delta x_9\left[m\right]$	$r_0/m_0$
15H125-146	1.25	1.46	15H	0.0155	1.1752	0.2485	1200	7823	(84,102,117,138,149,169)	(244,199,155)
125H125-146	1.25	1.46	125H	0.0155	1.1752	0.2485	858	7323	(78,95,110,129,140,158)	(244,199,155)
H125-146	1.25	1.46	H	0.0155	1.1752	0.2485	605	6824	(73,89,102,121,130,147)	(244,199,155)
HB125-146	1.25	1.46	HB	0.0155	1.1752	0.2485	423	6491	(69,84,97,115,124,140)	(244,199,155)
B125-146	1.25	1.46	B	0.0155	1.1752	0.2485	290	5992	(64,78,90,106,114,129)	(244,199,155)
15H118-155	1.18	1.55	15H	0.0155	1.1752	0.2455	1194	7889	(84,102,118,139,150,170)	(242,198,154)
125H118-155	1.18	1.55	125H	0.0155	1.1752	0.2455	855	7390	(79,96,111,131,141,159)	(242,198,154)
H118-155	1.18	1.55	H	0.0155	1.1752	0.2455	606	6990	(75,91,105,124,133,151)	(242,198,154)
HB118-155	1.18	1.55	HB	0.0155	1.1752	0.2455	423	6491	(69,84,97,115,124,140)	(242,198,154)
B118-155	1.18	1.55	B	0.0155	1.1752	0.2455	292	5992	(64,78,90,106,114,129)	(242,198,154)
15H117-156	1.17	1.56	15H	0.0155	1.1752	0.2450	1170	7889	(84,102,118,139,150,170)	(242,198,154)
125H117-156	1.17	1.56	125H	0.0155	1.1752	0.2450	837	7323	(78,95,110,129,140,158)	(242,198,154)
H117-156	1.17	1.56	H	0.0155	1.1752	0.2450	592	6990	(75,91,105,124,133,151)	(242,198,154)
HB117-156	1.17	1.56	HB	0.0155	1.1752	0.2450	414	6491	(69,84,97,115,124,141)	(242,198,154)
B117-156	1.17	1.56	B	0.0155	1.1752	0.2450	285	6058	(65,79,91,107,115,131)	(242,198,154)
15H112-140	1.12	1.40	15H	0.0150	1.0882	0.2470	1842	7989	(85,104,120,141,152,172)	(262,214,167)
125H112-140	1.12	1.40	125H	0.0150	1.0882	0.2470	1332	7490	(80,97,112,132,143,162)	(262,214,167)
H112-140	1.12	1.40	H	0.0150	1.0882	0.2470	955	6990	(75,91,105,124,133,151)	(262,214,167)
HB112-140	1.12	1.40	HB	0.0150	1.0882	0.2470	677	6491	(69,84,97,115,124,140)	(262,214,167)
B112-140	1.12	1.40	B	0.0150	1.0882	0.2470	475	6092	(65,79,91,108,116,131)	(262,214,167)
15H107-146	1.07	1.46	15H	0.0150	1.0882	0.2440	1845	7989	(85,104,120,141,152,172)	(261,213,166)
125H107-146	1.07	1.46	125H	0.0150	1.0882	0.2440	1335	7490	(80,97,112,132,143,162)	(261,213,166)
H107-146	1.07	1.46	H	0.0150	1.0882	0.2440	957	6990	(75,91,105,124,133,151)	(261,213,166)
HB107-146	1.07	1.46	HB	0.0150	1.0882	0.2440	684	6591	(71,86,99,117,126,142)	(261,213,166)
B107-146	1.07	1.46	B	0.0150	1.0882	0.2440	481	6091	(65,79,91,108,116,131)	(261,213,166)
SFHo135-135	1.35	1.35	SFHo	0.0155	1.1752	0.2500	460	6491	(69,84,97,115,124,140)	(244,200,156)

Table 2. Table 2: List of ρ 1 subscript 𝜌 1 \rho_{1} in two-piecewise polytropic EOSs.

EOS	${\rho }_1\left[g{\mathrm{cm}}^{-3}\right]$
15H	$9.3108\times {10}^{13}$
125H	$1.0711\times {10}^{14}$
H	$1.2323\times {10}^{14}$
HB	$1.4177\times {10}^{14}$
B	$1.6309\times {10}^{14}$

Table 3. Table 3: The radius, R M subscript 𝑅 𝑀 R_{M} , and the dimensionless tidal deformability, Λ M subscript Λ 𝑀 \Lambda_{M} , for spherical NSs with gravitational mass M = 1.07 𝑀 1.07 M=1.07 , 1.12 1.12 1.12 , 1.17 1.17 1.17 , 1.18 1.18 1.18 , 1.25 1.25 1.25 , 1.35 1.35 1.35 , 1.40 1.40 1.40 , 1.46 1.46 1.46 , 1.55 1.55 1.55 and 1.56 ​ M ⊙ 1.56 subscript 𝑀 direct-product 1.56\,M_{\odot} for the given EOS. R M subscript 𝑅 𝑀 R_{M} is listed in units of km km {\rm km} . For SFHo (tabulated) EOS, the quantities for the spherical star with M = 1.35 ​ M ⊙ 𝑀 1.35 subscript 𝑀 direct-product M=1.35\,M_{\odot} are listed. The last column in the upper table shows the maximum mass of the spherical NS in units of M ⊙ subscript 𝑀 direct-product M_{\odot} .

EOS	$R_{1.07}$	$R_{1.12}$	$R_{1.17}$	$R_{1.18}$	$R_{1.25}$	$R_{1.35}$	$R_{1.40}$	$R_{1.46}$	$R_{1.55}$	$R_{1.56}$	$M_{\text{max}}$
15H	13.54	13.58	13.61	13.62	13.65	13.69	13.71	13.72	13.74	13.74	2.53
125H	12.86	12.89	12.91	12.92	12.94	12.97	12.98	12.99	12.98	12.98	2.38
H	12.22	12.23	12.24	12.24	12.26	12.27	12.28	12.18	12.26	12.25	2.25
HB	11.60	11.59	11.60	11.60	11.61	11.61	11.60	11.59	11.55	11.55	2.12
B	10.97	10.97	10.98	10.98	10.98	10.96	10.95	10.92	10.87	10.86	2.00
SFHo	–	–	–	–	–	11.91	–	–	–	–	2.06
EOS	${\Lambda }_{1.07}$	${\Lambda }_{1.12}$	${\Lambda }_{1.17}$	${\Lambda }_{1.18}$	${\Lambda }_{1.25}$	${\Lambda }_{1.35}$	${\Lambda }_{1.40}$	${\Lambda }_{1.46}$	${\Lambda }_{1.55}$	${\Lambda }_{1.56}$
15H	4361	3411	2692	2575	1871	1211	975	760	530	509
125H	3196	2490	1963	1875	1351	863	693	535	366	350
H	2329	1812	1415	1354	966	607	484	369	249	238
HB	1695	1304	1013	966	684	422	333	252	165	157
B	1216	933	719	681	477	289	225	168	107	101
SFHo	–	–	–	–	–	460	–	–	–	–

Table 4. Table 4: Residual phase error (rad) and order of the convergence of the gravitational-wave phase at the peak time calculated by Eq. ( 16 ) for N max = 182 , 150 subscript 𝑁 max 182 150 N_{\text{max}}=182,150 , and 130 130 130 .

System	$(\Delta {\Phi }_{\text{peak}}^{2,2}\left(182\right),p)$	$(\Delta {\Phi }_{\text{peak}}^{2,2}\left(150\right),p)$	$(\Delta {\Phi }_{\text{peak}}^{2,2}\left(130\right),p)$
15H125-146	(0.11, 4.1)	(0.58, 2.7)	(5.44, 0.7)
125H125-146	(0.31, 2.6)	(0.15, 4.5)	(0.45, 3.6)
H125-146	(0.17, 3.4)	(0.78, 2.2)	(0.73, 2.8)
HB125-146	(0.13, 3.7)	(1.10, 1.7)	(1.00, 2.2)
B125-146	(0.12, 3.8)	(0.28, 3.7)	(0.45, 3.8)
15H118-155	(0.22, 3.1)	(0.75, 2.2)	(0.47, 3.5)
125H118-155	(0.26, 2.9)	(0.83, 2.1)	(1.44, 1.7)
H118-155	(0.23, 3.1)	(0.48, 3.0)	(0.56, 3.4)
HB118-155	(0.44, 2.3)	(1.21, 1.6)	(0.79, 2.5)
B118-155	(0.29, 2.7)	(0.69, 2.2)	(0.47, 3.3)
15H117-156	(0.26, 2.9)	(0.36, 3.2)	(0.39, 4.0)
125H117-156	(0.28, 2.8)	(0.38, 2.8)	(0.92, 2.4)
H117-156	(0.24, 3.0)	(0.31, 3.5)	(0.74, 2.9)
HB117-156	(0.22, 3.0)	(0.84, 2.0)	(1.42, 1.7)
B117-156	(0.42, 2.3)	(0.43, 2.8)	(0.23, 4.8)
15H112-140	(0.19, 3.4)	(0.70, 2.5)	(0.66, 3.2)
125H112-140	(0.21, 3.4)	(0.53, 3.0)	(0.66, 3.3)
H112-140	(0.17, 3.5)	(0.92, 2.1)	(1.00, 2.4)
HB112-140	(0.42, 2.5)	(0.48, 3.0)	(0.21, 5.5)
B112-140	(0.19, 3.6)	(0.34, 3.7)	(39.59, 0.13)
15H107-146	(0.38, 2.6)	(0.86, 2.2)	(0.43, 3.9)
125H107-146	(0.54, 2.2)	(2.93, 1.0)	(0.61, 3.2)
H107-146	(0.41, 2.4)	(0.60, 2.5)	(1.03, 2.3)
HB107-146	(0.35, 2.8)	(0.44, 3.3)	(0.43, 4.2)
B107-146	(0.33, 2.8)	(0.73, 2.4)	(1.05, 2.4)
SFHo135-135	(0.43, 2.3)	(0.76, 2.2)	(0.33, 4.2)

Table 5. Table 5: Mismatch between the inspiral waveform model and hybrid waveforms.

System	$\overline{F}(\times {10}^{-5})$
15H125-146	0.83
125H125-146	0.36
H125-146	0.29
HB125-146	0.28
B125-146	0.22
15H118-155	0.82
125H118-155	0.26
H118-155	0.30
HB118-155	0.32
B118-155	0.31
15H117-156	0.97
125H117-156	0.31
H117-156	0.25
HB117-156	0.30
B117-156	0.17
15H112-140	0.88
125H112-140	0.24
H112-140	0.37
HB112-140	0.71
B112-140	0.91
15H107-146	1.82
125H107-146	0.45
H107-146	0.30
HB107-146	0.79
B107-146	1.12
SFHo135-135	0.45

Table 6. Table 6: Binary tidal deformability Λ ~ ~ Λ \tilde{\Lambda} , f peak subscript 𝑓 peak f_{\text{peak}} , h peak subscript ℎ peak h_{\text{peak}} , f 2 subscript 𝑓 2 f_{2} , E GW , i 2 , 2 subscript superscript 𝐸 2 2 GW i E^{2,2}_{\rm GW,i} , E GW , p 2 , 2 subscript superscript 𝐸 2 2 GW p E^{2,2}_{\rm GW,p} , J GW , p 2 , 2 subscript superscript 𝐽 2 2 GW p J^{2,2}_{\rm GW,p} , J rem subscript 𝐽 rem J_{\text{rem}} , and m 0 − M ADM,0 subscript 𝑚 0 subscript 𝑀 ADM,0 m_{0}-M_{\text{ADM,0}} . M ADM,0 subscript 𝑀 ADM,0 M_{\text{ADM,0}} is the Arnowitt-Deser-Misner mass of the initial condition of the simulations. We adopt 2 % percent \% relative error for f peak subscript 𝑓 peak f_{\rm peak} and h peak subscript ℎ peak h_{\text{peak}} and 5 % percent 5 5\% relative error for f 2 subscript 𝑓 2 f_{2} , respectively, as a typical value. For f 2 subscript 𝑓 2 f_{2} , we exclude binary systems which collapse to a black hole within a few ms after the merger. For E GW , i 2 , 2 subscript superscript 𝐸 2 2 GW i E^{2,2}_{\rm GW,i} and J rem subscript 𝐽 rem J_{\text{rem}} , we adopt 2 % percent 2 2\% and 1 % percent 1 1\% relative error, respectively. E GW subscript 𝐸 GW E_{\text{GW}} and m 0 − M ADM,0 subscript 𝑚 0 subscript 𝑀 ADM,0 m_{0}-M_{\text{ADM,0}} are given in units of M ⊙ subscript 𝑀 direct-product M_{\odot} . J GW subscript 𝐽 GW J_{\text{GW}} and J rem subscript 𝐽 rem J_{\text{rem}} are in units of M ⊙ 2 superscript subscript 𝑀 direct-product 2 M_{\odot}^{2} .

System	${\tilde{\Lambda }}^{1/5}$	$f_{\mathrm{peak}}$ [Hz]	$D{h}_{\text{peak}}/m_0$	$f_2\left[\mathrm{Hz}\right]$	$E_{\text{GW},i}^{2,2}$	$E_{\text{GW},p}^{2,2}$	$J_{\text{GW},p}^{2,2}$	$J_{\text{rem}}$	$m_0-M_{\text{ADM,0}}$
15H135-135	4.14	1503$\pm$30	0.226$\pm$0.005	2321$\pm$116	(7.90$\pm 0.16)\times 10{}^{-3}$	1.35$\times {10}^{-2}$	0.40	6.64$\pm 0.07$	$1.65\times {10}^{-2}$
125H135-135	3.87	1652$\pm$33	0.236$\pm$0.005	2517$\pm$126	(9.04$\pm 0.18)\times 10{}^{-3}$	1.76$\times {10}^{-2}$	0.48	6.54$\pm 0.07$	$1.64\times {10}^{-2}$
H135-135	3.60	1820$\pm$36	0.249$\pm$0.005	2790$\pm$139	(1.03$\pm 0.02)\times 10{}^{-2}$	2.32$\times {10}^{-2}$	0.56	6.46$\pm 0.06$	$1.63\times {10}^{-2}$
HB135-135	3.35	1986$\pm$40	0.261$\pm$0.005	3243$\pm$162	(1.17$\pm 0.02)\times 10{}^{-2}$	2.89$\times {10}^{-2}$	0.59	6.39$\pm 0.06$	$1.64\times {10}^{-2}$
B135-135	3.11	2133$\pm$43	0.274$\pm$0.005	–	(1.30$\pm 0.03)\times 10{}^{-2}$	7.39$\times {10}^{-3}$	0.13	6.33$\pm 0.06$	$1.65\times {10}^{-2}$
15H121-151	4.13	1356$\pm$27	0.212$\pm$0.004	2261$\pm$163	(7.47$\pm 0.15)\times 10{}^{-3}$	5.47$\times {10}^{-3}$	0.17	6.66$\pm 0.07$	$1.66\times {10}^{-2}$
125H121-151	3.86	1490$\pm$30	0.224$\pm$0.004	2379$\pm$119	(8.53$\pm 0.17)\times 10{}^{-3}$	8.24$\times {10}^{-3}$	0.23	6.57$\pm 0.07$	$1.66\times {10}^{-2}$
H121-151	3.60	1637$\pm$33	0.236$\pm$0.005	2749$\pm$137	(9.70$\pm 0.19)\times 10{}^{-3}$	1.05$\times {10}^{-2}$	0.26	6.49$\pm 0.06$	$1.66\times {10}^{-2}$
HB121-151	3.35	1809$\pm$36	0.249$\pm$0.005	3268$\pm$161	(1.10$\pm 0.02)\times 10{}^{-2}$	2.26$\times {10}^{-2}$	0.48	6.41$\pm 0.06$	$1.66\times {10}^{-2}$
B121-151	3.11	1994$\pm$40	0.263$\pm$0.005	–	(1.23$\pm 0.02)\times 10{}^{-2}$	6.85$\times {10}^{-3}$	0.13	6.35$\pm 0.06$	$1.66\times {10}^{-2}$
15H125-125	4.51	1450$\pm$29	0.211$\pm$0.004	2159$\pm$108	(6.26$\pm 0.13)\times 10{}^{-3}$	7.98$\times {10}^{-3}$	0.25	5.95$\pm 0.06$	$1.53\times {10}^{-2}$
125H125-125	4.23	1568$\pm$31	0.222$\pm$0.004	2350$\pm$118	(7.19$\pm 0.14)\times 10{}^{-3}$	9.29$\times {10}^{-3}$	0.27	5.87$\pm 0.06$	$1.53\times {10}^{-2}$
H125-125	3.95	1710$\pm$34	0.234$\pm$0.005	2749$\pm$137	(8.15$\pm 0.16)\times 10{}^{-3}$	1.67$\times {10}^{-2}$	0.42	5.80$\pm 0.06$	$1.52\times {10}^{-2}$
HB125-125	3.69	1900$\pm$38	0.245$\pm$0.005	2873$\pm$144	(9.35$\pm 0.19)\times 10{}^{-3}$	1.66$\times {10}^{-2}$	0.39	5.74$\pm 0.06$	$1.53\times {10}^{-2}$
B125-125	3.43	2099$\pm$42	0.257$\pm$0.005	3353$\pm$168	(1.06$\pm 0.02))\times 10{}^{-2}$	2.19$\times {10}^{-2}$	0.44	5.69$\pm 0.06$	$1.53\times {10}^{-2}$
15H116-158	4.12	1273$\pm$26	0.205$\pm$0.004	2148$\pm$107	(7.19$\pm 0.14)\times 10{}^{-3}$	4.63$\times {10}^{-3}$	0.15	6.84$\pm 0.07$	$1.65\times {10}^{-2}$
125H116-158	3.85	1406$\pm$28	0.214$\pm$0.004	2276$\pm$124	(8.20$\pm 0.16)\times 10{}^{-3}$	1.01$\times {10}^{-2}$	0.28	6.76$\pm 0.07$	$1.65\times {10}^{-2}$
H116-158	3.60	1540$\pm$31	0.227$\pm$0.005	2767$\pm$138	(9.30$\pm 0.19)\times 10{}^{-3}$	1.23$\times {10}^{-2}$	0.31	6.69$\pm 0.07$	$1.66\times {10}^{-2}$
HB116-158	3.35	1709$\pm$34	0.240$\pm$0.005	3242$\pm$162	(1.05$\pm 0.02)\times 10{}^{-2}$	1.40$\times {10}^{-2}$	0.30	6.63$\pm 0.06$	$1.65\times {10}^{-2}$
B116-158	3.11	1885$\pm$37	0.254$\pm$0.005	–	(1.18$\pm 0.02)\times 10{}^{-2}$	4.64$\times {10}^{-3}$	0.10	6.58$\pm 0.07$	$1.65\times {10}^{-2}$
15H125-146	4.13	1401$\pm$28	0.214$\pm$0.004	2336$\pm$117	(7.62$\pm 0.02)\times 10{}^{-3}$	1.01$\times {10}^{-2}$	0.30	6.81$\pm 0.07$	$1.66\times {10}^{-2}$
125H125-146	3.86	1560$\pm$31	0.226$\pm$0.005	2576$\pm$129	(8.77$\pm 0.18)\times 10{}^{-3}$	1.26$\times {10}^{-2}$	0.34	6.73$\pm 0.07$	$1.66\times {10}^{-2}$
H125-146	3.60	1691$\pm$34	0.238$\pm$0.003	2827$\pm$141	(9.91$\pm 0.20)\times 10{}^{-3}$	1.89$\times {10}^{-2}$	0.45	6.66$\pm 0.07$	$1.66\times {10}^{-2}$
HB125-146	3.35	1856$\pm$37	0.252$\pm$0.005	3251$\pm$163	(1.12$\pm 0.20)\times 10{}^{-2}$	2.50$\times {10}^{-2}$	0.52	6.60$\pm 0.07$	$1.66\times {10}^{-2}$
B125-146	3.11	2039$\pm$41	0.265$\pm$0.005	–	(1.26$\pm 0.25)\times 10{}^{-2}$	7.99$\times {10}^{-3}$	0.14	6.56$\pm 0.06$	$1.66\times {10}^{-2}$
15H118-155	4.12	1308$\pm$26	0.206$\pm$0.004	2161$\pm$108	(7.31$\pm 0.15)\times 10{}^{-3}$	5.72$\times {10}^{-3}$	0.18	6.83$\pm 0.07$	$1.66\times {10}^{-2}$
125H118-155	3.86	1441$\pm$29	0.218$\pm$0.004	2358$\pm$118	(8.35$\pm 0.17)\times 10{}^{-3}$	7.12$\times {10}^{-3}$	0.21	6.75$\pm 0.07$	$1.67\times {10}^{-2}$
H118-155	3.60	1590$\pm$32	0.230$\pm$0.005	2782$\pm$139	(9.49$\pm 0.19)\times 10{}^{-3}$	1.59$\times {10}^{-2}$	0.39	6.68$\pm 0.07$	$1.66\times {10}^{-2}$
HB118-155	3.35	1759$\pm$35	0.243$\pm$0.005	3259$\pm$163	(1.08$\pm 0.02)\times 10{}^{-2}$	2.03$\times {10}^{-2}$	0.43	6.62$\pm 0.07$	$1.66\times {10}^{-2}$
B118-155	3.11	1942$\pm$39	0.257$\pm$0.005	–	(1.20$\pm 0.02)\times 10{}^{-2}$	5.54$\times {10}^{-3}$	0.11	6.66$\pm 0.07$	$1.66\times {10}^{-2}$
15H117-156	4.11	1293$\pm$26	0.204$\pm$0.004	2161$\pm$108	(7.26$\pm 0.15)\times 10{}^{-3}$	5.09$\times {10}^{-3}$	0.17	6.83$\pm 0.07$	$1.66\times {10}^{-2}$
125H117-156	3.84	1425$\pm$29	0.216$\pm$0.004	2416$\pm$121	(8.30$\pm 0.17)\times 10{}^{-3}$	8.09$\times {10}^{-3}$	0.23	6.76$\pm 0.07$	$1.66\times {10}^{-2}$
H117-156	3.58	1574$\pm$32	0.229$\pm$0.005	2775$\pm$139	(9.43$\pm 0.19)\times 10{}^{-3}$	1.39$\times {10}^{-2}$	0.34	6.69$\pm 0.07$	$1.66\times {10}^{-2}$
HB117-156	3.34	1724$\pm$35	0.242$\pm$0.005	3201$\pm$160	(1.06$\pm 0.02)\times 10{}^{-2}$	1.61$\times {10}^{-2}$	0.35	6.62$\pm 0.07$	$1.66\times {10}^{-2}$
B117-156	3.10	1933$\pm$38	0.256$\pm$0.005	–	(1.20$\pm 0.02)\times 10{}^{-2}$	5.26$\times {10}^{-3}$	0.11	6.58$\pm 0.06$	$1.64\times {10}^{-2}$
15H112-140	4.50	1281$\pm$26	0.197$\pm$0.004	2188$\pm$109	(5.91$\pm 0.12)\times 10{}^{-3}$	5.37$\times {10}^{-3}$	0.17	5.97$\pm 0.06$	$1.49\times {10}^{-2}$
125H112-140	4.21	1412$\pm$28	0.208$\pm$0.004	2269$\pm$113	(6.80$\pm 0.14)\times 10{}^{-3}$	4.80$\times {10}^{-3}$	0.15	5.89$\pm 0.06$	$1.49\times {10}^{-2}$
H112-140	3.94	1558$\pm$31	0.220$\pm$0.004	2470$\pm$123	(7.78$\pm 0.16)\times 10{}^{-3}$	6.18$\times {10}^{-3}$	0.17	5.82$\pm 0.06$	$1.50\times {10}^{-2}$
HB112-140	3.68	1717$\pm$34	0.231$\pm$0.005	2791$\pm$140	(8.84$\pm 0.18)\times 10{}^{-3}$	9.52$\times {10}^{-3}$	0.23	5.76$\pm 0.06$	$1.50\times {10}^{-2}$
B112-140	3.43	1890$\pm$38	0.244$\pm$0.005	3271$\pm$164	(9.98$\pm 0.20)\times 10{}^{-3}$	1.59$\times {10}^{-2}$	0.33	5.71$\pm 0.06$	$1.52\times {10}^{-2}$
15H107-146	4.50	1203$\pm$24	0.189$\pm$0.004	2054$\pm$103	(5.70$\pm 0.11)\times 10{}^{-3}$	3.63$\times {10}^{-3}$	0.13	5.99$\pm 0.06$	$1.51\times {10}^{-2}$
125H107-146	4.22	1328$\pm$27	0.200$\pm$0.004	2291$\pm$115	(6.57$\pm 0.13)\times 10{}^{-3}$	4.56$\times {10}^{-3}$	0.14	5.91$\pm 0.06$	$1.50\times {10}^{-2}$
H107-146	3.94	1475$\pm$30	0.212$\pm$0.004	2546$\pm$127	(7.49$\pm 0.15)\times 10{}^{-3}$	7.82$\times {10}^{-3}$	0.21	5.84$\pm 0.06$	$1.49\times {10}^{-2}$
HB107-146	3.69	1620$\pm$32	0.224$\pm$0.004	2870$\pm$143	(8.51$\pm 0.17)\times 10{}^{-3}$	1.02$\times {10}^{-2}$	0.25	5.78$\pm 0.06$	$1.50\times {10}^{-2}$
B107-146	3.44	1786$\pm$36	0.237$\pm$0.005	3298$\pm$165	(9.60$\pm 0.19)\times 10{}^{-3}$	1.29$\times {10}^{-2}$	0.27	5.73$\pm 0.06$	$1.51\times {10}^{-2}$
SFHo135-135	3.41	1987$\pm$40	0.261$\pm$0.005	3250$\pm$163	(1.17$\pm 0.02)\times 10{}^{-2}$	2.91$\times {10}^{-3}$	0.61	6.60$\pm 0.07$	$1.68\times {10}^{-2}$

Table 7. Table 7: Summary of the assessment of the universal relations for the non-spinning and non-magnetized binary systems. Neutrino radiation is not taken into account. We show the maximum relative errors produced by the original relation (upper row) and by the improved relation derived in this paper (lower row). For f 1 subscript 𝑓 1 f_{1} , the error is unable to be estimated because of the absence of f 1 subscript 𝑓 1 f_{1} peak in the asymmetric binary systems. Therefore, we conclude there is no universal relation between f 1 subscript 𝑓 1 f_{1} and Λ ~ ~ Λ \tilde{\Lambda} . For f 2 subscript 𝑓 2 f_{2} – R 1.6 subscript 𝑅 1.6 R_{1.6} relation, we do not propose an improved relation and sym. (asym.) in the parenthesis means the symmetric (asymmetric) binary. For E GW , p 2 , 2 subscript superscript 𝐸 2 2 GW p E^{2,2}_{\rm GW,p} and J GW , p 2 , 2 subscript superscript 𝐽 2 2 GW p J^{2,2}_{\rm GW,p} , we do not propose an improved relation because uncertainties of the life time of the merger remnant NSs are large.

$m_0{f}_{\text{peak}}$–${\tilde{\Lambda }}^{1/5}$	$D{h}_{\text{peak}}/m_0$–${\tilde{\Lambda }}^{1/5}$	$f_1$–${\tilde{\Lambda }}^{1/5}$	$m_0{f}_2$–${\tilde{\Lambda }}^{1/5}$	$f_2$–$R_{1.6}$	$E_{\mathrm{GW},i}^{\text{tot}}/\left(m_0\eta \right)$–${\tilde{\Lambda }}^{1/5}$	$J_{\text{rem}}/\left(m_0^2\eta \right)$–${\tilde{\Lambda }}^{1/5}$
$\approx 17\%$	N/A	–	$\approx 14\%$	$\approx 6\%\text{(sym.)}$ and $\approx 10\%\text{(asym.)}$	N/A	N/A
$\approx 3\%$	$\approx 4\%$	–	$\approx 9\%$	–	$\approx 4\%$	$\approx 3\%$

Equations

$$P_{\text{cold}}(\rho )={\kappa }_i{\rho }^{{\Gamma }_i},$$

$${ϵ}_{\text{cold}}(\rho )=\frac{{\kappa }_i}{{\Gamma }_i-1}{\rho }^{{\Gamma }_i-1}+\Delta {ϵ}_i({\rho }_i\le \rho <{\rho }_{i+1}),$$

$$P=P_{\text{cold}}(\rho )+({\Gamma }_{\text{th}}-1)\rho (ϵ-{ϵ}_{\text{cold}}(\rho )),$$

$$\tilde{\Lambda }$$

$$\displaystyle-\sqrt{1-4\eta}(1+9\eta-11\eta^{2})(\Lambda_{1}-\Lambda_{2})\Big{]},$$

$${\Psi }_4(t_{\text{ret}},r_0,\theta ,\varphi )=\sum _{l,m}{\Psi }_4^{l,m}(t_{\text{ret}},r_0{)}_{-2}{Y}_{lm}(\theta ,\varphi ),$$

$$t_{\text{ret}}\equiv t-\left[D+2m_0\ln \left(\frac{D}{2m_0}-1\right)\right],$$

$$D{\Psi }_4^{l,m,\infty }(t_{\text{ret}})$$

$$\displaystyle-\frac{(l-1)(l+2)}{2}\int^{t_{\text{ret}}}\Psi_{4}^{l,m}(t^{\prime},r_{0})dt^{\prime}\Big{]},$$

$$h^{l,m,\infty }(t_{\text{ret}})$$

$$=-{\int }^{t_{\text{ret}}}d{t}^{\prime }{\int }^{t^{\prime }}{\Psi }_4^{l,m,\infty }(t^{\prime \prime })d{t}^{\prime \prime }.$$

$$h^{l,m,\infty }(t_{\text{ret}})=\int d{f}^{\prime }\frac{{\tilde{\Psi }}_4^{l,m,\infty }(f^{\prime })}{(2\pi \max [f^{\prime },f_{\text{cut}}]{)}^2}\exp (2\pi i{f}^{\prime }{t}_{\text{ret}}),$$

$$h^{l,m,\infty }(t_{\text{ret}})=A^{l,m,\infty }(t_{\text{ret}})e^{-i{\Phi }^{l,m}(t_{\text{ret}})},$$

$$f_{\text{GW}}=\frac{1}{2\pi }\mathrm{Im}\left(\frac{h^{\ast 2,2,\infty }{\dot{h}}^{2,2,\infty }}{|h^{2,2,\infty }{|}^2}\right),$$

$$\frac{d{E}_{\text{GW}}^{l,m}}{dt}$$

$$\frac{d{J}_{\text{GW}}^{l,m}}{dt}$$

$$\displaystyle\times\int^{t}dt^{\prime}\int^{t^{\prime}}dt^{\prime\prime}\Psi_{4}^{l,m,\infty}(t^{\prime\prime})\Big{]}.$$

$$E_{\text{GW}}^{l,m}$$

$$J_{\text{GW}}^{l,m}$$

$$\delta {\Phi }^{\text{shift}}(t_{\text{ret}};\text{EOS1},\text{EOS2},N)$$

$$={\Phi }^{2,2}(t_{\text{ret}};EOS1,N)-{\Phi }^{2,2}(t_{\text{ret}};EOS2,N)$$

$$\delta {\Phi }^{\text{error}}(t_{\text{ret}};\text{EOS},N_1,N_2)$$

$$={\Phi }^{2,2}(t_{\text{ret}};\text{EOS},N_1)-{\Phi }^{2,2}(t_{\text{ret}};\text{EOS},N_2),$$

$${\Phi }^{2,2}(t_{\text{peak}};\text{EOS},N)$$

$$={\Phi }_{\text{peak}}^{2,2,\infty }(N_{\text{max}})-\Delta {\Phi }_{\text{peak}}^{2,2}(N_{\text{max}}){\left(\frac{N_{\text{max}}}{N}\right)}^p,$$

$$\delta {\Phi }^{\text{error}}(t_{\text{peak}};\text{EOS},N_{\text{max}},N)$$

$$=\Delta {\Phi }_{\text{peak}}^{2,2}(N_{\text{max}})\left[{\left(\frac{N_{\text{max}}}{N}\right)}^p-1\right].$$

$${\tilde{h}}_{+,\times }(f)={\int }_{t_{\mathrm{i}}}^{t_{\mathrm{f}}}{h}_{+,\times }^{2,2,\infty }(t)e^{-2\pi ift}dt,$$

$${\tilde{h}}_{+}\left(f\right)=A\left(f\right){\mathrm{e}}^{-i\Psi \left(f\right)}.$$

$$A^{\mathrm{tidal}}\left(f\right)=A\left(f\right)-A_{\mathrm{BBH}}\left(f\right)$$

$${\Psi }^{\mathrm{tidal}}\left(f\right)=\Psi \left(f\right)-{\Psi }_{\mathrm{BBH}}\left(f\right),$$

$${\Psi }_{\mathrm{model}}^{\mathrm{tidal}}=\frac{3}{128\eta }\left[-\frac{39}{2}\tilde{\Lambda }\left(1+a{\tilde{\Lambda }}^{2/3}{x}^p\right)\right]x^{5/2}$$

$$\times \left(1+\frac{3115}{1248}x-\pi x^{3/2}+\frac{28024205}{3302208}{x}^2-\frac{4283}{1092}\pi x^{5/2}\right)$$

$$A_{\mathrm{model}}^{\mathrm{tidal}}$$

$$\times \left(-\frac{27}{16}{x}^5-\frac{449}{64}{x}^6-b{x}^q\right)$$

$$\Delta \Psi (f)={\Psi }_{\mathrm{Hybrid}}^{\mathrm{tidal}}(f)-{\Psi }_{\mathrm{model}}^{\mathrm{tidal}}(f)-2\pi f{t}_0+{\varphi }_0,$$

$$\Delta A(f)/A(f)=(A_{\mathrm{Hybrid}}^{\mathrm{tidal}}(f)-A_{\mathrm{model}}^{\mathrm{tidal}}(f))/A_{\mathrm{model}}(f),$$

Full text

I Introduction

On August 17, 2017, advanced LIGO TheLIGOScientific:2014jea and advanced Virgo TheVirgo:2014hva detected gravitational waves from a binary neutron star (BNS) merger, GW170817, for the first time TheLIGOScientific:2017qsa . In this event, not only gravitational waves but also the electromagnetic signals in the gamma-ray Monitor:2017mdv ; Goldstein:2017mmi ; Savchenko:2017ffs , ultraviolet-optical-infrared Evans:2017mmy ; Drout:2017ijr ; Kilpatrick:2017mhz ; Kasliwal:2017ngb ; Nicholl:2017ahq ; Utsumi:2017cti ; Tominaga:2017cgo ; Chornock:2017sdf ; Arcavi:2017vbi ; Diaz:2017uch ; Shappee:2017zly ; Coulter:2017wya ; Soares-Santos:2017lru ; Valenti:2017ngx ; Pian:2017gtc ; Smartt:2017fuw , X-ray Haggard:2017qne ; Margutti:2017cjl ; Troja:2017nqp , and radio bands Alexander:2017aly ; Hallinan:2017woc ; Margutti:2018xqd ; Dobie:2018zno ; Mooley:2017enz ; Mooley:2018dlz were detected. This monumental event GW170817, GRB170817A, and AT2017gfo heralded the opening of the multi-messenger astrophysics. Furthermore, advanced LIGO and advanced Virgo have started a new observation run, O3, from April 2019 and a new BNS merger event, GW190425, was reported Abbott:2020uma and 7 candidates of a BNS merger as of Feb. 17, 2020, have been detected GCN .

One noteworthy finding in GW170817 is that tidal deformability of the neutron star (NS) was constrained for the first time. Due to a tidal field generated by a companion, NSs in a binary system could be deformed significantly in the late inspiral stage Flanagan:2007ix . The response to the tidal field, the tidal deformability, is imprinted as a phase shift in gravitational waves and its measurement gives a constraint on the equation of state (EOS) of NSs because the tidal deformability depends on EOSs. GW170817 constrained the binary tidal deformability in the range of $100\lesssim\tilde{\Lambda}\lesssim 800$ with the binary total mass of $2.73^{+0.04}_{-0.01}M_{\odot}$ TheLIGOScientific:2017qsa ; Abbott:2018exr ; De:2018uhw ; Abbott:2018wiz where the precise value depends on the analysis methods.

To extract information of the tidal deformability from observed gravitational wave data, a high precision template for gravitational waveforms plays an essential role. Numerical relativity simulation is the unique tool to derive high-precision gravitational waveforms in the late inspiral stage during which the gravitational-wave phase shift due to the tidal deformation becomes prominent. During this stage, any analytic techniques break down. Dietrich and his collaborators constructed a gravitational wave template for the inspiral stage based on the numerical relativity simulations in a series of papers Dietrich:2015pxa ; Dietrich:2017feu ; Dietrich:2017aum ; Dietrich:2018uni ; Dietrich:2018phi ; Dietrich:2019kaq and their template was used in gravitational wave data analysis by LIGO Scientific and Virgo Collaborations to infer the tidal deformability from GW170817 Abbott:2018exr . However, the residual phase error caused mainly by the finite grid resolution in their simulations is $\approx 0.5$–$2.3$ rad Dietrich:2019kaq . The phase error of $O(1)$ rad could be an obstacle to construct a high-quality inspiral gravitational waveform template (see also Refs. Haas:2016cop ; Foucart:2018lhe ).

In Ref. Kiuchi:2017pte , we tackled this problem by using our numerical relativity code SACRA-MPI and performed long-term simulations with the highest grid resolution to date (see also Refs. Hotokezaka:2013mm ; Hotokezaka:2015xka ; Hotokezaka:2016bzh ; Shibata:2005xz for our effort in the early stage of this project). In our numerical results, the gravitational-wave phase error caused by the finite grid resolution is less than $0.5$ rad for $31$–$32$ inspiral gravitational wave cycles. On the basis of these high-precision gravitational waveforms, Ref. Kawaguchi:2018gvj presented a waveform template, the SACRA inspiral gravitational waveform template, of BNS mergers. Specifically, we multiply the tidal-part phase of the $2.5$ Post-Newtonian (PN) order derived in Ref. Damour:2012yf by a correction term composed of the PN parameter and the binary tidal deformability. Then, we validated it by confirming that it reproduces the high-precision gravitational waveforms derived in Ref. Kiuchi:2017pte . We also validated a correction term in the tidal-part amplitude of the $1$ PN order derived in Refs. Damour:2012yf ; Vines:2011ud .

In Refs. Kiuchi:2017pte ; Kawaguchi:2018gvj , we performed simulations for a limited class of BNS systems, i.e., two equal-mass and two unequal-mass systems. Thus, the applicable range of the SACRA inspiral gravitational waveform template has not quantified precisely yet. In this paper, we derive a number of gravitational waveforms from BNS mergers by performing numerical-relativity simulations in a wider parameter space for EOSs, binary total mass, and mass ratio than that in the previous papers Kiuchi:2017pte ; Kawaguchi:2018gvj . For each binary parameter, we perform an in-depth resolution study to assess the accuracy of our waveforms. On the basis of newly derived high-precision gravitational waveforms, we validate the template.

In addition, we analyze post-merger gravitational wave signals derived in this paper. The post-merger signal in GW170817 has not been detected Abbott:2017dke , but a post-merger signal could be detected in near future for the nearby events or in the third generation detectors such as Einstein Telescope or Cosmic Explorer Punturo:2010zz ; Evans:2016mbw . The signal could bring us information of the EOS complementary to that imprinted in the late inspiral signal. To extract such information, we should explore a heuristic relation between post-merger signals and the tidal deformability/NS radius in numerical relativity simulations. In several previous papers, such an attempt has been made Rezzolla:2016nxn ; Read:2013zra ; Zappa:2017xba ; Bauswein:2011tp ; Bauswein:2012ya ; Bernuzzi:2014owa ; Bernuzzi:2015rla . However, systematics contained in these relations are unclear because of the lack of resolution study, the approximate treatment of relativistic gravity, the lack of the estimation for the systematics with the uncertainty of the NS EOS, and the narrow range of the BNS parameter space explored. In this paper, we assess to what extent the proposed universal relations between the post-merger gravitational wave signal and tidal deformability/NS radius Rezzolla:2016nxn ; Read:2013zra ; Zappa:2017xba ; Bauswein:2011tp ; Bauswein:2012ya ; Bernuzzi:2014owa ; Bernuzzi:2015rla hold.

To stimulate an independent attempt by other researchers for constructing a gravitational waveform template based on the numerical relativity simulations and/or to stimulate a comparison to numerical relativity waveforms derived by other groups, we release our simulation data on a website SACRA Gravitational Waveform Data Bank DB .

This paper is organized as follows. Section II describes our method, grid setup, and initial condition of the simulations. Section III is devoted to describing the accuracy of inspiral gravitational waveforms. Section IV presents validation of the SACRA inspiral gravitational waveform template. Section V describes the assessment of the universal relations of the post-merger signals. This section also presents the energy and angular momentum carried by gravitational waves. We summarize this paper in Sec. VI. Throughout this paper, we employ the geometrical unis of $c=G=1$ where $c$ and $G$ are the speed of light and the gravitational constant, respectively.

II Method, grid setup, and initial models

II.1 Method and grid setup

We use our numerical relativity code, SACRA-MPI Yamamoto:2008js ; Kiuchi:2017pte , to simulate a long-term inspiral stage of BNS up to early post-merger. SACRA-MPI implements the Baumgarte-Shapiro-Shibata-Nakamura-puncture formulation SN ; BS ; Capaneli ; Baker , locally incorporating a Z4c-type constraint propagation prescription Hilditch:2012fp , to solve Einstein’s equation. We discretize the field equation with the 4th-order accuracy in both the space and time. We also apply the 4th-order lop-sided finite difference scheme for the advection term Bruegmann:2006at .

In SACRA-MPI, a conservation form of general relativistic hydrodynamics equations is employed and we implement a high-resolution shock capturing scheme proposed by Kurganov and Tadmor Kurganov together with the 3rd-order accurate cell reconstruction Colella:1982ee .

We also implement the Berger-Oliger type adaptive mesh refinement (AMR) algorithm BergerOliger to enlarge a simulation domain to a local wave zone of gravitational waves while guaranteeing a high spatial grid resolution around NSs. A simulation domain consists of two sets of the 4 Cartesian AMR domains which follow orbital motion of each of NSs and the 6 Cartesian AMR domains whose center is fixed to the coordinate origin throughout all the simulations. The grid spacing of a coarser refinement level is twice as large as that of its finer refinement level. Thus, the grid spacing of a refinement level $l$ is given by $\Delta x_{l}=L/(2^{l}N)$ with $l=0,1,\cdots 9$. $L$ denotes the distance from the coordinate origin to the outer boundary along each coordinates axis. $N$ is an even number and each of the AMR domains possesses the grid point $(2N+1,2N+1,N+1)$ in the $(x,y,z)$ directions where we assumed the orbital plane symmetry.

In this work, we performed simulations with $N=182,150,130,110,102,$ and $90$ for all the systems to check the convergence of gravitational waveforms with respect to the grid resolution. The values of $L$ and $\Delta x_{9}$ are summarized in Table 1.

II.2 Binary system parameters and gravitational wave extraction

Table 1 shows the list of the binary systems as well as the grid setup for the simulations.

II.2.1 Equation of state

Following the previous papers Kiuchi:2017pte ; Kawaguchi:2018gvj , we employ a parameterized piecewise polytropic EOS to describe the NS matter rlof2009 . Specifically, we assume that the pressure and specific internal energy consist of two segments with respect to the rest-mass density:

[TABLE]

with $i=0,1$, $\rho_{0}=0{\rm~{}g~{}cm^{-3}}$, and $\rho_{2}=\infty$. $\rho_{1}$ is the rest-mass density which divides the pressure and specific internal energy into the two segments. Given the adiabatic indices $\Gamma_{0},\Gamma_{1}$ and one of the polytropic constants $\kappa_{0}$, the other polytropic constant $\kappa_{1}$ is calculated from the continuity of the pressure at $\rho=\rho_{1}$ by $\kappa_{0}\rho_{1}^{\Gamma_{0}}=\kappa_{1}\rho_{1}^{\Gamma_{1}}$. $\Delta\epsilon_{1}$ is also calculated from the continuity of the specific internal energy at $\rho=\rho_{1}$ by $\kappa_{0}\rho_{1}^{\Gamma_{0}-1}/(\Gamma_{0}-1)=\kappa_{1}\rho_{1}^{\Gamma_{1}-1}/(\Gamma_{1}-1)+\Delta\epsilon_{1}$. Note that $\Delta\epsilon_{0}=0$. Following Ref. rlof2009 , we fix $\Gamma_{0}=1.3562395$, $\Gamma_{1}=3$, and $\kappa_{0}=3.594\times 10^{13}$ in cgs units. By varying the remaining parameter $\rho_{1}$ for a wide range as shown in Table 2, we can derive plausible NS models with a variety of the radii and tidal deformability (see Table 3).

In addition to the piecewise polytropic EOS, we employ one tabulated EOS, SFHo Steiner:2012rk . To model an EOS for cold NS, we simply set $T=0.1$ MeV which is the minimum temperature in the table of SFHo EOS. We also impose the neutrinoless low-temperature $\beta$-equilibrium condition to set the value of $Y_{e}$. Then, the original tabulated EOS is reduced to a one dimensional SFHo (tabulated) EOS, i.e., $P_{\text{cold}}(\rho)$ and $\epsilon_{\text{cold}}(\rho)$ (see also Table 3 for the NS radius and tidal deformability).

During simulations (in particular for the post-merger stage), we employ a hybrid EOS to capture the shock heating effect. Specifically, we assume that the pressure consists of the cold and thermal parts:

[TABLE]

where $\epsilon$ is the specific internal energy and we assumed that the thermal part is described by the $\Gamma$-law EOS with the index $\Gamma_{\text{th}}$. Following Refs. Kiuchi:2017pte ; Kawaguchi:2018gvj , we fix $\Gamma_{\text{th}}=1.8$. We note that gravitational waveforms for the post-merger stage depend on the value of $\Gamma_{\rm th}$ Shibata:2005ss , although inspiraling waveforms do not. Since the major purpose of the present paper is to derive the accurate inspiraling waveforms, the choice of $\Gamma_{\rm th}$ does not have any essential importance. On the other hand, it has been long known that the post-merger waveform depends strongly on this value (see, e.g., Ref. Shibata:2005ss ). Thus, we have to keep in mind that the systematics exist due to the uncertainty of this value Carbone:2019pkr .

II.2.2 Binary systems

In this paper, we consider 6 irrotational binary systems assuming that NSs have no spin before merger. We fix a chirp mass, ${\cal M}_{c}$, and symmetric mass ratio, $\eta$, to be $({\cal M}_{c},\eta)=(1.1752M_{\odot},0.2500)$, $(1.1752M_{\odot},0.2485)$, $(1.1752M_{\odot},0.2455)$, $(1.1752M_{\odot},0.2450)$, $(1.0882M_{\odot},0.2470)$, and $(1.0882M_{\odot},0.2440)$. With this setting, gravitational masses of a less massive and massive components for the infinite orbital separation is $(m_{1},m_{2})=(1.35M_{\odot},1.35M_{\odot})$, $(1.25M_{\odot},1.46M_{\odot})$, $(1.18M_{\odot},1.55M_{\odot})$, $(1.17M_{\odot},1.56M_{\odot})$, $(1.12M_{\odot},1.40M_{\odot})$, and $(1.07M_{\odot},1.46M_{\odot})$ (see Table 1). For the SFHo (tabulated) EOS, we only consider the equal-mass binary system with $m_{1}=1.35M_{\odot}$ and $m_{2}=1.35M_{\odot}$.

Table 1 also shows the binary tidal deformability for all the binary systems Wade:2014vqa ; Favata:2013rwa :

[TABLE]

where $\Lambda_{1}(\Lambda_{2})$ is the tidal deformability of the less massive (massive) component. The value of the tidal deformability in this paper covers a wide range of $\approx 300$–$1800$.

Figure 1 plots the BNS systems simulated for long durations by our group to date. For the SFHo (tabulated) EOS case, an interpolation of the thermodynamic variables is necessary in the simulations. Because we implement the linear interpolation scheme for this purpose, the associated truncation error can be a non-negligible error source for generating high-precision gravitational waveforms. This system is used to assess the error budget possibly caused by employing tabulated EOS (see also Ref. Foucart:2019yzo for the gravitational-wave phase error stemming from different analytical descriptions of the EOSs).

We name all the systems according to the EOS, the mass of the less massive component, and that of the massive component. For example, 15H125-146 refers to the system with 15H EOS, $m_{1}=1.25M_{\odot}$, and $m_{2}=1.46M_{\odot}$. We set the initial orbital angular velocity to be $m_{0}\Omega_{0}=0.0150$–$0.0155$ with $m_{0}=m_{1}+m_{2}$. With this, the BNSs experience $15$–$16$ orbits before the onset of merger for all the systems.

To generate a high-precision inspiral waveform from BNS inspirals by a numerical relativity simulation, initial data with low orbital eccentricity are necessary because the orbital motion of a BNS in the late inspiral stage is circularized due to the gravitational-wave emission. We numerically obtain quasi-equilibrium sequences of the BNSs by a spectral-method library, LORENE LORENE ; Taniguchi . Then, we reduce orbital eccentricity by using the prescription in Ref. Kyutoku:2014yba . With this method, we confirm that the initial orbital eccentricity is reduced typically to $\approx 10^{-3}$ which is low enough to generate a high-precision inspiral waveform (see also Appendix in Refs. Kiuchi:2017pte ; Kawaguchi:2018gvj ).

II.3 Gravitational wave extraction

We calculate a complex Weyl scalar $\Psi_{4}$ from simulation data to derive gravitational waveforms Yamamoto:2008js . Given an extraction radius $r_{0}$, the Weyl scalar $\Psi_{4}$ is decomposed into $(l,m)$ modes with the spin-weighted spherical harmonics by

[TABLE]

where $t_{\text{ret}}$ is a retarded time defined by

[TABLE]

with $D=\sqrt{A/4\pi}$. $A$ is a proper area of the extraction sphere. We apply the Nakano’s method Nakano:2015 to extrapolate $\Psi_{4}^{l,m}$ to infinity by

[TABLE]

where $C(r_{0})$ is a function of $r_{0}$. Following Ref. Kiuchi:2017pte , we choose $D\approx r_{0}[1+m_{0}/(2r_{0})]^{2}$ and $C(r_{0})=1-2m_{0}/D$ because our coordinates are similar to isotropic coordinates of non-rotating black holes in the wave zone.

Gravitational waves of each harmonic mode are calculated by integrating $\Psi_{4}^{l,m,\infty}$ twice in time:

[TABLE]

For the time integration, we employ the fixed frequency method Reisswig:2010di by

[TABLE]

where $\tilde{\Psi}_{4}^{l,m,\infty}(f)$ is the Fourier component of $\Psi_{4}^{l,m,\infty}(t)$ and $f_{\text{cut}}$ is set to be $0.8m\Omega_{0}/(2\pi)$.

To check the convergence with respect to the extraction radius $r_{0}$, we repeat this analysis for $r_{0}$ $=244\,m_{0}$, $199\,m_{0},$ and $155\,m_{0}$ for ${\cal M}_{c}=1.1752M_{\odot}$ and $r_{0}=262\,m_{0}$, $213\,m_{0}$, and $156\,m_{0}$ for ${\cal M}_{c}=1.0882M_{\odot}$ (see Table 1).

In general, gravitational waves for each $(l,m)$ mode are decomposed into the amplitude and phase as

[TABLE]

and instantaneous gravitational-wave frequency is defined by $d\Phi^{l,m}/dt_{\text{ret}}$. In Sec. III, we explore the accuracy of the gravitational-wave phase of the $(l,m)=(2,2)$ mode, and simply refer to $\Phi^{2,2}$ as the gravitational-wave phase. With Eq. (8), the instantaneous frequency of the $(l,m)=(2,2)$ mode is calculated by

[TABLE]

where the asterisk symbol denotes the complex conjugate of $h^{2,2,\infty}$.

We also calculate the energy and angular momentum flux due to gravitational-wave emission by Shibata:text

[TABLE]

Thus, the energy and angular momentum carried by gravitational waves are calculated by

[TABLE]

where $t_{\text{sim}}$ denotes the time we terminate the simulations.

III Accuracy of waveforms

To date, we have simulated for long durations 46 binary systems with 6 grid resolutions for each model. 26 binary systems are newly reported in this paper and 20 binary systems have been reported in Refs. Kiuchi:2017pte ; Kawaguchi:2018gvj . Our waveform data are publicly available on the website:

SACRA Gravitational Waveform Data Bank DB .

On the website, the waveform data are tabulated according to the system name, dimensionless initial orbital angular velocity, and grid resolution. For example, 15H_135_135_00155_182 refers to the employed EOS as 15H, $m_{1}=1.35M_{\odot}$, $m_{2}=1.35M_{\odot}$, $m_{0}\Omega_{0}=0.0155$, and $N=182$ (see also Table 1). A user can download the data for $\Psi^{2,2}_{4}(t_{\text{ret}},r_{0})$ extracted at several values of $r_{0}$ and $h^{2,2,\infty}_{+,\times}(t_{\text{ret}})$ from the link on the system name.

III.1 Overview of physical and numerical phase shifts

First we briefly illustrate that the waveforms depend on EOSs and each mass of binary systems. The top panel of Fig. 2 shows the dependence of the gravitational waveforms on the EOSs for the binary systems with $m_{1}=1.12M_{\odot}$, $m_{2}=1.40M_{\odot}$ and $N=182$. It shows that the systems with the larger values of $\tilde{\Lambda}$ merge earlier than those with the smaller values of $\tilde{\Lambda}$ because the tidal force due to its companion induces the quadrupole moment and the resultant attractive force accelerates the orbital shrinkage. The bottom panel of Fig. 2 shows the dependence of the gravitational waveforms on the symmetric mass ratio for the binary systems with 15H125-125, 15H112-140, and 15H107-146 with $N=182$. It shows that the systems with the larger values of $\eta$ merge earlier than those with the smaller values of $\eta$ because the emissivity of gravitational waves decreases as the symmetric mass ratio decreases Blanchet:2013haa .

The top panel of Fig. 3 shows the dependence of the gravitational waveforms on the grid resolutions for 15H112-140 with $N=182,110$, and $N=90$. Errors in the amplitude and phase caused by the finite grid resolution become prominent for the late inspiral and post-merger stages. The bottom panel of Fig. 3 plots the phase shift among the systems of different EOSs for $m_{1}=1.12M_{\odot}$, $m_{2}=1.40M_{\odot}$, and $N=182$. The phase shift is defined by

[TABLE]

where $\Phi^{2,2}(t_{\text{ret}};\text{EOS},N)$ is the gravitational-wave phase for $l=|m|=2$ mode derived from a simulation with employing EOS and the grid number $N$. Because we compare the phase among models with common masses of components, we omit the masses from the argument. The shaded region shows the evolution of the phase error defined by

[TABLE]

where $N_{1}$ and $N_{2}$ denote the employed grid numbers. The red shaded region shows $\delta\Phi^{\text{error}}(t_{\text{ret}};15{\rm H},150,182)$ and the blue shaded region shows $\delta\Phi^{\text{error}}(t_{\text{ret}};15{\rm H},90,182)$, respectively, for 15H112-140. The overlapped region has purple color. The vertical dashed line denotes the peak time, $t_{\text{peak}}$, at which the gravitational-wave amplitude becomes maximal for 15H112-140 with $N=182$. Just after the peak time, burst-type gravitational waves are emitted for a short time as shown in the upper panel of Fig. 3, i.e., for $58~{}{\rm ms}\lesssim t_{\text{ret}}\lesssim 59~{}{\rm ms}$. These waves cause very rapid increase in phase during this short-term interval and consequently the phase shift shows very rapid increase. This feature can be also seen in the phase error and the very rapid increase appears later in $\delta\Phi^{\text{error}}(t_{\text{ret}};15{\rm H},150,182)$ than in $\delta\Phi^{\text{error}}(t_{\text{ret}};15{\rm H},90,182)$ because the peak time becomes later with improving the grid resolution.

The phase shift and the phase error up to the peak time are comparable, in particular, for the case with the coarser grid resolution. Therefore, unless a convergence study is sufficiently carried out, a capability of inspiral waveform models to measure the tidal deformability is unclear. This is also the case for the the post-merger stage. In particular, the phase evolution loses the convergence as found in the bottom panel of Fig. 3, i.e., $\delta\Phi^{\text{error}}(t_{\text{ret}};{\rm 15H},150,182)$ (red shaded region) is larger than $\delta\Phi^{\text{error}}(t_{\text{ret}};{\rm 15H},90,182)$ (blue shaded region). Therefore, time-domain post-merger gravitational waves derived in numerical-relativity simulations are not very reliable. Instead, we will discuss the post-merger signal in terms of the energy and angular momentum carried by gravitational waves and their spectrum amplitude. These quantities are calculated by a time integration of the gravitational waveforms and the convergence in the phase could be subdominant as discussed in Sec. V.

III.2 Estimation of the residual phase error in gravitational waves

Following Refs. Kiuchi:2017pte ; Kawaguchi:2018gvj , we estimate a residual gravitational-wave phase error at the peak time in the simulations. The left panel of Fig. 4 plots evolution of the phase error, $\delta\Phi^{\text{error}}(t_{\text{ret}};{\rm B},N,182)$, with $N=150,130,110,102$, and $90$ for B107-146. The vertical dashed line denotes the peak time for B107-146 with $N=182$. Although the phase error is accumulated with time, its value at the peak time decreases as improving the grid resolution. We estimate the residual phase error by assuming that the gravitational-wave phase at the peak time obeys the following functional form,

[TABLE]

where $\Phi_{\text{peak}}^{2,2,\infty}(N_{\text{max}})$ and $p$ denote the gravitational-wave phase at the peak time in the continuum limit of the finite difference $(N\to\infty)$ and an order of the convergence, respectively. $\Delta\Phi^{2,2}_{\text{peak}}(N_{\text{max}})$ should be recognized as the residual phase error for the simulation with $N=N_{\text{max}}$. $N_{\text{max}}$ denotes a reference value of $N$ to estimate unknown quantities $\Phi^{2,2,\infty}_{\text{peak}}(N_{\text{max}})$, $\Delta\Phi^{2,2}_{\text{peak}}(N_{\text{max}})$, and $p$. For example, with $N_{\text{max}}=182$, these unknowns are obtained by fitting the simulation results of $N=150,130,110,102,$ and $90$ with Eq. (16) given an EOS, a chirp mass, and a symmetric mass ratio.

The right panel of Fig. 4 plots the gravitational-wave phase error at the peak time, $\delta\Phi^{\text{error}}(t_{\text{peak}};\text{B},N_{\text{max}},N)$, as a function of $1/N^{p}$ with a reference grid number $N_{\text{max}}$ and $N=90,102,\cdots,N_{\text{max}}$. Assuming Eq. (16), the phase error at the peak time in a binary system is given as

[TABLE]

The values of $\Delta\Phi^{2,2}_{\text{peak}}(N_{\text{max}})$ and $p$ are shown in the legend of this plot. It is clear that the order of the convergence $p$ is improved and the residual gravitational-wave phase error is reduced as increasing $N_{\text{max}}$.

Table 4 summarizes the residual phase error and the order of the convergence of the gravitational-wave phase at the peak time for all the systems. We estimate the residual phase error with respect to three reference values of $N_{\text{max}}$ as $182,150,$ and $130$. In some systems, the residual phase error and the order of the convergence show an irregular behavior. That is, the residual phase error (the order of convergence) for $N_{\text{max}}=130$ happens to be smaller (higher) than that for $N_{\text{max}}=150$. Nonetheless, the residual phase error (the order of convergence) for $N_{\text{max}}=182$ is smaller (higher) than that for $N_{\text{max}}=150$ except for 125H125-146. Thus, we adopt the values for $N_{\text{max}}=182$ as the residual phase error in our waveforms and it is in the range of $\approx 0.1$–$0.5$ rad.

For the SFHo (tabulated) EOS, we find that the residual phase error still remains within sub-radian accuracy. Because SFHo135-135 and HB135-135 have nearly identical values of $\tilde{\Lambda}$ Kiuchi:2017pte , The phase error due to the tabulated EOS is estimated by comparing the results for them. For HB135-135, the residual phase error and the order of the convergence are $(\Delta\Phi^{2,2}_{\text{peak}}(182),p)=(0.17~{}\text{rad},3.6)$, $(\Delta\Phi^{2,2}_{\text{peak}}(150),p)=(0.48~{}\text{rad},3.2)$, and $(\Delta\Phi^{2,2}_{\text{peak}}(130),p)=(2.0~{}\text{rad},1.7)$ Kiuchi:2017pte . For SFHo135-135, the residual phase error and the order of the convergence are $(\Delta\Phi^{2,2}_{\text{peak}}(182),p)=(0.43~{}\text{rad},2.3)$, $(\Delta\Phi^{2,2}_{\text{peak}}(150),p)=(0.76~{}\text{rad},2.2)$, and $(\Delta\Phi^{2,2}_{\text{peak}}(130),p)=(0.33~{}\text{rad},4.2)$, respectively. Thus, the system with the SFHo (tabulated) EOS has slightly larger residual phase error than with the piecewise polytropic EOS. This indicates that the linear interpolation of the thermodynamic quantities could cause a phase error of $\approx 0.2$–$0.3$ rad. Nonetheless, it is encouraging that our waveforms have the sub-radian accuracy even for the SFHo (tabulated) EOS. For more detailed estimate of the error budget due to tabulated EOSs, we need to perform BNS simulations with a wide class of tabulated EOSs. In particular, we speculate that the phase error when using a tabulated EOS with a phase transition could be even larger.

IV Inspiral gravitational waveform modeling

IV.1 SACRA inspiral gravitational waveform template

In the previous paper Kawaguchi:2018gvj , we developed a frequency-domain gravitational waveform model for inspiralling BNSs (with $l=|m|=2$) based on high-precision numerical-relativity data. In this section, we extend the examination of the inspiral waveform model to a parameter space wider than the previous papers Kiuchi:2017pte ; Kawaguchi:2018gvj by employing new waveforms obtained in this paper.

Before moving on to the comparison, we briefly review our inspiral waveform model. First we calculate the Fourier component for the quadrupole mode of gravitational waves for all the systems by

[TABLE]

where $t_{\rm i}$ and $t_{\rm f}$ are the initial and final time of the waveform data, respectively. Then, we decompose $\tilde{h}_{+}(f)$ in Eq. (18) into the frequency-domain amplitude, $A\left(f\right)$, and phase, $\Psi\left(f\right)$, (with an ambiguity in the origin of the phase) by

[TABLE]

We only use $h^{2,2,\infty}_{+}$ for modeling the inspiral gravitational waveforms because the difference between $h^{2,2,\infty}_{+}$ and $h^{2,2,\infty}_{\times}$ is approximately only the phase difference of $\pi/2$. We define the corrections due to the NS tidal deformation to the gravitational-wave amplitude and phase by

[TABLE]

and

[TABLE]

respectively. Here, $A_{\rm BBH}\left(f\right)$ and $\Psi_{\rm BBH}\left(f\right)$ are the gravitational-wave amplitude and phase of a binary black hole with the same mass as the BNS, respectively (hereafter referred to as the point-particle parts: see Ref. Kawaguchi:2018gvj for details).

Our numerical-relativity waveforms only contain the waveforms for the frequency higher than $\approx 400\,{\rm Hz}$. Thus, we employ the effective-one-body waveforms of Refs. Hinderer:2016eia ; Steinhoff:2016rfi ; Lackey:2018zvw ; Taracchini:2013rva (SEOBNRv2T) to model the low-frequency part waveforms, in which the effect of dynamical tides is taken into account, and construct hybrid waveforms combining them with the numerical-relativity waveforms. The hybridization of the waveforms is performed in the time-domain by the procedure described in Refs. Hotokezaka:2016bzh ; Kawaguchi:2018gvj and we set the matching region to be from $t_{\text{ret}}\approx 7.38$ ms to $14.78$ ms. After the hybridization, the waveforms are transformed into the frequency domain employing Eq. (18), and the tidal-part amplitude and phase are extracted by Eqs. (20) and (21).

For modeling the tidal-part phase and amplitude, we employ the following functional forms motivated by the 2.5 PN order formula Damour:2012yf :

[TABLE]

for the phase correction and

[TABLE]

for the amplitude correction where $D_{\rm eff}$ is the effective distance to the binary Hotokezaka:2016bzh and $x\equiv(\pi m_{0}f)^{2/3}$. $a$, $p$, $b$, and $q$ are the free parameters of the models. To focus on the inspiral waveform and to avoid the contamination from the post-merger waveforms of high frequency, which would have large uncertainties, we restrict the gravitational-wave frequency range in $10$–$1000\,{\rm Hz}$. The fitting parameters were determined by employing the hybrid waveforms of 15H125-125, which has the largest value of binary tidal deformability in the systems studied in the previous study Kawaguchi:2018gvj . By performing the least square fit with respect to the phase shift and relative difference of the amplitude, we obtained $a=12.55$, $p=4.240$, $b=4251$, and $q=7.890$.

In Ref. Kawaguchi:2018gvj , the validity of the inspiral waveform model was examined employing hybrid waveforms which were not used for the parameter determination. We should stress again that the parameters $a,p,b$, and $q$ in Eqs. (22) and (23) were determined by the particular system 15H125-125. We found that the tidal-part waveform model always reproduced the tidal-part phase and amplitude of the hybrid waveforms within $\sim 0.1\,{\rm rad}$ and $15\%$, respectively, for the equal-mass and unequal-mass cases with ${\cal M}_{\rm chirp}=1.1752\,M_{\odot}$ and the equal-mass cases with ${\cal M}_{\rm chirp}=1.0882\,M_{\odot}$, covering the parameter space of $0.244\leq\eta\leq 0.250$ and $300\lesssim\tilde{\Lambda}\lesssim 1800$.

IV.2 Validation of SACRA inspiral gravitational waveform template

While the validity of our inspiral waveform model was already examined in the most interesting part of the parameter space of BNSs Kawaguchi:2018gvj , there still remain some important cases which were not examined in the previous study Kawaguchi:2018gvj . First, the dependence of the error of the tidal correction on the mass ratio has to be checked for less massive BNSs. While unequal-mass cases with total mass of $\approx\,2.7\,M_{\odot}$ were checked in the previous study Kawaguchi:2018gvj , it is important to check whether our inspiral waveform models are also applicable to unequal-mass cases with smaller total mass, for which tidal effect is enhanced due to increase of tidal deformability. Second, the systematics due to simplification on the high-density part of the EOS should be checked. For the inspiral waveforms, we expect that the high-density part of the EOS has a minor effect, and thus, we employ simplified two-piecewise polytropic EOS models. However, we should confirm that this assumption is indeed valid.

To check the points listed above, we compare our inspiral waveform model with hybrid waveforms employing the numerical-relativity waveforms obtained in this paper. Hybrid waveforms are constructed in the same manner as in the previous study Kawaguchi:2018gvj employing the SEOBNRv2T waveforms as the low-frequency part waveforms. In particular, we focus on the validity of the tidal correction model to the waveform, comparing it with the tidal-part phase and amplitude of the hybrid waveforms computed based on Eqs. (20) and (21) using the SEOBNRv2 waveforms with no-tides as the point-particle parts.

Figures 5 and 6 show the difference of the tidal-part phase and amplitude between our inspiral waveform model (22) and (23) and the hybrid waveforms for the models with ${\cal M}_{c}=1.1752\,M_{\odot}$ and ${\cal M}_{c}=1.0882\,M_{\odot}$. Here, the phase difference between the tidal-part phase of hybrid waveforms, $\Psi_{\rm Hybrid}^{\rm tidal}$, and that of our inspiral waveform model, $\Psi_{\rm model}^{\rm tidal}$, is computed by

[TABLE]

where $t_{0}$ and $\phi_{0}$ are the free parameters which correspond to the degrees of freedom in choosing the origins of time and phase, respectively, and are determined by minimizing $\int|\Delta\Psi(f)|^{2}df$ integrated in the range of $f=10$–$1000\,{\rm Hz}$. For the comparison of the tidal-part amplitude, relative difference of the amplitude,

[TABLE]

is computed, where $A_{\rm Hybrid}^{\rm tidal}$ and $A_{\rm model}=A_{\rm model}^{\rm tidal}+A_{\rm BBH}$ are the tidal-part amplitude of hybrid waveforms and the amplitude of the model waveforms including the point-particle part, respectively. Again, we employ the amplitude of the SEOBNRv2 waveforms with no-tides for $A_{\rm BBH}$.

The systems of mass $1.25M_{\odot}$–$1.46M_{\odot}$, $1,18M_{\odot}$–$1.55M_{\odot}$, and $1.17M_{\odot}$–$1.56M_{\odot}$ are within the parameter space which we studied in the previous study Kawaguchi:2018gvj , and thus, we expect that those waveforms are well reproduced by our inspiral waveform model. Indeed Fig. 5 shows that differences in both phase and amplitude are within the error which we observed in the previous study Kawaguchi:2018gvj . Figure 5 also shows that tidal-part phase and amplitude for system SFHo135-135 are well reproduced by our inspiral waveform model. This confirms that, at least for the frequency range and $m_{0}$ we focus on, employing an EOS whose high-density part is simplified has only a minor effect on the systematics of the model. Figure 6 shows the results in the unequal-mass cases with ${\cal M}_{c}=1.0882\,M_{\odot}$. The difference in the tidal-part phase is larger than the cases with ${\cal M}_{c}=1.1752\,M_{\odot}$. This is reasonable because we found that the error of tidal-part model becomes relatively large for a small mass ratio or a large value of tidal deformability in the previous study Kawaguchi:2018gvj . Nevertheless, the phase error is always smaller than $\approx 0.1\,{\rm rad}$, which is smaller than the systematics in the waveforms stemming from the finite difference as shown in the previous section. The deviation for the amplitude model is also the same level as for the models with ${\cal M}_{c}=1.1752\,M_{\odot}$.

To quantify the deviation of our inspiral waveform model from the new sets of hybrid waveforms, we calculate the mismatch between those waveforms, ${\bar{F}}$, defined by

[TABLE]

where $(\cdot|\cdot)$ and $||\cdot||$ are defined by

[TABLE]

where $f_{\rm min}=10$ Hz and $f_{\rm max}=1000$ Hz and

[TABLE]

Here, $h_{1}$ and $h_{2}$ denote the hybrid waveforms and our inspiral waveform models, respectively. The inspiral waveform model employs Eqs. (22) and (23) as the tidal part and the SEOBNRv2 waveforms with no-tides as the point-particle baseline. $S_{\rm n}$ denotes the one-sided noise spectrum density of the detector, and we employ the noise spectrum density of the ZERO_DETUNED_HIGH_POWER configuration of advanced LIGO aLIGOnoise for it.

We summarize the values of mismatch between our inspiral waveform model and hybrid waveforms in Table 5. For all the cases, the value of mismatch is smaller than $\approx 2\times 10^{-5}$. According to our previous results Kawaguchi:2018gvj , these results indicate that the the signal to noise ratio of the difference between our inspiral waveform model and hybrid waveforms are as small as $1$ even for the case in which the total signal to noise ratio is as large as $200$.

V Assessment of universal relation for late inspiral and post-merger gravitational waves

V.1 frequency and amplitude

Instantaneous gravitational-wave frequency defined by Eq. (9) at some characteristic time in the late inspiral or post-merger stage is reported to be correlated with the tidal deformability or the tidal coupling constant Read:2013zra ; Rezzolla:2016nxn ; Bernuzzi:2015rla ; Bernuzzi:2014owa . In addition, characteristic peak frequencies imprinted in the spectrum amplitude of post-merger gravitational waves are reported to be correlated with the tidal coupling constant or NS radius Rezzolla:2016nxn ; Shibata:2005xz ; Hotokezaka:2013iia ; Bauswein:2011tp . We assess these proposed universal relations using our waveform data, for which the systematic study has been conducted in a wide range of the binary parameters with a wide range of the grid resolution of the simulations. We also propose new relations in terms of the binary tidal deformability.

V.1.1 Peak frequency and binary tidal deformability relation

Reference Read:2013zra reported that the instantaneous gravitational-wave frequency (of $l=|m|=2$ mode) at the peak time $(t_{\text{peak}})$, $f_{\text{peak}}$, has a tight correlation with the binary tidal deformability $\tilde{\Lambda}$ (see also Refs. Rezzolla:2016nxn ; Bernuzzi:2015rla ; Bernuzzi:2014owa for the relation with the tidal coupling constant: In Ref. Rezzolla:2016nxn , they referred to it as $f_{\text{max}}$). Figure 7 plots the dependence of $f_{\text{peak}}$ on the grid resolution where $f_{\text{peak,ave}}$ is the average of $f_{\text{peak}}$ over the results with different grid resolutions. $f_{\text{peak}}$ does not converge perfectly with respect to the grid resolution, but the fluctuation around the averaged value is less than 2$\%$ for a wide range of the grid resolution. This is also the case for all the binary systems. Thus, we estimate a relative error due to the finite grid resolution in $f_{\text{peak}}$ to be 2$\%$ and tabulate the values of $f_{\text{peak}}$ in Table 6.

The right panel of Fig. 7 plots $m_{0}f_{\text{peak}}$ as a function of $\tilde{\Lambda}^{1/5}$. The error bar shows the systematics associated with the finite grid resolution in $f_{\text{peak}}$. We also plot the universal relations reported in Refs. Read:2013zra (black dashed line) and Rezzolla:2016nxn (black dotted line). We find that the universal relation in Ref. Rezzolla:2016nxn holds only for the symmetric binary systems with ${\cal M}_{c}=1.1752M_{\odot}$ and ${\cal M}_{c}=1.0882M_{\odot}$ (see also Table 6). Given an EOS and a chirp mass, $f_{\text{peak}}$ shifts to a lower value as the symmetric mass ratio decreases. This is attributed to following three facts. First, given the total mass $m_{0}$ and $f_{\text{GW}}$, $df_{\text{GW}}/dt$ decreases as the symmetric mass ratio decreases because the gravitational-wave luminosity is proportional to $\eta^{2}$ Blanchet:2013haa . Second, the time at which the two NSs come into contact becomes earlier as the symmetric mass ratio decreases because the less massive companion is more subject to the tidal elongation and the resultant mass accretion on the massive component starts earlier than for the symmetric binary. Third, the difference between the peak time and the contact time becomes small as the symmetric mass ratio decreases because the peak time corresponds to the moment when a dumbbell-like density structure with double dense cores formed after the contact disappears as discussed in Ref. Kiuchi:2017pte and the dumbbell-like density structure becomes less prominent in the asymmetric binary systems. Due to these effects, $f_{\text{peak}}$ becomes lower as the symmetric mass ratio decreases.

In a short summary, the $m_{0}f_{\text{peak}}$–$\tilde{\Lambda}^{1/5}$ relation depends strongly on the symmetric mass ratio and the universal relations reported in Refs. Read:2013zra and Rezzolla:2016nxn suffer from this systematics (see also Ref. Kiuchi:2017pte ). This finding is consistent with a discussion in Ref. Rezzolla:2016nxn . They mentioned that the mass asymmetry could break the universality in the $m_{0}f_{\text{peak}}$–$\tilde{\Lambda}^{1/5}$ relation for a possibly unrealistic mass ratio. We find that the realistic value of the mass ratio breaks the universality as the symmetric mass ratio adopted in this paper is consistent with that in GW170817 TheLIGOScientific:2017qsa . The scatter from the proposed universal relation in Ref. Rezzolla:2016nxn is as large as $\approx$ 18–19$\%$ at the maximum for $0.244\leq\eta\leq 0.250$.

We propose an improved fitting formula:

[TABLE]

With $\eta=0.2500$, $a_{0}(\eta)$ and $a_{1}(\eta)$ approximately reduce to be $a_{0}$ and $a_{1}$ footnote1 reported in Ref. Rezzolla:2016nxn . Figure 8 plots the improved relation with the simulation data and we confirm that the relative error between the data and the fitting formula (29) is smaller than $3\%$.

We should keep in mind that this relation could still suffer from systematics associated with physical effects that are not taken into the simulation. Because of the spin-orbit coupling, high NS spin could change $f_{\text{peak}}$ compared to the non-spinning case. NS magnetic fields also could produce systematics in Eq. (29) because at the contact of the two NSs, which occurs before the peak time, the magnetic field could be exponentially amplified by the Kelvein-Helmholtz instability within a very short timescale $\ll 1$ ms Kiuchi:2014hja ; Kiuchi:2015sga and the magnetic pressure could reach near the equipartition of the pressure locally, affecting the value of $f_{\text{peak}}$. These points should be explored in future work.

V.1.2 Peak amplitude and binary tidal deformability relation

References Read:2013zra ; Kiuchi:2017pte reported that the gravitational-wave amplitude at the peak time, $h_{\text{peak}}$, correlates with $f_{\text{peak}}$, i.e., with $\tilde{\Lambda}^{1/5}$. Because we do not find perfectly convergent result for $h_{\text{peak}}$ with respect to the grid resolution, first, we assess deviation of $h_{\text{peak}}$ relative to the averaged value of $h_{\text{peak}}$ (average of the results with different grid resolutions) in the left panel of Fig. 9 for the binary systems with $m_{1}=1.07M_{\odot}$ and $m_{2}=1.46M_{\odot}$. It is found that fluctuation around the averaged value is $\approx 1$–$2\%$. This is also the case for all the binary systems. Thus, we adopt $2\%$ as the systematics associated with the finite grid resolution in $h_{\text{peak}}$ and summarize the values of $h_{\text{peak}}$ in Table 6.

The right panel of Fig. 9 plots $Dh_{\text{peak}}/m_{0}$ as a function of $\tilde{\Lambda}^{1/5}$. The error bar shows the systematics associated with the finite grid resolution in $h_{\text{peak}}$. This figure shows that the relation depends strongly on the symmetric mass ratio. That is, the relation proposed in Refs. Read:2013zra ; Kiuchi:2017pte is not in general satisfied.

We propose a fitting formula for $Dh_{\text{peak}}/m_{0}$:

[TABLE]

Figure 10 plots the improved relation with the simulation data. We find that the relative error between the data and the fitting formula (30) is within $4\%$. Again note that this relation is calibrated in a limited class of the binary systems, i.e., non-magnetized non-spinning binary systems. We should keep in mind this point in using this relation to infer the tidal deformability from observational data.

V.1.3 $f_{1},f_{2}$ and binary tidal deformability relation

Reference Rezzolla:2016nxn reported that several gravitational-wave frequencies associated with the main peaks in the spectrum amplitude for post-merger gravitational waves correlate with the tidal coupling constant. Figures 11–13 show the spectrum amplitudes for the quadrupole mode of gravitational waves for all the systems defined by

[TABLE]

with $\tilde{h}_{+}(f)$ and $\tilde{h}_{\times}(f)$ in Eq. (18). In Figs. 11–13, the vertical dashed lines indicate the so-called $f_{1}$ frequency for the fitting formula in Ref. Rezzolla:2016nxn . This peak is a side-band peak of the main peak of $f=f_{2}$, and it is naturally understood as a result of the modulation of the main peak. According to Ref. Takami:2014tva , the remnant might be represented by a mechanical toy model composed of a rotating disk with two spheres. In this model, the two spheres, which mimic the double dense cores appearing after merger, are connected with a spring and oscillate freely (see their Fig. 17). $f_{1}$ frequency corresponds to the spin frequency when the separation between the two spheres becomes largest if we assume the angular momentum conservation. They claimed this scenario for the interpretation of $f_{1}$ frequency.

In Ref. Rezzolla:2016nxn , $f_{1}$ frequency is determined by identifying one of the main peaks in the spectrum amplitude and the spectrogram of post-merger gravitational waves. For the symmetric binary systems, $f_{1}$ peak could be identified in our numerical results using the same methods. However, the structure of the spectrum amplitude around $f=f_{1}$ depends highly on the grid resolution (see 125H135-135 and H135-135 systems for example). For a sequence with the fixed EOS and chirp mass, e.g., 15H135-135, 15H125-146, 15H121-151, 15H118-155, 15H117-156, and 15H116-158, we find it more difficult to identify $f_{1}$ peak as the symmetric mass ratio decreases. This was also pointed out in Ref. Dietrich:2015iva although their grid resolution was much lower than those in our present study and the resolution study on the spectrum amplitude of gravitational waves is not performed (see their Fig. 13). As demonstrated in Fig. 11, $f_{1}$ peak cannot be clearly identified for the asymmetric binary systems. Figure 12 shows that this is also the case for binary systems of relatively small mass $\sim 2.5M_{\odot}$ as discussed in Refs. Foucart:2015gaa ; Bauswein:2015yca ; Rezzolla:2016nxn .

We also analyze the spectrogram of post-merger gravitational waves and confirm that there is no prominent peak around $f_{\text{GW}}=f_{1}$ for the asymmetric binary systems. Therefore, we conclude that the universal relation for $f_{1}$ could be only applicable to nearly symmetric binary systems: essentially no universal relation is present. We speculate that for the asymmetric binary systems, the mechanical toy model proposed in Ref. Takami:2014tva could not describe the merger remnant because the less massive NS is tidally disrupted before merger and there is no prominent double dense cores. We also note that the method for constraining the EOS proposed in Ref. Takami:2014zpa could not be applied unless the symmetric mass ratio is measured precisely to be $0.25$ because this method relies on $f_{1}$ universal relation.

In Ref. Rezzolla:2016nxn , the peak frequency, $f_{2}$, in the spectrum amplitude footnote2 is reported to have a correlation with the tidal coupling constant. This peak frequency approximately corresponds to the f–mode oscillation of the remnant massive NS (see also Refs. Shibata:2005xz ; Hotokezaka:2013iia ; Bauswein:2011tp ; Shibata:2005ss ). The left panel of Fig. 14 plots fluctuation around the averaged value of $f_{2}$ (average of the results with different grid resolutions) for the binary systems with $m_{1}=1.12M_{\odot}$ and $m_{2}=1.40M_{\odot}$. We measure $f_{2}$ in the spectrum amplitude as a prominent peak for $f\geq 2$ kHz. The fluctuation is within $\approx 4$–$5\%$ and we find that this is also the case for all the binary systems. Thus, we adopt $5\%$ as a relative error of $f_{2}$ (see also Table 6). The right panel of Fig. 14 shows $f_{2}$ as a function of $\tilde{\Lambda}^{1/5}$. We exclude the systems which collapse to a black hole within a few ms after merger because the peak associated with $f_{2}$ is not prominent or absent in the spectrum amplitude. We also overplot the fitting formula proposed in Ref. Rezzolla:2016nxn . It is found that with this fitting formula, the scatter is $\approx 14\%$ at the maximum. Thus, we propose an improved fitting formula for $m_{0}f_{2}$;

[TABLE]

Even with this formula, the relative error is as large as $9\%$ (see also Fig. 15). This implies that even if the value of $f_{2}$ is determined precisely in the data analysis of gravitational waves, $\tilde{\Lambda}^{1/5}$ will be constrained with the error of $\approx\pm 0.1$.

V.1.4 $f_{2}$ and NS radius with $1.6M_{\odot}$ relation

References Bauswein:2011tp ; Bauswein:2012ya reported that $f_{2}$ frequency has a tight correlation with the NS radius of $1.6M_{\odot}$ (see Eq. (3) in Ref. Bauswein:2012ya ). In Ref. Hotokezaka:2013iia , we assessed their relation by using our numerical-relativity results and found that the scatter in the relation is larger than that reported in Ref. Bauswein:2012ya . We revisit this assessment because the initial orbital eccentricity reduction was not implemented in Ref. Hotokezaka:2013iia . In addition, the grid resolution in Ref. Hotokezaka:2013iia is much lower than that in this paper. These ingredients could modify the post-merger dynamics and the resulting gravitational waveforms.

Because the relation in Ref. Bauswein:2012ya holds only for symmetric binary systems of $m_{0}=2.7M_{\odot}$, we first assess this relation by employing binary systems of $(m_{1},m_{2})=(1.35M_{\odot},1.35M_{\odot})$ and found that the error is $\approx 6\%$ crust_comment . Second, we assess the relation by employing binary systems of $(m_{1},m_{2})=(1.25M_{\odot},1.46M_{\odot})$, $(1.21M_{\odot},1.51M_{\odot})$, $(1.18M_{\odot},1.55M_{\odot})$, $(1.17M_{\odot},1.56M_{\odot})$, and $(1.16M_{\odot},1.58M_{\odot})$. We found that the scatter from their fitting formula is $\approx 10\%$. Therefore, the scatter larger than that reported in Ref. Bauswein:2012ya stems from the mass asymmetry of the binary. Our numerical results suggest that the fitting formula in Ref. Bauswein:2012ya could infer the radius of the $1.6M_{\odot}$ NS within the $1$ km accuracy only if the symmetric mass ratio is well constrained to be $0.25$. Otherwise, we constrain the radius of the $1.6M_{\odot}$ NS with the accuracy of $\approx\pm 1$ km if the value of $f_{2}$ is determined precisely,

In Table 7, we summarize to what extent the so-called universal relations hold.

V.2 Energy and angular momentum

Using Eqs. (10)–(13), we calculate the energy and angular momentum carried by gravitational waves. We define $E_{\rm GW,i}^{\text{tot}}$ and $E_{\rm GW,p}~{}(J_{\rm GW,p})$ as the energy (angular momentum) emitted in the inspiral stage and in the post-merger stage, respectively. The subscripts $\rm i$ and $\rm p$ in these quantities denote the inspiral and the post-merger stage, respectively. The peak time introduced in Sec. III.1 defines the boundary between the inspiral and post-merger stages. In the following we summarize the energy and angular momentum emitted in each stage for all the systems. Their values are presented in Table 6.

V.2.1 inspiral stage

Table 6 and Fig. 16 show the energy, $E_{\rm GW,i}^{2,2}$, carried by gravitational waves with $(l,m)=(2,2)$ mode during the inspiral stage. We measure the relative error with respect to the averaged value in the left panel of Fig. 16 and find that the error relative to its averaged value of $E^{2,2}_{\rm GW,i}$ (average of the results with different grid resolutions) never exceeds $2\%$ for a wide range of the grid resolution. This is also the case for all the binary systems. Thus, we adopt this fluctuation as an error in $E_{\rm GW,i}^{2,2}$. Note that the other modes such as $(l,m)=(2,1)$ and $(3,3)$ are $\lesssim 0.1\%$ and $\lesssim 0.5\%$, respectively, of $E^{2.2}_{\rm GW,i}$.

The right panel of Fig. 16 plots $E^{\text{tot}}_{\rm GW,i}/(m_{0}\eta)$ as a function of $\tilde{\Lambda}^{1/5}$. We include the contribution due to the gravitational-wave emission during evolution from infinite separation to the initial orbital separation of the simulation, $m_{0}-M_{\text{ADM,0}}$ in Table 6, by $E^{\text{tot}}_{\rm GW,i}\approx 2E_{\rm GW,i}^{2,2}+m_{0}-M_{\text{ADM,0}}$. $M_{\text{ADM,0}}$ is the Arnowitt-Deser-Misner mass of the initial condition of the simulations. As proposed in Ref. Zappa:2017xba , this quantity correlates with the tidal coupling constant. We explicitly derive a fitting formula with the binary tidal deformability as

[TABLE]

It is reasonable that $E^{\text{tot}}_{\rm GW,i}$ decreases as $\tilde{\Lambda}$ increases because the binary systems with larger values of $\tilde{\Lambda}$ merge earlier than those with smaller values of $\tilde{\Lambda}$. This fitting formula reproduces the simulation data of $E^{\text{tot}}_{\rm GW,i}$ within an error of $\approx 4\%$. In the limit to a binary black hole merger $(\tilde{\Lambda}\to 0)$, the fitting formula predicts $E^{\text{tot}}_{\text{GW,i}}\approx 0.034m_{0}$ for $\eta=0.250$ and $E^{\text{tot}}_{\text{GW,i}}\approx 0.033m_{0}$ for $\eta=0.244$, respectively. On the other hand, high-precision binary black hole merger simulations for non-spinning system suggests $E^{\text{tot}}_{\text{GW,i}}\approx 0.03m_{0}$ for $0.247\leq\eta\leq 0.250$ Blackman:2017dfb ; Boyle:2019kee . We conclude that the fitting formula Eq. (33) reproduces the BBH result with $\approx 10\%$ error.

V.2.2 Post-merger stage

We estimate angular momentum of the remnant, $J_{\text{rem}}$ at the peak time of the gravitational-wave amplitude in the retarded time (4) by performing a surface integral on the sphere of $r=r_{0}$;

[TABLE]

$K_{ij}$, $K$, ${\delta^{i}}_{j}$, and $dS_{l}$ are the extrinsic curvature, its trace part, the Kronecker delta, and an element of the surface integral, respectively. We typically integrate it on the sphere of $r_{0}=200m_{0}$ and $214m_{0}$ for the binary systems with ${\cal M}_{c}=1.1752M_{\odot}$ and $1.0882M_{\odot}$, respectively. Table 6 and Fig. 17 show the result. In the left panel of Fig. 17, we estimate the residual error in $J_{\text{rem}}$ for HB$118$–$155$. We again assume that the numerical result obeys the following form;

[TABLE]

where $J_{\text{rem}}^{\infty}(N_{\text{max}})$ is the angular momentum of the remnant in the continuum limit of the finite difference. We estimate three unknowns, $J_{\text{rem}}^{\infty}(N_{\text{max}})$, $\Delta J_{\text{rem}}(N_{\text{max}})$, and $p$ by fitting the numerical data with $N=90,102,\cdots,$ and $N_{\text{max}}$ with Eq. (35). By comparing $N_{\text{max}}=150$ and $182$ cases, we confirm that adding a result of the higher resolution simulation reduces the residual error (see the legend of Fig. 17 for $p$ and $\Delta J_{\text{rem}}(N_{\text{max}})$). We find that $\Delta J_{\text{rem}}(N_{\text{max}})$ is $\lesssim 1\%$ of the continuum limit, $J_{\text{rem}}^{\infty}(N_{\text{max}})$, for $N_{\text{max}}=182$. This is also the case for all the binary systems. Thus, we adopt $1\%$ as a systematics associated with the finite grid resolution in $J_{\text{rem}}$.

Because $J_{\text{rem}}$ could correlate with $\tilde{\Lambda}^{1/5}$, we propose a fitting formula of $J_{\text{rem}}/(m_{0}^{2}\eta)$:

[TABLE]

The right panel of Fig. 17 plots this relation and we confirm that it is accurate within $3\%$ error.

Figures 18 and 19 plot $E_{\rm GW,p}^{2,2}$ and $J_{\rm GW,p}^{2,2}$ emitted in the post-merger stage. It is worth noting that energy and angular momentum radiated by gravitational waves in $(l,m)=(2,1)$ and $(3,3)$ modes are $\lesssim 2.5\%$ of $E^{2,2}_{\rm GW,p}$ and $\lesssim 2.4\%$ of $J^{2,2}_{\rm GW,p}$, respectively, even for the highly asymmetric binary systems, e.g., 15H107-146 (see also the upper panel of Fig. 23). The left panels in these figures show that it is hard to achieve a perfect convergence and the scatter is rather large compared to $E_{\rm GW,i}^{2,2}$, although the scatter never exceeds $50\%$ in $E_{\rm GW,p}^{2,2}$ and $J_{\rm GW,p}^{2,2}$. This is also the case for all the binary systems. The right panels in Figs. 18 and 19 show $E_{\rm GW,p}^{2,2}/(m_{0}\eta)$ and $J_{\rm GW,p}^{2,2}/(m_{0}^{2}\eta)$ as a function of $\tilde{\Lambda}^{1/5}$. As discussed in Ref. Zappa:2017xba , the energy and angular momentum radiated in the post-merger stage peak around $\tilde{\Lambda}\approx 400$ because the binary systems with $\tilde{\Lambda}\lesssim 350$ collapse to a black hole within a few ms after the peak time. However, $\tilde{\Lambda}$ at the peak in $E_{\rm GW,p}^{2,2}$ and $J_{\rm GW,p}^{2,2}$ could decrease for general EOSs because as discussed in Ref. Kiuchi:2019lls the remnant would survive for more than $20$ ms after the peak time even for the binary systems with $\tilde{\Lambda}\lesssim 300$. For $\tilde{\Lambda}\gtrsim 400$, correlation between $E_{\rm GW,p}^{2,2}$ and the binary tidal deformability is not as tight as that in $E^{\text{tot}}_{\rm GW,i}/(m_{0}\eta)$–$\tilde{\Lambda}^{1/5}$. For $J_{\rm GW,p}^{2,2}$, the correlation with the binary tidal deformability is also not very tight.

Note that $E_{\rm GW,p}^{2,2}$ and $J_{\rm GW,p}^{2,2}$ could increase from the values listed in Table 6 because we artificially terminated the simulations at $10$–$15$ ms after the peak time. At that moment, the gravitational-wave amplitude is still comparable to that in the late inspiral stage except for the systems which collapse to a black hole within a few ms after the peak time.

We also should keep in mind that we might miss relevant physics such as effective turbulent viscosity generated by the magneto-hydrodynamical instabilities during the merger Kiuchi:2014hja ; Kiuchi:2015sga ; Kiuchi:2017pte and/or the neutrino cooling Sekiguchi:2011zd ; Foucart:2015gaa for modeling the post-merger signal. Reference Shibata:2017xht suggests that the post-merger signal could be significantly suppressed in the presence of efficient angular momentum transport by the viscous effect inside the remnant NS.

As already mentioned, the post-merger gravitational wave signal is dominated by the f–mode oscillation with $(l,m)=(2,2)$ of the remnant massive NS Bauswein:2011tp ; Hotokezaka:2013iia . Thus, it is natural to expect that a relation holds between the energy emission rate and angular momentum emission rate (10)–(11) with instantaneous gravitational-wave frequency (9);

[TABLE]

where $dE_{\rm GW}^{\text{post}}/dt=\sum_{l,m}dE_{\rm GW}^{l,m}/dt$ and $dJ_{\rm GW}^{\text{post}}/dt=\sum_{l,m}dJ_{\rm GW}^{l,m}/dt$ for $t\geq t_{\text{peak}}$ in Eqs. (10) and (11). To investigate to what extent this relation is satisfied, we generate Figs. 20–22. In these figures, the solid curve is the left hand side of Eq. (37) and the dashed curve is the right hand side of Eq. (37). We find that they agree with each other with a relative error $\lesssim 8\%$ for any time. Because the emissivity reduces quickly to zero at $t_{\text{ret}}-t_{\text{peak}}\approx 0.5$ ms as shown in Figs. 20–22, we estimate the error for $t_{\text{ret}}-t_{\text{peak}}\gtrsim 1$ ms. We also find that the time integrated values of Eq. (37) agree with each other with a relative error $\lesssim 1\%$. This is also the case for the relation of $E_{\rm GW,p}\approx\pi f_{2}J_{\rm GW,p}$.

We also confirm that a contribution from the one-arm spiral instability in the post-merger stage Paschalidis:2015mla ; Radice:2016gym is negligible because the energy flux for $(l,m)=(2,1)$ mode is $\lesssim 0.5\%$ of that for $(l,m)=(2,2)$ mode even for the symmetric binary systems as shown in the bottom panel of Fig. 23. Thus, we conclude that Eq. (37) is well satisfied and confirm that the main gravitational-wave emission mechanism during the post-merger stage is the f–mode oscillation of the remnant massive NS, i.e, $f_{\text{GW}}\approx f_{2}$ (see also Figs. 11–13). These findings encourage us to build a model for the post-merger gravitational-wave emission (see Ref. Shibata:2019ctb ).

In Table 7, we summarize to what extent $E^{\text{tot}}_{\rm GW,i}/(m_{0}\eta)$–$\tilde{\Lambda}^{1/5}$ and $J_{\text{rem}}/(m_{0}^{2}\eta)$–$\tilde{\Lambda}^{1/5}$ relations of Eqs. (33) and (34) hold.

VI Summary

We performed long-term simulations for new 26 systems of the non-spinning BNS mergers in numerical relativity. To derive high-precision gravitational waveforms in a large parameter space, we systematically vary the EOSs of NS, the chirp mass, and the mass ratio. To assess gravitational-wave phase error stemming from a finite grid resolution, we change the grid spacing by a factor of two for simulating each binary system.

First, we found that the residual gravitational-wave phase error at the peak time of gravitational-wave amplitude is $\lesssim 0.5$ rad irrespective of the binary mass and NS EOS. By comparing the results for the piecewise polytropic and SFHo (tabulated) EOS systems, we also found that the interpolation of the thermodynamic quantities during the simulations could generate the phase error of $\approx 0.2$–$0.3$ rad. However the gravitational-wave phase error for the SFHo (tabulated) EOS system still remains within the sub-radian accuracy level.

Second, we validated our SACRA inspiral gravitational waveform template Kawaguchi:2018gvj by comparing with the high-precision gravitational waveforms derived in this paper. We found that for a variety of BNS the error in our inspiral waveform model is less than $0.1$ rad in the gravitational-wave phase and less than $20\%$ in the amplitude up to $f_{\text{GW}}=1000$ Hz. This template can be used for a new gravitational wave data analysis for extracting tidal deformability from GW170817 Narikawa:2019xng and for future event of BNS merger.

Third, we assessed the universal relations between the gravitational-wave related quantities and the binary tidal deformability/NS radius proposed in the literature Rezzolla:2016nxn ; Read:2013zra ; Zappa:2017xba ; Bauswein:2011tp ; Bauswein:2012ya ; Bernuzzi:2014owa ; Bernuzzi:2015rla . We found that the gravitational-wave frequency at the peak time $f_{\text{peak}}$, the gravitational-wave amplitude at the peak time $h_{\text{peak}}$, and the peak frequency $f_{2}$ associated with the f–mode oscillation of the remnant massive NS in the spectrum amplitude of post-merger gravitational waves depend strongly on the symmetric mass ratio and/or the grid resolution. This clearly illustrates that the universal relations proposed in the literature Rezzolla:2016nxn ; Read:2013zra ; Zappa:2017xba ; Bauswein:2011tp ; Bauswein:2012ya ; Bernuzzi:2014owa ; Bernuzzi:2015rla are not as universal as proposed.

We proposed improved fitting formulae (29) for $m_{0}f_{\text{peak}}$–$\tilde{\Lambda}^{1/5}$, (30) for $Dh_{\text{peak}}/m_{0}$–$\tilde{\Lambda}^{1/5}$, and (32) and for $m_{0}f_{2}$–$\tilde{\Lambda}^{1/5}$. However these fitting formulae may still suffer from systematics such as NS spin, NS magnetic fields, and the neutrino radiation, which are not taken into account in our simulations. In addition, the EOS of NS, in particular, for a high-density part of the NS, is still uncertain, and hence, the systematics due to this uncertainty should be kept in mind. We also note that we assessed the errors of these formulae only with our simulation data. A close comparison among the results of the independent BNS simulations with the existing numerical relativity codes is necessary to better understand the systematic error in these formulae. This should be done as a future project. We also found that $f_{1}$ frequency in the spectrum amplitude could be extracted only for the nearly symmetric binary systems. Unless we can determine the symmetric mass ratio accurately, using the universal relation for $f_{1}$ could lead to a misleading result in the gravitational-wave data analysis.

Finally, we assessed the energy, $E_{\text{GW}}$, and angular momentum, $J_{\text{GW}}$, carried by gravitational waves in the inspiral and post-merger stages. As proposed in Ref. Zappa:2017xba , the correlation between $E_{\rm GW,i}^{\text{tot}}$ and the binary tidal deformability is tight and it does not depend significantly on the symmetric mass ratio. We found that the relation $E_{\text{GW}}\approx\pi f_{2}J_{\text{GW}}$ is well satisfied in the post-merger gravitational wave signal irrespective of the binary mass and NS EOS because the signal from the remnant NSs is approximately monochromatically emitted by the f–mode oscillation. The angular momentum of the remnant massive NS, $J_{\text{rem}}$, correlates with the binary tidal deformability. This quantity is relevant to build a model of post-merger evolution of merger remnants Shibata:2019ctb .

Acknowledgements.

Numerical computation was performed on K computer at AICS (project numbers hp160211, hp170230, hp170313, hp180179, hp190160), on Cray XC50 at cfca of National Astronomical Observatory of Japan, Oakforest-PACS at Information Technology Center of the University of Tokyo, and on Cray XC40 at Yukawa Institute for Theoretical Physics, Kyoto University. This work was supported by Grant-in-Aid for Scientific Research (16H02183, 16H06342, 16H06341, 16K17706, 17H01131, 17H06361, 17H06363, 18H01213, 18H04595, 18H05236, 18K03642, 19H14720) of JSPS and by a post-K computer project (Priority issue No. 9) of Japanese MEXT. Our waveform data is publicly available on the web page.