Reanalysis of the binary neutron star mergers GW170817 and GW190425 using numerical-relativity calibrated waveform models

Tatsuya Narikawa1,2,∗ Nami Uchikata3,2,† Kyohei Kawaguchi2,4,‡ Kenta Kiuchi4,5,§ Koutarou Kyutoku1,6,7,5,¶ Masaru Shibata4,5,∗∗ and Hideyuki Tagoshi2†† Department of Physics, Kyoto University, Kyoto 606-8502, Japan Institute for Cosmic Ray Research, The University of Tokyo, Chiba 277-8582, Japan Graduate School of Science and Technology, Niigata University, Niigata 950-2181, Japan Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am Mühlenberg 1, Potsdam-Golm 14476, Germany Center for Gravitational Physics, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan Theory Center, Institute of Particle and Nuclear Studies, KEK, Tsukuba 305-0801, Japan Interdisciplinary Theoretical and Mathematical Science (iTHEMS) Research Group, RIKEN, Wako, Saitama 351-0198, Japan


I. INTRODUCTION
Binary-neutron-star (BNS) mergers are valuable laboratories for nuclear astrophysics. Matter effects influence the orbital evolution and gravitational radiation through the tidal interaction between the neutron stars (NSs) in the late inspiral phase. Additionally, the presence of material gives rise to electromagnetic emission approximately coincident with gravitational radiation. Because these signatures depend on the properties of nuclear matter, their observations allow us to study various nuclear properties such as the equation of state (EOS) for NS matter.
GW170817 [1] and associated electromagnetic counterparts are used to derive various constraints on NS properties and the underlying EOS. The existence of a blue component in the kilonova/macronova AT 2017gfo [2] might suggest that the merger remnant did not collapse promptly to a black hole. Thus, the maximum mass of the NS should not be as small as ∼ 2M [3] and also the radii of high-mass NS may not be very small, e.g., the radius of the maximum-mass configuration is likely to be larger than 9.60 km [4] (but see also Ref. [5]). At the same time, the short gamma-ray burst GRB 170817A [6] and the absence of magnetar-powered emission in AT 2017gfo suggest that the remnant NS collapsed early in the postmerger phase (but see also Refs. [7][8][9][10]). Accordingly, a maximum mass of 2.3M is also unlikely [3,[11][12][13][14].
Tidal deformability extracted via cross-correlating gravitational-wave (GW) data of GW170817 with theoretical waveforms gives us more concrete information about the NS than electromagnetic counterparts. The LIGO-Virgo collaborations (LVC) initially put an upper limit on the binary tidal deformabilityΛ of the binary asΛ 800 with the prior on the dimensionless NS spin being chosen to be |χ| ≤ 0.05 [1]. This limit is later corrected to beΛ 900 in Ref. [15], where the result of arXiv:1910.08971v3 [gr-qc] 18 Jun 2020 updated analysis is also reported as, e.g.,Λ = 300 +420 −230 for a particular set of assumptions. The constraint can be further improved by assuming the EOS to be common for both NS [16,17] (but see also Ref. [18]) as is also done in an independent analysis [19,20]. These constraints are used to investigate the NS EOS [21][22][23] as well as those for quark and hybrid stars [24][25][26]. While it has been claimed based on a limited number of numericalrelativity (NR) simulations thatΛ 400 is necessary to account for the ejecta mass of ≈ 0.05M required to explain AT 2017gfo [27] (see Ref. [28] for an updated analyses stating the bound on the tidal deformability in 323 ≤Λ ≤ 776 and Ref. [29], stating the similar bound to be in 302 ≤Λ ≤ 860), later systematic investigations reveal that this argument is premature [5].
Recently, the discovery of the second BNS merger event, GW190425, was reported [30]. This binary system is massive and it is intrinsically difficult to measure the tidal effect. LVC have reported that GW190425 con-strainΛ to below 650 for the low-spin prior (|χ| ≤ 0.05), after reweighing the posterior to derive approximately the result corresponding to a flat prior inΛ [30]. While GW190425 does not carry novel information on the NS properties, multi-messenger constraints on the NS EOS have been studied [31,32].
An accurate theoretical waveform template is crucial to extract accurately the tidal deformability of NSs from the observed GW data. For the early stage of the inspiral, the waveforms including the linear-order tidal effects derived by post-Newtonian (PN) calculation are useful [33,34]. However, the PN expansion becomes invalid as the orbit becomes relativistic, and thus, the error of the waveform becomes large in the late stage [35][36][37][38]. Such errors would cause the systematic bias in the parameter estimation, and it would be in particular problematic for estimating the tidal deformability because the tidal effects on the waveform become most significant just before the merger [39,40]. The effective-one-body (EOB) formalism can solve this problem by incorporating the higher-order PN correction by re-summation techniques and calibrating them to NR waveforms [39,[41][42][43][44][45]. However, such calibration is performed only focusing on binary black holes (BBHs), and the calibration of the tidal correction employing NR simulation data of BNSs is also required.
Dietrich et al. have derived a gravitational waveform model, NRTidal, for BNSs based on high-precision NR simulations [46]. Improved reanalyses of GW170817 with more sophisticated waveform models calibrated by NR simulation of BNS merger have been performed employing such a model [15]. Indeed, it is pointed out that the value of the tidal deformability tends to be overestimated if the PN models are employed for the parameter estimation [1]. Recently, its upgrade, the NRTidalv2 model, which is calibrated by more precise NR waveforms, has been derived [47]. Kawaguchi et al. have also developed a model (hereafter the KyotoTidal model) for frequencydomain gravitational waveforms of inspiraling BNSs [48].
In particular, this model is derived independently from the NRTidal model employing different NR waveforms. Since the NRTidal model is so far the only NR calibrated waveform model that is used for parameter estimation of GWs from BNS mergers, the analysis comparing these three NR-calibrated waveform models would help us to understand the systematic biases in resulting constraints on tidal deformability.
In this paper, we reanalyze the data around GW170817 and GW190425 against a NR calibrated waveform model, the TF2+ KyotoTidal model and present constraints on the binary tidal deformability. We also reanalyze the events with other waveform models: two PN (TF2 PNTidal and TF2+ PNTidal), TF2+ NRTidal, and TF2+ NRTidalv2 models. Here, TF2 is the abbreviation of TaylorF2, which is the PN waveform model for a pointparticle part [54,55] and TF2+ [48] is a phenomenologically extended model of TF2.
The remainder of this paper is organized as follows. In Sec. II, we explain the methods for parameter estimation including waveform models used to reanalyze GW170817 and GW190425. In Sec. III, we present results of our analysis of GW170817 and a comparison of our analysis with the LVC analysis. In sec. IV, we present results of GW190425 and discuss constraints on NS EOS by combining information obtained from GW170817 and GW190425. Section V is devoted to a summary. In Appendix, we present an in-depth study of our results for GW170817 by separate analysis for the LIGO twin detectors to interpret the origin of the complex structure at the high-Λ region for the posterior probability density function (PDF) ofΛ (see also Ref. [49]). Unless otherwise stated, we employ the units c = G = 1, where c and G are the speed of light and the gravitational constant, respectively.

A. Data and Bayesian inference
We use Bayesian inference to reanalyze GW170817 and GW190425 with various waveform models that incorporate tidal effects in a different manner. Our analysis follows the one performed in our recent work [49], and uses the public data by LVC 1 . We calculate the posterior PDF, p( θ| s(t), H), for the binary parameters θ for the gravitational waveform model, H, given the LIGO Hanford, LIGO Livingston, and Virgo data s(t) via p( θ| s(t), H) ∝ p( θ|H)p( s(t)| θ, H).
(1) p( θ|H) is the prior for the binary parameters. The likelihood p( s(t)| θ, H) is evaluated by assuming stationarity and Gaussianity for the detector noise using the noise power spectrum density derived with BayesLine 2 . We compute PDFs by using stochastic sampling engine based on nested sampling [50,51]. Specifically, we use the parameter estimation software, LALInference [52,53], which is one of the software of LIGO Algorithm Library (LAL) software suite. We take the frequency range from 23  , as the point-particle and spin parts, and NR calibrated tidal effects. The TF2 approximant employs the 3.5PN-and 3PN-order formulas for the phase and amplitude, respectively as the point-particle part, and treats aligned spins and incorporates 3.5PN-order formula in spin-orbit interactions, 2PN-order formula in spin-spin, and self-spin interactions. TF2+ is the TF2 approximant supplemented with phenomenological higher-order PN terms calibrated by SEOBNRv2 for the point-particle part. The TF2+ NRTidal model is another model whose tidal effects are calibrated by NR. The TF2+ NRTidalv2 model is the upgrade of the TF2+ NRTidal model. The TF2 PNTidal and TF2+ PNTidal models employ the PN tidal-part phase formula.

B. Waveform models for inspiraling BNSs
We use various analytic frequency-domain waveform models for the inspiral phase. The features of each waveform model are summarized in Table I. The Fourier transform of the gravitational waveform can be written as where the amplitude A(f ) and the phase Ψ(f ) can be decomposed into the point-particle evolution, the spin effects, and the tidal effects as and We use TF2 [54,55] and phenomenologically extended model of TF2, called TF2+ (see Ref. [48] and below) as BBH baseline, which consists of point-particle and spin parts. Here, the 3.5PN-order formula for the phase and 3PN-order formulas for the amplitude are employed as the point-paticle part of TF2 [56]. For TF2+, both the phase and amplitude of the point-particle part are extended to the 6PN-order by fitting the SEOBNRv2 model [57,58].
All waveform models used in our parameter estimation analyses assume that the spins of component stars are aligned with the orbital angular momentum, and incorporate 3.5PN-order formula in couplings between the orbital angular momentum and the component spins [59], 2PN-order formula in point-mass spin-spin, and self-spin interactions [60,61].
During the BNS inspiral, at the leading order, the induced quadrupole moment tensor Q ij is proportional to the external tidal field tensor E ij as Q ij = −λE ij . The information about the NS EOS can be quantified by the tidal deformability parameter λ [33,62]. The leading order tidal contribution to the GW phase evolution (relative 5PN-order) is governed by the symmetric contribution of NS tidal deformation, characterized by the binary tidal deformability [33] Λ = 16 13 which is a mass-weighted linear combination of the tidal deformability of the both components, where m 1,2 is the component mass and Λ 1,2 is the dimensionless tidal deformability parameter of each star Λ = λ/m 5 . The an- Tidal phase in the frequency domain normalized by the leading, Newtonian (relative 5PN-order) tidal phase formula. Here, we use (m1, m2) = (1.35M , 1.35M ). We showΛ = 1000 (dot-dashed, blue), 400 (dashed, blue), and 100 (dotted, blue) for the KyotoTidal model. The NRTidal model (solid, red), the NRTidalv2 model (solid, cyan), and the 5+2.5PN-order tidal-part phase formula, PNTidal (solid, green), are also presented, which are independent ofΛ when normalized by the leading tidal phase.
tisymmetric contribution δΛ terms are always subdominant on the tidal effects to the GW phase and the symmetric contributionΛ terms dominate [35,38]. In this paper, we ignore the δΛ contribution.
The TF2 PNTidal and TF2+ PNTidal models denote the waveform models employing TF2 and TF2+ as the BBH baseline, respectively. Both the TF2 PNTidal and the TF2+ PNTidal models employ the 2.5PN-order (relative 5+2.5PN-order) tidal-part phase formula [39] Ψ PNTidal is the symmetric mass ratio, and z is the source redshift. The tidal-part amplitude for both TF2 PNTidal and TF2+ PNTidal models employ the 5+1PN-order amplitude formula given by [39] A PNTidal tidal = 5πη 24 where d L is the luminosity distance to the source. The TF2+ KyotoTidal model is a NR calibrated waveform model for the inspiral phase of BNS mergers [48,63]. The TF2+ KyotoTidal model employs TF2+ as the BBH baseline and extends the 2.5PN-order (relative 5+2.5PNorder) tidal-part phase formula [39] by multiplyingΛ by a nonlinear correction to model the tidal part of the GW phase. The functional forms of the tidal-part phase is where a = 12.55 and p = 4.240. The tidal-part amplitude is extended by adding the higher-order PN tidal effects to Eq. (7) as where b = 4251 and r = 7.890. In the KyotoTidal model, the hybrid waveforms constructed from high-precision NR waveforms and the SEOBNRv2T waveforms [45,57,58,64,65] are used for model calibration in the frequency range of 10-1000 Hz. The phase difference between the TF2+ KyotoTidal model and the hybrid waveforms is smaller than 0.1 rad up to 1000 Hz for 300 Λ 1900 and for the mass ratio q = m 2 /m 1 ≤ 1 between 0.73 and 1 [48]. In [48], it is shown that the mismatch between the TF2+ KyotoTidal model and the hybrid waveforms is always smaller than 1.1 × 10 −5 in the frequency range of 10-1000 Hz.
The NRTidal model is another approach to describe tidal effects calibrated by NR waveforms [46]. The TF2+ NRTidal model employs TF2+ as the BBH baseline. For the tidal effects, this model extends the linear-order effects by effectively adding the higher-order PN terms of the tidal contribution to the GW phase. As shown in Ref. [46], the expression of the tidal phase is given by the form of a rational function:  [46]. The TF2+ NRTidalv2 model is an upgrade of the TF2+ NRTidal model [47]. Specifically, they derive a new expression for the tidal phase which is calibrated to more accurate NR waveforms, where the known coefficients arec 1 = 3115/1248,c 3/2 = −π,c 2 = 28024205/3302208,c 5/2 = −4283π/1092, and the fitting coefficients areñ 5/2 = 90.550822,ñ 3 = −60.253578,d 1 = −15.111208,d 2 = 8.0641096. They also introduce the tidal amplitude, where d = 13477.8. Although the new phase model, NRTidalv2, introduced in Ref. [47], includes higher order spin-squared and spin-cubed terms with their associated spin-induced moments, we do not add them in this work. In Fig. 1, we show differences in the phase evolution of tidal part among the KyotoTidal, NRTidal, NRTidalv2, and PNTidal models. A difference in the treatment of the tidal effects makes differentΛ-dependence. The tidal phase normalized by the leading (relative 5PN-order) tidal phase formula for the KyotoTidal model depends on the binary tidal deformabilityΛ due to the nonlinear correction. Since the NRTidal, NRTidalv2, and PNTidal models employ the linear-order effects of the tidal deformability, they are independent ofΛ when normalized by the leading tidal effect. Figure 1 shows good agreement between the TF2+ KyotoTidal model and the TF2+ NRTidalv2 model forΛ 1000 below 1000 Hz as suggested in Ref. [47]. The NRTidal model gives the largest phase shift, the second is the NRTidalv2 model, the third is the KyotoTidal model, and the PNTidal model gives the smallest, forΛ ≤ 1000, up to ∼1000 Hz. The TF2+ KyotoTidal model is calibrated only up to 1000 Hz and overestimates tidal effects at frequencies above 1000 Hz. The KyotoTidal model gives the largest phase shift at frequency above 1200 Hz forΛ = 1000, and larger phase shift than the one for the NRTidalv2 model at frequency above about 1000 Hz (1400 Hz) for Λ = 1000 (400).

C. Source parameters
The source parameters and their prior probability distributions are chosen to follow those adopted in our recent work [49], and we mention specific choices made in this work.
For GW170817, we fix the sky location to the position of AT 2017gfo, which is an electromagnetic counterpart of GW170817 [66], and estimates of the remaining source parameters. Specifically, we estimate the luminosity distance to the source d L , the binary inclination θ JN , which is the angle between the total angular momentum and the line of sight, the polarization angle ψ, the coalescence time t c , the phase at the coalescence time φ c , component masses m 1,2 , where we assume m 1 ≥ m 2 , the orbit-aligned dimensionless spin components of the stars χ 1,2 where χ 1,2 = cS 1,2 /(Gm 2 1,2 ) is the orbit-aligned dimensionless spin components of the stars with S 1,2 are the magnitudes of the spin angular momenta of the components, and the binary tidal deformabilityΛ.
For GW170817, we assume a uniform distribution as the detector-frame component mass prior m 1,2 ∼ U [0.83, 7.7]M with an additional constraint on the detector-frame chirp mass M det : where the chirp mass is the best estimated mass parameter defined as M = (m 1 m 2 ) 3/5 (m 1 + m 2 ) −1/5 . The prior range for M det is the same as that used for LVC analysis [15]. The impact of wider prior range for M det on parameter estimation is negligible. We assume a uniform prior on the spin magnitudes and we enforce χ 1,2 ∼ U [−0.05, 0.05]. This prior range of spin is consistent with the observed population of known BNSs that will merge within the Hubble time [67,68], and is referred to as the low-spin prior for the LVC analysis [15]. We assume a uniform prior on the binary tidal deformability, withΛ ∼ U [0, 3000].
For GW190425, we also estimate the sky location of the source with an isotropic prior. We assume the detectorframe component mass prior m 1,2 ∼ U [1.0, 5.0]M and the spin and the binary tidal deformability priors are the same as the ones for GW170817.

III. RESULTS OF GW170817
A. Source properties other than the tidal deformability In this subsection, we show validity of our analysis as a sanity check by comparison with the LVC results. Figure 2 shows the marginalized posterior PDFs of parameters other than the tidal deformability for various waveform models for f max = 1000 Hz. Table II presents the 90% credible intervals of the luminosity distance d L , the binary inclination θ JN , mass parameters (the component masses m 1,2 , the detector-frame chirp mass M det , the source-frame chirp mass M, the total mass M tot , and the mass ratio q), and the effective spin parameter χ eff = (m 1 χ 1 + m 2 χ 2 )/M tot , which is the most measurable combination of spin components, estimated using various waveform models. The source-frame chirp mass is derived by assuming a value of the Hubble constant H 0 = 69 km s −1 Mpc −1 (a default value in LAL).
For comparison of our analysis with the results of the previous LVC analysis [15,69], we also analyze GW170817 by using the restricted TF2 approximant as the waveform model with 5+1PN-order tidal-part phase formula. This model has the BBH baseline whose amplitude is constructed only from the Newtonian-order pointparticle evolution [54,55,[59][60][61] and is implemented in LALInference. We checked that estimates of parameters other than the tidal deformability we obtained by using the restricted TF2 model are broadly consistent with the LVC results presented in [15,69]. The estimates of parameters other than the tidal deformability presented in Fig. 2 and Table II show almost no systematic bias associated with a difference among waveform models for both BBH baseline and tidal parts. The posterior PDFs of these parameters for f max = 2048 Hz are almost the same as the ones for f max = 1000 Hz as illustrated for the TF2 PNTidal model in Fig. 2. This is due to the fact that the parameters other than the tidal deformability are mainly measured from information at low frequency region [39] and terms up to 3.5PN-order of the point-particle part for the phase are the same among different waveforms. On the other hand, the tidal deformability is mainly measured from information at high frequency region as discussed in the next section 3 .   III. 90% credible interval of the binary tidal deformability,Λ, for GW170817 for various waveform models. We report both the symmetric 90% credible interval (Symmetric) and the 90% highest-posterior-density intervals (HPD), for both fmax = 1000 Hz (left side) and 2048 Hz (right side), where the median is shown as a representative value. 3 We note that the spin-induced quadrupole moments can affect largely estimates of the component spins and mass ratio for large NSs' spins, and combination of the effects of the spin-induced quadrupole moments and the tidal deformability is important to investigate NS EOS as shown in [70].

B. Posterior of binary tidal deformability
Before presenting our results obtained with various waveform models, we first compare our results obtained by using the restricted TF2 model that incorporates the 5+1PN-order tidal-part phase with those from the LVC analysis [15] as a sanity check. The restricted Marginalized posterior PDFs of various parameters for GW170817 derived by various waveform models. The blue, cyan, red, green, and orange curves correspond to the TF2+ KyotoTidal, TF2+ NRTidalv2, TF2+ NRTidal, TF2+ PNTidal, and TF2 PNTidal models, respectively. The top-left, top-middle, top-right, middle-left, center, middle-right, bottom-left, bottommiddle, and bottom-right panels show the mass ratio q, the primary mass m1, the secondary mass m2, the source-frame chirp mass M, the detector-frame chirp mass M det , the total mass Mtot, the luminosity distance to the source dL, the inclination angle θJN, and the effective spin parameter χ eff , respectively. Here, we show the distribution for fmax = 1000 Hz, except for the TF2 PNTidal model, for which the intervals for both fmax = 1000 Hz and fmax = 2048 Hz are given.
sult of 90% credible symmetric (highest posterior density (HPD)) interval onΛ is 347 +564 −243 (347 +453 −295 ) for restricted TF2 with 5+1PN-order tidal-part phase, low-spin prior (|χ 1,2 | ≤ 0.05), and f max = 2048 Hz, the LVC report Λ = 340 +580 −240 (340 +490 −290 ) in [15]. Here, uniform priors in Λ 1 and Λ 2 are adopted in both analyses, and the posterior ofΛ is divided by its prior determined by those of other parameters following [15] to derive approximate re-sults for the case of a uniform prior onΛ. The closeness of the inferred credible ranges indicates that our analysis successfully reproduces the results derived by the LVC. If we assume a uniform prior onΛ from the beginning, 90% credible symmetric (HPD) interval onΛ is 316 +504 −224 (316 +367 −291 ) for restricted TF2 with 5+1PN-order tidal-part phase.  Table III. the binary tidal deformabilityΛ for various waveform models for both f max = 1000 Hz (left panel) and 2048 Hz (right panel). The corresponding 90% credible intervals are presented in Table III. We caution that the TF2+ KyotoTidal model is calibrated only up to 1000 Hz and can overestimate tidal effects at frequencies above 1000 Hz. Thus, the results for f max = 2048 Hz should be regarded as only a reference.
For f max = 1000 Hz (left panel of Fig. 3), the peak values ofΛ are 400-500 and the 90% credible intervals do not extend to 900 for NR calibrated waveform models: the TF2+ KyotoTidal, TF2+ NRTidalv2, and TF2+ NRTidal models. Our results show that the posterior of binary tidal deformability for GW170817 depend on waveform models. The TF2+ KyotoTidal, TF2+ NRTidal, TF2+ NRTidalv2, and TF2+ PNTidal models are constructed from the same BBH baseline, TF2+, but with different tidal descriptions.
Therefore, a difference of estimates among these waveform models reflects directly their different tidal description. The TF2+ NRTidal model gives the smallest median value onΛ of 403, the second is the TF2+ NRTidalv2 model of 445, the third is the TF2+ KyotoTidal model of 481, and the TF2+ PNTidal model gives the largest one of 569. This order is derived from the order of the phase shift of different waveform models for a given value ofΛ = 400, up to about 1400 Hz as shown in Fig. 1. The tendency to give smaller estimated values for NR calibrated waveform models than for PN waveform models are consistent with previous results derived in Ref. [71] (see also Ref. [72] for the detailed study of systematic biases associated with spin effects). The TF2+ PNTidal and TF2 PNTidal models are constructed from the same tidal part and the different point-particle part. A difference in the posterior PDFs of estimatedΛ between these models are very small for f max = 1000 Hz. This result shows that the higherorder point-particle terms do not significantly affect the estimate of the binary tidal deformability of GW170817 for f max = 1000 Hz. (See Ref. [73] for systematic study on the binary tidal deformability by injection of signals with incomplete baselines) For f max = 2048 Hz (right panel of Fig. 3), the peak values ofΛ are 250-400 and the 90% credible intervals do not extend to 850 for NR calibrated waveform models. The width of symmetric 90% credible intervals for f max = 2048 Hz are narrower than those for f max = 1000 Hz, by about 7% for the TF2+ KyotoTidal model, 4% for the TF2+ NRTidal model, 5% for the TF2+ NRTidalv2 model, 13% for the TF2+ PNTidal model, and about 5% for the TF2 PNTidal model, as shown in Table III. These decrease in the width of the interval are consistent with the fact that higher-frequency data are more informative to measureΛ [39]. The peak values of the posterior PDFs ofΛ tend to decrease as f max increases for all waveform models as shown in Fig. 3. The order of peak values of Λ for the different waveform models that incorporate the same BBH baseline, TF2+, is not affected by varying f max as shown in Fig. 3. This is explained by the same reason as that for f max = 1000 Hz. We note that 1400 Hz approximately corresponds to f ISCO for estimated mass range. The TF2 PNTidal model gives a slightly smaller peak value than the TF2+ KyotoTidal model. This cannot be explained only by the feature of the tidal part as shown in Fig. 1. This might be due to the effects of the higher-order point-particle terms or the fact that the data at frequencies above 1000 Hz are dominated by the detector's noise. The difference in the posterior PDFs of estimatedΛ between the TF2+ PNTidal and TF2 PNTidal models for f max = 2048 Hz is larger than that for f max = 1000 Hz (see Fig. 3 and Table III). This is due to the effects of higher-order point-particle terms as discussed in [36,73].
We present marginalized posterior PDFs for source parameters for GW190425 in Fig. 4 and corresponding 90% credible interval in Table  IV. The estimates of parameters other thanΛ presented in Fig. 4 and Table IV are broadly consistent with the LVC results presented in [30] and show almost no systematic bias among different waveform models. The posterior PDFs of these parameters for f max = 2048 Hz are almost the same as the ones for f max = 1000 Hz as illustrated for the TF2+ PNTidal model in Fig. 4. Figure 5 shows marginalized posterior PDFs ofΛ for GW190425 for three waveform models. While this figure indicates that there is a small difference between PN and NR calibrated models, there is almost no difference between two NR calibrated models. The posterior PDF ofΛ have a large value aroundΛ = 0 and this fact means that no significant tidal effect is detected as found in [30]. HPD upper limit on the binary tidal deformability isΛ ≤ 610 for the TF2+ KyotoTidal model for f max = 1000 Hz. The posterior PDF ofΛ for the TF2+ KyotoTidal model for f max = 2048 Hz is bimodal. Investigation of the secondary peak's origin remains as a future work, but it may result from the nonlinear tidal terms ∝ x p of this model, which increase rapidly at 1000 Hz where the calibration by the hybrid waveforms is not performed. In order to discuss constraints on NS EOS by combining information obtained from GW170817 and GW190425, we plot prediction of various NS EOS on posterior of the binary tidal deformability and binary's chirp radius, which is a conveniently scaled dimensionful radius-like parameter [38]. Figure 6 shows 50% and 90% credible regions in theΛ-M plane (left panel) and the M-R plane (right panel), for GW170817 and GW190425, where R = 2MΛ 1/5 is the binary's chirp radius. Five colored curves are posteriors predicted by various NS EOS models; MS1 [74], H4 [75], MPA1 [76], APR4 [77], and WFF1 [78]. For these plots, we use the masses uniformly distributed in the mass ratio range 0.7 ≤ q ≤ 1, which include the 90% credible regions of mass posteriors for GW170817 and GW190425. Our results using TF2+ KyotoTidal model show that softer NS EOS models are preferred, which is consistent with the LVC results presented in [1,[15][16][17]. In particular, the MS1 and H4 models lie outside the 90% credible region for GW170817, while they are not disfavored from GW190425.

V. SUMMARY
We reanalyze GW170817 and GW190425 with a NR calibrated waveform model, the TF2+ KyotoTidal model, which has been developed independently from the one used in previous studies by LVC. The TF2+ KyotoTidal model is calibrated in the frequency range of 10-1000 Hz by hybrid waveforms composed of high-precision NR waveforms and the Marginalized posterior PDFs of source parameters for GW190425 using TF2+ PNTidal, TF2+ KyotoTidal, and TF2+ NRTidalv2 models for the low-spin prior (|χ1,2| ≤ 0.05). Here, we show the distribution for fmax = 1000 Hz, except for the TF2+ PNTidal model, for which the intervals for both fmax = 1000 Hz and fmax = 2048 Hz are given.
SEOBNRv2T waveforms, and reproduces the phase of the hybrid waveforms within 0.1 rad error up to 1000 Hz. In the TF2+ KyotoTidal model, the nonlinear effects of the tidal deformability is incorporated. We also reanalyze the events with other waveform models: two PN waveform models (TF2 PNTidal and TF2+ PNTidal), the TF2+ NRTidal model that is another NR calibrated waveform model, and its upgrade, the TF2+ NRTidalv2 model.
We compare parameter estimation results with different tidal waveform models. For GW170817, there seems to be almost no systematic biases for extraction of source parameters other than the binary tidal deformability using different waveform models. We find that the PN model tends to overestimateΛ compared to the NR calibrated waveform models, while the estimates of Λ also depend on NR calibrated waveform models for f max = 1000 Hz. But the difference is smaller than the statistical uncertainties.
Our results for GW170817 indeed indicate thatΛ is constrained more tightly for f max = 2048 Hz than for f max = 1000 Hz. For the TF2+ KyotoTidal model, the 90% symmetric interval ofΛ for f max = 2048 Hz is about 7% narrower than that for f max = 1000 Hz. Though the credible interval ofΛ becomes narrower as the f max increases, the TF2+ KyotoTidal model is calibrated only up to 1000 Hz. Since higher frequency data are more informative forΛ [39], it is important to improve current waveform models at high-frequencies above 1000 Hz to accurately determineΛ from the GW data, toward third generation detector era.
For the second BNS merger event, GW190425, we use three waveform models; TF2+ KyotoTidal, TF2+ NRTidalv2, and TF2+ PNTidal models. Similarly to GW170817, there are almost no systematic biases for extraction of source parameters other than the bi-nary tidal deformability among different waveform models. This binary system is massive and it is intrinsically difficult to measure the tidal effect. While our results show that the 90% credible interval ofΛ for the PN waveform model is slightly wider than for NR models, there is almost no difference between NR calibrated waveform models.
We discuss constraints on NS EOS models by combining information obtained from GW170817 and GW190425. Our results using TF2+ KyotoTidal model show that softer NS EOS models are preferred, which is consistent with the LVC results presented in [1,[15][16][17]. By using independent waveform model (TF2+ KyotoTidal model) and independent analysis, we obtain consistent results with the LVC's one that a low SNR event from a massive binary like GW190425 cannot contribute very much to constrain NS EOS as shown in [30]. As the number of BNS merger events increases and sensitivities of detectors are improved, the systematic errors will become significant. ACKNOWLEDGMENT We thank John Veitch for very helpful explanation of LALInference and Chris van den Broeck for useful discussions. This work is supported by Japanese Society for the Promotion of Science (JSPS) KAKENHI Grant Numbers  JP15K05081, JP16H02183, JP16H06341, JP16H06342,  JP17H01131, JP17H06133, JP17H06358, JP17H06361,  JP17H06364, JP18H01213, JP18H04595, JP18H05236, and JP19K14720, and by a post-K project hp180179. This work is also supported by JSPS Core-to-Core Program A. Advanced Research Networks and by the joint research program of the Institute for Cosmic Ray Research, University of Tokyo, Computing Infrastructure Project of KISTI-GSDC in Korea, and Computing Infrastructure ORION in Osaka City University. T. Narikawa is supported in part by a Grant-in-Aid for JSPS Research Fellows, and he also thanks hospitality of Chris's group during his stay at Nikhef. K. Kawaguchi was supported in part by JSPS overseas research fellowships. We are also grateful to the LIGO-Virgo collaborations for the public release of GW data of GW170817. This research has made use of data, software, and web tools obtained from the Gravitational Wave Open Science Center (https://www.gw-openscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. LIGO is funded by the U. S. National Science Foundation. Virgo is funded by the French Centre National de la Recherche Scientifique (CNRS), the Italian Istituto Nazionale di Fisica Nucleare (INFN), and the Dutch Nikhef, with contributions by Polish and Hungarian institutes. Marginalized posterior PDFs of binary tidal deformability,Λ, for GW170817 derived by data of different detector combinations with both fmax = 1000 Hz (solid) and 2048 Hz (dashed) for the TF2+ KyotoTidal (left panel) and the TF2+ NRTidalv2 (right panel) models. The distribution derived by the Hanford-only data (blue), that by the Livingston-only data (orange), and that by combined data of Advanced LIGO twin detectors and Advanced Virgo (green, denoted by HLV) are presented. For fmax = 2048 Hz, a multimodal (bump) structure at high-Λ for the TF2+ KyotoTidal (TF2+ NRTidalv2) model appear due to Livingston data. There is a multimodal structure at the high-Λ region in the posterior PDF ofΛ for GW170817 for the TF2+ KyotoTidal model and a bump structure for the TF2+ NRTidal and TF2+ NRTidalv2 models for f max = 2048 Hz as shown in the right panel of Fig. 3. In this appendix, we present an in-depth study to interpret these features by separate analysis for the LIGO twin detectors for GW170817. Figure 7 shows marginalized posterior of Λ derived by separate analysis for the Hanford and Livingston detectors with both f max = 1000 Hz and 2048 Hz for the TF2+ KyotoTidal model (left panel) and the TF2+ NRTidalv2 model (right panel). Table V shows corresponding 90% credible intervals ofΛ.
In the case of the TF2+ KyotoTidal model, the left panel in Fig. 7 suggests that the origin of the bump at the high-Λ region for f max = 2048 Hz for the HLV combined data is as follows. On the one hand, for the Livingston data, the unimodal distribution for f max = 1000 Hz, whose peak is at about 600, is separated into a bimodal distribution for f max = 2048 Hz that is constructed from twin peaks, a low-Λ bump, and a few high-Λ bumps. On the other hand, for the Hanford data, the unimodal distribution for f max = 1000 Hz, whose peak is at low-Λ region, shrinks for f max = 2048 Hz. As a result, for f max = 2048 Hz, the remaining high-Λ peak for the Livingston data produces the bump for the HLV combined data. Moreover, a few high-Λ bumps in the case of HLV combined data for f max = 2048 Hz are inherited from the bumps of the Livingston-only data, which are associated with the high-frequency data. The location of the low-Λ bump derived by the Livingston-only data is close to the peak ofΛ of about 250 derived by the Hanford-only data.
In the case of the TF2+ NRTidalv2 model, as shown in the right panel of Fig. 7, a bump at the high-Λ region in the case of HLV combined data for f max = 2048 Hz are inherited from the peak of the Livingston-only data,Λ ∼ 750. While a bimodal distribution appears in the posterior PDF ofΛ with the SEOBNRv4 ROM NRTidal model in the case of LVC analysis as shown in Fig. 11 in [15], a small high-Λ bump atΛ ∼ 600 appears in that with the TF2+ NRTidal model presented for f max = 2048 Hz in the right panel of Fig. 3.
Here, the SEOBNRv4 ROM NRTidal model is constructed from the SEOBNRv4 model [79,80] as the BBH baseline and the NRTidal model as the tidal part. Supplementary analysis with the TF2+ NRTidal model as shown in Fig. 8 demonstrates that the different priors inΛ (one uniform and one non-uniform) make such different distribution between our analysis and the LVC analysis. The LVC used "Weighted" prior. In this prior, they assume uniform priors in Λ 1 and Λ 2 , and weight the posterior ofΛ by dividing by its prior determined by those of other parameters [15]. "Weighted" prior approximately corresponds to imposing a uniform prior onΛ. Figure 8 shows the dependence of the results on different priors inΛ, "Λ 1,2 -flat", "Weighted", and "Λ-flat" for the TF2+ NRTidal model with f max = 2048 Hz. This figure demonstrate that the distribution for "Λ 1,2 -flat" and "Weighted" prior tends to be a bimodal rather than a high-Λ bump.
In Ref. [49], it is found that there is a discrepancy in the estimates of binary tidal deformability of GW170817 between the Hanford and Livingston detectors of Advanced LIGO by using the restricted TaylorF2 waveform model. Figure 7 shows that the discrepancy is enhanced with sophisticated waveform models (the TF2+ KyotoTidal and TF2+ NRTidalv2 models). While the two distributions in the cases of the Hanford-only and Livingston-only data seem to be consistent with each other and also consistent with what we would expect from noise realization (see e.g., Ref. [38]), the results that the width of the 90% credible interval for the Livingston-only data does not shrink as f max increases indicate that the Livingston's high-frequency data are not very useful to determine the tidal deformability for GW170817. Dependence of the marginalized posterior PDFs ofΛ on different priors inΛ for GW170817 for the TF2+ NRTidal model with fmax = 2048 Hz. In addition to PDF ofΛ for uniform priors in Λ1 and Λ2 (dotted, cyan), we show the PDF for "Weighted"-prior (dashed, magenta), which is weighted by dividing the original prior (also shown by solid yellow curve) and the PDF for a uniform prior inΛ (solid, green).