  The mass function of neutron stars (NSs) contains information about the late
evolution of massive stars, the supernova explosion mechanism, and the
equation-of-state of cold, nuclear matter beyond the nuclear saturation
density. A number of recent NS mass measurements in binary millisecond pulsar
(MSP) systems increase the fraction of massive NSs (with $M > 1.8$ M$_{\odot}$)
to $\sim 20\% $ of the observed population. In light of these results, we
employ a Bayesian framework to revisit the MSP mass distribution. We find that
a single Gaussian model does not sufficiently describe the observed population.
We test alternative empirical models and infer that the MSP mass distribution
is strongly asymmetric. The diversity in spin and orbital properties of
high-mass NSs suggests that this is most likely not a result of the recycling
process, but rather reflects differences in the NS birth masses. The asymmetry
is best accounted for by a bimodal distribution with a low mass component
centred at $1.393_{-0.029}^{+0.031}$ M$_{\odot}$ and dispersed by
$0.064_{-0.025}^{+0.064}$ M$_{\odot}$, and a high-mass component with a mean of
$1.807_{-0.132}^{+0.081}$ and a dispersion of $0.177_{-0.072}^{+0.115}$
M$_{\odot}$. We also establish a lower limit of $M_{max} \ge 2.018$ M$_{\odot}$
at 98% C.L. for the maximum NS mass, from the absence of a high-mass truncation
in the observed masses. Using our inferred model, we find that the measurement
of 350 MSP masses, expected after the conclusion of pulsar surveys with the
Square-Kilometre Array, can result in a precise localization of a maximum mass
up to 2.15 M$_{\odot}$, with a 5% accuracy. Finally, we identify possible
massive NSs within the known pulsar population and discuss birth masses of
MSPs.
