Discrepancy in tidal deformability of GW170817 between the Advanced LIGO twin detectors

We find that the Hanford and Livingston detectors of Advanced LIGO derive a distinct posterior probability distribution of binary tidal deformability tilde{Lambda} of the first binary-neutron-star merger GW170817. By analyzing public data of GW170817 with a nested-sampling engine and the default TaylorF2 waveform provided by the LALInference package, the probability distribution of the binary tidal deformability derived by the LIGO-Virgo detector network turns out to be determined dominantly by the Hanford detector. Specifically, by imposing the flat prior on tidal deformability of individual stars, symmetric 90% credible intervals of tilde{Lambda} are estimated to be 527^{+619}_{-345} with the Hanford detector, 927^{+522}_{-619} with the Livingston detector, and 455^{+668}_{-281} with the LIGO-Virgo detector network. Furthermore, the distribution derived by the Livingston detector changes irregularly when we vary the maximum frequency of the data used in the analysis. This feature is not observed for the Hanford detector. While they are all consistent, the discrepancy and irregular behavior suggest that an in-depth study of noise properties might improve our understanding of GW170817 and future events.

We find that the Hanford and Livingston detectors of Advanced LIGO derive a distinct posterior probability distribution of binary tidal deformabilityΛ of the first binary-neutron-star merger GW170817. By analyzing public data of GW170817 with a nested-sampling engine and the default TaylorF2 waveform provided by the LALInference package, the probability distribution of the binary tidal deformability derived by the LIGO-Virgo detector network turns out to be determined dominantly by the Hanford detector. Specifically, by imposing the flat prior on tidal deformability of individual stars, symmetric 90% credible intervals ofΛ are estimated to be 527 +619 −345 with the Hanford detector, 927 +522 −619 with the Livingston detector, and 455 +668 −281 with the LIGO-Virgo detector network. Furthermore, the distribution derived by the Livingston detector changes irregularly when we vary the maximum frequency of the data used in the analysis. This feature is not observed for the Hanford detector. While they are all consistent, the discrepancy and irregular behavior suggest that an in-depth study of noise properties might improve our understanding of GW170817 and future events.

I. INTRODUCTION
Tidal deformability of neutron stars can be a key quantity to understand the hitherto-unknown nature of supranuclear density matter (see Ref. [1] for reviews). The relation between the mass and tidal deformability is uniquely determined by the neutron-star equation of state [2,3] as is the mass-radius relation [4]. Thus, simultaneous measurements of the mass and tidal deformability are eagerly desired, and gravitational waves from binary-neutron-star mergers give us a perfect opportunity. Once the mass-tidal deformability relation is understood accurately, binary neutron stars can be used as standard sirens to explore the expansion of the universe even in the absence of electromagnetic counterparts [5]. Motivated by these facts, the influence of tidal deformability on gravitational waves from binary neutron stars has been studied vigorously in this decade [6][7][8][9][10][11][12].
The direct detection of gravitational waves from a binary-neutron-star merger, GW170817, enabled us to measure the tidal deformability of a neutron star for the first time [13]. The LIGO-Virgo collaboration (LVC) reported an upper bound on the most influential combination of tidal deformability parameters of two neutron stars, the so-called binary tidal deformabilityΛ, to be 800 (all the values in this paper refer to 90% credibility) in their discovery paper [13] under the reasonable assumption of small neutron-star spins (later corrected to 900 [14]). Independent analysis in Ref. [15] reported, e.g.,Λ = 222 +420 −138 with the flat prior on the mass of neu-tron stars and the reasonable assumption of a common, causal equation of state for both neutron stars. LVC also reported an updated highest-posterior-density interval, Λ = 300 +420 −230 [14] using sophisticated waveform models [16,17] (see also Ref. [18] for an update), and this is further restricted to 190 +390 −120 if a common equation of state is assumed [19].
All these inferences are made by combining the output of Advanced LIGO twin detectors, i.e., the Hanford and Livingston detectors (and Advanced Virgo). It should be important to examine the extent to which results derived by individual detectors agree, particularly in the presence of a glitch near merger [13]. A study of p-g instability presented the posterior probability distribution ofΛ derived by individual detectors [20], but this is estimated only by incorporating this effect and by allowing high spins, which broaden the distribution ofΛ. Neither consistency nor discrepancy of the derived distribution is discussed.
In this paper, we present our independent analysis of GW170817 to show that the Advanced LIGO twin detectors derive a distinct posterior probability distribution ofΛ (and only for this quantity; see the Appendix). Although the 90% credible intervals ofΛ are nominally consistent between the twin detectors, the distribution derived by the Livingston detector tends to prefer larger values ofΛ than those reported in the literature. Close inspection of the distribution suggests that the difference between the twin detectors might not be purely statistical. Specifically, the distribution derived by the Livingston detector does not behave smoothly with respect to the variation of the maximum frequency of the data used for parameter estimation. This behavior is not expected from physics of tidal deformation and should be contrasted with that of the distribution derived by the Hanford detector. The discrepancy between the twin detectors presages a challenge for determining tidal deformability accurately in future detections.

