We combine equation of state of dense matter up to twice nuclear saturation density ($n_{rm sat}=0.16, text{fm}^{-3}$) obtained using chiral effective field theory ($chi$EFT), and recent observations of neutron stars to gain insights about the high-density matter encountered in their cores. A key element in our study is the recent Bayesian analysis of correlated EFT truncation errors based on order-by-order calculations up to next-to-next-to-next-to-leading order in the $chi$EFT expansion. We refine the bounds on the maximum mass imposed by causality at high densities, and provide stringent limits on the maximum and minimum radii of $sim1.4,{rm M}_{odot}$ and $sim2.0,{rm M}_{odot}$ stars. Including $chi$EFT predictions from $n_{rm sat}$ to $2,n_{rm sat}$ reduces the permitted ranges of the radius of a $1.4,{rm M}_{odot}$ star, $R_{1.4}$, by $sim3.5, text{km}$. If observations indicate $R_{1.4}<11.2, text{km}$, our study implies that either the squared speed of sound $c^2_{s}>1/2$ for densities above $2,n_{rm sat}$, or that $chi$EFT breaks down below $2,n_{rm sat}$. We also comment on the nature of the secondary compact object in GW190814 with mass $simeq 2.6,{rm M}_{odot}$, and discuss the implications of massive neutron stars $>2.1 ,{rm M}_{odot},(2.6,{rm M}_{odot})$ in future radio and gravitational-wave searches. Some form of strongly interacting matter with $c^2_{s}>0.35, (0.55)$ must be realized in the cores of such massive neutron stars. In the absence of phase transitions below $2,n_{rm sat}$, the small tidal deformability inferred from GW170817 lends support for the relatively small pressure predicted by $chi$EFT for the baryon density $n_{rm B}$ in the range $1-2,n_{rm sat}$. Together they imply that the rapid stiffening required to support a high maximum mass should occur only when $n_{rm B} gtrsim 1.5-1.8,n_{rm sat}$.