  The sample of neutron stars with a measured mass is growing quickly. With the
latest sample, we adopt both a flexible Gaussian mixture model and a Gaussian
plus Cauchy-Lorentz component model to infer the mass distribution of neutron
stars and use the Bayesian model selection to explore evidence for
multimodality and a sharp cutoff in the mass distribution. The two models yield
rather similar results. Consistent with previous studies, we find evidence for
a bimodal distribution together with a cutoff at a mass of $M_{\rm
max}=2.26_{-0.05}^{+0.12}M_\odot$ (68% credible interval; for the Gaussian
mixture model). If such a cutoff is interpreted as the maximum gravitational
mass of nonrotating cold neutron stars, the prospect of forming supramassive
remnants is found to be quite promising for the double neutron star mergers
with a total gravitational mass less than or equal to 2.7$M_\odot$ unless the
thermal pions could substantially soften the equation of state for the very hot
neutron star matter. These supramassive remnants have a typical kinetic
rotational energy of approximately $1-2\times 10^{53}$ ergs. Together with a
high neutron star merger rate approximately $10^{3}~{\rm Gpc^{-3}~yr^{-3}}$,
the neutron star mergers are expected to be significant sources of
EeV($10^{18}$eV) cosmic-ray protons.
