Resolving phase transition properties of dense matter through tidal-excited g-mode from inspiralling neutron stars

The investigation of the phase state of dense matter is hindered by complications of first-principle nonperturbative quantum chromodynamics. By performing the first consistent general-relativistic calculations of tidal-excited g-mode of neutron stars with a first-order strong interaction phase transition in the high-density core, we demonstrate that gravitational wave signal during binary neutron star inspiral probes their innermost hadron-quark transition and provides potent constraints from present and future gravitational-wave detectors.


INTRODUCTION
Despite decades of exploration, the phase state of cold, dense matter is still a fundamental problem that remains unresolved.This problem pertains to low-energy nonperturbative quantum chromodynamics (QCD) (Kurkela et al. 2010) and is one of the key scientific goals in current and anticipated relativistic heavyion collision experiments (Huth et al. 2022).It is widely believed that this type of matter exists in the core of neutron stars (NSs) (Weber 1999;Heiselberg & Pandharipande 2000;Lattimer & Prakash 2001).The equation of state (EOS) of dense stellar matter has drawn much attention in the era of gravitational wave (GW) astronomy.This is because the observations of GW170817 led to the first measurement of NS tidal deformability (Abbott et al. 2017a), which offered valuable insight into the phase state of dense matter (Abbott et al. 2018).
The tidal deformability accounts for the quadrupole deformations of inspiralling NSs, which are induced by the tidal force from each other, when the quasiequilibrium approximation can be utilized (i.e., the tidal field changes slowly).While in the late stage of the inspiral phase, once the binary's orbital frequency approaches the frequency of a particular normal mode of NSs, a resonance can occur.This resonance will transfer Corresponding author: Ang Li liang@xmu.edu.cn the orbital energy to the stellar oscillation via a mechanism referred to as the dynamical tide.This event will accelerate the orbital decay and result in a phase shift in the GW waveform, which will provide insights into the NS inner structure with the help of GW detection from inspiralling NSs.
Previous studies have focused on the interfacial (i) modes, which are excited at the interface between the solid crust and fluid core (Tsang et al. 2012;Pan et al. 2020;Zhu et al. 2023).These modes typically have frequencies of ∼ 100 Hz and contribute ∼ O(0.1) phase shifts during resonance, making them promising candidates for detection using Advanced LIGO (aLIGO) (LIGO Scientific Collaboration et al. 2015) and 3rd generation GW detectors such as Einstein Telescope (ET) (Maggiore et al. 2020) and Cosmic Explorer (CE) (Reitze et al. 2019).Attempts have also been made to infer NS EOS through tide-induced f-modes and g-modes, which arise from the composition gradients within the NS's interior.Studies (see e.g., Kokkotas & Schmidt 1999;Ho & Lai 1999;Reisenegger & Goldreich 1994;Lai 1994;Yu & Weinberg 2017a;Xu & Lai 2017) have shown that f-modes are challenging to excite due to their high resonant frequencies, whereas g-modes demonstrate negligible phase shifts because their couplings to tidal fields are meager.
In this work, we investigate a unique type of g-modes, known as discontinuity g-modes, which are expected to occur in NSs that possibly have a strong phase transi-tion in their high-density inner cores (Finn 1987) 1 .Our study showcases how discontinuity g-modes' dynamical tide can be leveraged to probe the dense matter's phase transition properties using the data obtained from GW170817, along with possible observations from the next-generation GW detectors.

