Abstract Motivated by the recent discovery of a compact object with mass in the range 2.5–2.67 M ⊙ in the binary merger GW190814, we revisit the question of the maximum mass of neutron stars (NSs). We use a Markov Chain Monte Carlo approach to generate about two million phenomenological equations of state with and without first-order phase transitions. We fix the crust equation of state and only assume causality at higher densities. We show how a strict upper bound on the maximum NS mass can be inferred from upcoming observation of NS radii and masses. The derived upper bounds depend only on relativity and causality, so it is not affected by nuclear physics uncertainties. We show how a lower limit on the maximum mass of NSs, in combination with upcoming measurements of NS radii by LIGO/Virgo and NICER, would constrain the equation of state of dense matter. Finally, we discuss the implications for GW190814.