  Bekenstein's Tensor-Vector-Scalar (TeVeS) theory has had considerable success
in explaining various phenomena without the need for dark matter. However, it
is difficult to observationally discern the differences between TeVeS and
predictions made within the Lambda-cold dark matter concordance model. This
implies that alternative tests are required that independently verify which
theory is correct. For this we turn to the strong-field regime of TeVeS. In
particular, we solve the spherically symmetric equations of hydrostatic
equilibrium for a perfect fluid with a realistic equation of state to build
models of neutron stars in TeVeS. We show that causality within the neutron
star is only maintained for certain cosmological values of the scalar field,
which allows us to put constraints on this value independently of cosmological
observations. We also discuss in detail the internal structure of neutron stars
and how each of the free parameters in the theory effects the overall size and
mass of the neutron stars. In particular, the radii of neutron stars in TeVeS
can significantly differ from those in General Relativity for certain values of
the vector field coupling, which allows us to also place extra constraints on
this parameter. Finally, we discuss future observations of neutron stars using
both the electromagnetic and gravitational wave spectrums that will allow for
tests of the appropriate theory of gravity.