Relativistic g-modes and tidal coupling
The existence of a quark core inside NSs (also called hybrid stars) has been discussed with Newtonian formulations of perturbations and tidal couplings, focusing on the i-mode at the interface of a crystalline quark core and a fluid hadronic envelope (Lau & Yagi 2021).With the solid quark core condition relaxed, we here employ a full general relativistic (GR) formalism as follows: 1) Constructing the background spherical hybrid stars in hydrostatic equilibrium using Einstein's field equations, i.e., solving the Tolman-Oppenheimer-Volkoff (TOV) equations, and 2) Computing their discontinuity g-modes by solving the linear pulsation equations in general relativity (Lindblom & Detweiler 1983;Detweiler & Lindblom 1985).Our study focuses solely on the l = 2 g-modes, as they have the strongest coupling to the tidal gravitational field.
Our approach prioritizes maintaining consistency in how we treat the oscillation mode and the tidal overlap integral Q within the framework of GR (see Appendix .2 for a derivation of this form): Here, g represents the determinant of the metric of the background spacetime and e 2Φ = g tt is the (t, t) component of the metric.Y 2m is the l = 2 spherical harmonic.
, where e 2λ = g rr is the (r, r) component of the metric and ⟨ξ µ |ψ µ ⟩ = √ −ge −2Φ d 3 x(p+ε)ξ * µ ψ µ is the relativistic inner product.f is the frequency of the mode and ω = 2πf is the angular frequency.The present self-consistent approach effectively mitigates the issue of "f-mode contamination" (Reisenegger 1994) in the calculation of Q encountered in the hybrid (GR background + Newtonian perturbation) approach commonly used previously.We provide the full description of this relativistic formula in Appendix .2,where we've also validated that our relativistic Q satisfies both orthogonality and the sum rule.
We determine the mass-radius (M -R) relationship of NSs through the solution of the TOV equations, using the EOS -i.e., pressure (p) as a function of energy density (ε) -as our basic input.To describe the fiducial soft and stiff hadronic matter (HM), we employ the QMF (Zhu et al. 2018) and NL3ωρ (Horowitz & Piekarewicz 2001) models, ensuring their consistency with results from laboratory nuclear experiments.For modeling the sharp first-order hadron-quark phase transition and representing quark matter within the highdensity cores of NSs, we utilize the well-established generic constant-speed-of-sound (CSS) parameterization (Alford et al. 2013).Previous studies have thoroughly discussed the weak density dependence of the speed of sound in quark matter (see e.g., Alford et al. 2013Alford et al. , 2015)).This approach enables us to establish the EOS from the start of the phase transition up to the maximum central density of the star by using three phase-transition parameters: the transition pressure p t (or the transition density ε t ≡ ε HM (p t )), the transition strength ∆ε and the speed of sound in quark matter c s .The transition density and transition strength affect the tidal overlap integral and its typical value remains at ∼ O(0.01).
It is worth noting that since the discontinuity g-modes are a result of the discontinuity in the density profile of the star, we expect the frequency of these modes to depend on both the transition strength and the transition density.Based on our calculations, we observe that the g-mode frequency typically falls within the range of 100 Hz to 1500 Hz, which is strongly dependent on the transition strength.The transition strength varies between ∼ 5 MeV/fm 3 and ∼ 300 MeV/fm 3 .Our study also indicates that the speed of sound (c s ) has a subdominant effect on mode frequency.Specifically, the difference in mode frequency between the c s = 1/ √ 3 case (corresponding to the conformal limit) and the c s = 1 case (corresponding to the casual limit) is found to be less than ∼ 10%.Thus, for the purpose of our analysis, we solely focus on extremely stiff quark matter with a speed of sound c s = 1.

Waveform correction
Throughout the inspiral process, the binary's individual components will be subject to external tidal forces from their respective companions.Upon the tidal driving frequency nearing the g-mode frequency of the NSs, it becomes possible to trigger internal stellar oscillations through resonance, ultimately causing a phase shift in ) denotes the GW frequency corresponding to the innermost stable circular orbit, at which we assume the signal of inspiral phase shut off (Cutler & Flanagan 1994).Also plotted are the detectability thresholds of aLIGO, ET and CE, which are calculated from their designed sensitivities by assuming the system is located at DL = 100 Mpc.
the GW waveform.With the mode frequency and mode eigenfunction obtained in our GR framework, we then are able to estimate the phase shift of the waveform as (Lai 1994): where q = M ′ /M with M ′ representing the companion star mass.M 1.4 = M/1.4M ⊙ and R 12 = R/12 km.The total phase correction for a binary NS waveform h(f ) = Ae iΨ(f ) , which takes into account the mode resonances of both components, is expressed as (Yu & Weinberg 2017b;Pan et al. 2020): where f i is the g-mode frequency of the i-th star and δ φi represents the corresponding phase shift due to resonance.Here, Θ(•) is the Heaviside step function.
The second line of Eq. ( 3) introduces the total phase shift δ φ = i δϕ i and the resonant frequency f = δ φ/ i (δϕ i /f i ) to reduce the number of parameters (Pan et al. 2020).For simplicity, |δ φ| will be used as the waveform parameter from here on.

