The equation of state (EOS) of cold dense matter is a central open problem in nuclear astrophysics. Its inference is hindered by the lack of ab initio control above about twice nuclear saturation density, requiring extrapolation. Parametric schemes such as piecewise polytropes (PP) are efficient but restrictive, while nonparametric approaches like Gaussian processes (GP) allow more flexibility at the cost of larger prior volumes. We extend our hybrid EOS framework by replacing the high-density polytropic extension with a GP representation of the squared sound speed, anchored at low densities by the SLy crust EOS and a nuclear meta-model constrained by χEFT and laboratory data. Using hierarchical Bayesian analysis, we jointly constrain the EOS and neutron star mass distribution with multi-messenger observations, including NICER radii, GW170817 and GW190425 tidal deformabilities, pulsar masses, and neutron skin experiments. We examine four scenarios defined by high-density extrapolation (PP vs. \ GP) and hotspot geometry in the NICER modeling of PSR~J0030+0451 (ST+PDT vs. \ PDT-U). GP extrapolations generally yield softer EOS posteriors with broader uncertainties. Hotspot assumptions also play an important role, shifting inferred mass--radius relations. Bayesian evidence strongly favors the ST+PDT geometry over PDT-U under both extrapolations, while GP is mildly preferred over PP. These results underscore the impact of observational modeling and EOS extrapolation on neutron star inferences, and show that a GP-based extension offers a robust way to quantify systematic uncertainties in high-density matter.
B. Biswas (Sun,) studied this question.