  We perform a Bayesian analysis of the maximum mass $M_{\rm 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 \nicer. 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" (CSS) parameterization. The mass distribution measured
for the $2.14 \,{\rm M}_{\odot}$ pulsar, MSP J0740+6620, is used as the lower
limit on $M_{\rm TOV}$. We find the most probable values of the hybrid star
maximum mass are $M_{\rm TOV}=2.36^{+0.49}_{-0.26}\,{\rm M}_{\odot}$
($2.39^{+0.47}_{-0.28}\,{\rm M}_{\odot}$) for QMF (DD2), with an absolute upper
bound around $2.85\,{\rm M}_{\odot}$, 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 post-merger product of GW170817, short gamma-ray bursts, and other likely
binary neutron star mergers.
