Bayesian analysis of multimessenger M-R data with interpolated hybrid EoS

We introduce a family of equations of state (EoS) for hybrid neutron star (NS) matter that is obtained by a two-zone parabolic interpolation between a soft hadronic EoS at low densities and a set of stiff quark matter EoS at high densities within a finite region of chemical potentials $\mu_H<\mu<\mu_Q$. Fixing the hadronic EoS as the APR one and choosing the color-superconducting, nonlocal NJL model with two free parameters for the quark phase, we perform Bayesian analyses with this two-parameter family of hybrid EoS. Using three different sets of observational constraints that include the mass of PSR J0740+6620, the tidal deformability for GW170817, and the mass-radius relation for PSR J0030+0451 from NICER as obligatory (set 1), while set 2 uses the possible upper limit on the maximum mass from GW170817 as an additional constraint and set 3 instead of the possibility that the lighter object in the asymmetric binary merger GW190814 is a neutron star. We confirm that in any case, the quark matter phase has to be color superconducting with the dimensionless diquark coupling approximately fulfilling the Fierz relation $\eta_D=0.75$ and the most probable solutions exhibiting a proportionality between $\eta_D$ and $\eta_V$, the coupling of the repulsive vector interaction that is required for a sufficiently large maximum mass. We used the Bayesian analysis to investigate with the method of fictitious measurements the consequences of anticipating different radii for the massive $2~M_\odot$ PSR J0740+6220 for the most likely equation of state. With the actual outcome of the NICER radius measurement on PSR J0740+6220 we could conclude that for the most likely hybrid star EoS would not support a maximum mass as large as $2.5~M_\odot$ so that the event GW190814 was a binary black hole merger.


