Parameter estimation for strong phase transitions in supranuclear matter using gravitational-wave astronomy

At supranuclear densities, explored in the core of neutron stars, a strong phase transition from hadronic matter to more exotic forms of matter might be present. To test this hypothesis, binary neutron-star mergers offer a unique possibility to probe matter at densities that we can not create in any existing terrestrial experiment. In this work, we show that, if present, strong phase transitions can have a measurable imprint on the binary neutron-star coalescence and the emitted gravitational-wave signal. We construct a new parameterization of the supranuclear equation of state that allows us to test for the existence of a strong phase transition and extract its characteristic properties purely from the gravitational-wave signal of the inspiraling neutron stars. We test our approach using a Bayesian inference study simulating 600 signals with three different equations of state and find that for current gravitational-wave detector networks already twelve events might be sufficient to verify the presence of a strong phase transition. Finally, we use our methodology to analyze GW170817 and GW190425, but do not find any indication that a strong phase transition is present at densities probed during the inspiral.


I. INTRODUCTION
Neutron Stars (NSs) are remnants of core-collapse supernovae and contain matter at the highest densities that we can observe in the Universe, up to several times nuclear saturation density, n sat = 0.16 fm −3 , which corresponds to a mass density of 2.7×10 14 g cm −3 .Hence, NSs are perfect laboratories to determine the unknown equation of state (EOS) of dense matter.The EOS relates the pressure with the energy density in the NS interior and is determined by the fundamental degrees of freedom inside the NS and their interactions among each other.Each possible EOS determines the global structure of NSs, i.e., their masses and radii, in a unique way.Therefore, detailed astronomical observations of NSs, in particular of binary NS (BNS) coalescences, are of extreme importance to nuclear physicists and allow us to constrain the densematter EOS.To date, most EOS constraints stem from NS mass measurement [1][2][3] or radius extractions from X-ray observations [4][5][6].The latter however, suffer from relatively large statistical [5] or systematic [7] uncertainties.In addition to electromagnetic (EM) observations, the remarkable observation of gravitational waves (GWs) from a BNS merger in 2017, GW170817, by Advanced LIGO [8] and Advanced Virgo [9] provided another avenue to determine NS properties and the EOS [10][11][12].Numerous efforts have been made to extract information on the EOS from the GW signal of BNS mergers, see e.g., Refs.[13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32]; cf.Ref. [33] for a recent review and further references.* thopang@nikhef.nlConstraints from GW170817 arise either purely from the analysis of the GW signal, e.g., Refs.[10,[34][35][36][37][38], from a combination of GW and EM information, e.g, Refs.[32,[39][40][41][42][43][44], or from analyses of large sets of possible EOSs constrained by nuclear-physics theory at low densities, e.g., Refs.[20,22,24,26,30,32,45].Furthermore, the NS mass distribution, and therefore also the maximum density in the NS core, is bounded from above by the NS maximum-mass configuration.This is the highest NS mass that can be supported against gravitational collapse by the dense matter in the NS interior and depends on the EOS.While this maximum mass can in principle be as high as 3−4 M (see e.g.Ref. [24]), the EM observation of the kilonova associated with GW170817 has constrained the maximum mass to be much smaller, of the order of 2.2 − 2.3 M [41,46,47].Most NS observations so far have explored NSs in a mass range of 1.4 − 2.1 M , and hence below the maximal possible density.In contrast, the coalescence of two typical NSs of approximately 1.4 M creates an object that is likely above the maximum mass, and truly explores the EOS at the highest densities in the Universe.
In the context of the dense-matter EOS, an important problem is to determine the nature of matter inside of NSs.For example, at very high energy densities the fundamental theory of strong interactions, Quantum Chromodynamics, predicts that matter undergoes a phase transition to quark matter but it is unknown at what densities such a transition occurs.A long-standing question is whether NSs explore such phase transitions to new and exotic forms of matter in their cores or whether they solely consist of nucleonic matter [48][49][50], and whether these transitions are observable [51,52].Among NSs that explore a phase transition, "twin stars" are a partic-ular family.For these stars, the phase transition in the EOS is sufficiently strong so that the mass-radius curve around the phase transition is disconnected and contains two or more branches [53,54].
For EOSs where a phase transition is present, two scenarios can be distinguished based on its onset density.First, it is possible that the transition happens at very high densities, i.e., only in heavy stars.In this case, the phase transition is probed only after the collision of the two individual stars; see, e.g., Refs.[55][56][57][58][59][60].A second possibility is that the onset density of the phase transition is at lower densities explored in typical NS around 1.4 M , that are probed already during the inspiral phase of the NS merger, e.g., Refs.[61,62].In the latter case, one could imagine scenarios in which a mass asymmetry between the two stars in the BNS leads to one (lighter) star containing only nuclear matter, while in the other (more massive) star a quark core is already present.In such a case the two individual stars could have very different radii and tidal deformabilities while their masses are comparable 1 .This is of particular importance since a number of existing GW analyses, e.g., Refs.[35,37], and multi-messenger constraints on the EOS, e.g., Refs.[39,40,[42][43][44], rely on the assumption that both inspiraling objects are NSs following some given quasi-universal relations [35,66,67].In the presence of a phase transitions those quasi-universal relations might be violated, in which case their employment in GW analyses are likely to lead to biases of the determined binary properties and the EOS.
While the current analysis of GW170817 seems to disfavor NSs with too large radii and tidal deformabilities, consistent with the appearance of phase transitions, the data from this single event is insufficient to conclusively answer this question; see, e.g., Ref. [27].Many recent works have addressed the question whether GWs allow to constrain the existence of hybrid stars, i.e., NSs that explores strong phase transitions to exotic forms of matter in their cores, and in particular twin stars [29,30,55,56,59,61,68].For example, Ref. [29] searched for the presence of a phase transition by applying quasi-universal relations.In particular, the presence of a strong phase transition was probed via observing the breakdown of quasi-universal relations.Ref. [29] found that the mass at which the phase transition occurs, M t , can be measured with 50 − 100 detections and the corresponding microscopic parameters can be estimated via quasi-universal relations.Furthermore, Refs.[27,28] looked for indications of phase transitions in GW data using a non-parametric inference approach to the EOS.By combining heavy pulsar observations, GW170817, and the recent NICER observation [5], a Bayes factor in favor of the presence of multiple stable branches of 1.8±0.2[28] was found.
When looking for the imprint of a phase transition in the GW170817 data, these previous works have mainly searched for the presence of multiple stable branches in the mass-radius relation.However, this is only one among various scenarios.In this work, we aim at quantifying whether we are capable of determining the presence of a strong phase transition from GW data even when only one stable branch is present.In particular, we ask the question how many GW observations are necessary to observe a phase transition and recover the parameters of an injected EOS from GW data.We focus on three different EOS that experience a phase transition in the typical mass range explored in BNS systems and which show 3 different behaviors in the mass-radius relation.
For this purpose, we introduce a novel method, based on a new parameterization for EOSs at supranuclear densities, of testing GW data from the inspiral phase of a BNS merger for the appearance of a strong phase transition.This new approach is based on Bayesian inference methods, and can be used with current GW detectors.Simulating 600 signals for three different EOSs, we find that already 12 events might be sufficient to confidently find the presence of a phase transition.However, when analyzing the signals GW170817 and GW190425 with our method, we do not find a hint of a strong phase transition.
The major differences between previous studies and our work are that we simultaneously (i) analyze EOSs with different phase-transition signatures, i.e., one EOS with a twin-star solution which is commonly searched for, but also two EOSs with phase transitions leading to single-branch solutions, (ii) analyze both simulated data and actual events with state-of-the-art Bayesian GW data analysis techniques, which allows for hypothesis testing and parameter estimation at once, and (iii) explicitly demonstrate that our method is able to measure the microscopic characteristics of strong phase transitions, by comparing injected with recovered parameters.Hence, our method allows us to make statistically robust statements on the presence of strong phase transitions.
The paper is structured as follows.We describe our methods and our mock data setup in Sec.II and III, respectively.Main results are shown in Sec.IV, and in Sec.V we apply our method to GW170817 and GW190425.We conclude in Sec.VI.

A. The equation of state of NS matter
The structure of NSs is completely determined by solving the Tolman-Oppenheimer-Volkoff (TOV) equations.The only necessary input is the EOS, a relation between the pressure, energy density, temperature, and composi- tion inside the NS.The EOS is determined by the microscopic degrees of freedom in the NS interior and their interactions.At lower densities, these are mostly neutrons with a few percent protons interacting via nuclear forces, but at higher densities new degrees of freedom might appear.With typical radii of the order of 12 km and masses of 1 − 2 M , the densities inside NSs are so large that thermal energies are much smaller than typical Fermi energies, except in the most violent astrophysical scenarios.Hence, for isolated NSs and NSs during the inspiral phase of a BNS merger, finite temperature effects can be neglected and the EOS is simply a relation of pressure and energy density for a given composition.
From the theoretical side, the EOS of cold dense matter at the densities explored in the NS core is very uncertain.There exists a multitude of models for the EOS which explore a wide range of pressures at densities beyond n sat , leading to large uncertainties for the radii of typical NSs.The models can differ both in the degrees of freedom that they assume and in the effective interactions among them.At low densities explored in the NS crust and outermost core, where experimental input on, e.g., saturation properties or extractions of the symmetry energy are available to constrain models, the EOS can instead be constrained rather reliably.In this density regime, approaches ranging from density functionals [69] or relativistic mean-field models [70,71] to ab initio calculations using a variety of models for the nuclear interactions, e.g., Refs.[72,73], lead to consistent results.
In recent years, important constraints on the EOS of NS matter at low densities have been obtained from microscopic calculations of neutron matter using systematic interactions from chiral effective field theory (EFT) [73][74][75][76][77][78].Chiral EFT [79,80] represents a systematic lowmomentum expansion of nuclear forces, which is connected to the symmetries of the fundamental theory of strong interactions, Quantum Chromodynamics.Instead of using quark and gluon degrees of freedom, it uses the fact that quarks are confined to hadrons at low densities and models interactions in terms of effective degrees of freedom, nucleons and pions.Chiral EFT naturally includes two-body and many-body interactions among nucleons, and provides an order-by-order scheme for the interactions among nucleons in terms of contact interactions, whose couplings are fit to experimental data, and long-range pion-exchange interactions.By going to higher orders in this expansion, calculations increase in difficulty, but results become more precise and accurate.In addition, the systematic order-by-order scheme can be used to obtain theoretical uncertainty estimates [81,82].Hence, calculations of neutron matter with chiral EFT interactions provide constraints on the EOS with reliable uncertainties.
However, chiral EFT calculations are valid only within the radius of convergence of the theory.Typically, the momentum expansion breaks down at momenta of the order of 500 − 600 MeV; see, e.g., Ref. [82] for a recent analysis.Hence, in neutron matter the chiral EFT approach might be reliable only up to (1 − 2)n sat , but most likely fails beyond that [30,78].Therefore, at densities beyond two times saturation density, currently no reliable statement about the EOS can be made from microscopic nuclear theory.At these densities, models suffer from the absence of available experimental data and our ignorance of strong interactions in this regime.In particular, while we know the relevant degrees of freedom at lower densities to be nucleons, it is not clear which de-grees of freedom appear at larger densities.While many astrophysical EOSs assume nucleonic degrees of freedom to be valid in the whole NS, a phase transition to new degrees of freedom, e.g., quark matter or exotic condensates, might occur [68].This phase transition might be strong and of first order, in which case it can lead to interesting features in the mass-radius relation, like kinks or disconnected branches [52].
To extend microscopic low-density results for the EOS to higher densities, the uncertainty in the degrees of freedom must be taken into account.This is typically done by applying general extension schemes, e.g., by using sets of polytropes for the energy density and pressure [20,73,83] or an expansion in the speed of sound [24,25,29].An extension of this approach are nonparametric inference schemes which have prior support for all possible EOS curves [84] and have recently been combined with chiral EFT calculations [30].Such extension schemes abandon explicit assumption about the degrees of freedom at higher densities but instead model all EOS curves permitted by the chosen functional form in the case of parametric extensions 2 and general physics considerations such as causality.By sampling all allowed functions, uncertainties at low densities can be systematically extended to high densities.
These general EOS sets contain both smooth EOSs, i.e., EOSs for which the change of pressure with energy density is continuous at all densities, as well as EOSs with drastic changes in the pressure.While EOSs of the first type might be obtained by using a purely nucleonic description of NSs, the latter type contains EOSs with strong first-order phase transitions.Such transitions can be modeled within a Maxwell or Gibbs construction, depending on the properties of the considered phases.
In a Maxwell construction, no mixed phases appear and the phase transitions can be modeled by an EOS segment where the speed of sound, c 2 S = ∂p/∂ = 0, vanishes.This EOS segment, and hence the phase transition, can be described by its onset density, where the speed of sound becomes 0, and its width, i.e., the density jump between the two phases with nonvanishing speed of sound.Depending on these properties, different features might be observed in the M -R relation (or the M -Λ relation).In Fig. 1, we show two examples of such phase transitions, which we have selected from the EOS set of Ref. [24].This EOS set was constrained by microscopic chiral EFT calculations below nuclear saturation density, including a consistent NS crust [85], and extended to larger densities by using a speed-of-sound extension scheme.Hence, it ensures NS stability (c S > 0) and causality (c S < c, 2 Because no particular choice for the functional form of the EOS is made in nonparametric inference schemes, they are less limited in this sense.However, in this work we choose a parameterized approach for the EOS because it is straightforward to implement a c 2 s = 0 segment that appears in case of a strong first-order phase transition. with c the speed of light) by design.For one of the two chosen EOS (labeled KINK, orange line), an intermediate width is chosen for the phase transition which leads to a visible kink in the mass-radius curve.For the other EOS (labeled TWIN, blue line), a larger width leads to a stronger phase transition which results in the appearance of two disconnected branches in the mass-radius relation, a so-called twin-star solution.For both EOS, the onset density is chosen such that the interesting feature appears already in typical NSs.
In the Gibbs construction, on the other hand, a mixed phase appears and smoothens the resulting EOS around the phase transition.In that case an EOS with a phase transition might be indistinguishable from a purely nucleonic EOS, which is known as the masquerade problem [51].We compare the two EOSs TWIN and KINK with the model ALF2 [51], which is a hybrid EOS with a phase transition to quark matter that leads to the formation of a mixed phase.For all three EOS models, the maximum mass is greater than 1.93 M and, hence, consistent with observed masses of heavy NSs [2,3].

B. Imprint of phase transitions on the GW signal
A possible phase transition can imprint itself upon the GW signal in different ways during the inspiral and during the postmerger phase.
Inspiral: During the inspiral, the GW signal depends on the properties of the two binary stars (masses, spins, and tidal deformabilities), as well as on the source location and orientation.Of particular importance for the description of tidal effects are the parameters that captures the leading-order contribution at fifth Post-Newtonian (PN) order, and with the symmetric mass ratio which captures additional contributions at sixth PN order.The presence of a phase transition in the EOS might lead to a significant change in the radii and tidal deformabilities for almost equal mass systems, if in any of the two stars the onset density for the transition is already reached in the core, cf.Fig. 1.In Fig. 2, we show the expected GW signal for a non-spinning BNS system with component masses of 1.50 M and 1.45 M for the three EOSs that we have chosen in this work.For the The GW waveforms for a non-spinning BNS system with component masses of 1.50M and 1.45M for the three EOS used in this work.The dashed lines are the waveforms assuming that the Binary-Love relation [67,87] holds.
two EOSs with strong first-order phase transitions, due to the different tidal deformabilities of both stars, we find a dephasing of the waves compared to the ALF2 EOS.Since the leading-order tidal contribution enters at the fifth PN order [63][64][65]86], the dephasing is most prominent in the late-inspiral phase.
In addition to the GWs for the EOSs described above, we also present as dashed lines the waveforms when we assume that the quasi-universal relation of Refs.[67,87] holds for TWIN, KINK, and ALF2.To apply this relation, we fix the tidal deformability of the lower-mass star and compute the tidal deformability of the primary component by using the quasi-universal relation.We find that the resulting waveform significantly deviates from the full waveform for EOSs with a very strong phase transitions, i.e., TWIN, but approximates the waveform well in the other cases.These differences suggest the failure of the quasi-universal relation with respect to EOSs with strong phase transitions like TWIN, while the relation holds approximately for the KINK EOS.For smooth EOSs, like ALF2, there is no observable difference and the quasi-universal relation seems to be valid.This suggests that detecting a phase transition that does not result in a twin-star solution will be challenging if a methodology purely based on quasi-universal relations is employed.
Postmerger: In cases in which the phase transition happens at densities beyond the ones probed during the inspiral, the postmerger signature can change if a phase transition is present, as outlined in, e.g., Refs.[29,55,56,[59][60][61]68]. In most cases, the presence of a phase transition will lead to a different (often shorter) lifetime of the remnant and a shift of the main postmerger GW emission mode.However, due to the missing sensitivity of existing GW detectors in the highfrequency range [9,88] and the absence of high-quality GW models describing the postmerger evolution of BNS mergers -see Refs.[89][90][91][92][93][94] for some first attempts -it seems natural to investigate, at the current stage, pos-sible phase transition effects that can be extracted from the GW signal during the inspiral.

C. EOS parameterization for phase transitions
When analyzing NS observations, one needs to assume an EOS describing the relation between pressure and energy density.For the "true" NS EOS realized in nature, however, the functional form is unknown.Hence, an EOS parameterization needs to be flexible enough to capture the various effects one might encounter in nature, in particular phase transitions.In this work, we consider the three EOSs of Fig. 1 as three possible "true" EOSs.In order to capture all features of these EOSs, in particular the first-order phase transitions, here we propose to use a 5-parameter piecewise-polytrope EOS parameterization scheme, which we refer to as Maxwell parameterization.This scheme is similar to the parameterization proposed in Ref. [83].
In our Maxwell parametrization, at low densities up to nuclear saturation density, we use an EOS constrained by the chiral EFT calculation of Ref. [78] (V E,1 parametrization).This EOS contains a consistent inner crust and uses the BPS model for the outer crust.For the highdensity part beyond n sat , we use a modified 5-parameter 4-piece polytrope.Each polytrope is characterized by the starting pressure p i and the adiabatic index Γ i .Therefore, our extension starts with 8 free parameters: To ensure continuity for the first polytrope, the starting pressure p 1 is chosen to be the pressure of the chiral EFT EOS at the nuclear saturation density, p CEFT (ρ 0 ).The adiabatic index of the second polytrope, Γ 2 , is set to be zero to represent a Maxwell construction for a phase transition extended across a density gap of ∆ρ.Therefore, p 2 = p 3 = p tr , where p tr is the phase transition pressure.Furthermore, we choose the transition pressure between the third and fourth polytropes to be 5 times the phase transition pressure, p 4 = 5p tr .Fixing p 4 reduces the numbers of free parameters and, therefore, helps during the recovery.The particular value of 5p tr is chosen ad-hoc by comparison against various different EOSs.We have also explored leaving this parameter free, but this only improved the fitting performance marginally.Therefore, to reduce dimensionality, we have fixed this parameter.A sketch of the parameterization is shown in Fig. 3 and its capability of representing the three EOSs is shown in Fig. 1.We find that the Maxwell parametrization works well for EOSs where phase transitions appear between 1 − 4n sat , corresponding to a NS in a typical mass range.As suggested in [95,96], we have found that our parametrization introduces systematic uncertainties due to imperfect fits to the "true" EOSs.Yet, as shown in the right panel of Fig. 1, the error is overall small and will be below the statistical uncertainty of an EOS measurement [97].However, it might be necessary to verify the generality of this assumption assuming different e e m a t + g 5 h j / y P n 8 A D p e Q t Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " M E J U G g l X g f s 3 n p 9 h z u 1 e e m a t + g 5 h j / y P n 8 A D p e Q t Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " M E J U G g l X g f s 3 n p 9 h z u 1 a z 6 a H 5 8 6 J U X p O P 1 a m I n S m 6 u + J j E i t x z I 0 n Z L g U C 9 6 E / E / r 5 N i / 9 L P e J S k y C I 6 W 9 R P h Y O x M 0 n B 6 X H F K I q x I Y Q q b m 5 1 6 J A o Q t F k V T I h e I s v L 5 P m W d V z q 9 6 t V 6 m d w w x F O I J j O A U P L q A G N 1 C H B l B I 4 R l e 4 c 1 6 t F 6 s d + t j 1 l q w 5 j M H 8 A f W 5 w 9 p / p T 0 < / l a t e x i t > < l a t e x i t s h a a z 6 a H 5 8 6 J U X p O P 1 a m I n S m 6 u + J j E i t x z I 0 n Z L g U C 9 6 E / E / r 5 N i / 9 L P e J S k y C I 6 W 9 R P h Y O x M 0 n B 6 X H F K I q x I Y Q q b m 5 1 6 J A o Q t F k V T I h e I s v L 5 P m W d V z q 9 6 t V 6 m d w w x F O I J j O A U P L q A G N 1 C H B l B I 4 R l e 4 c 1 6 t F 6 s d + t j 1 l q w 5 j M H 8 A f W 5 w 9 p / p T 0 < / l a t e x i t > < l a t e x i t s h a r s 1 6 s d + t j 1 l q w 5 j O H 8 A f W 5 w / i j J U z < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 5

a t e x i t s h a 1 _ b a s e 6 4 = " D x 1 Z W I V 4 R 5 H n d T L T n W o z g f N n y p 8 = " >
B S 2 M Q P I U Z w e U Y 8 e I x g l k g G Y a e T i d p 0 r P Q X R M M w 9 z 8 D C 8 e F P H q n 3 j z b + w s B 0 1 8 U P B 4 r 4 q q e k E i h r s 1 6 s d + t j 1 l q w 5 j O H 8 A f W 5 w / i j J U z < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 5     parametrizations.Because of the large computational cost, this is not part of this work.Furthermore, we do not expect systematics induced by our parametrization to significantly affect the present study and its main goal, namely to identify strong phase transitions and their parameters.

3 m h j i / F v M v Y 2 O E I c t o c a k l T y I = " >
In total, our parametrization is described by five parameters Γ 1 , Γ 3 , Γ 4 , the phase transition onset pressure p tr , and the transition density jump ∆ρ.In practice, in case of an EOS with a twin-star solution, it is possible to have NSs with the same mass but different radii, i.e., NSs that live on different branches of the M -R relation.These NSs have central densities around the onset density of the phase transition.Because in this case the NS mass cannot be used to distinguish the individual stars -see Fig. 1 (right panel) -we need to add two extra parameters B 1 and B 2 for each of the two NSs to indicate on which branch the star lives.Therefore, the EOS parameters E are given by where E c denotes the common parameters of all stars, assuming they follow the same EOS.

A. Bayesian Analysis
According to Bayes' theorem, the posterior p( θ|d, H) ≡ P( θ) on the parameters θ under hypothesis H and with data d is given by where L( θ), π( θ), and Z(d) are the likelihood, prior, and evidence, respectively.The prior describes our knowledge of the source or model parameters prior to the experiment or observation.The likelihood and evidence quantify how well the hypothesis describes the data with the given set of parameters and over the whole parameter space, respectively.By assuming Gaussian noise, the likelihood L( θ) that the data d is a sum of noise and a GW signal h with parameters θ is given by [98] where the inner product a|b is defined by Here, ã(f ) is the Fourier transform of a(t), * denotes complex conjugation, and S n (f ) is the one-sided power spectral density of the noise.In our study we will set f low and f high to 20 Hz and 2048 Hz, respectively.The evidence Z is given by which is the normalization constant for the posterior distribution.Moreover, we can compare the plausibilities of two hypotheses, H 1 and H 2 , by using the odd ratio, which is given by where B 1 2 and Π 1 2 are the Bayes factor and prior odds, respectively.If O 1 2 > 1, H 1 is more plausible than H 2 , and vice versa.Throughout this study, we have the prior odds set to 1, in which case the Bayes factor is the same as the odd ratio.
Within the Bayesian framework we can combine the information from multiple detections.For parameters that are expected to be the same across several detections (e.g.EOS parameters), the combined posterior for common parameters P c ( θ c ) can be obtained as where θ c are the common parameters and P i ( θ c ) is the posterior including the i-th detection.We can also combine the odds ratios into a catalog odds ratio O , which is given by where B 1 2,i is the Bayes factor for the i-th detection and B 1 (cat) 2 is the catalog Bayes factor.As we have the prior odds set to 1, the catalog odd ratio O is the same as the catalog Bayes factor B 1 (cat) 2 .

B. Waveform approximants
In this paper, we restrict our studies to the PN model TaylorF2.This model was also employed in the analysis of GW170817 and GW190425 by the LIGO and Virgo Collaborations; see e.g.Refs.[36,38,99].
The version of TaylorF2 we use is based on on a 3.5 PN order point-particle description [100][101][102][103][104] that includes spin-orbit effects [105] and spin-spin effects [106][107][108][109]. Tidal effects are added up to 7.5 PN following Refs.[86,110].We note that we also incorporate the EOS-dependency of the 2PN and 3PN spin-spin contributions.For this purpose, we use quasi-universal relations outlined in [66] to connect the spin-induced quadrupole moments to the tidal deformability of the NS stars.This approach is commonly used for GW data analysis and was by default used here, but a phase transition could also affect the employed quasi-universal relation.This might introduce additional biases.Sampling over the individual quadrupole moments of the NSs (see, e.g., Ref. [111]), would cause an increase of the computational costs and it seemed more appropriate to employ the quasi-universal relations of [66] than simply neglecting the EOS imprint on the quadrupole moment.However, since all simulated signals are non-spinning, we do not expect that any significant biases appear during our analysis and refer to the study of Ref. [112], where it was shown that potential biases only arise for high spins.

C. Injection setup
We choose to use an astrophysically motivated distribution for the parameters of the simulated sources.We distribute the sources uniformly in a co-moving volume with the optimal network SNR range ρ ∈ [30,100].Thus, a relatively high lower bound on SNR is assumed.Indeed, to probe the phase transition, an accurate measurement of Λ i is needed, which can not be achieved with BNS signals that have low or medium SNR [97].
The orientation (ι, ψ) and the sky location (α, δ) of the sources are placed uniformly on a sphere.Since NS spins are expected to be small [113], we set them to zero for all simulated sources.The component masses of the binaries are sampled from the uniform distribution [1 M , M TOV ], where M TOV is the maximum allowed mass of a NS with the given EOS.
For the EOSs, we want to investigate under what circumstances one can distinguish between the existence and the absence of a strong phase transition.Therefore, we choose two EOSs that have phase transitions with different onset pressure and density jumps but that lead to an observable and distinct feature within the mass range of [1,2]M in the mass-radius curve, labeled as TWIN and KINK; cf.Sec.II.In addition, we choose another EOS with phase transition but a smooth density dependence of the pressure.As the phase transition masquerades, this model is indistinguishable from a purely nucleonic EOS model, see Fig. 1.Here, we choose the ALF2 EOS described in Ref [50] because of its plausibility based on the multi-messenger analysis of GW170817 [39,44].The three EOSs are all able to support the observed heavy NSs [1][2][3].
The simulated GW signals are injected coherently into the data of the Advanced LIGO and Advanced Virgo detectors.The detector noise is simulated as stationary Gaussian noise with the power spectral density to be that of the design sensitivities of each detectors [9,88].

D. Implementation
Our analysis follows a similar approach as previous works [18,27,31].The analysis consist of a two-stage process: I. Estimation of the posterior of the macroscopic parameters P( θ macro ) based upon a GW analysis.
II. Estimation of the posterior of the microscopic parameters (i.e. the EOS parameters) P( E) with P( θ macro ) given.
For stage I, the posterior P( θ macro ) is estimated with the Nested Sampling algorithm [114] implemented in LALInference [115] with a prior of m i ∈ [0.5, 3.0] M and Λ i ∈ [0, 5000].For stage II, the posterior P( E) is given by where π( θ macro | E) and π( θ macro |I) are the priors on θ macro with and without the EOS given, respectively.For our study, the macroscopic parameters of interest are the component masses m 1,2 and the corresponding tidal deformabilities Λ 1,2 .Therefore the likelihood L( E) is given by where Λ(m, E) is the tidal deformability as a function of mass with an EOS given.We have also chosen the priors π(m i | E) and π(m i |I) to be the same.
Therefore, the joint posterior P( E, m i ) is given by The joint posterior is estimated with the Nested Sampling algorithm Multinest [116] implemented in PyMultinest [117,118].The posterior P( E) is then obtained via marginalizing P( E, m i ).
For the stage II process, we choose the prior for the parameters to be m i ∈ [0.5, 3.0] M , Γ i ∈ (1, 10], log 10 p tr [dyne cm −2 ] ∈ [33.7, 38.0] and log 10 ∆ρ[g cm −3 ] ∈ [13.85, 16].We also impose the constraints of M TOV ≥ 1.93 M as part of the prior.To increase efficiency, sampling over masses is done in terms of the chirp mass M and ln ∆ν rather than individual masses.The chirp mass M and ln ∆ν are given by IV. LOCATING PHASE TRANSITIONS FROM GW SIGNALS

A. Method description
Because the pressure p within a compact star is monotonically decreasing from p c in the center to p = 0 at the surface, only the part of the EOS with pressures below the central pressure p c is observable.With this in mind, we define the hypotheses to be tested as follows: • H PT : The phase transition pressure p tr is below p c for one or both of the stars and the phase transition experiences a density jump ∆ρ > 0; • H NPT : The phase transition density jump ∆ρ is zero or p tr is larger than p c , and therefore the transition is not observable.
For H NPT , we found that it is sufficient to test for the condition ∆ρ = 0. Within our parameterization, the condition p tr > p c is equivalent to fitting the whole observable EOS with a single polytrope (instead of 3-4 polytropes in the case of p tr < p c ), which is penalized by the fit.Moreover, a significant decrease of the number of fit degrees of freedom complicates our interpretation of the evidence and, correspondingly, the Bayes factors.
The evidences for the two hypotheses are given by where the central pressure p c,i is estimated via interpolation of m i -p c,i for given EOS parameters E. The Bayes factor B PT NPT between H PT and H NPT is given by By examining the Bayes factor B PT NPT , one can deduce if a phase transition is observed.for injections with the TWIN, KINK, and ALF2 EOSs.Each catalog consists of 5 events (left), 9 events (center ), or 12 events (right).The presence of a strong phase transition can be recognized from the distributions for 9 or more events.The 5σ threshold, indicated by the black dashed line, is estimated with respect to the log Bayes factor distribution for the ALF2 EOS using a Gaussian kernel density estimation (the smooth grey curve).
of the three EOSs are shown in Fig. 4. In our simulations, 85%, 98% and 70% of the injections lead to a positive ln B PT NPT for the TWIN, KINK, and ALF2 EOSs, respectively.Even though the EOSs with a strong phase transition, KINK and TWIN, do shift the distribution of ln B PT NPT towards larger values, the shift is not pronounced enough to draw a statistically robust conclusion.Furthermore, while there is no strong first-order phase transition in the ALF2 EOS, H PT is favored over H NPT for the ALF2 injections.However, in H PT , the EOS below p c is fitted with 3 − 4 polytropes while in H NPT it is fitted with only 1 − 3 polytropes.The additional degrees of freedom for H PT improve the fit to the M -Λ curve and lead to a higher evidence.
To improve on the situation, we follow the catalog technique described in Sec.III A, and estimate the catalog log Bayes factors ln B PT (cat) NPT .In Fig. 5, we show the distributions of ln B PT (cat) NPT with 5, 9, and 12 events per catalog for the three EOSs.We find that the presence of a strong phase transition in the EOS can be clearly recognized from the distributions for 9 or more events per catalog.For both the TWIN and KINK EOSs, with 12 events per catalog, all the catalogs result in a higher than 5σ statistical significance with respect to the catalogs estimated for the ALF2 EOS.
Before continuing, we note that by combining the Bayes factors with simple multiplication we are not making use of the knowledge that all BNSs share the same EOS.
The inclusion of this additional constraint would require us to do the analysis on all BNSs in a catalog simultaneously, which would be computationally very demanding.While our procedure is sub-optimal, it is a conservative one.Most importantly, we see that by using our catalog Bayes factor as a detection statistic for strong phase transitions, with 12 sources we can already draw significant conclusions, regardless of the interpretation of the Bayes factor.
Returning to Fig. 5, we see that the phase transition in the KINK EOS is easier to identify than the TWIN EOS, even though the TWIN EOS has the more pronounced feature in the M -R relation.Moreover, in the case of TWIN, the parameters B 1 and B 2 , which indicate on which branch the stars live, are not redundant, and one would think that this should boost the evidence in that case.However, the individual Λ i for each component in a BNS are not well measured4 ; as a result, the inclusion of the {B i } does not necessarily lead to the evidence being elevated.In addition, as seen in Fig. 1, our parameterization can fit the KINK EOS better than the TWIN EOS, which contributes to the higher evidence of KINK compared to TWIN.
In addition to observing the presence of a phase transition with statistical significance, we can also probe its characteristics.In particular, the phase transition onset pressure p tr , and the phase transition density jump ∆ρ, can be measured.In Fig. 6, we show the maximum a posteriori (MAP) and 95% credible interval evolution for log 10 p tr and log 10 ∆ρ with an increasing number of events, where we only include events with positive ln B PT NPT for TWIN and KINK.Indeed, we expect that events with positive ln B PT NPT will tend to have p c > p tr , and hence be informative for the estimation of log 10 p tr and log 10 ∆ρ; this is something we will return to momentarily.Looking at the results, for both EOSs the true phase transition parameters are recovered within the 95%  credible interval.For log 10 ∆ρ, the true value is recovered after ∼ 10 events are included.The statistical errors decrease as more events are combined, though the rate of decrease appears to slow down after ∼ 10 events.
In Fig. 7, we show the joint posteriors for log 10 p tr and log 10 ∆ρ with 25 combined events for the TWIN and KINK EOSs.For both EOSs, the phase transition parameters (log 10 p tr , log 10 ∆ρ) are measured with 10% statistical uncertainty, with the true values lying within the 95% credible interval.
Finally, for the same 25 events, we show the distribution of the central pressure p c for both stars in Fig. 8.For the majority of the events, both stars have a central pressure above the phase transition onset pressure.This matches our expectation from the positive ln B PT NPT and provides a valuable crosscheck for our analysis.In particular, none of the events have both components' central pressure below phase transition pressure.
Based on our findings, we conclude that: (i) It is possible to establish the presence of a strong first-order phase transition with twelve BNS observations; (ii) With ∼ 10 BNSs, the phase transition parameters (log 10 p tr , log 10 ∆ρ) can be measured with 10% statistical uncertainty.
Because we are imposing a strict bound on the M TOV during the stage II analysis, systematics might be induced as suggested in Ref. [28,119].However, since we are not interested in recovering the full EOS or M TOV but the parameters of the phase transition, the simple hard cut on M TOV is sufficient for the purposes of this analysis.Also, no significant systematics are observed with respect to the simulation.

C. Limitations of our analysis
As we have shown in the previous section, it is possible to confirm the existence of a phase transition and extract its parameters.However, our approach is limited to such EOSs where the phase transition is Maxwell-like, i.e., it is described by a segment with c S = 0.As the comparison with the ALF2 EOS clearly shows, our approach cannot establish the existence of a phase transition in case a mixed phase appears, that smears out an observable EOS feature.Due to the masquerade problem [50], macroscopic structure properties (M -R or M -Λ relations) for such EOS cannot be distinguished from purely nucleonic EOS.As only such properties affect the GW waveforms, the inspiral phase does not provide information on the phase transition in this case.Should such a case be realized in nature, information might only be obtained from the postmerger GW signal.
Furthermore, we can only observe a Maxwell-like transition that is strong enough to leave an observable feature, e.g., at least a kink, in the M -Λ curve.Should the density jump be too small, the M -Λ curve would be smooth and, again, the inspiral phase would not provide information on the phase transition.
Finally, our method only works if inspiraling NSs have central pressures above the onset of the phase transition.If NS masses in binaries are limited to be around 1.4 M , exploring lower central pressures, but phase transitions appear at much higher pressures in heavier stars, it will only be probed by the postmerger GW signal.However, the observation of GW190425 shows that also BNS mergers of heavy NS might be observed by GW interferometers (note that GW190425 could potentially have been a neutron-star-black-hole merger [120][121][122]).Should the phase transition appear at low pressures/densities, in NS below 1 M , such that both NS in a binary are hybrid stars, our method might also not be able to identify its presence.In this case, while the observed macroscopic NS observables represent integrals over the EOS at all densities in the star, and, hence, in principle include information of the phase transition, observations might not be able to distinguish between the "true" EOS and an EOS that connects smoothly between the low-density nuclear-physics constraints and the observed part of the M -R curve.However, in this case the phase transition onset would likely be between (1 − 2)n sat , at low energy densities, which might be identified using other analysis techniques [30].
These shortcomings highlight that GW observations alone might not be able to answer the question of whether phase transitions exists in NSs or not.Hence, interdisciplinary studies including both nuclear physics and GW astrophysics are crucial.

V. ANALYSIS OF GW170817 AND GW190425
Having shown with simulations that our methodology is in principle capable of uncovering and characterizing strong phase transitions, we now apply it to the real signals GW170817 [10] and GW190425 [123].The former can confidently be assumed to have come from a BNS inspiral.We will also assume that the latter was emitted by a BNS.For both events, we take the publicly available posterior samples [124,125] as the input for the stage II analysis as described in Sec.III D.
Table I shows the log Bayes factors ln B PT NPT and the Kullback-Leibler (KL) divergence for the two events separately and combined.The KL divergence is estimated from the posterior and prior for the EOS parameters, i.e., KL divergence = P( E c ) ln and it quantifies to what extent the posterior distribution is different from the prior distribution.A KL divergence of zero indicates that the two distributions are identical.
Based on the values of KL divergence for the two events, it would seem that GW170817 is carrying more information regarding the EOS, while GW190425 is not very informative, similar to the findings of Ref. [28].
In order to make a statistically robust statement, we would need to have a reliable distribution of ln B PT (cat) NPT in the absence of strong phase transitions, which we could use as "background" to estimate the significance of the "foreground" values in Table I.This would require (i) an accurate or at least representative simulated BNS population, and (ii) a justifiable representation of EOSs without a strong phase transition.One can systematically generate a representative ensemble of EOSs by using parameteric or non-parameteric methods and used them to calculate a background distribution for ln B

PT (cat) NPT
, but the required computational resources will be significant.Due to the uncertainty in the BNS population and the limitation of available computational resources, such an estimation is currently not achievable.Instead, we follow the interpretation of the Bayes factor described in Ref. [126], and find no strong evidence for or against the presence of a strong phase transition.
In addition to the Bayes factors, a small value for the KL divergence indicates similarity of posterior and prior, and the extracted phase transition parameters are strongly influenced by the priors.The combined posteriors for the phase transition onset density, n tr , and the corresponding density jump in terms of n sat , ∆n, are shown together with the prior in Fig. 9. 5 The posteriors distributions for log 10 p tr , log 10 ∆ρ, m 1 , m 2 , and Λ for GW170817 and GW190425 are shown in Figs. 10 and 11, respectively.We conclude that our measurements of n tr and n sat are not very informative, as could be expected based on the KL divergence. 5Initially we have priors ntr ∈ [1,4]

VI. CONCLUSION
In this paper, we have presented a reliable way of searching for phase transitions in supranuclear matter using the inspiral waveform of binary NS mergers, which is successful both if the resulting mass-radius relation has only one or multiple stable branches.In contrast to previous works, which calculated the preference for multiple stable branches to search for strong phase transitions [27,28,30], our approach searches for an extended segment with c S = 0, i.e., we do not explicitly assume a multiple-branch feature in the M -R curve.
As long as there is some observable imprint of the phase transition in the mass-radius relation, our method can recover injected phase-transition parameters, and, hence, represents an important step forward in the search for a possible phase transition.We have explicitly demonstrated this by injecting simulated BNS mergers with different equations of state into synthetic stationary Gaussian noise.We have shown that our method can detect the presence of phase transitions at 5σ confidence with 12 signals.Moreover, the phase-transition onset pressure and the corresponding density jump (log 10 p tr , log 10 ∆ρ) were recovered with 10% statistical uncertainty with ∼10 events.Finally, we have applied the method to GW170817 and GW190425, but found no strong evidence for or against the presence of strong phase transitions.
FIG.2.The GW waveforms for a non-spinning BNS system with component masses of 1.50M and 1.45M for the three EOS used in this work.The dashed lines are the waveforms assuming that the Binary-Love relation[67,87] holds.
z o P z 7 L w 7 H / P S n L P o O Y I / c j 5 / A B M j k L g = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " y u s I D O 7 6 b z G / 8 7 Z I + h A + w V s c W t 4 = " > A A A B 7 3 i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g F e 5 E 1 D J g o W U E E w P J E e Y 2 m 2 T J 7 t 6 5 u y e E I 3 / C x k I R W / + O n f / G T X K F J j 4 Y e L w 3 w 8 y 8 K B H c W N / / 9 g o r q 2 v r G 8 X N 0 t b 2 z u 5 e e f + g a e J U U 9 a g s Y h 1 K 0 L D B F e s Y b k y v 8 O Y 9 e i / e u / c x b y 1 4 + c w h / I H 3 + Q N + 9 Y + N < / l a t e x i t > ⇢ < l a t e x i t s h a 1 _ b a s e 6 4 = " w e 4 l o i s l k n z P k E u 9 + E u a 7 Q D D 4 6 M = " > A A A B 8 X i c b Z C 7 S g N B F I b P x l u M t 0 R L m 8 E g W I V d C 7 U M a G E Z w V w w u 4 T Z y d l k y O z s M j M r h C W P Y G d j o Y i t b 2 P n 2 z i 5 F J r 4 w 8 D H / 5 / D n H P C V H B t X P f b K a y t b 2 x u F b d L O 7 t 7 + w f l y m F L J 5 l i 2 G S J S F Q n p B o F l 9 g 0 3 A j s p A p p H A p s h 6 P r a d 5 + R K V 5 I u / N O M U g p g P J I 8 6 o s d a D f 4 P C U F 8 N k 1 6 5 6 t b c m c g q e A u o 1 i t P F b B q 9 M p f f j 9 h W Y z S M E G 1 7 n p u a o K c K s O Z w E n J z z S m l I 3 o A L s W J Y 1 R B / l s 4 g k 5 t U 6 f R I m y T x o y c 3 9 3 5 D T W e h y H t j K m Z q i X s 6 n 5 X 9 b N T H Q V 5 F y m m U H J 5 h 9 F m S A m I d P 1 S Z 8 r Z E a M L V C m u J 2 V s C F V l B l 7 p J I 9 g r e 8 8 i q 0 z m u e W / P u v G r 9 A u Y q w j G c w B l 4 c A l 1 u I U G N I G B h G d 4 h T d H O y / O u / M x L y 0 4 i 5 4 j + C P n 8 w e 9 E p G p < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " n p E g v 7 Q K j U U c Y r 3 6 H 6 z D k 5 b E e P g = " e 1 M C I 7 M s j c T / / N 6 G c b X Y c 5 l m i G T d L E o z o S L y p 2 9 7 w 6 4 Z h T F x B J C N b e 3 u n R E N K F o Q 6 r Y E P z l l 1 d J + 6 L u e 3 X / 3 q 8 1 L o s 4 y n A C p 3 A O P l x B A + 6 g C S 2 g I O E Z X u H N M c 6 L 8 + 5 8 L F p L T j F z D H / g f P 4 A e X 2 Q u g = = < / l a t e x i t > ⇢ < l a t e x i t s h a 1 _ b a s e 6 4 = " w e 4 l o i s l k n z P k E u 9 + E u a 7 Q D D 4 6 M = " > A A A B 8 X i c b Z C 7 S g N B F I b P x l u M t 0 R L m 8 E g W I V d C 7 U M a G E Z w V w w u 4 T Z y d l k y O z s M j M r h C W P Y G d j o Y i t b 2 P n 2 z i 5 F J r 4 w 8 D H / 5 / D n H P C V H B t X P f b K a y t b 2 x u F b d L O 7 t 7 + w f l y m F L J 5 l i 2 G S J S F Q n p B o F l 9 g 0 3 A j s p A p p H A p s h 6 P r a d 5 + R K V 5 I u / N O M U g p g P J I 8 6 o s d a D f 4 P C U F 8 N k 1 6 5 6 t b c m c g q e A u o 1 i t P F b B q 9 M p f f j 9 h W Y z S M E G 1 7 n p u a o K c K s O Z w E n J z z S m l I 3 o A L s W J Y 1 R B / l s 4 g k 5 t U 6 f R I m y T x o y c 3 9 3 5 D T W e h y H t j K m Z q i X s 6 n 5 X 9 b N T H Q V 5 F y m m U H J 5 h 9 F m S A m I d P 1 S Z 8 r Z E a M L V C m u J 2 V s C F V l B l 7 p J I 9 g r e 8 8 i q 0 z m u e W / P u v G r 9 A u Y q w j G c w B l 4 c A l 1 u I U G N I G B h G d 4 h T d H O y / O u / M x L y 0 4 i 5 4 j + C P n 8 w e 9 E p G p < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " n p E g v 7 Q K j U U c Y r 3 6 H 6 z D k 5 b E e P g = " e 1 M C I 7 M s j c T / / N 6 G c b X Y c 5 l m i G T d L E o z o S L y p 2 9 7 w 6 4 Z h T F x B J C N b e 3 u n R E N K F o Q 6 r Y E P z l l 1 d J + 6 L u e 3 X / 3 q 8 1 L o s 4 y n A C p 3 A O P l x B A + 6 g C S 2 g I O E Z X u H N M c 6 L 8 + 5 8 L F p L T j F z D H / g f P 4 A e X 2 Q u g = = < / l a t e x i t > p tr < l a t e x i t s h a 1 _ b a s e 6 4 = " T u C j R 3 z 7 I 7 r u t f u 7 O S r q p t Q S 2 T w = " > A A A B + H i c b V B N S 8 N A E J 3 U r 1 o / G v U i e A k W w V N J P K j H i h e P F e w H t C F s t t t 2 6 W 4 S d i d i D b n 5 L 7 x 4 U M S r P 8 W b / 8 b t x 0 F b H w w 8 3 p t h Z l 6 Y C K 7 R d b + t w s r q 2 v p G c b O 0 t b 2 z W 7 b 3 9 p s 6 T h 1 _ b a s e 6 4 = " D + p b P 1 1 W 8 V 4 D T G v A O Z P L p n I x 6 7 k = " > A A A B + H i c b V B N S 8 N A E N 3 4 W e t H o x 6 9 L B b B U 0 k 8 q M e C F 4 8 V 7 A e 0 I W y 2 m 3 b p 7 i b s T s Q a 8 k u 8 e F D E q z / F m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 K B X c g O d 9 O 2 v r G 5 t b 2 5 W d 6 u 7 e / k H N P T z q m 6 M 1 5 c l 6 c d + d j 0 b r m l D P H 6 A + c z x / o m p P X < / l a t e x i t > p tr < l a t e x i t s h a 1 _ b a s e 6 4 = " T u C j R 3 z 7 I 7 r u t f u 7 O S r q p t Q S 2 T w = " > A A A B + H i c b V B N S 8 N A E J 3 U r 1 o / G v U i e A k W w V N J P K j H i h e P F e w H t C F s t t t 2 6 W 4 S d i d i D b n 5 L 7 x 4 U M S r P 8 W b / 8 b t x 0 F b H w w 8 3 p t h Z l 6 Y C K 7 R d b + t w s r q 2 v p G c b O 0 t b 2 z W 7 b 3 9 p s 6 T h 1 _ b a s e 6 4 = " D + p b P 1 1 W 8 V 4 D T G v A O Z P L p n I x 6 7 k = " > A A A B + H i c b V B N S 8 N A E N 3 4 W e t H o x 6 9 L B b B U 0 k 8 q M e C F 4 8 V 7 A e 0 I W y 2 m 3 b p 7 i b s T s Q a 8 k u 8 e F D E q z / F m / / G b Z u D t j 4 Y e L w 3 w 8 y 8 K B X c g O d 9 O 2 v r G 5 t b 2 5 W d 6 u 7 e / k H N P T z q m 6 M 1 5 c l 6 c d + d j 0 b r m l D P H 6 A + c z x / o m p P X < / l a t e x i t > 5p tr < l a t e x i t s h a 1 _ b a s e 6 4 = " D x 1 Z W I V 4 R 5 H n d T L T n W o z g f N n y p 8 = " > A A A B + X i c b V D J S g N B E K 2 J W 4 z b q B f B S 2 M Q P I U Z w e U Y 8 e I x g l k g G Y a e T i d p 0 r P Q X R M M w 9 z 8 D C 8 e F P H q n 3 j z b + w s B 0 1 8 U P B 4 r 4 q q r x 3 7 2 N e u u Y t e o 7 h j 7 z P H 4 T / j 8 o = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " e 9N R Y W a N Z Y r a p m C d r L E K J y L 9 p U 4 = " > A A A B 7 X i c b V A 9 S w N B E J 3 z M 8 a v q K X N Y h C s w p 2 F W g Z s L C O Y D 0 i O s L f Z S9 b s 7 R 6 7 c 0 I 4 8 h 9 s L B S x 9 f / Y + W / c J F d o 4 o O B x 3 s z z M y L U i k s + v 6 3 t 7 a + s b m 1 X d o p 7 + 7 t H x x W j o 5 b V m e G 8 S b T U p t O R C 2 X Q v E m C p S 8 k x p O k 0 j y d j S + n f n t J 2 6 s 0 O o B J y k P E z p U I h a M o p N a P T P S f b 9 A 1 1 u I M G N I H B I z z D K 7 x 5 2 n v x 3 r 2 P R e u a V 8 y c w B 9 4 n z 9 B a o 7 b < / l a t e x i t > ⇢ 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " A M j J N v o j c 0 l j L p n D I F s F 4 q z x Y I M = " > A A A B 7 X i c b Z D L S g M x F I b P e K 3 1 1 u r S T b A I r s q M r x 3 7 2 N e u u Y t e o 7 h j 7 z P H 4 T / j 8 o = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " e 9N R Y W a N Z Y r a p m C d r L E K J y L 9 p U 4 = " > A A A B 7 X i c b V A 9 S w N B E J 3 z M 8 a v q K X N Y h C s w p 2 F W g Z s L C O Y D 0 i O s L f Z S9 b s 7 R 6 7 c 0 I 4 8 h 9 s L B S x 9 f / Y + W / c J F d o 4 o O B x 3 s z z M y L U i k s + v 6 3 t 7 a + s b m 1 X d o p 7 + 7 t H x x W j o 5 b V m e G 8 S b T U p t O R C 2 X Q v E m C p S 8 k x p O k 0 j y d j S + n f n t J 2 6 s 0 O o B J y k P E z p U I h a M o p N a P T P S f b 9

FIG. 4 .
FIG. 4. The distribution of ln B PTNPT for injections with the TWIN, KINK, and ALF2 EOSs.The presence of a strong phase transition does shift the distribution of ln B PT NPT towards larger values.

FIG. 6 .
FIG. 6. Maximum a posteriori (MAP) and 95% credible interval evolution for log 10 ptr (left panel ) and log 10 ∆ρ (right panel ) for the TWIN and KINK EOSs.The dashed lines indicate the true values, which are found to lie within the extracted 95% credible interval for both EOS.

3 4 . 4 3 5 FIG. 7 .
FIG. 7. The joint posteriors for log 10 ptr and log 10 ∆ρ with 25 combined events for the TWIN EOS (left panel ) and the KINK EOS (right panel ).The dashed lines indicate the 95% credible interval and the solid lines the true value.For both EOSs, the true values of the phase transition parameters (log 10 ptr, log 10 ∆ρ) are constrained with 10% statistical uncertainty.

34.50 34 .FIG. 8 .
FIG.8.The distribution of the central pressure pc for the 25 events included in the joint parameter-estimation analysis with the phase transition pressure for the two EOSs indicated by the dashed line.Most of the events have both of the components' central pressure above the phase transition pressure.In particular, none of the events have both components' central pressure below phase transition pressure.

FIG. 9 .
FIG.9.Corner plots showing the posterior distribution of ntr and ∆n in terms of nsat with GW170817 and GW190425 combined (blue) on top of the prior with heavy pulsars constraint (grey).The dashed lines mark the 95% credible intervals.The posterior does not significantly deviate from the prior.
. 1. Density-pressure relations (left panel ) and mass-tidal-deformability relations of the three EOSs (solid) used in this work, and the least-squares Maxwell fits (dashed) to the the EOSs.The Maxwell parameterization successfully captures all features, including the phase transition, for both the EOS and the mass-radius curve.Moreover, both the EOSs and the EOSs' best-fits support heavy NSs. FIG

TABLE I .
The log Bayes factor ln B PT NPT and the KL divergence estimated with the two BNS event.