  Neutron stars provide an excellent laboratory for physics under the most
extreme conditions. Up to now, models of axisymmetric, stationary,
differentially rotating neutron stars were constructed under the strong
assumption of barotropicity, where a one-to-one relation between all
thermodynamic quantities exists. This implies that the specific angular
momentum of a matter element depends only on its angular velocity. The physical
conditions in the early stages of neutron stars, however, are determined by
their violent birth processes, typically a supernova or in some cases the
merger of two neutron stars, and detailed numerical models show that the
resulting stars are by no means barotropic. Here, we construct models for
stationary, differentially rotating, non-barotropic neutron stars, where the
equation of state and the specific angular momentum depend on more than one
independent variable. We show that the potential formulation of the
relativistic Euler equation can be extended to the non-barotropic case, which,
to the best of our knowledge, is a new result even for the Newtonian case. We
implement the new method into the XNS code and construct equilibrium
configurations for non-barotropic equations of state. We scrutinize the
resulting configurations by evolving them dynamically with the numerical
relativity code BAM, thereby demonstrating that the new method indeed produces
stationary, differentially rotating, non-barotropic neutron star
configurations.
