Was GW170817 a canonical neutron star merger? Bayesian analysis with a third family of compact stars

We investigate the possibility that GW170817 has not been the merger of two conventional neutron stars (NS) but that at least one of them was a hybrid star with a quark matter core, possibly belonging to a third family of compact stars. We present a Bayesian analysis method for selecting the most probable equation of state (EoS) under a set of constraints, which include besides the maximum mass also the tidal deformability from GW170817 and the first mass and radius measurement by the NICER Collaboration for PSR J0030+0451. We apply this method for the first time to a two-parameter family of hybrid EoS based on the DD2 model with nucleonic excluded volume for hadronic matter and the color superconducting generalized nlNJL model for quark matter. The model has a variable onset density for deconfinement and can mimic effects of pasta phases with the possibility of producing a third family of hybrid stars in the mass-radius ($M-R$) diagram. The main findings of this study are that: 1) the presence of multiple configurations for a given mass corresponds to a set of disconnected lines in the $\Lambda_1-\Lambda_2$ diagram of tidal deformabilities for binary mergers, so that merger events from the same mass range may result in a probability landscape with different peak positions; 2) the Bayesian analysis with the above observational constraints favors an early onset of the deconfinement transition, at masses of $M_{\rm onset}\approx 0.5~M_\odot$ with a $M-R$ relationship that in the range of observed neutron star masses is almost indistinguishable from that of a soft hadronic APR EoS; 3) a yet fictitious measurement of the NICER experiment with a $1\sigma$ range that is half of the present value and a larger mass (within the present $1\sigma$ range) would change the posterior likelihood so that a phase transition onset at $M_{\rm onset} = 1.6~M_\odot $ would be favored.


Introduction
The detection of gravitational waves from the merger of two compact stars by the state-of-the-art interferometers on Earth has opened new possibilities of studying the equation of state (EoS) of matter under extreme conditions [1]. The process of fusion of compact stars releases gravitational wave signals, which for the event GW170817 could be detected from the inspiralling phase when they probe the tidal deformability of the merging stars and allow to derive constraints on their radii [2]. The gravitational waves from the postmerger phase could not be detected for this event. But once they will become accessible to observation, the typical peak frequency in the Fourier spectrum of their ringdown signal will reveal information about the compactness of the hypermassive star that forms as a result of the fusion. It has been suggested that both signals together would allow to identify a first-order phase transition in the neutron star merger [3]. If deconfined quark matter can occur in NS cores, then the fascinating question to be settled with NS observations is that for the onset mass M onset of the deconfinement phase transition and whether the quark-core hybrid stars could form a separate family of compact stars. It is known since half a century that a third family of compact stars in the M − R diagram would serve as an indicator of a strong phase transition in neutron-star matter [4].
In this work we consider a class of hybrid star EoS based on state-of-the-art microphysical approaches to hadronic and quark matter phases which are joined by an interpolating mixed phase construction that mimics the effects of pasta phases in the transition. This theoretical approach is capable of describing hybrid stars as a third family of compact stars with varying M onset and with the possibility to remove the gap in the M − R diagram between this third family and the second family of canonical NS. In order to address the above-mentioned goal of NS merger research, we compute the mass-radius and tidal deformability relations for our class of hybrid EoS and perform a Bayesian analysis that takes into consideration the recent constraints from astrophysical observations.
The main important constraint for the EoS is the lower limit on the maximum mass of a compact star which comes from the precise measurement of the mass for the most massive pulsars, like PSR J0740+6220 [5] and PSR J0348+0432 [6]. Fortunately, up to now there are published data from two binary neutron star merger events denoted as GW170817 [2] and GW190425 [7]. Since for the latter no measurement of the tidal deformability has been reported, we shall not include it into our study. Furthermore, the recent X-ray timing observation of the isolated millisecond pulsar PSR J0030+0451 performed by the NICER experiment has provided a simultaneous measurement of both mass and radius for that object [8,9].
Before describing our approach and its results in detail, we want to give a short introduction that relates our work to the actual status of research in this field.

