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

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.

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 [34]. 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 Λ 800 with the binary total mass of 2.73 +0.04 −0.01 M [3, [35][36][37] 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 [38][39][40][41][42][43] and their template was used in gravitational wave data analysis by LIGO Scientific and Virgo Collaborations to infer the tidal deformability from GW170817 [35]. However, the residual phase error caused mainly by the finite grid resolution in their simulations is ≈ 0.5-2.3 rad [43]. The phase error of O(1) rad could be an obstacle to construct a high-quality inspiral gravitational waveform template (see also Refs. [44,45]).
In Ref. [46], we tackled this problem by using our numerical relativity code SACRA-MPI and performed longterm simulations with the highest grid resolution to date (see also Refs. [47][48][49][50] for our effort in the early stage of this project). In our numerical results, the gravitationalwave 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. [51] presented a waveform template, the SACRA inspiral gravitational waveform template, of BNS mergers. Specifically, we multiply the tidalpart phase of the 2.5 Post-Newtonian (PN) order derived in Ref. [52] by a correction term composed of the PN parameter and the binary tidal deformability. Then, we validated it by confirming that it reproduces the highprecision gravitational waveforms derived in Ref. [46]. We also validated a correction term in the tidal-part amplitude of the 1 PN order derived in Refs. [52,53].
In Refs. [46,51], 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 [46,51]. 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 [54], 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 [55,56]. 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 [57][58][59][60][61][62][63]. 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 [57][58][59][60][61][62][63] 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 [64].
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 A. Method and grid setup
We use our numerical relativity code, SACRA-MPI [46,65], to simulate a long-term inspiral stage of BNS up to early post-merger. SACRA-MPI implements the Baumgarte-Shapiro-Shibata-Nakamura-puncture formulation [66][67][68][69], locally incorporating a Z4c-type constraint propagation prescription [70], 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 [71].
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 [72] together with the 3rdorder accurate cell reconstruction [73].
We also implement the Berger-Oliger type adaptive mesh refinement (AMR) algorithm [74] 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 ∆x l = L/(2 l N ) with l = 0, 1, · · · 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 ∆x 9 are summarized in Table I. B. Binary system parameters and gravitational wave extraction Table I shows the list of the binary systems as well as the grid setup for the simulations.

Equation of state
Following the previous papers [46,51], we employ a parameterized piecewise polytropic EOS to describe the NS matter [75]. Specifically, we assume that the pressure and specific internal energy consist of two segments with respect to the rest-mass density: with i = 0, 1, ρ 0 = 0 g cm −3 , and ρ 2 = ∞. ρ 1 is the rest-mass density which divides the pressure and specific internal energy into the two segments. Given the adiabatic indices Γ 0 , Γ 1 and one of the polytropic constants κ 0 , the other polytropic constant κ 1 is calculated from the continuity of the pressure at ρ = ρ 1 by κ 0 ρ Γ0 1 = κ 1 ρ Γ1 1 . ∆ 1 is also calculated from the continuity of the specific internal energy at ρ = ρ 1 by [75], we fix Γ 0 = 1.3562395, Γ 1 = 3, and κ 0 = 3.594 × 10 13 in cgs units. By varying the remaining parameter ρ 1 for a wide range as shown in Table II, we can derive plausible NS models with a variety of the radii and tidal deformability (see Table III).
In addition to the piecewise polytropic EOS, we employ one tabulated EOS, SFHo [76]. 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 β-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 cold (ρ) and cold (ρ) (see also Table III 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: where is the specific internal energy and we assumed that the thermal part is described by the Γ-law EOS with the index Γ th . Following Refs. [46,51], we fix Γ th = 1.8. We note that gravitational waveforms for the postmerger stage depend on the value of Γ th [77], although inspiraling waveforms do not. Since the major purpose of the present paper is to derive the accurate inspiraling waveforms, the choice of Γ 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. [77]). Thus, we have to keep in mind that the systematics exist due to the uncertainty of this value [78].

Binary systems
In this paper, we consider 6 irrotational binary systems assuming that NSs have no spin before merger. We fix a chirp mass, M c , and symmetric mass ratio, η, to be (M c , η) = (1.  Table I). For the SFHo (tabulated) EOS, we only consider the equal-mass binary system with m 1 = 1.35M and m 2 = 1.35M . Table I also shows the binary tidal deformability for all the binary systems [79,80]: where Λ 1 (Λ 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 ≈ 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. [81] 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 , and m 2 = 1.46M . We set the initial orbital angular velocity to be m 0 Ω 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 [82,83].
Then, we reduce orbital eccentricity by using the prescription in Ref. [84]. With this method, we confirm that the initial orbital eccentricity is reduced typically to ≈ 10 −3 which is low enough to generate a high-precision inspiral waveform (see also Appendix in Refs. [46,51]).

C. Gravitational wave extraction
We calculate a complex Weyl scalar Ψ 4 from simulation data to derive gravitational waveforms [65]. Given an extraction radius r 0 , the Weyl scalar Ψ 4 is decomposed into (l, m) modes with the spin-weighted spherical harmonics by where t ret is a retarded time defined by with D = A/4π. A is a proper area of the extraction sphere. We apply the Nakano's method [85] to extrapolate Ψ l,m 4 to infinity by where C(r 0 ) is a function of r 0 . Following Ref. [46], we choose D ≈ 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 Ψ l,m,∞ 4 twice in time: For the time integration, we employ the fixed frequency method [86] by (t) and f cut is set to be 0.8mΩ 0 /(2π).
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 M c = 1.1752M and r 0 = 262 m 0 , 213 m 0 , and 156 m 0 for M c = 1.0882M (see Table I).
In general, gravitational waves for each (l, m) mode are decomposed into the amplitude and phase as h l,m,∞ (t ret ) = A l,m,∞ (t ret )e −iΦ l,m (tret) , (2.8) and instantaneous gravitational-wave frequency is defined by dΦ l,m /dt ret . In Sec. III, we explore the accuracy of the gravitational-wave phase of the (l, m) = (2, 2) mode, and simply refer to Φ 2,2 as the gravitational-wave phase. With Eq. (2.8), the instantaneous frequency of the (l, m) = (2, 2) mode is calculated by where the asterisk symbol denotes the complex conjugate of h 2,2,∞ . We also calculate the energy and angular momentum flux due to gravitational-wave emission by [87] dE l,m where t 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. [46,51]. Our waveform data are publicly available on the website: On the website, the waveform data are tabulated according to the system name, dimensionless initial orbital angular velocity, and grid resolution.  Table I). A user can download the data for Ψ 2,2 4 (t ret , r 0 ) extracted at several values of r 0 and h 2,2,∞ +,× (t ret ) from the link on the system name. 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 , m 2 = 1.40M and N = 182. It shows that the systems with the larger values ofΛ merge earlier than those with the smaller values ofΛ 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 η merge earlier than those with the smaller values of η because the emissivity of gravitational waves decreases as the symmetric mass ratio decreases [88].
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 where Φ 2,2 (t ret ; 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 where N 1 and N 2 denote the employed grid numbers.
The red shaded region shows δΦ error (t ret ; 15H, 150, 182) and the blue shaded region shows δΦ error (t ret ; 15H, 90, 182), respectively, for 15H112-140. The overlapped region has purple color. The vertical dashed line denotes the peak time, t 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 ms t ret 59 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 δΦ error (t ret ; 15H, 150, 182) than in δΦ error (t ret ; 15H, 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  is larger than δΦ error (t ret ; 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. 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, where Φ 2,2,∞ peak (N max ) and p denote the gravitational-wave phase at the peak time in the continuum limit of the finite difference (N → ∞) and an order of the convergence, respectively. ∆Φ 2,2 peak (N max ) should be recognized as the residual phase error for the simulation with N = N max . N max denotes a reference value of N to estimate unknown quantities Φ 2,2,∞ peak (N max ), ∆Φ 2,2 peak (N max ), and p. For example, with N max = 182, these unknowns are obtained by fitting the simulation results of N = 150, 130, 110, 102, and 90 with Eq. (3.3) 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, δΦ error (t peak ; B, N max , N ), as a function of 1/N p with a reference grid number N max and N = 90, 102, · · · , N max . Assuming Eq. (3.3), the phase error at the peak time in a binary system is given as The values of ∆Φ 2,2 peak (N 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 max . Table IV 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 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 max = 130 happens to be smaller (higher) than that for N max = 150. Nonetheless, the residual phase error (the order of convergence) for N max = 182 is smaller (higher) than that for N max = 150 except for 125H125-146. Thus, we adopt the values for N max = 182 as the residual phase error in our waveforms and it is in the range of ≈ 0.1-0.5 rad.

IV. INSPIRAL GRAVITATIONAL WAVEFORM MODELING
A. SACRA inspiral gravitational waveform template In the previous paper [51], we developed a frequencydomain 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 [46,51] 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 bỹ where t i and t f are the initial and final time of the waveform data, respectively. Then, we decomposeh + (f ) in Eq. (4.1) into the frequency-domain amplitude, A (f ), and phase, Ψ (f ), (with an ambiguity in the origin of the phase) byh We only use h 2,2,∞ + for modeling the inspiral gravitational waveforms because the difference between h 2,2,∞ + and h 2,2,∞ × is approximately only the phase difference of π/2. We define the corrections due to the NS tidal deformation to the gravitational-wave amplitude and phase by and respectively. Here, A BBH (f ) and Ψ BBH (f ) 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. [51] for details).
Our numerical-relativity waveforms only contain the waveforms for the frequency higher than ≈ 400 Hz. Thus, we employ the effective-one-body waveforms of Refs. [89][90][91][92] (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 timedomain by the procedure described in Refs. [50,51] and we set the matching region to be from t ret ≈ 7.38 ms to 14.78 ms. After the hybridization, the waveforms are transformed into the frequency domain employing Eq. (4.1), and the tidal-part amplitude and phase are extracted by Eqs. (4.3) and (4.4).
For modeling the tidal-part phase and amplitude, we employ the following functional forms motivated by the 2.5 PN order formula [52]: for the phase correction and for the amplitude correction where D eff is the effective distance to the binary [50] and x ≡ (π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 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 [51]. 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. [51], 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. (4.5) and (4.6) 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 ∼ 0.1 rad and 15%, respectively, for the equal-mass and unequal-mass cases with M chirp = 1.1752 M and the equal-mass cases with M chirp = 1.0882 M , covering the parameter space of 0.244 ≤ η ≤ 0.250 and 300 Λ 1800.

B. 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 [51], there still remain some important cases which were not examined in the previous study [51]. 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 ≈ 2.7 M were checked in the previous study [51], 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 [51] employing the SEOB-NRv2T 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 tidalpart phase and amplitude of the hybrid waveforms computed based on Eqs. (4.3) and (4.4) using the SEOBNRv2 waveforms with no-tides as the point-particle parts. Figures 5 and 6 show the difference of the tidalpart phase and amplitude between our inspiral waveform model (4.5) and (4.6) and the hybrid waveforms for the models with M c = 1.1752 M and M c = 1.0882 M . Here, the phase difference between the tidal-part phase of hybrid waveforms, Ψ tidal Hybrid , and that of our inspiral waveform model, Ψ tidal model , is computed by where t 0 and φ 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 |∆Ψ(f )| 2 df integrated in the range of f = 10-1000 Hz. For the comparison of the tidal-part amplitude, relative difference of the amplitude, is computed, where A tidal Hybrid and A model = A tidal model +A BBH are the tidal-part amplitude of hybrid waveforms and the amplitude of the model waveforms including the pointparticle part, respectively. Again, we employ the amplitude of the SEOBNRv2 waveforms with no-tides for A BBH .
The systems of mass 1.25M -1.46M , 1, 18M -1.55M , and 1.17M -1.56M are within the parameter space which we studied in the previous study [51], 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 [51]. 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 M c = 1.0882 M . The difference in the tidal-part phase is larger than the cases with M c = 1.1752 M . 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 [51]. Nevertheless, the phase error is always smaller than ≈ 0.1 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 M c = 1.1752 M .
To quantify the deviation of our inspiral waveform model from the new sets of hybrid waveforms, we calculate the mismatch between those waveforms,F , defined byF   where (·|·) and || · || are defined by  Here, h 1 and h 2 denote the hybrid waveforms and our inspiral waveform models, respectively. The inspiral waveform model employs Eqs. (4.5) and (4.6) as the tidal part and the SEOBNRv2 waveforms with no-tides as the point-particle baseline. S 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 [93] for it. We summarize the values of mismatch between our inspiral waveform model and hybrid waveforms in Table V. For all the cases, the value of mismatch is smaller than ≈ 2 × 10 −5 . According to our previous results [51], 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.

A. frequency and amplitude
Instantaneous gravitational-wave frequency defined by Eq. (2.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 [57,58,60,61]. 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 [47,58,63,94]. 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.

Peak frequency and binary tidal deformability relation
Reference [57] reported that the instantaneous gravitational-wave frequency (of l = |m| = 2 mode) at the peak time (t peak ), f peak , has a tight correlation with the binary tidal deformabilityΛ (see also Refs. [58,60,61] for the relation with the tidal coupling constant: In Ref. [58], they referred to it as f max ). Figure 7 plots the dependence of f peak on the grid resolution where f peak,ave is the average of f peak over the results with different grid resolutions. f 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 peak to be 2% and tabulate the values of f peak in Table VI. The right panel of Fig. 7 plots m 0 f peak as a function ofΛ 1/5 . The error bar shows the systematics associated with the finite grid resolution in f peak . We also plot the universal relations reported in Refs. [57] (black dashed line) and [58] (black dotted line). We find that the universal relation in Ref. [58] holds only for the symmetric binary systems with M c = 1.1752M and M c = 1.0882M (see also Table VI). Given an EOS and a chirp mass, f 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 GW , df GW /dt decreases as the symmetric mass ratio decreases because the gravitational-wave luminosity is proportional to η 2 [88]. 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 dumbbelllike density structure with double dense cores formed after the contact disappears as discussed in Ref. [46] and the dumbbell-like density structure becomes less prominent in the asymmetric binary systems. Due to these effects, f peak becomes lower as the symmetric mass ratio decreases.
In a short summary, the m 0 f peak -Λ 1/5 relation depends strongly on the symmetric mass ratio and the universal relations reported in Refs. [57] and [58] suffer from this systematics (see also Ref. [46]). This finding is consistent with a discussion in Ref. [58]. They mentioned that the mass asymmetry could break the universality in the m 0 f peak -Λ 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 [3]. The scatter from the proposed universal relation in Ref. [58] is as large as ≈ 18-19% at the maximum for 0.244 ≤ η ≤ 0.250.
We propose an improved fitting formula: With η = 0.2500, a 0 (η) and a 1 (η) approximately reduce to be a 0 and a 1 [95] reported in Ref. [58]. Figure 8 plots the improved relation with the simulation data and we confirm that the relative error between the data and the fitting formula (5.1) 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 peak compared to the non-spinning case. NS magnetic fields also could produce systematics in Eq. (5.1) 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 1 ms [96, 97] and the magnetic pressure could reach near the equipartition of the pressure locally, affecting the value of f peak . These points should be explored in future work.

Peak amplitude and binary tidal deformability relation
References [46,57] reported that the gravitationalwave amplitude at the peak time, h peak , correlates with f peak , i.e., withΛ 1/5 . Because we do not find perfectly convergent result for h peak with respect to the grid resolution, first, we assess deviation of h peak relative to the averaged value of h 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 and m 2 = 1.46M .
It is found that fluctuation around the averaged value is ≈ 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 peak and summarize the values of h peak in Table VI.
The right panel of Fig. 9 plots Dh peak /m 0 as a function ofΛ 1/5 . The error bar shows the systematics associated with the finite grid resolution in h peak . This figure shows that the relation depends strongly on the symmetric mass ratio. That is, the relation proposed in Refs. [46,57] is not in general satisfied.
We propose a fitting formula for Dh peak /m 0 : Figure 10 plots the improved relation with the simulation data. We find that the relative error between the data and the fitting formula (5.2) 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.

f1, f2 and binary tidal deformability relation
Reference [58] 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 withh + (f ) andh × (f ) in Eq. (4.1). In Figs. 11-13, the vertical dashed lines indicate the so-called f 1 frequency for the fitting formula in Ref. [58]. This peak is a sideband 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. [98], 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. [58], 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   7. (Left) A deviation of instantaneous gravitational-wave frequency at the peak time f peak relative to f peak,ave as a function of 1/N for the binary systems with m1 = 1.17M and m2 = 1.56M . f peak,ave is an average of f peak over the results with different grid resolutions. (Right) m0f peak as a function ofΛ 1/5 . Meaning of the color and symbols is the same as that in Fig. 1. The error bar of ±2% comes from the systematics associated with the finite grid resolution in f peak . The proposed universal relations in Refs. [57,58] are shown. 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. [99] 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 ∼ 2.5M as discussed in Refs. [58,100,101]. We also analyze the spectrogram of post-merger grav-itational waves and confirm that there is no prominent peak around f 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. [98] 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. [102] 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. [58], the peak frequency, f 2 , in the spectrum amplitude [103] 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. [47,63,77,94]). 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 and m 2 = 1.40M . We measure f 2 in the spectrum amplitude as a prominent peak for f ≥ 2 kHz. The fluctuation is within ≈ 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 VI). The right panel of Fig. 14 shows f 2 as a function ofΛ 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. [58]. It is found that with this fitting formula, the scatter is ≈ 14% at the maximum. Thus, we propose an improved fitting   9. (Left) A deviation of the gravitational-wave amplitude at the peak time, h peak , relative to h peak,ave as a function of 1/N for the binary systems with m1 = 1.07M and m2 = 1.46M . h peak,ave is an average of h peak over the results with different grid resolutions. (Right) Dh peak /m0 as a function ofΛ 1/5 . Meaning of the color and symbols is the same as Fig. 1. The error bar of ±2% comes from the uncertainty associated with the finite grid resolution in h peak . 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,Λ 1/5 will be constrained with the error of ≈ ±0.1.

f2 and NS radius with 1.6M relation
References [62,63] reported that f 2 frequency has a tight correlation with the NS radius of 1.6M (see Eq. (3) in Ref. [62]). In Ref. [94], 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. [62]. We revisit this assessment because the initial orbital eccentricity reduction was not implemented in Ref. [94]. In addition, the grid resolution in Ref. [94] 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. [62] holds only for symmetric binary systems of m 0 = 2.7M , we first assess this relation by employing binary systems of (m 1 , m 2 ) = (1.35M , 1.35M ) and found that the error is ≈ 6% [104].
Second, we assess the relation by employing binary systems of (m 1 , m 2 ) = (1.25M , 1.46M ), (1.21M , 1.51M ), (1.18M , 1.55M ), (1.17M , 1.56M ), and (1.16M , 1.58M ). We found that the scatter from their fitting formula is ≈ 10%. Therefore, the scatter larger than that reported in Ref. [62] stems from the mass asymmetry of the binary. Our numerical results suggest that the fitting formula in Ref. [62] could infer the radius of the 1.6M 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 NS with the accuracy of ≈ ±1 km if the value of f 2 is determined precisely, In Table VII, we summarize to what extent the socalled universal relations hold.

B. Energy and angular momentum
Using Eqs. (2.10)-(2.13), we calculate the energy and angular momentum carried by gravitational waves. We define E tot GW,i and E GW,p (J GW,p ) as the energy (angular momentum) emitted in the inspiral stage and in the post-merger stage, respectively. The subscripts i and 1752M . The number attached in the right-hand side vertical axis is the symmetric mass ratio η. We also show f1 frequency proposed in Ref. [58] with vertical dashed lines. For completeness, we also show the systems reported in Refs. [46,51].
p in these quantities denote the inspiral and the postmerger stage, respectively. The peak time introduced in Sec. III A 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 VI. 1. inspiral stage Table VI and Fig. 16 show the energy, E 2,2 GW,i , 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 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 2,2 GW,i . Note that the other modes such as (l, m) = (2, 1) and (3,3) are 0.1% and 0.5%, respectively, of E 2.2 GW,i .
The right panel of Fig. 16 plots E tot GW,i /(m 0 η) as a function ofΛ 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 ADM,0 in Table VI, by E tot GW,i ≈ 2E 2,2 GW,i + m 0 − M ADM,0 . M ADM,0 is the Arnowitt-Deser-Misner mass of the initial condition of the simulations. As proposed in Ref. [59], this quantity correlates with the tidal coupling constant. We explicitly derive a fitting formula with the binary tidal deformability as It is reasonable that E tot GW,i decreases asΛ increases be-cause the binary systems with larger values ofΛ merge earlier than those with smaller values ofΛ. This fitting formula reproduces the simulation data of E tot GW,i within an error of ≈ 4%. In the limit to a binary black hole merger (Λ → 0), the fitting formula predicts E tot GW,i ≈ 0.034m 0 for η = 0.250 and E tot GW,i ≈ 0.033m 0 for η = 0.244, respectively. On the other hand, highprecision binary black hole merger simulations for nonspinning system suggests E tot GW,i ≈ 0.03m 0 for 0.247 ≤ η ≤ 0.250 [106,107]. We conclude that the fitting formula Eq. (5.5) reproduces the BBH result with ≈ 10% error.

Post-merger stage
We estimate angular momentum of the remnant, J rem at the peak time of the gravitational-wave amplitude in the retarded time (2.4) by performing a surface integral on the sphere of r = r 0 ; K ij , K, δ 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 M c = 1.1752M and 1.0882M , respectively. Table VI and Fig. 17 show the result. In the left panel of Fig. 17, we estimate the residual error in J rem for HB118-155. We again assume that the numerical result obeys   the following form; where J ∞ rem (N max ) is the angular momentum of the remnant in the continuum limit of the finite difference. We estimate three unknowns, J ∞ rem (N max ), ∆J rem (N max ), and p by fitting the numerical data with N = 90, 102, · · · , and N max with Eq. (5.7). By comparing N 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 ∆J rem (N max )). We find that ∆J rem (N max ) is 1% of the continuum limit, J ∞ rem (N max ), for N 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 rem .
Because J rem could correlate withΛ 1/5 , we propose a fitting formula of J rem /(m 2 0 η): 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 2,2 GW,p and J 2,2 GW,p 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 2.5% of E 2,2 GW,p and 2.4% of J 2,2 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 2,2 GW,i , although the scatter never exceeds 50% in E 2,2 GW,p and J 2,2 GW,p . This is also the case for all the binary systems. The right panels in Figs. 18 and 19 show E 2,2 GW,p /(m 0 η) and J 2,2 GW,p /(m 2 0 η) as a function ofΛ 1/5 . As discussed in Ref. [59], the energy and angular momentum radiated in the post-merger stage peak aroundΛ ≈ 400 because the binary systems withΛ 350 collapse to a black hole within a few ms after the peak time. However,Λ at the peak in E 2,2 GW,p and J 2,2 GW,p could decrease for general EOSs because as discussed in Ref.
[105] the remnant would survive for more than 20 ms after the peak time even for the binary systems withΛ 300. ForΛ 400, correlation between E 2,2 GW,p and the binary tidal deformability is not as tight as that in E tot GW,i /(m 0 η)-Λ 1/5 . For J 2,2 GW,p , the correlation with the binary tidal deformability is also not very tight.
Note that E 2,2 GW,p and J 2,2 GW,p could increase from the values listed in Table VI 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 [46,96,97] and/or the neutrino cooling [100,108] for modeling the post-merger signal. Reference [109] 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 [63,94]. Thus, it is natural to expect that a relation holds between the energy emission rate and angular momentum emission rate (2.10)-(2.11) with instantaneous gravitational-wave frequency (2.9);  20-22. In these figures, the solid curve is the left hand side of Eq. (5.9) and the dashed curve is the right hand side of Eq. (5.9). We find that they agree with each other with a relative error 8% for any time. Because the emissivity reduces quickly to zero at t ret − t peak ≈ 0.5 ms as shown in Figs. 20-22, we estimate the error for t ret − t peak 1 ms. We also find that the time integrated values of Eq. (5.9) agree with each other with a relative error 1%. This is also the case for the relation of E GW,p ≈ πf 2 J GW,p .
We also confirm that a contribution from the one-arm spiral instability in the post-merger stage [110,111] is negligible because the energy flux for (l, m) = (2, 1) mode is 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. (5.9) is well satisfied and confirm that the main gravitationalwave emission mechanism during the post-merger stage is the f-mode oscillation of the remnant massive NS, i.e, f GW ≈ f 2 (see also Figs. 11-13). These findings encourage us to build a model for the post-merger gravitationalwave emission (see Ref. [112]).

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 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 ≈ 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 [51] 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 GW = 1000 Hz. This template can be used for a new gravitational wave data analysis for extracting tidal deformability from GW170817 [113] 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 [57][58][59][60][61][62][63]. We found that the gravitational-wave frequency at the peak time f peak , the gravitational-wave amplitude at the peak time h 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 [57][58][59][60][61][62][63] are not as universal as proposed.
We proposed improved fitting formulae (5.1) for m 0 f peak -Λ 1/5 , (5.2) for Dh peak /m 0 -Λ 1/5 , and (5.4) and for m 0 f 2 -Λ 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   (5.5). In the right panel, the error bar of ±2% comes from the systematics associated with the finite grid resolution in E 2,2 GW,i .    18. (Left) A deviation of E 2,2 GW,p relative to E 2,2 GW,p,ave as a function of 1/N for binary systems with m1 = 1.25M and m2 = 1.46M . E 2,2 GW,p,ave is an average of E 2,2 GW,p over the results with different grid resolutions. (Right) E 2,2 GW,p /(m0η)-Λ 1/5 relation. In the right panel, the error bar of ±50% comes from the systematics associated with the finite grid resolution in E 2,2 GW,p . 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 GW , and angular momentum, J GW , carried by gravitational waves in the inspiral and post-merger stages. As proposed in Ref. [59], the correlation between E tot GW,i and the binary tidal deformability is tight and it does not depend significantly on the symmetric mass ratio. We found that the relation E GW ≈ πf 2 J 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 fmode oscillation. The angular momentum of the remnant massive NS, J rem , correlates with the binary tidal deformability. This quantity is relevant to build a model of post-merger evolution of merger remnants [112].   VII. 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 f1, the error is unable to be estimated because of the absence of f1 peak in the asymmetric binary systems. Therefore, we conclude there is no universal relation between f1 and Λ. For f2-R1.6 relation, we do not propose an improved relation and sym. (asym.) in the parenthesis means the symmetric (asymmetric) binary. For E 2,2 GW,p and J 2,2 GW,p , we do not propose an improved relation because uncertainties of the life time of the merger remnant NSs are large.