  The high-density behavior of nuclear matter is analyzed within a relativistic
mean-field description with non-linear meson interactions. To assess the model
parameters and their output, a Bayesian inference technique is used. The
Bayesian setup is limited only by a few nuclear saturation properties, the
neutron star maximum mass larger than 2 M$_\odot$, and the low-density pure
neutron matter equation of state (EOS) produced by an accurate N$^3$LO
calculation in chiral effective field theory. Depending on the strength of the
non-linear scalar vector field contribution, we have found three distinct
classes of EOSs, each one correlated to different star properties
distributions. If the non-linear vector field contribution is absent, the
gravitational maximum mass and the sound velocity at high densities are the
greatest. However, it also gives the smallest speed of sound at densities below
three times saturation density. On the other hand, models with the strongest
non-linear vector field contribution predict the largest radii and tidal
deformabilities for 1.4 M$_\odot$ stars, together with the smallest mass for
the onset of the nucleonic direct Urca processes and the smallest central
baryonic densities for the maximum mass configuration. These models have the
largest speed of sound below three times saturation density, but the smallest
at high densities, in particular, above four times saturation density the speed
of sound decreases approaching approximately $\sqrt{0.4}c$ at the center of the
maximum mass star. On the contrary, a weak non-linear vector contribution gives
a monotonically increasing speed of sound. A 2.75 M$_\odot$ NS maximum mass was
obtained in the tail of the posterior with a weak non-linear vector field
interaction. We found that pQCD favors models with a large non-linear vector
field contribution or hyperons.
