We combine equation of state of dense matter up to twice nuclear saturation density (n_sat=0.16fm⁻³) obtained using chiral effective field theory (χ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 χ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 ∼1.4M_⊙ and ∼2.0M_⊙ stars. Including χEFT predictions from n_sat to 2n_sat reduces the permitted ranges of the radius of a 1.4M_⊙ star, R1.4, by ∼3.5km. If observations indicate R1.4<11.2km, our study implies that either the squared speed of sound for densities above 2n_sat, or that χEFT breaks down below 2n_sat. We also comment on the nature of the secondary compact object in GW190814 with mass ≃2.6M_⊙, and discuss the implications of massive neutron stars >2.1M_⊙(2.6M_⊙) in future radio and gravitational-wave searches. Some form of strongly interacting matter with >0.35(0.55) must be realized in the cores of such massive neutron stars. In the absence of phase transitions below 2n_sat, the small tidal deformability inferred from GW170817 lends support for the relatively small pressure predicted by χEFT for the baryon density nB in the range 1−2n_sat. Together they imply that the rapid stiffening required to support a high maximum mass should occur only when n_B≳1.5−1.8n_sat.