Detectability
Fig. 1 displays the total phase shift |δ φ| against the resonant frequency f for two equal-mass binary systems, along with the detection thresholds of aLIGO and 3rd generation detectors.These thresholds are computed using a Fisher matrix method (Cutler & Flanagan 1994) and are based on the designed sensitivities (Abbott et al. 2017b) assuming that the binary systems are located at D L = 100 Mpc (details in Appendix .5).Generally, |δ φ| decreases with the increase of f since δϕ ∝ f −2 [cf.Eq. ( 2)] (Lai 1994).For a significant number of hybrid star models, the total phase shifts are greater than the detectability threshold of 3rd generation detectors.It is worth noting that there is a threshold frequency f th , which is independent of the hadronic EOS.Binaries with f ≲ f th may be detectable by 3rd generation detectors, while those with f ≳ f th may not be detectable.For a 1.4 M ⊙ − 1.4 M ⊙ system, f th ≃ 1200 Hz, while for a 1.8 M ⊙ − 1.8 M ⊙ system, f th ≃ 1050 Hz.These threshold frequencies lead to constraints on the phase transition parameters, which are illustrated in the shaded regions (orange for 1.4 M ⊙ − 1.4 M ⊙ system and blue for 1.8 M ⊙ −1.8 M ⊙ system) in Fig. 2. It is evident that the detection of a signal of mode resonance provides a new opportunity to probe the properties of the hadron-quark phase transition, which can serve as a complementary constraint to other observations, such as the ∼ 2 M ⊙ maximum mass constraint from the most massive pulsars (grey regions in Fig. 2).

GW170817 constraints
We proceed by making use of the available GW170817 data and conduct a search for mode resonance signals in full range of plausible resonant frequencies [50, 1600] Hz (the upper limit corresponds to the shut-off frequency of the waveform, 4400 Hz/(M + M ′ ) ≈ 1600 Hz, for the GW170817 system)2 , with the application of waveform correction as described in Eq. ( 3).This is similar to a previous search conducted in Pan et al. (2020) for signals in the frequency range of [30,300] Hz.We employ a Bayesian parameter estimation approach, following the data analysis framework presented in Abbott et al. (2019).Our analysis utilizes PyMultiNest (Buchner et al. 2014) implemented in Bilby (Ashton et al. 2019), enabling us to simultaneously obtain the posterior distribution of parameters and the model evidence Z.By denoting H 1 as the hypothesis with mode resonance and H 0 as the one without, we are able to compute the Bayes factor from the evidence, B 1 0 = Z 1 /Z 0 .We obtain a Bayes factor of B 1 0 = 0.72, which implies that the current data does not strongly favor one hypothesis over the other (see Appendix .6 for more discussions on the Bayesian parameter estimation).
Though the search yields no results, the posterior distribution of ∆| φ| and f can be utilized to set constraints on the phase transition parameters.In Fig. 3 we present the calculated ∆| φ| and f for a GW170817-like system with a chirp mass of M = 1.186M ⊙ and a mass ratio range of q = 0.7-1.The posterior distribution for ∆| φ| and f together with their 95% and 68% contours are also plotted.It can be seen that for a certain transition density, the hybrid star models with small resonant frequencies should be excluded.For example, those models with p t /ε t = 0.11 and f ≤ 485 Hz are ruled out at 95% confidence level for QMF+CSS.Converting these constraints on f to phase transition properties, we conclusively exclude the possibility of a weak phase transition occurring at low densities with a confidence level of 95%.This exclusion is represented in Fig. 2 by the red region.Our conclusion here is in agreement with our previous findings based on static tidal deformability constraints from GW170817 (Miao et al. 2020) and the study of proton and lambda production during relativistic heavy-ion collisions at several GeV/nucleon incident beam energies (Li et al. 2023).

