  We study current bounds on strong first-order phase transitions (PTs) along
the equation of state (EOS) of dense strongly interacting matter in neutron
stars, under the simplifying assumption that on either side of the PT the EOS
can be approximated by a simple polytropic form. We construct a large ensemble
of possible EOSs of this form, anchor them to chiral effective field theory
calculations at nuclear density and perturbative QCD at high densities, and
subject them to astrophysical constraints from high-mass pulsars and
gravitational-wave observations. Within this setup, we find that a PT permits
neutron-star solutions with larger radii, but only if the transition begins
below twice nuclear saturation density. We also identify a large parameter
space of allowed PTs currently unexplored by numerical-relativity studies.
Additionally, we locate a small region of parameter space allowing twin-star
solutions, though we find them to only marginally pass the current
astrophysical constraints. Finally, we find that sizeable cores of high-density
matter beyond the PT may be located in the centers of some stable neutron
stars, primarily those with larger masses.