II. PARAMETER ESTIMATION
We perform a Bayesian parameter estimation of GW170817 for (1) the Hanford data (Hanford-only), (2) the Livingston data (Livingston-only), and (3) combined data of the twin detectors and Advanced Virgo (HLV) as in previous work [13][14][15]18]. The data of GW170817 are made public by LVC. 1 The calibration error is taken into account by using the calibration uncertainty envelope release for GWTC-1. 2 As far as we tested, results derived by the HLV data change very little when the data from Virgo are discarded, as expected from the small signalto-noise ratio [13].
Presuming that gravitational waves are detected in the relevant data, s(t), we compute the posterior probability distribution of binary parameters θ via where p(s(t)|θ) is the likelihood and p(θ) is the prior probability distribution. The parameters, θ, consist of two masses, (aligned components of) two spins, two tidal deformability parameters, the luminosity distance, the sky position, the binary inclination, the polarization angle, the coalescence time, and the coalescence phase [13]. We employ the nested sampling [21,22] for the practical analysis using an engine implemented in the public LALInference package [23], a part of LSC Algorithmic Library Suite. We checked that all the independent realizations of the nested-sampling chains derive consistent results. We also checked that Markov Chain Monte Carlo parameter estimation derive consistent results.
We evaluate the likelihood following the standard procedure (see, e.g., Refs. [24,25]) using the noise power spectrum derived by the relevant data using BayesWave (see Appendix A and Appendix B of Ref. [18]). 3 The noise is assumed to be stationary and Gaussian during the parameter estimation. Following previous work [13,15], we adopt the restricted post-Newtonian Tay-lorF2 approximant as the waveform model (see Ref. [14] and references therein). This choice facilitates comparisons with previous results. Because this approximant is implemented in LALInference, results presented in this work should not be affected by our own analysis method and should be easy to reproduce. The minimum frequency of the data used in the analysis is fixed to 23 Hz, and the maximum frequency f max is varied to investigate its influence on estimation ofΛ. Because we truncate the TaylorF2 approximant above the frequency at the innermost stable circular orbit around a nonspinning black hole with its mass being equal to the total mass of the binary, the results become identical for f max 1600 Hz. In this work, we represent them by f max = 2048 Hz for simplicity.
We caution that the inferred value ofΛ entails systematic errors associated with inaccuracy of the TaylorF2 approximant. While the systematic error is subdominant compared to the statistical error for GW170817 [14], we are also conducting further analysis employing a sophisticated waveform model developed based on numericalrelativity simulations by the Kyoto group [26,27]. Regarding the topic of this paper, preliminary results suggest that this model only enhances the discrepancy between the twin detectors. In particular, the posterior probability distribution derived by the HLV data begins to exhibit a multiply-peaked structure consistently with the LVC analysis performed employing other sophisticated waveform models [14]. These results will be presented in a separate publication focusing on the comparison among waveform models [28].
The prior probability distribution is chosen to follow those adopted in the LVC analysis [14], and we mention specific choices made in this work. The sky position is fixed to the location determined by optical followup observations [29] to save the computational cost. We checked that this has a negligible impact on estimation ofΛ. This is expected, becauseΛ is determined entirely by the phase of gravitational waves, while the sky position affects only the amplitude (see also Refs. [15,30]). It should be cautioned that the sky position cannot be determined to any accuracy by a single detector [31] in the absence of electromagnetic information. The lowspin prior (see Ref. [14]) is adopted for the neutron-star spins for simplicity. The prior of the tidal deformability is chosen to be flat in [0 : 5000] for individual components. This choice neglects the underlying equation of state, and its appropriate incorporation will tighten the constraint onΛ [15,19]. If we impose the flat prior onΛ, the discrepancy between the twin detectors is alleviated but remains. This alleviation is reasonable, because the flat prior onΛ gives weight to the low-Λ region where the discrepancy is mild as we see below.
In this study, we focus primarily on the marginalized posterior probability distribution ofΛ. For completeness, we present estimates of other parameters in the Appendix. The discrepancy between the Advanced LIGO twin detectors is not observed significantly for parameters other than binary tidal deformability. −345 andΛ = 927 +522 −619 , respectively. The distribution obtained by combined data of Advanced LIGO twin detectors and Advanced Virgo (HLV, green) is closer to that derived by the Hanford-only data with its symmetric 90% credible interval beingΛ = 455 +668 −281 . Figure 1 shows the marginalized posterior probability distribution of binary tidal deformability,Λ, derived by the Hanford-only data, Livingston-only data, and combined HLV data with f max = 2048 Hz. The corresponding 90% credible intervals are presented in Table I. The HLV distribution exhibits a peak atΛ ≈ 370 with a tail extending to the high-Λ region consistently with previous work [14,15]. Estimates of other parameters are also broadly consistent with those derived in previous work (see the Appendix).