EoS constraints from M − R measurements before GW170817
Bayesian methods have proven to be a powerful tool for parameter estimation of a selected model using empirical data as well as for statistical inference. In the context of the determination of the neutron star equation of state, seminal works before the multi-messenger era include the analysis of X-ray bursters [10,11] which unfortunately had to deal with uncertainties related to the atmospheric composition of the associated neutron star. Recent progress has been made using Bayesian methods for the mass and radius determination data of the Rossi X-ray Timing Explorer for some accretion-powered millisecond X-ray pulsars like, e.g., 4U 1702-429 [12,13], SAX J1808.4-3658 [14] and others.
On the other hand, other Bayesian studies include the ones of Raithel et al. [15] which introduced a popular piecewise multipolytrope EoS parametrization together with mass and radius measurements of several objects, whereas [16] included realistic EoS based on the hadronic DD2 and NJL quark matter approaches to estimate the model parameters not only on already performed mass measurements but also considering fictitious radius measurements on high-mass pulsars which could provide evidence for mass-twin stars. In addition, Lackey and Wade [17] performed a Bayesian inference study based on future detections of gravitation radiation from populations of inspiralling neutron stars using a piecewise polytrope EoS model with four polytropes to conclude that the EOS determination above nuclear density could be significantly improved.

New constraints on M − R relations in the era of multi-messenger astronomy
Following both the detection of gravitational radiation from the GW170817 event [1] and the mass and radius measurement of the millisecond pulsar PSR J0030+0451 by NICER [8,9], a few works employing the Bayesian techniques appeared. Since these new measurements are more precise and thus more reliable than previous ones, they allow for tighter constraints on the compact star EoS. Consequently, we base our present work on a study carried out in [18] where a Bayesian analysis considers EoS models featuring effects of geometrical structures (pasta phases) at the hadron-quark interface on mass twins by using as inputs the tidal deformability estimates from GW170817 besides mass-radius measurements. This Bayesian analysis has recently been extended by including the NICER mass and radius measurement for PSR J0030+0451 [8,9] in a study that used a purely hadronic EoS [19]. Regarding mass twins, a few other works have also recently studied their properties by varying model parameters with the result of finding disconnected branch structures in their deformability plots for compact star mergers [20,21]. Christian and Schaffner-Bielich [22] use a hadronic relativistic mean field EoS together with the simple constant speed of sound (CSS) model for the quark matter EoS to rule out strong phase transitions below 1.7n 0 within their study. This strong claim should, however, be relaxed because they based it on the 1σ confidence level only and neglected the fact that the mass and radius measurements are correlated. K. Chatziioannou and S. Han [23] use two different hadronic EoS, namely DBHF and SFHo together with the same CSS model EoS to find out by Bayesian methods that the onset mass and the strength of a strong first order phase transition can be constrained with 50 -100 neutron star merger detections.
Miller et al. [24] base their study on a multi-polytrope EoS to asses the impact of the improvement of the precision on the determination of the nuclear the symmetry energy in the laboratory for the determination of the compact star EoS, particularly below nuclear saturation densities. An important aspect that was first implemented in [16] and is stressed in [24] is the need to use the full posteriors of neutron star masses. The implementation of a hard lower bound on the maximum star mass leads to the loss of important information since this value changes with every new observation albeit the uncertainties inherent of the measurements themselves, see [25] which is also based on a multi-polytrope parameterization. Moreover, a follow up work by Miller et al. [8] extends the EoS to include first order phase transitions in the multi-polytrope approach.
Found among the several state-of-the-art studies is the work of Most et al.
[26] that considers a set of more than two million of EoS models with and without phase transitions fitted to the multi-polytrope parametrization. They impose constraints on the upper bound for the maximum compact star mass to be less than 2.16 M as well on the tidal deformability to be less than 800. The resulting R value for a 1.4M pure hadronic star lies in the range 12.00 km < R 1.4 < 13.45 km and for a hybrid twin star within 8.53 km < R 1.4 < 13.74 km at the 2σ confidence level. Interestingly, other recent works find consistent results. The follow up study of Raaijmakers et al. [27] features two EoS classes: a five-polytrope piecewise EoS and a piecewiese constant speed of sound EoS together with the measurements of the GW170817 event, the NICER data for PSR J0030+0451, the most massive NS of about 2.14M , and chiral EFT constraints. Capano et al. [28] used a speed of sound parameterization based in chiral EFT constraints as well but different priors to the ones used by Raaijmakers et al. resulting in skewed radius results in lower values for the radii 9.2 km < R 1.4 < 12.3 km, similar to the value found by Lim et al. [29] who base their studies on empirical nuclear observables. Essick et al. [30] used a non-parametric EoS Bayesian inference under different priors, also finding consistent results. The LIGO-Virgo collaboration has presented a study that incorporates realistic compact star equations of state under the framework of long-lived rotation remnants following the compact star merger [31]. The result is that in the case of a slowly rotating remnant (less than 1.9 kHz), the baryonic mass of the static compact star turns out to be 3.5M , allowing for ruling out three of the EoS they have considered. Last but not least, in Ref. [32] for the first time the results from X-ray observations of accretion-powered pulsars in low-mass X-ray binaries are included in the Bayesian analysis together with the LVC results for GW170817 and the NICER results for PSR J0030+0451 to find based on a relativistic mean field EoS that R 1.4 ∼ 12 km. These authors report that in a version of the model that includes Λ hyperons the radius should be R 1.4 ∼ 14 km, at tension with the constraints from GW170817.