DISCUSSION
The investigation of phase transitions in QCD matter is a central area of interest in physics and has relevance for a range of open problems, such as early Universe (Hogan 1983;Cai et al. 2022), core-collapse supernovae (Fischer et al. 2018), cosmological gamma-ray bursts (Cheng & Dai 1996), as well as the nonperturbative properties of QCD itself (Huber 2020).We propose to probe the phase transition properties possible in dense NS matter through the tidal resonant excitation of g-modes.
Our analysis demonstrates that the tidal excitation of g-modes from inspiralling neutron stars can induce significant phase shifts in the gravitational waveform, on the order of ∼ O(0.1) − O(1).These shifts are detectable by both aLIGO and 3rd generation detectors.Notably, these low-frequency signatures complement previous efforts to identify first-order phase transitions using high-frequency post-merger gravitational waves (Bauswein et al. 2019), expanding the range of density regimes that can be probed for potential phase transitions.We provide two analytical formulas to clarify how the phase shift in the gravitational waveform relates to the properties of the phase transition.These formulas will serve as valuable tools for future attempts to directly estimate phase transition parameters using GW data.Our analysis suggests that the estimation er-rors are moderate, with a maximum of ∼ 15% for the mode frequencies and 25% for the tidal overlap integrals.This implies a potential underestimation or overestimation of the phase shift by no more than a factor of 2 to 3, as detailed in the Appendix .3-.4.Even in cases where no mode resonance signal is detected, the search results can effectively rule out a significant parameter space, as demonstrated earlier with GW170817.
To calculate the correction to the gravitational waveform induced by the g-mode resonance, we use the form presented in Eq. ( 3).However, it is important to note that this form is derived under the assumption that the energy transfer of the resonance is instantaneous (Yu & Weinberg 2017b).
Therefore, it is only valid when the duration of the resonance, t res ≃ 0.009s (M/1.2M ⊙ ) −5/6 (f /600 Hz) −11/6 , is much shorter than the orbital decay time-scale, t D ≃ 0.08s (M/1.2M ⊙ ) −5/3 (f /600 Hz) −8/3 (Lai 1994).For resonant frequencies that are very high (e.g., ≳ 1000 Hz), the waveform correction ∆Ψ(f ) should be modified to account for the continuous process of energy transfer during the resonance.It is also important to consider the relativistic corrections to the Newtonian orbit as the binary approaches closer, especially for high resonant frequencies.As demonstrated by Bildsten & Cutler (1992), tidal interaction is not sufficient to spin up the binary NS system before the merger and an irrotational merger is expected in most cases due to the spinning-down by magnetic dipole radiation before the merger.At present, we do not take into consideration the star rotation.Nevertheless, NS spin could be important in cases such as dynamic captures, which can shift the mode frequency to lower values and thus can contribute to larger GW phase shift (Ho & Lai 1999).These improvements will be the subject of future work.
In addition, it is crucial to identify the physical origin of the mode resonance signature once it is detected.This is because there could be other possible origins apart from the resonant excitation of hybrid star g-modes.
As suggested by Pan et al. (2020), other processes such as tidal-p-g coupling (Weinberg et al. 2013) or those in modified gravity (Palenzuela et al. 2014;Annulli et al. 2019;Krüger & Doneva 2021;Mendes et al. 2021) predict different corrections on the waveform.Therefore, one could perform a Bayesian model selection to compare these different origins.On the other hand, crustcore i-mode or core composition g-mode can produce a ∆Ψ(f ) similar to that studied in this work.For the former, its resonant frequency is typically around 100 Hz, while for the latter, its frequency can reach up to several hundred Hz according to recent studies using the Cowling approximation (see e.g.Jaikumar et al. (2021)).
Our Fisher analysis shows that the resonant frequencies can be measured with a precision of ≲10% using 3rd generation detectors (see Appendix .5).This precision would enable us to discriminate between the discontinuity g-mode resonance considered here and the crust-core i-mode or composition g-mode by comparing their frequencies.It is worth mentioning that the existence and properties of QCD phase transition are still very unclear in nuclear physics studies.Various scenarios, including cross-over/first-order transitions and fast/slow transitions, lack effective constraints.The analysis of this paper is based on the slow first-order phase transition scenarios.We might expect future GW observations of binary NS mergers could exclude or put constraints on such scenarios, hence contributing to our understanding of the QCD phase diagram.