Introduction
The observation of the first binary neutron star merger GW170817 in gravitational waves [1] and the subsequent electromagnetic signals from the gamma-ray burst to the light curve of the kilonova [2] have opened the era of multimessenger astronomy. This extends the available mass range for neutron star observations up to 2.6 M for the companion star of the 23 M black hole in the binary merger GW190814 [3], if that object was indeed the heaviest neutron star and not the lightest black hole, which is a currently disputed question. The observation of gravitational waves from the inspiral phase of the merger GW170817 did allow to extract for the first time a new constraint on the equation of state (EoS) of dense matter, the tidal deformability, to be in the range of 70 < Λ 1.4 < 580 [4] for a neutron star with the mass of 1.4 M . From this measurement, together with other constraints, the authors of [5] could constrain the radius of a neutron star in that mass range to the rather narrow limits of R 1.4 = 11.0 +0. 9 −0. 6 km. An open and controversially discussed question is the interior composition of neutron stars, in dependence of their mass.
It is very likely that the quark substructure of nucleons manifests itself at increasing densities first by a stiffening of the EoS due to quark Pauli blocking in nuclear matter [6] and at still higher densities by a delocalization of the quark wave function and the occurrence of deconfined quark matter. For a recent discussion of soft delocalization vs. hard deconfinement in the transition from nuclear to quark matter, see [7]. A crucial open question, to which the present work intends to contribute, concerns the onset mass of deconfinement and the character of the transition [8].
A standard approach to the hadron-to-quark-matter transition would start from separate EoS models for these two phases and obtain the phase transition from a Maxwell construction (for sufficiently large surface tension between these phases) or a Glendenning construction of a homogeneous mixed phase (for vanishing surface tension) [9]. Inbetween these limiting cases, the more realistic scenario of the first-order phase transition would consider structures of finite size formed by the balance between Coulomb interactions and surface tension (pasta phases), see [10] and references therein. This approach has been used recently for a Bayesian analysis with observational constraints for masses and radii of neutron stars [11] which reaches the conclusion that very likely the phase transition onset occurs in the center of neutron stars with masses around 1 M and would then match the observed compactness [5] in this way. For this scenario to work, it is customary to have a sufficient stiffness of nuclear matter at supersaturation densities so that the deconfinement transition is driven to relatively low densities.
It is worth noticing that the approach [11] is based on a strong first-order phase transition which entails the formation of hybrid stars as a third family of compact stars [12] featuring the mass twin [13] phenomenon.
In Fig. 1 we illustrate this situation. A soft hadronic EoS like that of Akmal, Pandharipande and Ravenhall (APR) [14] has either no crossing (Maxwell construction) with the color superconducting quark matter EoS (considering a nonlocal version of the NJL model, nlNJL, described below) for the weak diquark coupling strength (η D = 0.71) or, at slightly increased dimensionless diquark coupling (η D = 0.79) an unrealistically early transition that is followed by a "reconfinement" or there is again no transition, depending on the value of the dimensionless vector meson coupling η V .
In Fig. 2 we show how such pathologies of the phase transition construction (or its inapplicability) with too simple EoS which are not suitable for such a construction, could be cured. A stiffening of the hadronic EoS, here realized by an excluded nucleon volume, leads already to a reasonable transition at not too low densities and to circumvention of the reconfinement problem of a second (unphysical) crossing of hadronic and quark matter EoS at higher densities 1 . The situation would still be improved towards a more realistic description when confining effects 1 For a discussion of the reconfinement problem see, e.g., [16], for the related masquerade problem, see [17] and for their solution see, e.g., [18]. would be included to the quark matter description, e.g., by a (density-dependent) bag pressure that resembles the effect of a nonperturbative QCD vacuum surrounding color charges (quarks) and leads to their confinement in color singlet multiquark states (hadrons). We note that without such a negative pressure (and/or a confining force), by the larger number of quark and gluon degrees of freedom on the one hand and the larger masses of hadrons in the spectrum on the other, the quark-gluon plasma phase would be favorable over the hadronic matter phase at low temperatures and densities, see Fig. 1 and [19]. As we already noted above, within a first-order phase transition, the formation of structures such a bubbles, droplets, rods and plates (pasta phases) is likely, with their sizes defined by an interplay of surface tension and Coulomb interaction effects. The resulting pressure (green curve in Fig. 2 looks as if a direct interpolation between the hadronic and the quark matter EoS would have been performed and the underlying three main microphysical ingredients (quark Pauli blocking, quark confinement and pasta structures in the mixed phase) could be circumvented by a direct shortcut from the nuclear matter phase just above saturation density to the quark matter phase which would then appear as a crossover-like EoS.
Such crossovers have been invoked on physical grounds by symmetry arguments as a quark-hadron continuity in Refs. [20,21,22] and by the combined effects of chiral symmetry breaking and diquark condensation intertwined by the axial anomaly so that they result in a crossover at low temperatures which eventually entails a second critical endpoint in the QCD phase diagram [23,24,25]. The crossover behaviour has subsequently been realized in ef- fective interpolating constructions following [26,27,28,29] and further literature in this direction.
While in interpolating constructions the information about the composition of the matter in the mixed phase is lost, they allow for conclusions on the likely properties of the pure quark matter EoS that is used as input, once this phase is reached in neutron star interiors. In contrast, the dominant class of EoS used to extract EoS constraints from observations of masses and radii [30,31,32] as well as tidal deformabilities [33,34] and in near future also the moment of inertia [35] are the multi-polytrope extrapolations of the EoS at supersaturation densities. Here we referred only to some prominent examples. These constraints have also been used to perform Bayesian analyses of the most likely EoS, see for instance the recent work of Miller et al. [36]. We note that analyses based on the multi-polytrope ansatz for the high-density EoS are agnostic of the composition of matter as well as constraints on the microphysics of dense quark matter.
In this paper, we will perform a Bayesian analysis study with modern mass and radius constraints, as well as fictitious radius measurements, on the basis of a new, twozone interpolation construction for obtaining hybrid EoS that allows for conclusions on the most favorable parameter set for the lagrangian model of color superconducting quark matter with nonlocal chiral interactions of the Nambu-Jona-Lasinio type. The EoS model, including the new interpolation procedure is described in the next section. The results for masses, radii and tidal deformabilities of hybrid stars are outlined in Sect. 3 and astrophysical inputs for the Bayesian analyses are given in Section 4. The results of the Bayesian analysis are presented in Section 6 and in Section 7 we draw the conclusions from this study.
2 New class of quark-hadron hybrid EoS by two-zone interpolation The idea is to interpolate between hadron and quark EoS models from the trustable region of the hadronic EoS at the nuclear saturation density (n H = 1.0...1.5 n 0 ) to the trustable region for the quark matter model at n Q 3n 0 (see fig. 1). Before we outline in detail the new interpolation method in subsection 2.3, we specify in the following two subsections the hadronic and the quark matter EoS that we employ in the present study to describe the pure phases outside the region of "terry incognita" indicated in Figs. 1 and 2.

Hadronic EoS
Our choice of hadronic equation of state for this work is the well known APR model [14]. It is a non-relativistic model derived by means of variational chain summation methods which included Urbana potentials of two and three nucleon interactions and features a pion condensate. Moreover, it exhibits a causality breach in neutron star matter for massive stars, a problem that shall not appear for the hybrid star models build in this work. The APR EoS version we have chosen is A18 + δv+UIX * which is not extremely stiff, reaching the maximum neutron star mass right below 2M .
In addition, in order to complete the description of the neutron star matter EoS, we adjoin a low density region EoS corresponding to the crust of neutron stars, namely the SLy4 model [37].

Quark matter EoS
For the description of the quark matter phase we consider a nonlocal chiral quark model, as in Ref. [15], which includes scalar quark-antiquark interaction, anti-triplet scalar diquark interactions and vector quark-antiquark interactions. The grand canonical thermodynamic potential per unit volume at zero temperature and finite density in the mean field approximation (MFA) reads see Ref. [38] for details of the calculation. The input parameters of the model are determined as to reproduce meson properties in the vacuum, at vanishing temperature and densities, then, m c , p 0 (the cutoff) and G S can be determined under that conditions. The remaining coupling constants G D and G V are driving the terms that, after bosonization, give rise to the superconducting gap field and the vector field. Then, the ratios η D = G D /G S and η V = G V /G S are input parameters. For OGE interactions in the vacuum, Fierz transformation leads to η D = 3/4 and η V = 1/2. As the microscopic interaction is not derived directly from QCD then, the above coupling ratios have in principle no strong phenomenological constraint except for the fact that η D values larger than η * D = (3/2)m/(m − m c ) may lead to color symmetry breaking in the vacuum [39] (where m stands for the dressed mass and m c for the current quark mass). In the present work we consider η D and η V as free parameters to be varied in reasonable limits, as it has been done for the NJL model case in Ref. [40]. The mean field valuesσ,∆ andω satisfy the coupled equations As we are focused on describing the behaviour of quark matter in the core of NSs, we have to impose: equilibrium under weak interactions, chemical equilibrium, and color and electric charge neutrality. Then, the six different chemical potentials µ f c in Eqn. (1) (depending on the two quark flavors u and d and quark colors r, g and b), can be written in terms of three independent quantities: the baryonic chemical potential µ, the electron chemical potential µ e and a color chemical potential µ 8 . So basically, for each value of µ we self-consistently solve the gap equations (2), complemented with the conditions for β-equilibrium and electric charge and color charge neutrality (details of the calculation can be found in the Appendix of Ref. [38]).
In the present work, we consider a Gaussian form factor g(p) = exp −p 2 /p 2 0 in Euclidean 4-momentum space. The input parameters of the quark model are fixed to m c = 5.4869 MeV, p 0 = 782.16 MeV and G S p 2 0 = 19.804 so that the pion mass m π = 139 MeV, the pion decay constant f π = 92.4 MeV and the chiral condensate − qq 1/3 = 244 MeV are reproduced.As mentioned above, we perform our analysis considering η D and η V as free parameters.

Two-zone interpolation method
We introduce here a two-zone interpolation scheme that is inspired by the discussion in subsection V.D (Fig. 15) and V.G (Fig. 18) of Ref. [8], see also Fig. 2 and its discussion in the previous section. According to that discussion, the interpolated part of the hybrid EoS can be motivated as a result of hadronic interactions and many-body forces in the hadronic phase (leading to P H (µ) → P * H (µ)) and confining effects in the quark matter phase (resulting in P Q (µ) → P * Q (µ)). In general, one would have to expect that these effects lead to a "normal" crossing of hadronic and quark matter P (µ) curves, but with a "kink" at the location µ c of the crossing point which results in a jump the density ∆n(µ c ) = dP * Q /dµ| µc −dP * H /dµ| µc > 0. In this work we discuss the continuous interpolation as the limiting case for ∆n(µ c ) → 0. An advantage of the two-zone interpolation between n H and n Q over a simple single-zone interpolation is that one can obtain good results without employing higher-order polynomials [41,42], just by using two parabolic functions where µ H and µ Q correspond to n H and n Q respectively, and µ c is free parameter taking value between them: µ H < µ c < µ Q .
The four parameters b η , b ρ , c η and c ρ can be immediately defined from the boundary conditions The remaining parameters a η , a ρ are defined by the matching conditions at at the intermediate µ c where both functions should be sewed These conditions are equivalent to the following system of linear algebraic equations (SLAE) where The determinant of this SLAE shows that the system has always a solution when µ c = µ Q , µ c = µ H and µ H = µ Q . The solution to the SLAE is We note, that this two-zone interpolation allows for a generalization to describe a first-order phase transition. This could be achieved by choosing a nonzero value for the jump in the density ∆n(µ c ), which appears as additional parameter in the second equation of the system (8).

Constant speed-of-sound representation
For the nlNJL model EoS following from (1), a causality violation at high energy densities (which corresponds to very high chemical potentials) appears due to a backbending of the quark pressure as a function of the energy density. To circumvent that obstacle, we make use of the recently discovered fact that in the range of densities relevant for NS applications the nlNJL model can be represented with high accuracy by a constant-speed-of-sound (CSS) EoS [43]. This EoS is given by [44,45] Here µ x is a scale for the chemical potential which we set to µ x = 1 GeV in correspondence with [43,45]. The parameters P 0 , P 1 and c 2 s are obtained from the parameters η D and η V of the nlNJL model by a functional mapping that has been defined in Ref. [43]. We note that the values for the squared sound speed obtained by this fit are in the range c 2 s = 0.45 . . . 0.54.

Hybrid EoS from a two-zone interpolation
The interpolation method has been implemented for APR and NJL models with different η D = 0.71(0.02)0.79 and η V = 0.06(0.02)0.20 (see Figs. [4][5][6][7][8]. The onset density for the interpolation has been set to n H = n 0 , and the density where the interpolation matches the quark matter EoS has been varied depending on η V as n Q = 4.5 . . . 2.5 n 0 while simultaneously η V was incremented. The value of µ c has been fixed as It is interesting to observe in Figs. 7 and 8 the similarity in the behaviour of the squared speed of sound with that of recent models of quarkyonic compact star matter [46]. A fast rise (stiffening of nuclear matter) is followed by a dip (hadron-quark mixed phase) and then saturates at an asymptotic value. This value in our case is the value obtained for the CSS fit of the nlNJL model and ranges between 0.45 and 0.54. Since the present neutron star phenomenology requires a stiffness of the EoS that results in a maximum mass of at least 2 M , the corresponding energy densities that are probed by neutron star interiors do not exceed about 1 GeV/fm 3 , see Fig. 9. This means that with the present setting of the two-zone interpolation construction, the pure quark matter phase for n B > n Q is barely reached in stable neutron star configurations, if at all.

Masses, radii and tidal deformabilities for the hybrid EoS
Neutron stars are computed within the framework of general relativity by using the corresponding EoS in the form of p(ε) to solve the Tolman-Oppenheimer-Volkoff (TOV) equations [47,48] dp(r) dr = − (ε(r) + p(r)) m(r) + 4πr 3 p(r) r (r − 2m(r)) , dm(r) dr = 4πr 2 ε(r). which describe a static, spherical star. Radial mass m(r), energy density ε(r) and pressure p(r) stellar internal profiles help to determined the mass M and radius R of a star with central density ε c with boundary conditions m(r = 0) = 0 and p(r = R) = 0. In Fig. 9 and Fig. 10 we show the compact star mass as a function of the central energy density and baryon density, respectively, for the two-zone interpolation construction. The dashed blue line corresponds to the hadronic APR EoS and the remaining curves show the hybrid EoS for a range of input parameters for the quark matter phase (same patters/colors as in previous figures). The compact star mass-radius sequence is a benchmark for every EoS model commonly presented in a mass-radius diagram that includes pulsar measurement regions as well as excluded ones by other astrophysical observations, see for instance Fig. 11. The EoS sequences displayed in there are obtained by a systematic integration of the the TOV equations for increasing p c for each single star up to the value of the maximum mass for which the condition ∂M/∂ε c > 0 holds. In addition, deformation of the compact star is a feature of the EoS that is closely related to the last moments of the inspiral phase of compact star mergers. It is quantified by computing the tidal deformability Λ, for which estimated regions were derived by the observation of the GW170817 event [1,4], also displayed in Fig. 11. The corresponding equations are derived from perturbations of the spherical metric of the compact star supplemented with the stellar internal profiles for physical quantities derived from the TOV equations. The equation relates the dimensionless tidal deformability Λ with the tidal Love number k 2 and the total mass and radius of the star. Details on the derivation of the above formula can be found in [49,50,51,52,53].

Bayesian inference
In this section we introduce the Bayesian methodology and its application to the set of considered EoS characterized by the parameters η D and η V in order to find their best values that fulfill observational data. The a posteriori probability P (π q |E) is a conditional probability of the given vector of parameters π q (introduced below), where q denotes the indexes the values of parameters for each alliterative representation of the model of EoS. The condition E is the set of the observational data (events), its likelihood for the given model is represented as a product which is the the conjunction of all events E α (where α is an index of an event). The a posteriori probabilities,   likelihoods and a priori probabilities are connected to each other via the Bayes formula where the factor P (π q ) is the prior of a given model. Then, the vector of parameters π q will be an element of the set H D × H V ,  Fig. 12. Mass-radius relations for EoS models featuring an interpolation scheme [8] between the low density APR model for hadronic matter and high density nlNJL quark matter. A few compact star sequences are displayed for two fixed quark matter parameter values of ηD with varying ηV together with the hadronic APR EoS. Different color regions correspond to either pulsar measurements or forbidden regions that serve as constraints for the compact star EoS. The green band region above 2M corresponds to the updated mass measurement of PSR J0740+6620 [55], which was recently upgraded to a mass-radius measurement by NICER, shown by the blue ellipsoidal region for the result of the Maryland-Illinois team [56]. The other blue ellipse corresponds to the mass and radius measurement of PSR J0030+0451 by NICER [58] whereas the grey and light green regions correspond to the estimates of the components of the binary system labeled as M1 and M2 of the GW170817 merger [4]. Red bands correspond to excluded regions derived from GW170817 observations by Bauswein et al. [59], Annala et al. [33] and Rezzolla et al. [60]. The black dashed horizontal lines are the upper and lower limit for the mass 2.59 +0.08 −0.09 M of the lighter component in the binary merger event GW190814 [3].
where q = iN V + j and N V = 8, N D = 5. Therefore, q = 0..N − 1 and N = 40 is the full number of model representations. For the choice of the uniform distribution of the prior we have P (π q ) = 1/N.
In the next section, we discuss the specific astrophysical constraints that we will employ in our Bayesian analysis (BA). Tidal deformabilities diagram. Λ1 and Λ2 correspond to the dimensionless tidal deformability for each of the components of the binary in GW170817. Light and dark green regions correspond to the 50% and 90% credibility regions for the posteriors used in the LIGO-Virgo analysis [4]. All the EoS in this work result in curves that fall inside the lighter region as well as the APR EoS which is displayed in blue.

Lower limit of maximum mass
Recently, the mass of the PSR J0740+6620 was obtained recently by a Shapiro delay based measurement combining data from the North American Nanohertz Observatory for Gravitational Waves (NANOGrav) and data from the orbital-phase-specific observations using the Green Bank Telescope. The 68.3% of the credibility interval was given as 2.14 +0.10 −0.09 M in [54], but it has recently been updated to 2.08 ± 0.07 M in [55]. This value has been chosen as the lower limit of maximum mass (M low max ).
The likelihood for the lower limit of maximum mass constraint is given by where M q is the maximum mass of the sequence of neutron star configurations for the given π q , and Φ(M, µ, σ) is the cumulative distribution function (CDF) of the standard normal distribution. And µ l and σ l are the parameters of the uncertainty of a low limit maximum mass estimation. Additionally, the assumption that one of the component of the binary merger GW190814 is a neutron star gives an estimation of the lower limit of the maximum mass as 2.59 +0.08 −0.09 M [3]. This value has been used in the Bayesian Inference as an alternative scenario for the lower limit. In Fig. 12 we display the mass-radius relations for our hybrid configurations together with the APR model as the low density baseline for hadronic matter that becomes invalid at higher densities where it crosses over to the nlNJL quark matter model. The different colored regions correspond to either mass and/or radius measurements or to forbidden regions following from GW170817 phenomenology that serve as constraints for the compact star EoS. The horizontal black dashed lines show the mass range for the lighter object in the binary merger GW190814, that we employ as a possible lower limit on the maximum mass, in the case that this object was a neutron star.

Upper limit of maximum mass
There is an estimation of the upper limit of maximum mass of neutron star in the literature [60]. It was estimated with combination of the observation of gravitational waves (GW170817) and drawing from basic arguments on kilonova modeling of GRB 170817A, together with the quasiuniversal relation between the maximum masses of static neutron stars and the fastest stable star under uniform rotation [62]. The upper limit of the maximum mass is 2.16 +0. 17 −0.15 M as shown in Fig. 12. The likelihood for the upper limit of maximum mass constraint is given by where M q is the maximum mass of the sequence of neutron star configurations for the given π q , and µ u = 2.16 M and σ u = 0.17 M . However, in Ref. [61] a relationship between the onset mass of prompt collapse to a black hole, the tidal deformability at half this mass and the maximum TOV mass has been derived, according to which the fact that the merger GW170817 did not promptly collapse to a black hole implies a lower limit on the maximum TOV mass. Therefore, we will include the disputable maximum mass constraint of Ref. [60] only in one of the sets of our Bayesian analysis.

Gravitational wave constraint
The observation of the gravitational waves from binary NS-NS merger GW170817 allows to calculate relation between tidal deformabilities of the primary and secondary components [1,4]. In order to implement the tidal deformability constraint to Bayesian Inference the Gaussian kernel density estimation has been used to recover the probability distribution function with use of the data on the Λ 1 − Λ 2 publicly shared by LIGO collaboration 2 .
The likelihood of the gravitational wave constraint is introduced as where l is the length of the line on the Λ 1 -Λ 2 diagram produced by − → π q , and n c is the central baryon density of a star. β(Λ 1 , Λ 2 ) is the probability distribution function (PDF) that has been reconstructed (as previously done in [63,64]). In Fig. 13 we show the dimensionless tidal deformabilities of hybrid com-pact stars configurations. The line colors/patters are the same as in previous figures. We display as well the corresponding measurement from GW170817. Fig. 14 shows our results for the tidal deformabilities Λ 2 as a function of Λ 1 . Dark and light green regions correspond to the 50% and 90% credibility for the posteriors used in the LIGO-Virgo Collaboration analysis, respectively. All the hybrid EoS in this work result in curves that fall inside the 90% region as the APR EoS which is displayed as dashed blue line.

Mass-radius constraint
A simultaneous measurement of mass and radius has been performed use data collected by the Neutron Star Interior Composition Explorer (NICER) space observatory for the pulsar PSR J0030+0451. The results of the observation have been reported in a collection of publications 3 .
There were two independent analysis of the mass and equatorial radius based on mutually exclusive assumptions about the uniform-temperature emitting spots. The first result for radius and mass is M 1 = 1.44 +0. 15 −0.14 M and R 1 = 13.02 +1.24 −1.06 km [58] whereas the second one is M 2 = 1.34 +0. 15 −0.16 M and R 2 = 12.71 +1.14 −1.19 km [65]. A bivariate probability distribution function α(M, R) has been reconstructed by the method of the Gaussian kernel density estimation using the data [66]. A likelihood is formulated as where M (n c ) and R(n c ) are mass and radius of sequence of neutron star for a given q th equation of state, and n c is the central baryon density. Another simultaneous measurement of the mass and the equatorial radius has recently been published based on NICER observations combined with XMM Newton data of PSR J0740+6620. The two independent analysis teams reported R = 13.7 +2.6 −1.5 km [56] and R = 12.39 +1.3 −0.98 km [57]. To implement these results in a Bayesian analysis, the function α(M, R) has been reconstructed with a kernel density estimation using mass-radius data from [67]. The above constraints are shown in Fig. 12.

Fictitious radius measurements
Simultaneous measurements of masses and radii like the ones provided recently by NICER are very important for deriving constraints on the EoS. However, due to limited observation time, the difficult decision has to be made for which object with a known mass should a precise radius measurement (and thus a long-term observation campaign) bear the largest discovery potential. A BA can support such decisions by using fictitious results of a radius measurement on the object under scrutiny. An essential input for the BA would be precise radius measurements for pulsars with known masses. One example is the high-mass millisecond pulsar PSR J0740+6620, for which the mass is known from Ref. [54] with the recent update by [55], but for which by the time this paper was written the results of the NICER radius measurement were not yet available. Now, after their publication in Refs. [56,57], we can compare the these real results with the predictions of the fictitious radii and suggest the method of fictitious radii for supporting decisions on the selection of targets for future measurements. We consider here three different values for the radius of PSR J0740+6620, namely R = 11, 12 and 13 km with the the same design uncertainty σ = 0.5 km of the NICER experiment. The likelihood for fictitious measurements is the same as introduced above (23).

Sets of constraints
We suggest three sets of constraints: 1: the mass measurement for PSR J0740+6620 [55] as the lower limit for the maximum mass, the tidal deformability from GW170817 [4] and the mass-radius constraint from PSR J0030+0451 [58] (set 1); 2: in addition to the constraints of set 1, the constraint on the upper limit of the maximum mass from Ref. [60] is included; 3: as for set 1, but assuming that the lower mass companion of the black hole in the asymmetric binary merger GW190814 [3] was a neutron star, the lower limit for the maximum mass is replaced by the lower limit on its mass M 190814 = 2.59 +0.08 −0.09 M . Besides these "pure" sets of constraints, we investigate for each of them the possibility of an additional massradius constraint for the pulsar PSR J0740+6620, for which the mass 2.08 ± 0.07 M is rather precisely measured and a result for the radius measurement by the NICER experiment is anticipated to yield the values of 11, 12 and 13 km. We denote these subsets with the numbers 1, 2 and 3, respectively. The entirety of Bayesian constraint sets in this work is synoptically summarized in Tab. 1. A separate BA is performed using the recent measurement of the combined mass-radius constraint for PSR J0740+6620 on the basis of NICER and XMM Newton data, and compared to the predictions using fictitious radii.

Results of the Bayesian analysis
In this section we show and discuss our results for the BA of the hybrid EoS for compact stars under the constraints defined in the previous section and summarized in Tab. 1. The details of the considered families of hybrid EoS, the interpolation method, the constant speed of sound representation, the BA and astrophysical constraints were already presented in previous sections. Figure 15 shows the results of our BA for the η V and η D values of our EoS models under consideration. The difference between the three LEGO plots is that for their derivation the values of the maximum mass constraint have been changed. These three different cases comprise 1) the mass measurement of the object PSR J0740+6620, [55], 2) the previous measurement plus the upper limit of the maximum mass from Ref. [60], 3) the first measurement plus the mass for the lighter object in the binary merger GW190814 [3] assuming that it was a neutron star. We can see that whereas for the first to cases the posterior probability distributions peak at intermediate values of the η V parameter of quark matter, the third case favours its highest values. The probability distribution in η D direction remains flat for the three cases.
The situation becomes more interesting when we consider a new, fictitious mass and radius measurement which anticipates the outcome of the NICER observation of the heavy pulsar PSR J0740+6620. In figure 16 each column corresponds to the results for such a radius measurement with a value of R = 11 km, R = 12 km or R = 13 km, whereas the rows represent the same three aforementioned cases for the constraints. Just like before, case 3) displays a more selective effect on the η V parameter, however the inclusion of this new radius observation favours the highest η D values for most of the plots.
In addition, figures 17, 18 and 19 present the equivalent contour plot constraints for the sets 1, 2 and set 3, respectively, as well as the Bayesian most probable EoS together with their mass-radius curves and tidal deformabilities for the GW170817 case. While this paper was under revision in the refereeing process, the NICER experiment has published its results for the M −R measurement [56], shown in Fig. 12. We have performed a BA using the NICER M − R measurements for PSR J0740+6620 [56] and PSR J0030+0451 [58] as well as the LVC tidal deformability constraint from GW170817 [4]. The result is shown in Fig. 20. A comparison with the three columns for fictitious radii in Fig. 17 shows that the actual result of the NICER measurement is well compatible with the R = 12 km case. Therefore, the method of fictitious radii anticipation can be considered a reliable tool in predicting the implications of future M − R measurements, e.g., by the NICER collaboration. Moreover, for the most     15. Bayesian analysis using the mass measurement for PSR J0740+6620 [55] as the lower limit for the maximum mass, the tidal deformability from GW170817 [4] and the mass-radius constraint from PSR J0030+0451 for the class of hybrid EoS obtained with a two-zone interpolation between APR and nlNJL in the two-dimensional EoS parameter plane spanned by ηV and ηD (left panel; set 1). In the middle panel the constraint on the upper limit of the maximum mass from Ref. [60] has been added (set 2) and in the right panel this limit has been lifted again in favor of the new lower limit on the maximum mass from lower-mass companion of the black hole in the asymmetric binary merger GW190814 [3] that replaces the one from [55] (set 3), if the former object would not be a black hole.
probable parameter set, the corresponding EoS lead to the M − R sequences highlighted in Fig. 20 which fulfill the 2 M mass constraint but do not reach masses above 2.5 M which would be required if the lighter object in GW190814 was a (hybrid) neutron star.

