Bayesian Analysis for Extracting Properties of the Nuclear Equation of State from Observational Data including Tidal Deformability from GW170817

We develop a Bayesian analysis method for selecting the most probable equation of state under a set of constraints from compact star physics, which now include the tidal deformability from GW170817. We apply this method for the first time to a two-parameter family of hybrid equations of state that is based on realistic models for the hadronic phase (KVORcut02) and the quark matter phase (SFM$\alpha$) which produce a third family of hybrid stars in the mass-radius diagram. One parameter ($\alpha$) characterizes the screening of the string tension in the string-flip model of quark matter while the other ($\Delta_P$) belongs to the mixed phase construction that mimics the thermodynamics of pasta phases and includes the Maxwell construction as a limiting case for $\Delta_P=0$. We present the corresponding results for compact star properties like mass, radius and tidal deformabilities and use empirical data for them in the newly developed Bayesian analysis method to obtain the probabilities for the model parameters within their considered range.


Introduction
The determination of the equation of state (EoS) of compact star (CS) interiors is a subject of active research by means of both laboratory measurements and astronomical observations. One of the most urgent questions concerns the possible existence of deconfined quark matter in CS cores, where matter densities can exceed by several times the nuclear saturation value, n 0 = 0.15 fm −3 , the typical density in large atomic nuclei. Hybrid compact stars have an inner core composed of quark matter surrounded by an outer core of nuclear matter. Of particular interest is the CS mass twin (MT) phenomenon [1]: when a pair of stars has the same gravitational mass but different radii, with the larger star being composed only of hadronic matter and the smaller one being a hybrid star with a quark matter core. The presence of MTs is synonymous to the existence of a third family [2] of CS in the mass-radius (M − R) diagram. It requires the CS EoS to feature a strong first order phase transition from hadron to quark matter [3][4][5][6][7][8], where the latent heat of the transition fulfils the Seidov criterion [9].
Should it turn out that the mass twin phenomenon can be discovered in CS observations this would give indirect evidence for the existence of a critical endpoint in the QCD phase diagram [4,5,10] which is sought for in relativistic heavy-ion collisions, so far without being conclusive.
On the other hand, astronomical observations may provide measurements of masses and/or radii which are relevant for constraining the nuclear EoS via the one-to-one relationship provided by the Tolman-Oppenheimer-Volkoff (TOV) equations [11,12] which govern the general relativistic hydrodynamic stability of CS matter under its own gravity. Neutron star masses are accurately measured by monitoring binary pulsar dynamics whereas the radii were determined so far with big uncertainties. The situation has changed with the first observation of the gravitational wave signal from the inspiral phase of the binary CS merger GW170817 [13] which allowed to constrain the tidal deformability of both merging stars and thus their mass and radius [14][15][16][17].
Consequently, theoretical approaches are required to quantitatively assess the most probable EoS out of a whole class obtained by varying intrinsic model parameters. One of the most powerful methods is the Bayesian analysis (BA) or Bayesian interpretation of probability that allows for estimation of model parameters based on prior knowledge, in this case the already collected measurements. Several works in this direction include BA based on observations of X-ray bursters [18] that unfortunately suffer from uncertainties in the stellar atmosphere composition modeling required in the interpretation of the observed signal or the more general approach of [19] that includes a generic multipolytrope EoS with priors that include experimental nuclear matter pressure values and that is able to support the most massive CS of about 2M and reports an accuracy of about 30% in pressure values. The recent approach of [20] employs X-ray timing observations of accretion-powered millisecond pulsars, resulting in radius estimates of about 5% if the CS mass is known. Moreover, a more stringent analysis [21] that incorporates the hypothesis of the Direct Urca cooling constraint [22] in addition to the afore mentioned measurements concludes that the neutron star radius has a value of 12.7 ± 0.4 km for masses ranging from 1 up to 2 M .
Despite the fact that the recent detection of gravitational radiation from the inspiral phase of the binary CS merger GW170817 has led to constraints on the tidal deformability of CS matter and to the discussion of the possibility that one or even both of the CS in the binary could be hybrid stars with quark matter interior and there would be a possibility to discover the mass twin phenomenon and thus a strong first order phase transition in CS matter [23][24][25][26][27], these data have not yet been included in BA studies such as the ones already listed above. Thus, it is of great interest to update such BA studies by incorporating this new information to constrain the CS EoS.
In this work we consider the mass twin EoS and modifications to it in order to perform new Bayesian analysis calculations that, in addition to the previously used CS measurements in [28], will include the GW170817 data as priors to assess the probability of model parameters. Our model choice is the KVOR EoS for hadronic matter together with the String-Flip approach for the deconfined quark regime.

Hybrid EoS
The hybrid EoS has been produced with the one-parameter replacement interpolation method (for a second order polynomial ansatz for P(µ)) [29][30][31] which mimics the thermodynamic behavior of pasta phases in the coexistence region between the pure hadronic phase and the pure quark one.
The hadronic phase in this work is described by the well known KVOR equation of state [32] with a modification of stiffness as introduced in [33] and denoted as KVORcut2. This particular version of the KVOR EoS model is the stiffest one presented in that work and allows to fulfill the condition [9] for the latent heat of phase transition to quark matter, so that the disconnected hybrid star family at higher densities is possible.
The quark phase is based on the String-Flip Model with the so called density functional [34] including the available volume fraction Φ, where the parameter α describes the effective density-dependence of the confining interaction, and it is varied here in the range [0.1 . . . 0.3]. The value of this parameter is correlated with the critical density of the phase transition, as it is shown in the figures below. The larger α, the lower the critical density and therefore the critical mass for the onset of the phase transition in the hybrid star. This systematics allows to calibrate the range of the model parameters with the observational data without contradicting the known properties of nuclear matter at saturation density. For the construction of the hybrid star EoS from the hadronic and quark matter EoS we employ here the mixed phase construction that is described in detail in Refs. [29][30][31]. It introduces the free parameter ∆ P = P(µ c )/P c − 1, the pressure increment at the critical chemical potential µ c relative to the critical pressure P c of the corresponding Maxwell construction, which is then obtained as limiting case for ∆ P = 0. The resulting family of hybrid EoS corresponds to the one described in [29], see Fig. 1. This way of modelling the phase transition mimics the possibility of so-called pasta phases characterized by different geometric structures in the coexistence region of the transition sufficiently well, as has been shown in [35]. In that reference it has also been demonstrated that the parameter ∆ P of the transition construction can be related to the value of the surface tension at the hadron to quark matter interface.

