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 } \geq 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 .