  Neutron stars (NSs) which could contain exotic degrees of freedom in the core
and the self-bound quark stars (QSs) made purely of absolutely stable
deconfined quark matter are still two main candidates for the compact objects
observed in pulsars and gravitational wave (GW) events in binary star mergers.
We perform a Bayesian model-agnostic inference of the properties of NSs and QSs
by combining multi-messenger data of GW170817, GW190425, PSR J0030+0451, PSR
J0740+6620, PSR J1614-2230, PSR J0348+0432 as well as ab initio calculations
from perturbative quantum chromodynamics and chiral effective field theory. We
find the NS scenario is strongly favored against the QS scenario with a Bayes
factor of NS over QS $\mathcal{B}^\text{NS}_\text{QS} = 11.5$. In addition, the
peak of the squared sound velocity $c_s^2 \sim 0.5c^2$ around $3.5$ times
nuclear saturation density $n_0$ observed in the NS case disappears in the QS
case which suggests that the $c_s^2$ first increases and then saturates at
$c_s^2 \sim 0.5c^2$ above $\sim 4n_0$. The sound velocity and trace anomaly are
found to approach the conformal limit in the core of heavy NSs with mass $M
\gtrsim 2M_{\odot}$, but not in the core of QSs.
