The millisecond pulsar mass distribution: Evidence for bimodality and constraints on the maximum neutron star mass

Abstract: 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.