Maximal mass of neutron stars constrained by neutron star observations

We investigate constraints on the high-density equation of state (EOS) of neutron star matter by analyzing the probability distributions of the endpoints of mass-radius M(R) sequences within a Bayesian weighting framework. Starting from two representative hadronic baseline EOSs, SFHo and DD2, matched at higher densities to an extended linear sigma model description and constrained to approach perturbative QCD (pQCD) results, we construct families of causal hybrid EOSs spanning a broad range of stiffness at supranuclear densities. Observational constraints from the binary neutron-star merger GW170817, mass-radius measurements from the Neutron Star Interior Composition Explorer (NICER), and candidate low-mass and mass-gap compact objects are incorporated through Bayesian likelihood weighting. This approach allows us to determine probability distributions for the maximum neutron-star mass M$_{\rm TOV}$ and the corresponding radius R$_{\rm TOV}$, i.e., the endpoints of the M(R) sequences. We find that the maximum-mass distributions are largely determined by observational constraints and show only weak sensitivity to the choice of baseline EOS, favoring values around 2.2-2.3 M$_\odot$ when the most robust constraints are applied. In contrast, the corresponding radius distributions exhibit a stronger dependence on the underlying hadronic EOS, with typical preferred values near $12\pm 1$ km. Additional tidal-deformability constraints further restrict the allowed parameter space and disfavor very stiff EOS realizations when interpreted together with the possible mass-gap neutron-star candidate. Our results demonstrate that endpoint distributions of M(R) sequences provide a sensitive and complementary diagnostic for constraining the high-density behavior of the neutron-star EOS within a multimessenger Bayesian framework.