Abstract We perform a Bayesian analysis of the maximum mass M TOV of neutron stars with a quark core, incorporating the observational data from tidal deformability of the GW170817 binary neutron star merger as detected by LIGO/Virgo and the mass and radius of PSR J0030+0451 as detected by the Neutron Star Interior Composition Explorer. The analysis is performed under the assumption that the hadron–quark phase transition is of first order, where the low-density hadronic matter described in a unified manner by the soft QMF or the stiff DD2 equation of state (EOS) transforms into a high-density phase of quark matter modeled by the generic “constant-sound-speed” parameterization. The mass distribution measured for the 2.14 M ⊙ pulsar MSP J0740+6620 is used as the lower limit on M TOV . We find the most probable values of the hybrid star maximum mass are M TOV = 2.36 − 0.26 + 0.49 M ⊙ ( 2.39 − 0.28 + 0.47 M ⊙ ) for QMF (DD2), with an absolute upper bound around 2.85 M ⊙ , to the 90% posterior credible level. Such results appear robust with respect to the uncertainties in the hadronic EOS. We also discuss astrophysical implications of this result, especially on the postmerger product of GW170817, short gamma-ray bursts, and other likely binary neutron star mergers.