Goals of the present work
We investigate the question of how likely it may be that GW170817 has not been the merger of two canonical neutron stars but involved at least one if not two hybrid stars with a quark matter core. To this end we develop a Bayesian analysis method for selecting the most probable equation of state (EoS) under a set of constraints from compact star physics, which now include the tidal deformability from GW170817 and the first result for the mass and radius determination for PSR J0030+0451 by the NICER Collaboration. We apply this method for the first time to a two-parameter family of hybrid EoS that is based on realistic models for hadronic matter (DD2 with nucleonic excluded volume) and color superconducting quark matter (generalized nlNJL model) which produce a third family of hybrid stars in the M − R diagram. One parameter (µ < ) characterizes the chemical potential for the onset of quark deconfinement while the other (∆ P ) belongs to the mixed phase construction that mimics the thermodynamics of pasta structures and includes the Maxwell construction as a limiting case for ∆ P = 0. This paper is organised as follows. The next section contains a description of the equation of state followed by a section describing the methods of computing the neutron star configurations and tidal deformabilities. Then follows a section describing our Bayesian analysis that goes beyond our previous work [18] by now including the first results of the mass and radius determination for the millisecond pulsar PSR J0030+0451 of the NICER team [8], including the results of for the actual measurements as well as a fictitious result as one possible outcome of a forthcoming NICER analysis. In the final section we present our conclusions and an outline of ongoing research.

Quark-hadron hybrid equation of state
The hybrid equation of state is the basis of our work and consists of three ingredients. The Sly4 EoS [50][51][52] is employed for the crust in the subsaturation region and matched to the density-dependent relativistic mean field EoS "DD2" [53] describing the properties of atomic nuclei and nucleonic matter around the saturation density with stiffening effects of a nucleonic excluded volume of 4 fm 3 according to Typel [54] acting at supersaturation densities only. The nomenclature for this EoS is "DD2p40". For the high-density phase of matter we use the microscopic quark matter EoS which is matched to the hadronic phase by a mixed phase construction, both being described in the following subsections.