III. POSTERIOR OF TIDAL DEFORMABILITY
The separate analysis of the data obtained by the individual of twin detectors reveals a mild discrepancy between them. On the one hand, the posterior probability distribution derived by the Hanford-only data is similar to that derived by the HLV data. It exhibits a peak at Λ ≈ 440 with a tail at the high-Λ region. Because the tail structure is enhanced for the HLV data due to the Livingston detector as we discuss below, the 90% credible intervals are similar for these two distributions (see Table I). On the other hand, the distribution derived by the Livingston-only data peaks at a large value ofΛ ≈ 960, which is close to the edge of the 90% credible intervals of either the Hanford-only or HLV distribution. Conversely, the Livingston-only distribution has a small probability aroundΛ ≈ 400. We find that the HLV distribution is approximately reproduced by multiplying the Hanford-only distribution and the Livingston-only distribution if we appropriately incorporate the prior probability distribution ofΛ determined by that of other parameters. Specifically, we need to divide the multiplied distribution by the prior, because it is included in both the Hanford-only and Livingstononly distributions. The division by the prior reduces the probability at the high-Λ region, and thus the data of the Livingston detector have a minor influence on the combined result. Still, the HLV distribution shows a bump at Λ ≈ 1000 inherited from the peak of the Livingston-only distribution.
Detailed features of the data from individual detectors are clarified by examining changes of the posterior probability distribution with respect to the variation of the maximum frequency, f max , imposed in the data analysis. Figure 2 shows the results for f max = 800, 1100, 1400, and 2048 Hz (the same as Fig. 1). The Hanfordonly distribution shifts smoothly to the low-Λ region as f max increases, and the median value decreases approximately monotonically. This is reasonably expected, because the tidal deformability is primarily determined by the gravitational-wave data at high frequency [9,15]. This feature results from the nature of tidal interaction and does not rely on complicated relativistic effects.
Dependence of the marginalized posterior probability distribution on the maximum frequency, fmax, imposed in the data analysis. We adopt fmax = 800, 1100, 1400, and 2048 Hz. While the distribution derived by the Hanford-only data (blue, denoted by H) shifts smoothly to the low-Λ region as fmax increases, that by the Livingston-only data (orange, denoted by L) shows a turn at fmax ≈ 1100 Hz.
On another front, the posterior probability distribution derived by the Livingston-only data shows irregular behavior with respect to the increase of f max . At the beginning, the distribution shifts to the low-Λ region in a similar manner to the Hanford-only distribution when f max is decreased from 800 Hz. However, the shift turns around at f max ≈ 1100 Hz, and then the distribution moves back to the high-Λ region. The median values presented in Table I clearly exhibit the turn, and neither a systematic decrease nor increase is observed. This is not naturally anticipated from the viewpoint of measurability [9,15].
Taking the fact that the distribution ofΛ obtained by the combined HLV data is similar to that by the Hanfordonly data into account, the peculiar dependence on f max might indicate that the high-frequency data of the Livingston detector are not very helpful to determineΛ of GW170817. Note that the signal-to-noise ratio is larger for Livingston than for Hanford due to the higher sensitivity [13,14]. Specifically, the signal-to-noise ratios derived in our analysis are 18.8, 26.8, and 32.6 for Hanford, Livingston, and HLV, respectively. Thus, the dominance of Hanford data in the combined result is not reasonably understood from the signal-to-noise ratio. Actually, we find that the posterior probability distributions of other parameters are more strongly influenced by Livingston than Hanford (see the Appendix).
A careful examination of Fig. 2 reveals that a small bump appears for f max 1400 Hz at the low-Λ region of the posterior probability distribution derived by the Livingston-only data. The location of this bump is close to the peak ofΛ derived by the Hanford-only data. Thus, the Livingston-only distribution may consist of the main peak with a large value ofΛ associated with the lowfrequency data and the side peak with a small value ofΛ associated with the high-frequency data.