Conclusions
In the present work we have applied Bayesian analysis methods to investigate the most likely quark-hadron hybrid EoS among a family of models that are agnostic about the detailed microphysical scenario of the hadron-to-quark matter transition. The model uses the APR EoS for densities below 2.0 n 0 and above n = 4.0 n 0 a set of nonlocal, color superconducting chiral quark model parametriza-tions of the NJL type with a covariant Gaussian formfactor. The vector meson coupling η V and the scalar diquark coupling η D are varied as free parameters that determine the stiffness of high-density quark matter. The transition between both regimes is constructed by a new two-zone interpolation that realizes a smooth crossover behaviour, due to the assumption that the nature of the transition is a mixing of phases (crossover transition). Thus, in performing this interpolation we exclude the possibility of a first-order transition associated with a jump in the density, setting the corresponding parameter ∆n = 0. In future work, we plan to systematically investigate also the extension of the new two-zone interpolation construction to the cases with ∆n = 0. We note, however, that the interpolation construction may be viewed as a shortcut that replaces three microphysical effects: 1) stiffening of the nu-  16. Probabilities when an additional measurement of the radius R of PSR J0740+6220 (as provided recently by the NICER experiment [56]) is taken into account with anticipated values of R = 11 km (second column), 12 km (third column) or 13 km (fourth column) and with a standard deviation of σR = 0.5 km. The mass is taken from the measurement by Fonseca et al. [55]. The three rows correspond to set 1, set 2 and set 3, respectively.
clear EoS due to quark substructure effects (quark Pauli blocking modeled, e.g., by a baryon excluded volume), 2) softening of the quark matter EoS at low densities due to confining effects (modeled, e.g., by a medium-dependent  Fig. 17. Results of the BA for set 1 which includes the constraints (inf{Mmax} [55], Λ1.4 [4], (M, R)J0030+0451 and [58]) in the leftmost column and with an additional (yet fictitious) NICER radius measurement for PSR J0740+6620 of R = 11, 12 or 13 km with an estimated standard deviation of σR = 0.5 km in the other three columns. The highlighted most probable M-R sequences (2nd row), EoS (3rd row) and Λ1 − Λ2 (4th row) relationships correspond to the parameter sets with at least 75% of the maximum probability as shown in the first row.
bag pressure) and 3) mixed-phase effects due to the occurrence of finite-size structures (pasta phases). This possible microphysical underpinning of the interpolation construction, in particular within the two-zone version suggested in this paper, has been discussed and illustrated in the Introduction, but a detailed exploration is deferred to future work.
In our Bayesian study we have applied standard constraints for mass and radius measurements in set 1 and demonstrated their effect of narrowing the viable range of parameter values. We obtained the result that the most probable EoS lie along a line of proportionality between η V and η D , whereby the higher values of η V are favorable for obtaining larger maximum masses of hybrid neutron stars. This finding confirms similar results of earlier studies in [68,40] and the recent work [25] which do not employ Bayesian methods.
In the sets 2 and 3 we explored the nonstandard constraints of an upper limit on the maximum mass and the high mass of the lighter companion object of GW190814 as a lower limit on the maximum mass, respectively. While for the set 2 the narrowing of the parameter range due to the upper limit on the maximum mass leads to an exclusion of the higher values of the vector meson coupling, the high value of the lower limit on the maximum mass for set 3 allows only the highest possible vector couplings resulting in stiffer EoS and thus larger maximum masses and radii. It is a remarkable fact that within the present interpolation approach the 2.6 M companion object in GW190814 could be a hybrid star with quark matter core. This possibility has been pointed out before in several works, among them [45,70].
Finally, we have used the Bayesian approach to explore the consequences that the radius measurements on the 2 M pulsar PSR J0740+6620 by the NICER experiment may have for neutron star phenomenology.
Before the NICER radius measurements on J0740+6620 [56,57] were released, we have performed Bayesian anal-  Fig. 18. Results of the BA for set 2 which includes the constraints (inf{Mmax} [55], Λ1.4 [4], (M, R)J0030+0451 [58], sup{Mmax} [60]) in the leftmost column and with an additional (yet fictitious) NICER radius measurement for PSR J0740+6620 of R = 11, 12 or 13 km with an estimated standard deviation of σR = 0.5 km in the other three columns. The highlighted most probable M-R sequences (2nd row), EoS (3rd row) and Λ1 − Λ2 (4th row) relationships correspond to the parameter sets with at least 75% of the maximum probability as shown in the first row.
yses with fictitious results for that radius measurement, anticipating 11 km, 12 km and 13 km as a possible outcome. For a radius value as small as 11 km or less, the hybrid star scenario for the 2.6 M object in GW190814 could be excluded and this event was a binary black hole merger. In that case the present two-zone interpolation approach with continuous crossover is not suitable since it does not produce stable hybrid stars with small enough radii above 2 M .
For a radius of J0740+6620 as large as 13 km or more, the likely interpretation of GW190814 is that of a neutron star -black hole merger where the neutron star possibly had a color superconducting quark matter core.
We have performed a Bayesian analysis with the actual results for the NICER radius measurement on PSR J0740+6620 as reported by the Maryland-Illinois team [56] and compared the outcome for the probability contours in the model parameter space with those for the fictitious radius measurements. A striking similarity was obtained for R = 12 km, which proves the usefulness of the fictitious radius method for estimating the outcome of future NICER radius measurement campaigns and its impact for the dense matter EoS and neutron star phenomenology.
The outcome of the Bayesian analysis in our setting favors the scenario that GW190814 was a binary black hole merger and not a neutron star -black hole merger. The favored EoS which are highlighted in Fig. 20 favor M − R sequences with a maximum mass that lies well below 2.5 M , as shown also in that figure.
In future extensions of the Bayesian approach to neutron star phenomenology it is desirable to widen the class of EoS by extending the present direct smooth interpolation approach to a first-order phase transition scenario [11] where in the present two-zone interpolation method a nonzero density jump parameter at the matching point would be chosen, and to augment this with the possibility  Fig. 19. Results of the BA for set 3 of constraints (inf{Mmax} [3], Λ1.4 [4], (M, R)J0030+0451 [58]) in the leftmost column and with an additional (yet fictitious) NICER radius measurement for PSR J0740+6620 of R = 11, 12 or 13 km with an estimated standard deviation of σR = 0.5 km in the other three columns. The highlighted most probable M-R sequences (2nd row), EoS (3rd row) and Λ1 − Λ2 (4th row) relationships correspond to the parameter sets with at least 75% of the maximum probability as shown in the first row. of a subsequent smoothing of the transition by a pasta phase construction or its emulation.
In order to compare the results of two Bayesian analyzes obtained on the same astrophysical data, it is necessary to combine the factors of the normalization of the Bayesian formula (18). This allows to introduce the relative posterior probability of each analysis. The set of EoS models, whose analysis gives a greater value of the relative posterior probability can be considered to be more successful in describing the observational data. A comparison of the result presented in this paper with the previous result from Ref. [11], where constraints corresponding to set 1 have been employed, shows that the likelihood for models considered in [11] exceeds that of the present work by more than 5 times. Such a direct comparison is possible only for the case of set 1, because the other two sets have not been considered in [11].
The question has been raised whether one should generalize the Bayesian analysis by sytematically varying also the hadronic EoS at supersaturation densities like, e.g., within an excluded volume approach (see Ref. [71]). In the present setting of the two-zone interpolation model, with a choice of the onset of the interpolated part of the EoS at the saturation density (n H = n 0 ), there is no room for such a variation. Uncertainties in the knowledge of the nuclear EoS at supersaturation densities are covered by the variation of the parameters in the hadron-like zone of the interpolation.