Quark matter model
The quark matter approach to the EoS consists of the nonlocal NJL model [55] for which a generalization was introduced in [44] which consisted in adopting a chemical potential dependence of the vector mean field coupling η(µ) and a chemical potential-dependent bag pressure B(µ) by where the density-dependent vector coupling is given by (2) and the density-dependent bag function is Here, two smooth switch functions were introduced, one that changes from one to zero at a lower chemical potential µ < with a width Γ < , and one that switches off at a higher chemical potential µ with a width Γ , together with the complementary functions that serve to switch-on, With these switch functions the density-dependence is introduced in such a way that the additional parameters have an intuitively clear meaning. The vector coupling switches from a value η < at low chemical potentials to η > at high ones, whereby the switch is positioned at µ and extends over a region of the width Γ . For the introduction of the chemical-potential dependent vector coupling to obtain heavy hybrid stars, see also [33]. This density dependence of the vector coupling facilitates the density-dependent stiffening of quark matter as a prerequisite for obtaining a third family of compact stars and the directly related mass-twin phenomenon. For this, two conditions must be simultaneously fulfilled: 1) a strong deconfinement transition with a sufficiently large jump of the energy density so that it induces a gravitational instability (requiring a relatively soft quark matter EoS at low densities), and 2) a high maximum mass of the hybrid star sequence to fulfill the 2 M constraint (requiring a stiff quark matter EoS at higher densities). We note that an alternative to the introduction of a density dependent vector meson coupling would be the use of a higher order vector interaction like the 8-quark interaction of Ref. [56] with constant couplings for which it has been demonstrated in [34] that both above conditions can be fulfilled. The bag pressure switches in a region of width Γ < around µ = µ < from a maximal value at low chemical potentials to zero at high ones, with a slight modification by the factor f (µ). This bag function and its behavior are motivated in the following way. In dynamical chiral quark models like the nlNJL, there is a mean-field pressure associated to the dynamical chiral symmetry breaking by the appearance of a chiral (quark) condensate. It is well known that this contribution to the thermodynamics of the quark model can be interpreted as a bag pressure (see, e.g., Refs. [57,58]) with a chemical-potential dependence that has been extracted in [59] within a color superconducting nonlocal NJL model under neutron star conditions. While this effect of the melting of the quark condensate at increasing density is inherent in the dynamical quark model, one can expect that a similar contribution to the thermodynamics should come from the melting of the gluon condensate. As the gluon dynamics is not explicit in the nlNJL model it may therefore be in order to superimpose a density-dependent bag pressure that should resemble nonperturbative effects from the gluon sector which are to cease at the deconfinement transition [60]. In the present work we will use the parameter µ < in order to vary the onset density of the deconfinement transition, while the other parameters remain fixed at the values of the set 2 from Ref. [44].

Mixed phase construction
For the completion of the compact star EoS, the hadron and quark matter EoS are joined either by a Maxwell construction or an interpolation procedure aimed at mimicking the effect of pasta phases that may exist at the interface. As for the latter case, we employ the replacement interpolation method (RIM) [61,62] that replaces an EoS region around the Maxwell construction point by a polynomial function of the type with α 1 , α 2 , µ H and µ Q that are determined from the continuity conditions at the hadron and quark border points mu H and mu Q , The resulting EoS connects the three points P H (µ H ), P c + ∆P = P c (1 + ∆ P ), and P Q (µ Q ). The parameter ∆P corresponds to an extra pressure at the Maxwell construction point thus we use it as a free parameter to quantify the effect of the mixed phase of the hadron to quark transition.

Constant speed-of-sound extrapolation
Due to a backbending of the given quark curves on P-ε plot a causality violation at high energy densities appears. To avoid such violation the constant speed-of-sound (CSS) extrapolation has been implemented for the quark tables. The CSS EoS for this extrapolation is given by where µ x is the chemical potential where the matching to the CSS extrapolation starts. The parameter β is directly related to the squared speed of sound so that It is obvious that fulfilling the causality constraint c 2 s ≤ 1 entails that β ≥ 2. The coefficients P 0 and P 1 in Eq. (12) are defined as where The results below have been obtained for the choice of the following rule. For those EoS for which c 2 s has a peak (below c 2 s = 1) at an energy density beyond the upper limit ε Q of the deconfinement phase transition, the squared speed of sound stays constant at the value it has at ε x = 818 MeV/fm 3 . For the other EoS the squared speed of sound stays at c 2 s = 1 for ε > ε x , where ε x is the energy density where it first reached this value. In the latter case of the CSS extrapolation with c 2 s = 1 the parameter The resulting EoS P(ε) and the squared speed of sound c 2 s (ε) are shown in the left and right panels of Fig. 1, respectively.

Compact star configurations
In this section we describe how we compute the sequences of hybrid compact stars and their tidal deformabilities which became now accessible to observation with the gravitational wave detectors of the LIGO-Virgo Collaboration and the analysis of neutron star merger events like GW170817.

