We infer the mass distribution of neutron stars in binary systems using a flexible Gaussian mixture model and use Bayesian model selection to explore evidence for multi-modality and a sharp cut-off in the mass distribution . We find overwhelming evidence for a bimodal distribution , in agreement with previous literature , and report for the first time positive evidence for a sharp cut-off at a maximum neutron star mass . We measure the maximum mass to be 2.0 M _ { \odot } < m _ { \mathrm { max } } < 2.2 M _ { \odot } ( 68 % ) , 2.0 M _ { \odot } < m _ { \mathrm { max } } < 2.6 M _ { \odot } ( 90 % ) , and evidence for a cut-off is robust against the choice of model for the mass distribution and to removing the most extreme ( highest mass ) neutron stars from the dataset . If this sharp cut-off is interpreted as the maximum stable neutron star mass allowed by the equation of state of dense matter , our measurement puts constraints on the equation of state . For a set of realistic equations of state that support > 2 M _ { \odot } neutron stars , our inference of m _ { \mathrm { max } } is able to distinguish between models at odds ratios of up to 12 : 1 , whilst under a flexible piecewise polytropic equation of state model our maximum mass measurement improves constraints on the pressure at 3 - 7 \times the nuclear saturation density by \sim 30 - 50 \% compared to simply requiring m _ { \mathrm { max } } > 2 M _ { \odot } . We obtain a lower bound on the maximum sound speed attained inside the neutron star of c _ { s } ^ { \mathrm { max } } > 0.63 c ( 99.8 % ) , ruling out c _ { s } ^ { \mathrm { max } } < c / \sqrt { 3 } at high significance . Our constraints on the maximum neutron star mass strengthen the case for neutron star-neutron star mergers as the primary source of short gamma-ray bursts .