Bayesian inference favors quark matter in neutron star interiors

We perform a physics-informed Bayesian analyses of the equation of state of hybrid neutron stars that incorporates color-flavor-locked quark matter modeled by a three-flavor non-local Nambu-Jona-Lasinio framework with vector repulsion and diquark pairing. Contrary to the model-agnostic Bayesian analyses our scheme allows for distinguishing between the scenarios of neutron stars with quark cores and without them. The used quark model realizes asymptotic conformality at high densities in accordance with perturbative QCD. The hadronic sector is described by the density-dependent relativistic functional DD2Y-T, which satisfies chiral effective field theory constraints and includes hyperonic degrees of freedom. We construct a large set of candidate hybrid EOSs by varying the vector and diquark couplings and apply a Maxwell construction for the quark-hadron phase transition. Observational constraints from recent NICER pulsar mass-radius measurements and tidal deformability from GW170817 are incorporated into the likelihood. Depending on whether the observational data from the black widow pulsar PSR J0952-0607 and the HESS J1731-347 object are included to the analysis or not, the posterior distribution favors vector and diquark couplings around $(\eta_V,\eta_D)\simeq (0.82,0.40)$ or $(\eta_V,\eta_D)\simeq (0.64,0.36)$, respectively. This corresponds to equations of state that support two-solar-mass neutron stars with superconformal speed of sound and relatively low onset densities for deconfinement. Our findings indicate that the most probable hybrid EOSs are statistically preferred over the purely hadronic baseline. The corresponding probabilities of agreeing with the observational data differ by one or two orders of magnitude depending on the data set used. This suggests that quark cores may exist in all observed neutron stars.