  We provide a systematic study of hybrid neutron star equations of state (EoS)
consisting of a relativistic density functional for the hadronic phase and a
covariant nonlocal Nambu--Jona-Lasinio (nlNJL) model to describe the color
superconducting quark matter phase. Changing the values of the two free
parameters, the dimensionless vector and diquark coupling strengths $\eta_V$
and $\eta_D$ results in a set of EoS with varying stiffness and deconfinement
onset. The favorable parameters are obtained from a systematic Bayesian
analysis for which the multi-messenger constraint on the neutron star radius at
$1.4~$M$_\odot$ and the combined mass-radius constraint for PSR J0740+6620 from
NICER experiment are used as the constraints. Additionally, the transition from
hadronic matter to deconfined quark matter is constrained to occur above
nuclear saturation density. Hybrid stars modeled with these favorable
parameters are compatible with the NICER results for the radius of the highest
known mass neutron star, PSR J0740+6620. Three new observations interesting for
neutron star phenomenology are reported: 1) We show that the constant sound
speed (CSS) EoS provides an excellent fit to that of the nlNJL model which
implies the squared speed of sound at high densities to be about $0.5$ for the
optimized parameters; 2) we give a simple functional form for the mapping
between the parameter spaces of these two models valid for the whole range of
relevant chemical potentials and 3) we observe that the special point property
of hybrid EoS based on CSS quark matter generalizes to a set of lines
consisting of special points when two EoS parameters are varied instead of one.
A lower limit for the maximum mass of hybrid stars as a function of the vector
coupling strength is obtained.