Mass and Radius
The structure and global properties of compact stars are obtained by solving the TOV equations where P(r), ε(r) and M(r) are the profiles of pressure, energy density and enclosed mass as a function of the distance r from the center of the star. The radius R of the star is obtained from the condition of vanishing pressure at the surface P(R) = 0 and the gravitational mass of the star is M = M(R).

Tidal deformability
The recent observation of gravitational waves from the inspiral phase of the binary CS merger GW170817 gives an estimation of the tidal deformability of these stars, therefore we include also this information in the analyses. The tidal deformability of a star can be calculated from the Einstein equation for small elliptic deformation as described in [36]. For the results of the calculation see Fig. 2 [13] shows that if the hadronic EoS is as stiff as DD2p40 at least one of the stars in the binary has to be a hybrid star in order to avoid a violation of the Λ 1 -Λ 2 constraint. As can be seen from Fig. 17 of Ref. [23] the merger of two hybrid stars which is admissible when the onset mass for the deconfinement transition is lowered, e.g., by increasing α, would lead to the appearance of a new branch in this diagram mimicking the pattern of a merger of two neutron stars with a rather soft nuclear EoS (like APR or SLy4).

Vector of Parameters
The set of parameters of models could be represented in the parameter space with introduction of the vector of parameters, each vector is one fixed model from above described types of hadronic, quark phases and 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 correspondingly.

Likelihood of a model for the Λ 1 − Λ 2 constraint from GW170817
In order to implement the tidal deformability constraint on the compact stars EoS, reflected on the Λ 1 -Λ 2 diagram that includes probability regions from GW170817 event [13,37], we employ the definition of the probability as an integral over the probability distribution function (PDF) β(Λ 1 , Λ 2 ) by when both stars in the binary merger belong to the second family branch of neutron stars (and l 22 is the corresponding path in the Λ 1 − Λ 2 plane), or by in case the parameter vector π i corresponds to hybrid star equation of state with a third family of compact stars in the mass range relevant for the merger. Then l 23 is the path in the Λ 1 − Λ 2 plane corresponding to a binary merger composed of a neutron star and a hybrid star. The parameter τ defines the position of a point on the paths l 22 and l 23 in equations (5)- (6). Note that the PDF has been reconstructed by the method Gaussian kernel density estimation with Λ 1 − Λ 2 data given at LIGO web-page [38], see fig. 3.

Likelihood of a model for the mass constraint
The likelihood of the model is the conditional probability of the expected value of the possible maximum mass for the given model parameter vector: here M i is maximum mass of the given by π i , and µ A = 2.01 M and σ A = 0.04 M is the mass measurement [39].

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 The posterior distribution of models on parameter diagram is given by Bayes' theorem where P − → π j is a prior distribution of a models taken to be uniform: P − → π j = 1/N.   [29], the sequences with a lower onset mass for the deconfinement transition are not accessed for which an interpretation of GW170817 as a merger of two hybrid stars would be possible.

Results
The results of the Bayesian analysis are given in Fig. 4.The most likely EoS are those with a strong mixed phase effect described by ∆ P > 0.04 and with a large screening parameter α > 0.28, at the limit of the parameter range considered here. For these parameter sets the phase transition and therefore the compactification of the hybrid star configuration occurs within the mass range that is relevant for GW170817. These results indicate that it should be worthwile to repeat this exploratory calculation with a wider set of parameters, in particular a larger screening parameter α.

Conclusions
We have developed a Bayesian analysis method for selecting the most probable equation of state under a set of constraints from compact star physics, which now include the tidal deformability from GW170817. We have applied this method to a case study that employed the class of hybrid EoS introduced in [29] that allow for the existence of a third family of compact stars. We have investigated the requirements on future measurements to find the equivalent phenomenon of mass twins, i.e. two objects with the same mass but different, distinguishable radii. A mass -radius relation with two branches is mapped also on the Λ 1 − Λ 2 diagram as it is shown here in Fig. 3 and other recent publications [23][24][25][26][27]. Since binary compact star mergers are expected to be observed by the LIGO-Virgo Collaboration at a rate of 1-10 events per year, we expect that with the observation of next merger events a binodal structure of the PDF in the Λ 1 − Λ 2 plane could become apparent as a manifestation of the low-mass twin case with an onset of the third family branch below ∼ 1.3 M if such a branch exists in nature. A similar suggestion has been already proposed by Christian et al. [25]. Such an observation would support the existence of an early phase transition, around 2n 0 for strong in-medium screening of the string tension.

Acknowledgements
We acknowledge discussions with K. Maslov on the hybrid star EoS. A.A., D.B., and H.G. acknowledge support from the Russian Science Foundation under grant No. 17-12-01427 for the work described in sections 2 -6. D. A.-C. is grateful for partial support from the Ter-Antonian -Smorodinsky program for collaboration between JINR and Armenian scientific institutions and from the Bogoliubov-Infeld program for collaboration between JINR and Polish Institutions.