  Deviations from the predictions of general relativity due to energy-momentum
squared gravity (EMSG) are expected to become pronounced in the high density
cores of neutron stars. We derive the hydrostatic equilibrium equations in EMSG
and solve them numerically to obtain the neutron star mass-radius relations for
four different realistic equations of state. We use the existing observational
measurements of the masses and radii of neutron stars to constrain the free
parameter, $\alpha ,$ that characterizes the coupling between matter and
spacetime in EMSG. We show that $-10^{-38}\,\mathrm{cm^{3}/erg}<\alpha
<+10^{-37}\,\mathrm{cm^{3}/erg}$. Under this constraint, we discuss what
contributions EMSG can provide to the physics of neutron stars, in particular,
their relevance to the so called \textit{hyperon puzzle} in neutron stars. We
also discuss how EMSG alters the dynamics of the early universe from the
predictions of the standard cosmological model. We show that EMSG leaves the
standard cosmology safely unaltered back to $t\sim 10^{-4}$ seconds at which
the energy density of the universe is $\sim 10^{34}\,\mathrm{erg\,cm^{-3}}$.