APPENDIX .1. Solving pulsation equations in general relativity
For the background star, the static and spherically symmetric line element is given by where the metric components Φ and λ are functions of r only.We model the star as a perfect fluid, with the stressenergy-momentum tensor given by where u µ is the four-velocity.We also assume the metric perturbation as and the Lagrangian displacement as where Y lm (θ, ϕ) is the spherical harmonics.Note the Eulerian perturbations of energy density and pressure are given by δε = (p + ε)(∆n/n) − ξ i dε dx i , δp = γp(∆n/n) − ξ i dp dx i . ( where γ = [(p + ε)/p](∂p/∂ε) ad is the adiabatic index.The Lagrangian perturbation of baryon number density ∆n can be derived from the baryon number conservation ∇ µ (nu µ ) = 0, it reads with g (3) the determinant of the metric of the 3-geometry at constant time.Plugging these perturbations (Eqs.3-5) into linear Einstein equation δG µν = 8πδT µν , one finds (Lindblom & Detweiler 1983;Detweiler & Lindblom 1985) and Here the prime denotes the derivative of r.Quark matter is present from the first-order transition onset up to the maximum central pressure of a star (see Fig. 4).Our description of the EOS is written as follows: where 'HM' denotes hadronic matter.We then follow the method described in Lindblom & Detweiler (1983) to solve these equations and find the eigenvalue of ω for each mode.It's particularly important to note that when solving the aforementioned perturbations equations, we require the boundary condition that ∆p = −r l e −Φ XY lm e iωt to be continuous at the hadron-quark interface.This necessitates that at the transition radius R t , where the subscripts '+' and '−' denote the values at R t + 0 + and R t − 0 + , respectively.This implies that at both ends of the interface, although the radial displacement W is continuous, the tangential displacement V is discontinuous due to the presence of a density discontinuity.It is this particularity that gives rise to the discontinuity g-mode studied in this work.Fig. 5 provides a solved example, allowing us to observe the density jump and the corresponding variation of V at the interface. .

Relativistic tidal overlap integral
In some previous works, the hybrid (GR background + Newtonian perturbation) approach was typically used to calculate oscillation modes and the Newtonian formalism was employed to compute tidal overlap integral, i.e., Q N α ≡ ⟨ξ α |∇(r l Y lm )⟩/⟨ξ α |ξ α ⟩ with the inner product ⟨ξ|ψ⟩ = ρξ * • ψd 3 x.This hybrid approach is inconsistent and will result in computed eigenfunctions that are not orthogonal (Reisenegger 1994), means ⟨ξ α |ξ β ⟩ ̸ = 0 when α ̸ = β, see Table 1.In this case, the obtained Q N for g-modes would contain some f-mode contamination, leading to unreliable values of Q N .
In this work, we calculate the stellar oscillation modes within the GR framework.This naturally ensures that the eigenfunctions are orthogonal under a relativistic extension of inner product1 Based on this inner product, we find that the relativistic tidal overlap integral can still be defined as In the following, we give a brief derivation and test the accuracy of our numerical results.We begin from the orthogonal part of the Euler equation ∇ ν T µν = 0, which is expressed as where ⊥ µν = g µν + u µ u ν is the projection operator.Plugging the the Lagrangian displacement of a certain mode α, ξ µ α = (0, ξ α (r)e iωαt ), into Eq.( 17) one can derive the linear pulsation equation with During the inspiral, the linear perturbation of the tidal potential is governed by Plugging Eq. ( 18), u µ = (e −Φ , 0, 0, 0), δu i = e −Φ ∂ t ξ i and δΓ µ 00 = −1/2g µρ ∂ ρ h tid 00 into Eq.( 20), one obtains where U = −1/2h tid 00 = −M ′ /|r − D| is the tidal potential, D is the distance of the secondary component M ′ .Eq. ( 21) is similar to its Newtonian counterpart, as seen, for instance, in Eq.(2.5) of Lai (1994).Hence, similar to Lai (1994), we arrive to the relativistic definition of tidal overlap integral as To test the accuracy of our numerical results, we first calculate the inner product ⟨ξ g |ξ α ⟩, the exemplary results are shown in Table 1.It can be seen that for α ̸ = g, the inner products we obtained are very tiny, which verifies the orthogonality of the eigenfunctions.We also check the sum rule of , which is a general property of the tidal overlap integral (Reisenegger 1994).From Table 1 we see that our numerical Q obeys the sum rule.
.3.Approximated formulas of the mode frequency Previously, Finn (1987) pointed out that when the discontinuity near the surface (i.e., R − R t ≪ R t ), the frequency can be approximated as And this approximation has been improved by Zhao & Lattimer (2022) with where M t is the mass inside R t and D = 1.21 is a fitting coefficient.
However, both approximations of Finn (1987) and Zhao & Lattimer (2022) start from the analytical solution of the frequency of surface gravity waves between two incompressible fluids, obtained by Landau & Lifshitz (1959) assuming the slab configuration.Thus it should be noted that their approximations can not seamlessly cater to both the R − R t ≪ R and R t ≪ R cases, see Fig. 13 in Zhao & Lattimer (2022).In this work, we present a more accurate approximation that can accommodate both these two cases simultaneously, i.e., where ς = (1/5)(∆ε/ε t )[3(κ − 1) + (3 + 2κ)x 5 ] with x = R t /R and κ = (M t /M )(R/R t ) 3 .K = 0.59 is a fitting coefficient.Fig. 6 shows the approximation of Eq. ( 24).It can be seen that the estimation of Eq. ( 24) performs well for both R t < R/2 and R t > R/2 regions.The estimation errors of Eq. ( 24) are found to be not more than ∼ 15%.We should emphasize here that, apart from a factor K, Eq. ( 24) is the exact solution for a two-layer fluid in a spherical configuration within the Newtonian limit.In the following we present a brief derivation.24) versus the frequencies calculated numerically from relativistic pulsation equations.
We begin with the pulsation equation under the Cowling approximation (Cowling 1941) where ξ r and ξ h are the radial and tangential displacement, respectively.g is the gravitational acceleration.N 2 = g 2 (1/c 2 e −1/c 2 ad ) is the Brunt-Väisälä frequency with c e and c ad the equilibrium and adiabatic sound speed, respectively.For simplicity, we shall omit the gravity and composition gradient, i.e., we assume c ad = c e , gr/c 2 ad = 0 and ω 2 r 2 /c 2 ad = 0.The density for the lower layer (i.e., quark core) and upper layer (i.e., hadronic envelope) are denoted as ρ − and ρ + .Restrict to l = 2, one finds that Eq. ( 25) reduces to The general solutions can be obtained with and At the interface of two layers, ξ r and ∆p = ρg(ω 2 rξ h − ξ r ) should be continuous, which requires where ∆ρ ≡ ρ − − ρ + .At the surface, the Lagrangian perturbation of the pressure must vanish, i.e., ∆p = 0. Using Eq. (27-29) and replacing ρ + (ρ − ) with ε t (ε t + ∆ε), we finally obtain .4.Approximated formulas of the tidal overlap integral We find a reliable approximation for the tidal overlap integral, as shown in Fig. 7.We should mention that such an approximation is only valid for small transition strength, i.e., ∆ε ≲ 150 MeV/fm 3 .However, higher transition strength (∆ε ≳ 150 MeV/fm 3 ) corresponds to higher mode frequencies (≳ 1200 Hz) that we are not interested in, therefore we do not need to worry about them when using this approximation.By defining Q We fit the data points of the 1.4 M ⊙ stars and the resulting fitting coefficients are presented in .

