  This paper presents a new analysis of the thermal emission from the neutron
star surface to constrain the dense matter equation of state. It is based on
the use of a Markov-Chain Monte Carlo algorithm combined with an empirical
parametrization of the equation of state, as well as the consistent treatment
of seven neutron star quiescent low-mass X-ray binaries in globular clusters
with well-measured distances. Previous analyses have indicated that the thermal
emission of these neutron stars tends to prefer low neutron star radii,
questioning basic knowledge from nuclear physics. We show that it is possible
to reconcile the thermal emission analyses with nuclear physics knowledge, with
or without including a prior on the slope of the symmetry energy $L_{\rm sym}$.
We obtain radii of the order of about 12~km without worsening the fit
statistic. With an empirical parametrization of the equation of state, we
obtain the following values for the slope of the symmetry energy, its curvature
$K_{\rm sym}$, and the isoscalar skewness parameter $Q_{\rm sat}$: $L_{\rm
sym}=37.2^{+9.2}_{-8.9}$ MeV, $K_{\rm sym}=-85^{+82}_{-70}$ MeV, and $Q_{\rm
sat}=318^{+673}_{-366}$ MeV. For the first time, we measure the values of the
empirical parameters $K_{\rm sym}$ and $Q_{\rm sat}$. These values are only
weakly impacted by our assumptions, such as the distances or the number of free
empirical parameters, provided they are taken within a reasonable range. We
also study the weak sensitivity of our results to the set of sources analyzed,
and we identify a group of sources that dominates the constraints. The
resulting masses and radii obtained are also discussed in the context of the
independent constraints from GW 170817 and its electromagnetic counterpart, AT
2017gfo.
