  The merger remnant of a binary neutron star coalescence is initially strongly
differentially rotating. Some properties of these remnants can be accurately
modeled through building equilibrium neutron star models. In the present paper,
we study how a modification of general relativity, namely scalar-tensor theory
with a massive scalar field, will alter the picture. In contrast to previous
studies, we implement a realistic phenomenological differential rotational law
which allows for neutron star models to attain maximal angular velocity away
from the center. We find that solutions with much higher masses and angular
momenta exist in scalar-tensor theory compared to general relativity. They keep
their quasi-spherical energy-density distribution for significantly higher
values of the angular momentum before transitioning to quasi-toroidal models,
in contrast to pure general relativity. Constructing such neutron star
solutions is the first step to our final goal that is studying how
scalarization alters the stability and gravitational wave emission of
post-merger remnants.