Detectability analysis
Following Cutler & Flanagan (1994); Yu & Weinberg (2017b), we use the Fisher information matrix to estimate the detectability of the phase shift due to resonance, assuming that the GW signal have a large signal-to-noise ratio (SNR).The Fisher information matrix for a set of parameters {θ i } is defined by where the inner product is defined by with h and S n are the Fourier transform of GW strain data and the spectral noise density of the detector.The SNR is given by SNR = (h|h) 1/2 .The root-mean-square (rms) measurement error of θ i can be written as where Γ −1 is the inverse of the Fisher information matrix.If the rms error of phase shift ∆|δ φ|is smaller than |δ φ|, then the phase shift is detectable.We employ the IMRPhenomD NTidal waveform (Dietrich et al. 2017), together with the addition of phase correction given by Eq. ( 5), to calculate the GW signal.The waveform parameters are {θ i } = {M, q, Λ 1 , Λ 2 , d L , t c , ϕ c , |δ φ|, f }. (35) Here for simplicity, we omit the angular dependence and the stellar rotation.We also fix the tidal defromability as Λ 1.4 = 400 and Λ 1.8 = 300, because it is found that ∆|δ φ| is insensitive to the choice of tidal deformability.In Table 3 we list some quantitative results of rms measurement errors for a 1.4 M ⊙ − 1.4 M ⊙ system.
.6.Bayesian parameter estimation for GW170817 We perform a Bayesian parameter estimation to analysis the GW170817 data by incorporating the resonant parameter |δ φ| and f .We utilize the IMRPhenomD NTidal waveform for the basic hypothesis H 0 , whereas add the phase correction for the hypothesis H 1 .The waveform parameters are given by {θ i } = {M, q, Λ 1 , Λ 2 , χ 1z , χ 2z , θ jn , t c , ϕ c , Ψ, |δ φ|, f } (36) We fix the location of the source to the position determined by electromagnetic observations (Abbott et al. 2017c;Levan et al. 2017)  In Fig. 8 we show the posterior distribution of waveform parameters.Since we use the PyMultiNest to sample from the posterior distribution, we can obtain the evidence (Z) of each hypothesis directly.We then calculate their Bayes Factor with B 1 0 = Z 1 /Z 0 .We also calculate the Bayes Factor with the Savage-Dickey Density Ratio method (Dickey & Lientz 1970), by using the posterior distribution of hypothesis H 1 alone, as Both methods yield similar Bayes Factors.By evaluating the evidence, we obtain B 1 0 = 0.72, while using the Savage-Dickey Density Ratio method, we find B 1 0 = 1.11.According to the significance quantification approach in Kass & Raftery (1995), when the Bayes Factor is close to unity, it indicates that the current data does not strongly favor one hypothesis over the other.