Mass and radius
Using the full compact star EoS we obtain the mass-radius relations for stellar sequences by solving the Tolman-Oppenheimer-Volkoff (TOV) equations which describe a static, spherical star dP(r) dr = − (ε(r) + P(r)) m(r) + 4πr 3 P(r) r (r − 2m(r)) , where P(r) is the pressure profile, related to the energy density profile ε(r) by the EoS. The latter defines the profile function m(r) of the mass enclosed up to the distance r from the center of the spherical star configuration. The boundary conditions P(r = R) = 0, M = m(r = R) determine the total stellar mass M and radius R, respectively, in dependence on the initial value P c = P(r = 0) for the integration of the TOV equations. By varying this initial value, or equivalently the central energy density ε c = ε(r = 0) = ε(P c ), one obtains the M − R relationship which is uniquely characterizing the EoS. This equivalence between the EoS of compact star matter and the M − R relation for the compact star sequences it generates is the basis for constraining the EoS by mass and radius measurements on pulsars. One of the landmarks of a sequence of stars in the M − R diagram is the maximum mass which delimits a branch (family) of them that is stable against gravitational collapse and can thus be subject to comparison with observations. The observation of a pulsar mass that exceeds the maximum mass can be used to exclude the corresponding EoS. In Fig. 2  Mass-Radius diagram with solutions For comparison, the modern constraints from multi-messenger astronomy are shown: the mass of PSR J0740+6620 [5], the compactness constraint on tidal deformabilities [2] obtained by an analysis of the gravitational wave signal from the inspiral phase of the binary compact star merger GW170817, and one of both mass and radius measurement for PSR J0030+0451 by NICER [8]. Also shown for an orientation is the mass-radius constraint reported in Ref. [7] for the binary compact star merger GW190425, but it is not used in the Bayesian analysis presented here.
narrow range, there exist two possible star configurations with different radii on disconnected stable branches, belonging to the second family of ordinary neutron stars (here represented by the DD2p40 EoS) and to the third family of hybrid stars. The unstable branch between the second and third family is characterised by a positive slope. We did not remove it here in Fig. 2 in order to guide the eye from the 3rd family branch to the onset of the phase transition, where the unstable branch joins the 2nd family. When the phase transition proceeds through a mixed phase (when ∆ P = 0), then the mass M onset at the onset of the transition is lowered as compared to the Maxwell-construction case, and eventually the disconnected branches join so that there are no twins anymore. All hybrid star EoS in this study fulfill the maximum mass constraint that is represented here by the Shapiro-delay based mass measurement on PSR J0740+6620 [5]. Former observational maximum mass constraints from PSR J0348+0432 [6] and PSR J1614-2230 [63][64][65] are also shown, but no longer as constraining. The first joint observation of a gravitational wave signal from the inspiral phase of the merger GW170817 and the related electromagnetic signals in different bands of the spectrum (the kilonova) have resulted in additional constraints shown as magenta shaded excluded regions in Fig. 2. The constraint on the compactness of pulsars itself which was derived from the measurement of the tidal deformability by the LIGO-Virgo Collaboration [2] is shown as a shaded heart-shaped region in the mass range between 1.16 and 1.60 M . Also shown for an orientation is the mass-radius constraint reported in Ref. [7] for the binary compact star merger GW190425, but it is not used in the Bayesian analysis presented in the following section. The cigar-shaped region in Fig. 2 is one of both mass and radius measurements for PSR J0030+0451 by NICER [8]. The second one [9] is sufficiently similar and shall be explicitly considered in a subsequent work. For a comparison the soft hadronic EoS by Akmal, Pangharipande and Ravenhall (APR) [66] is also shown. We note that in the relevant mass range for GW170817 where tidal deformability constraints are available, the curves for the hybrid EoS with low M onset < M are almost indistinguishable from that for APR.

