  Motivated by the recent discovery of a compact object with mass in the range
$2.5{-}2.67\, \mathrm{M}_\odot$ 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 2 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.