Figure 1 .
Figure1.Total phase shift |δ φ| and resonant frequency f for a 1.4 M⊙ − 1.4 M⊙ system (left panel) and a 1.8 M⊙ − 1.8 M⊙ system (right panel), depending on the transition densities (represented by pt/εt).fISCO ≃ 4400 Hz/(M + M ′ ) denotes the GW frequency corresponding to the innermost stable circular orbit, at which we assume the signal of inspiral phase shut off(Cutler & Flanagan 1994).Also plotted are the detectability thresholds of aLIGO, ET and CE, which are calculated from their designed sensitivities by assuming the system is located at DL = 100 Mpc.

Figure 2 .
Figure2.Detectable regions of phase transition parameters for a 1.4 M⊙−1.4 M⊙ system (orange) and a 1.8 M⊙−1.8 M⊙ system (blue) on the 3rd generation detectors.The red region depicts the parameter space ruled out by the mode resonance analysis of GW170817 at 95% confidence level, and the grey regions represent the portions eliminated by the maximum mass of observed pulsars.The results of QMF+CSS and NL3ωρ+CSS are represented by dotted and dashed lines, respectively.The critical line corresponds to ∆ε/εt = 1/2 + 3pt/2εt, dividing the parameter space into two parts.Below the critical line, there is a hybrid star branch connected to the hadronic branch, while above the critical line, there is no connected hybrid star branch(Alford et al. 2013).

Figure 3 .
Figure3.Constraints on phase transition densities (represented by pt/εt) from GW170817 data for QMF+CSS (upper panel) and NL3ωρ+CSS (lower panel).Colored bands depict calculations for a GW170817-like system with chirp mass M = 1.186M⊙ and mass ratio q = 0.7-1.The posterior probability density for ∆| φ| and f inferred from the GW170817 data is shown in blue heatmap, with the 95% and 68% boundaries being overlaid by dashed lines.

Figure 4 .Figure 5 .
Figure 4.An illustration of our NS EOS with a first-order phase transition to quark matter.

Figure 8 .
Figure 8. Posterior distributions of the waveform parameters inferred from GW170817 strain data.The horizontal lines in last two diagonal panels denote the prior distributions of corresponding parameters.

Table 2 .
Fit parameters of the tidal overlap integral for 1.4 M⊙ stars in the functional form of Eq. (31).

Table 2 .
It is found the estimation errors are less than ∼ 25%. Figure 7. Scaled tidal overlap integral Q versus x, together with the fitting lines for 1.4 M⊙ stars.The associated fitting coefficients can be found in Table2.

Table 3 .
The rms measurement errors of |δ φ| and f for a 1.4 M⊙ − 1.4 M⊙ system.