Tidal deformabilities and GW170817
The tidal deformability (TD) for a given star configuration is computed following the approach based on a perturbation of the static metric of the star [67]. The resulting perturbation functions, together with the TOV mass and pressure profiles, are the ingredients for the determination of the Love number k 2 [68]. The dimensionless TD is defined as Λ = λ/M 5 in terms of the deformability λ related to the Love number In Fig. 3  are possible for the same M, lying on the second and third family branch of compact stars, labeled by (2) and (3), resp. There are four combinations of stars in the binary GW170817 possible which fulfill the constraint of the measured total mass M tot = 2.73 +0.04 −0.01 M and which we label with two numbers as (22), (23), (33) and (32), denoting the origin of the primary and secondary star, resp., stemming from the second (2) or third (3) family. With the twofold degeneracy in the case of mass twins there are four lines generated when the tidal deformabilities Λ 1 and Λ 2 of the two stars in a binary are plotted against each other. These lines have endpoints which form the corners of a square in the Λ 1 − Λ 2 diagram shown in the right panel of Fig. 4. These endpoints correspond to the equal-mass case. Of course, any specific merger event like GW170817 is represented by just one point in the Λ 1 − Λ 2 diagram so that different merger events which cover the same mass interval and have overlap with the twin mass range may result in different peaks in the probability landscape on the Λ 1 − Λ 2 plane.
We note that for the case of compact star mass triplets, which were introduced in [42] and discussed in the context of GW170817 in [46], the mass-degenerate case would correspond to nine points, lying on the corners of a 3x3 quadratic array in the Λ 1 − Λ 2 diagram. We would like to draw the attention to the fact that EoS with mass twins and triplets are examples for which the Λ 1 − Λ 2 diagram is not reflection symmetric to the diagonal Λ 1 = Λ 2 because of the existence of separate families of stars corresponding to one and the same EoS. This situation appears as if the two components in the binary which are lying on different star branches were described by different EoS and it should be taken into consideration in the analysis of the gravitational wave signals. If the true compact star EoS would allow mass twin (or triplet) stars in a mass range accessible to merger events then in a combined analysis of the inspiral GW signals from many events in that very mass interval should emerge a multi-peaked (4-or 9-peaked) probability landscape in the Λ 1 − Λ 2 plane.
In Fig. 5 we explore how the pattern in the Λ 1 − Λ 2 plane would evolve when the mass ranges covered by the twin phenomenon and by the merger would be shifted against each other. Here we consider the latter fixed at that of GW170817 and vary the former by changing µ < , the parameter for the onset of deconfinement in the EoS; in the left panel for the Maxwell transition case (∆ P = 0) and in the right panel for different mixed phase parameters ∆ P = 2, 4, 6 NS 8%. For the transition at low densities (µ < < 1100 MeV) and thus low M onset < M the merger probes the hybrid star branch only which is almost indistinguishable from that for the soft APR EoS. For high transition densities (µ < > 1250 MeV) the onset mass exceeds the merger mass range and thus the tidal deformabilities are those of the stiff hadronic EoS DD2p40. In the intermediate range 1100 < µ < [MeV] < 1250, as a consequence of the overlap between merger and twin mass ranges, multiple disconnected branches occur in the Λ 1 − Λ 2 diagram. Besides the binary merger of two second family stars (22) and that of two third family stars (33) also the mixed cases occur where a pure neutron star from the second DD2p40 APR Figure 5. Tidal deformability constraint from the binary compact star merger GW170817 in the Λ 1 − Λ 2 diagram. According to the LVC analysis [2], the light and dark green regions correspond to 90% and 50% confidence, respectively, see also [69]. In the left panel we display results for the hybrid EoS with selected values (all in MeV) of the parameter µ < for the onset of deconfinement in the Maxwell construction case, see also Figs. 3 and 4. The right panel shows the effect of the mixed phase construction is for these selected values of µ < and ∆ P = 2%, 4% , 6% and 8% highlighted in color.
family merges with a hybrid star from the third family (23) and (32). In the right panel of Fig. 5 one can follow how the separate branches of would merge to a single line when the mixed phase parameter increases and this would reflect the joining of second and third family branches for the twin EoS in the M − R diagram. An observational tomography of the Λ 1 − Λ 2 plane equivalent to the theoretical one presented here would require the detection of GW signals from a large number of sufficiently "loud" merger events that would allow for a precise measurement of tidal deformabilities and to map different mass windows in order to "detect" the emergence of disconnected maxima in the probability landscape as a signal for the mass twin phenomenon.

Vector of Parameters
The set of parameters defining a set of models could be represented in the parameter space by introducing the vector of parameters. Each vector then corresponds to a particular model defined by the value of µ < related to the onset density of deconfinement and the value ∆ P specifying the transition construction, where i = 0..N − 1 and i = N 2 × j + k and j = 0..N 1 − 1, k = 0..N 2 − 1 and N 1 and N 2 are number of values of model parameters µ < and ∆ P , respectively.

