Neutron stars provide a unique opportunity to study strongly interacting
matter under extreme density conditions. The intricacies of matter inside
neutron stars and their equation of state are not directly visible, but
determine bulk properties, such as mass and radius, which affect the star's
thermal X-ray emissions. However, the telescope spectra of these emissions are
also affected by the stellar distance, hydrogen column, and effective surface
temperature, which are not always well-constrained. Uncertainties on these
nuisance parameters must be accounted for when making a robust estimation of
the equation of state. In this study, we develop a novel methodology that, for
the first time, can infer the full posterior distribution of both the equation
of state and nuisance parameters directly from telescope observations. This
method relies on the use of neural likelihood estimation, in which normalizing
flows use samples of simulated telescope data to learn the likelihood of the
neutron star spectra as a function of these parameters, coupled with
Hamiltonian Monte Carlo methods to efficiently sample from the corresponding
posterior distribution. Our approach surpasses the accuracy of previous
methods, improves the interpretability of the results by providing access to
the full posterior distribution, and naturally scales to a growing number of
neutron star observations expected in the coming years.