  The interior of neutron stars contains matter at the highest densities
realized in our Universe. Interestingly, theoretical studies of dense matter,
in combination with the existence of two solar mass neutron stars, indicate
that the speed of sound $c_s$ has to increase to values well above the
conformal limit ($c_s^2\sim 1/3$) before decreasing again at higher densities.
The decrease could be explained by either a strong first-order phase transition
or a cross-over transition from hadronic to quark matter. The latter scenario
leads to a pronounced peak in the speed of sound reaching values above the
conformal limit, naturally explaining the inferred behavior. In this work, we
use the Nuclear-Physics Multi-Messenger Astrophysics framework \textsc{NMMA} to
compare predictions of the quarkyonic matter model with astrophysical
observations of neutron stars, with the goal of constraining model parameters.
Assuming quarkyonic matter to be realized within neutron stars, we find that
there can be a significant amount of quarks inside the core of neutron stars
with masses in the two solar mass range, amounting to up to $\sim 0.13M_\odot$,
contributing $\sim 5.9\%$ of the total mass. Furthermore, for the quarkyonic
matter model investigated here, the radius of a $1.4M_\odot$ neutron star would
be $13.44^{+1.69}_{-1.54} (13.54^{+1.02}_{-1.04})$ km, at $95\%$ credibility,
without (with) the inclusion of AT2017gfo.