Likelihood of a model under the Λ 1 -Λ 2 constraint from GW170817
In order to implement the tidal deformability constraint on the compact star EoS, reflected on the Λ 1 -Λ 2 diagram that includes probability regions from GW170817 event [1,2], we employ for the likelihood the formula where l ps is a line in the Λ 1 -Λ 2 plane formed by all admissible combinations of masses for the primary (p) and the secondary (s) star in the binary GW170817. The probability distribution function (PDF) β(Λ 1 , Λ 2 ) is given for GW170817 and the free parameter τ which, for instance, is the central density of the primary star, are included in equation (25). The PDF β(Λ 1 , Λ 2 ) (see Fig. 5) has been reconstructed (as previously in [18]) by the Gaussian kernel density estimation method with Λ 1 -Λ 2 data given at the LIGO web-page [69].

Likelihood of a model under the constraint on the lower limit of the maximum mass
The likelihood of the model is the conditional probability of the expected value of the possible maximum mass for the given model parameter vector: where M i is the maximum mass of the model given by π i . The values for µ A = 2.14 M and σ A = 0.095 M are taken from the mass measurement of PSR J0740+6620, M = 2.14 +0.10 −0.09 M [5]. We note that we combine this new measurement with the previously used mass measurement for the two solar-mass pulsar J0348+0432 with M = 2.01 +0.04 −0.04 M [6] as a product of independent likelihoods and thus obtain a strengthening of the selective power of the maximum mass constraint because in this way the models with a lower maximum mass get additionally suppressed.

Likelihood of a model under the combined M-R constraint of the NICER experiment
We have implemented in our Bayesian analysis the first simultaneous M-R measurement as provided by the NICER experiment for the millisecond pulsar PSR J0030+0451 [8,9]. We represent this measurement by a bivariate normal distribution in the M-R diagram where the parameters µ M = 1.44 M , σ M = 0.145 M , µ R = 13.02 km, σ R = 1.150 km, and the correlation coefficient ρ = 0.9217 have been calculated within the approximation (23) to the result of [8]. This corresponds to a tilted ellipse, see the red areas in Fig. 6. The likelihood for a model π i under the constraint from the NICER experiment is where the lines for the second and third families in M-R diagram are denoted as l 2 and l 3 , respectively. Besides the real NICER measurement one can assume a fictitious M − R measurement that could be a possible future result for either the same object or another one. We explore the consequences for the Bayesian analysis when such a measurement would have a standard deviation reduced by a factor 2 and a mean value of mass and radius that is shifted by their corresponding 1σvalues, i.e.

Posterior distribution
The full likelihood for the given π i can be calculated as a product of all likelihoods, since the considered constraints are independent of each other µ R = 13.02 km σ R = 1.15 km Figure 6. Probability measure for the mass-radius measurement of the NICER experiment for PSR J0030+0451 using values from the analysis of Ref. [8]. The blue lines correspond to a fictitious measurement with mass and radius shifted by 1σ to larger values and the error ellipse is narrowed by a factor two.
where m is the index of the constraints. The posterior distribution of models on parameter diagram is given by Bayes' theorem where P − → π j is the prior probability distribution in the space of N models that is taken to be uniform, P − → π j = 1/N.