IV. DISCUSSION
Our analysis suggests that the noise in the highfrequency region of the Livingston data somehow corrupted information about the tidal deformability of GW170817. Although this could simply be caused by the stationary and Gaussian noise, it could be worthwhile to look for possible peculiar features in the data of GW170817, e.g., a residual of the glitch in the Livingston data at about a half second before merger [13] (but see also Ref. [32]) as it is quite important to estimate tidal deformability accurately. By contrast, the Hanford data seem to be well-behaved during the reception of GW170817.
Having said that, it is ultimately impossible to judge which of the Advanced LIGO twin detectors provides us with a more reliable estimate of the binary tidal deformability than the other does, or whether simply their combination is the most reliable, without meaningful data obtained by a third detector. It should be emphasized that the 90% credible intervals are consistent between the twin detectors. What we may safely conclude is that the posterior probability distribution is exceptionally distinct for binary tidal deformability (see the Appendix for other parameters) and that the Livingston data are not very useful for constraining its value in the case of GW170817. Secure parameter estimation will be helped by unambiguous detection by other instruments such as Advanced Virgo or KAGRA [33]. However, if the irregular loss of information is typical for detections with a moderate signal-to-noise ratio, accurate determination of tidal deformability will remain challenging unless its origin is identified.
to-Core Program A, Advanced Research Networks and by the joint research program of the Institute for Cosmic Ray Research, University of Tokyo, and the Computing Infrastructure Project of KISTI-GSDC in the Republic of Korea. T. Narikawa was supported in part by a Grant-in-Aid for JSPS Research Fellows, and he is also grateful to the 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 collaboration for the public release of gravitational-wave 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.

Appendix A: Parameter other than tidal deformability
We present estimates of parameters other than binary tidal deformability with f max = 2048 Hz for completeness. Table II presents the 90% credible intervals of the luminosity distance, the binary inclination, mass parameters, and the effective spin parameter derived by different data. We recall that the sky location is fixed by the information from electromagnetic observations and that the low-spin prior is imposed. The cosmological redshift is not taken from the host galaxy NGC4993 [30] and is determined from the luminosity distance by assuming the Hubble constant H 0 = 69 km s −1 Mpc −1 (a default value in LAL) to derive the chirp mass in the source frame. The marginalized posterior distribution is presented in Fig. 3. The consistency of our results with previous work [14,15] serves as an important sanity check. They also clarify that the estimation of binary tidal deformability is exceptionally delicate for GW170817. Table II shows that the credible intervals agree remarkably between the Hanford-only and Livingston-only data. The marginalized posterior probability distribution depicted in Fig. 3 not only confirms this agreement but also shows that the distribution is approximately identical. The parameters shown here are estimated primarily from information at low frequency, where the gravitationalwave signal spends most of time with a large number of cycles [9,15]. The detector noise is also less severe at lower frequency. As a result, the situation is different from that of binary tidal deformability discussed in the main text.
The 90% credible intervals derived by combining the HLV data are very close to those of a single detector (Hanford-only or Livingston-only) except for the binary inclination, θ JN . The reason for the difference in θ JN is that the degeneracy of the reflection with respect to the orbital plane, i.e., face-on or face-off, is resolved when the HLV data are combined. The resolution is clearly shown in the middle panel of the bottom row of Fig. 3, where the bimodal structure for a single detector is changed to a single peak for the HLV data favoring the face-off orientation. The posterior probability distribution of other parameters shifts only moderately.
Close inspection reveals that the posterior probability distribution derived by the HLV data closely follows that of the Livingston-only data for the chirp mass in the detector frame, the mass ratio, and the effective spin. Because they are determined by the gravitational-wave phase at low frequency, the closeness of the distribution can be understood as a result of the large signal-to-noise ratio for the Livingston detector. Again, these features are different from what we observe for the binary tidal deformability, where the distribution derived by the HLV data is closer to that derived by the Hanford-only data than the Livingston-only data. II. 90% credible interval of the luminosity distance, the binary inclination, mass parameters, and the effective spin parameter derived by different data with fmax = 2048 Hz. We show 10%-100% regions of the mass ratio with the upper limit q = 1 imposed by the prior, and those of m1 and m2 are given accordingly. We give symmetric 90% credible intervals, i.e., 5%-95%, for the other parameters with the median as a representative value. The binary inclination derived by either the Hanford-only or Livingston-only data is undetermined up to the orbital-plane reflection, and thus the 90% credible interval is very large.
Hanford-only Livingston-only HLV   Table II for the definition of quantities). The distribution of θJN for a single-detector data exhibits a bimodal structure due to the degeneracy of the orbital-plane reflection, and this is resolved for the HLV data.