Results of the Bayesian analysis
In Fig. 7 we show the posterior probability distribution in the space of hybrid EoS models spanned by the values for the two model parameters µ < and ∆ P for four cases of constraints (in clockwise order): 1) the combined mass measurements of PSR J0740+6620 [5] and PSR J0348+0432 [6] as a lower limit on the maximum mass; 2) this mass constraint together with the compactness constraint from tidal deformabilities of the binary compact star merger GW170817 [2]; 3) additionally the constraint from one of the mass-radius measurements by NICER [8], and 4) these three constraints when the NICER result is replaced by a fictitious measurement with mass and radius shifted by 1σ to larger values and the error ellipse is narrowed by a factor two, see Fig. 6.
From this analysis we deduce the most likely parameter cases, i.e. for which the posterior probability is at least 75% of the maximal probability and highlight them in the M − R diagram of Fig. 8 and the EoS P(ε) of Fig. 9. We compare the case 3) of three measurements (left panels) with Figure 7. Bayesian analysis using only the mass measurements of PSR J0740+6620 [5] and PSR J0348+0432 [6] as a lower limit on the maximum mass (upper left panel); this mass constraint together with the compactness constraint from tidal deformabilities of the binary compact star merger GW170817 [2] (upper right panel) and the constraint from one of the mass-radius measurement by NICER [8] (lower left panel). The lower right panel shows the result of a Bayesian analysis with these three constraints when the NICER result is replaced by a fictitious measurement with mass and radius shifted by 1σ to larger values and the error ellipse is narrowed by a factor two, see Fig. 6. the case 4) when the real NICER measurement is replaced by a fictitious one (right panels). We observe a remarkable strong effect on the selectivity due to the shift in the mean value of the NICER mass-radius measurement and the improvement of its accuracy. Out of a set of 50 EoS with the same prior probability remain 10 models with a posterior probability higher than 75 % of the maximum one, corresponding to a selectivity of 20%, which all would favor a low onset of the deconfinement transition at M onset < 0.75 M , so that the resulting preferable EoS would be almost indistinguishable from the soft hadronic APR EoS in the relevant mass range of GW170817 which is representative for the typical neutron star mass range. Replacing the real NICER measurement by a fictitious one just one EoS remains under these sharpened constraints and it corresponds to a higher M onset = 1.60 M , see the right panel of Fig. 8. APR DD2p40 µ < = 950, ∆ P = 2 % µ < = 950, ∆ P = 0 % µ < = 950, ∆ P = 4 % µ < = 950, ∆ P = 6 % µ < = 950, ∆ P = 8 % µ < = 975, ∆ P = 0 % µ < = 975, ∆ P = 2 % µ < = 975, ∆ P = 4 % µ < = 975, ∆ P = 6 % µ < = 975, ∆ P = 8 %  For an orientation, we show also the preferable region deduced from mass-radius constraint in the pre-multi-messenger era by Hebeler et al. [70].

Conclusions
We have developed a Bayesian analysis method for selecting the most probable equation of state under a novel set of constraints from compact star observations, which now include, besides the actual lower limit on the maximum mass from PSR J0740+6620 and the tidal deformability from GW170817 also the first result for the mass and radius reported by the NICER Collaboration, for PSR J0030+0451. We have applied this method for the first time to a two-parameter family of hybrid equations of state that is based on realistic models for hadronic matter (DD2 with nucleonic excluded volume) and color superconducting quark matter (generalized nlNJL model) which produce a third family of hybrid stars in the mass-radius diagram. One parameter characterizes the chemical potential for the onset of quark deconfinement (µ < ) while the other (∆ P ) belongs to the mixed phase construction that mimics the thermodynamics of pasta phases and includes the Maxwell construction as a limiting case for ∆ P = 0.
It turns out that an early phase transition is favored and the results are rather independent of the mixed phase parameter and whether the hybrid stars form a third family branch or the hybrid star sequence is connected with the hadronic sequence for sufficiently large ∆ P > 4 %. The presence of multiple configurations for a given mass (twins or even triplets) corresponds to a set of disconnected lines in the Λ 1 − Λ 2 diagram of tidal deformabilities for binary mergers, so that merger events from the same (sufficiently narrow) mass range may result in probability landscapes in this diagram with different peak positions corresponding to different origins of the stars in the binary.
The hybrid equation of state favored under this set of present observables corresponds to a line in the M − R diagram which is very similar to that for the APR equation of state which it practically resembles in the range of observed neutron star masses (1.1 . . . 2.1 M ). This is a modern reappearance of the masquerade problem discussed first in [71]. The situation would change dramatically if NICER would improve the accuracy of its mass and radius determination for J0030+0451 so that it would have only half of the present variance σ and the mean value of the radius (as well as the mass) would be corrected upwards by 1σ. We have examined the effect that such a so far fictitious measurement would have on our Bayesian analysis and find that an early onset of the deconfinement transition at masses well below 1 M is no longer favored. There remains just one favorable equation of state with an onset mass of 1.6 M and a Maxwell construction with a disconnected third family branch. Ongoing work explores other phase transition constructions, such as a crossover type interpolation between a soft hadronic EoS and a stiff high-density quark matter EoS as the original nlNJL model. Such a construction allows to discuss softer hadronic EoS for the hybrid EoS construction as, e.g., the APR or DD2F equations of state. We are also investigating the hyperon puzzle in this context [72,73]. Finally, there are systematic studies ongoing to unify the hadronic and quark matter descriptions [74][75][76][77]. A review-type publication of our results is in preparation.