  We develop a formalism to model slowly rotating superfluid Newtonian neutron
stars. A two-fluid model is used to describe the matter, where one fluid is the
superfluid neutrons that are believed to exist in the inner crust and core of
mature neutron stars, while the other fluid is a charge neutral conglomerate of
the remaining constituents (crust nuclei, core superconducting protons,
electrons, etc). We include the entrainment effect, which is a non-dissipative
interaction between the two fluids whereby a momentum induced in one of the
fluids will cause part of the mass of the other fluid to be carried along. The
equations that describe rotational equilibria (i.e. axisymmetric and stationary
configurations) are approximated using the slow-rotation approximation; an
expansion in terms of the rotation rates of the two fluids where only terms up
to second-order are kept. Our formalism allows the neutrons to rotate at a rate
different from that of the charged constituents. For a particular equation of
state that is quadratic in the two mass-densities and relative velocities of
the fluids, we find an analytic solution to the slow-rotation equations. This
result provides an elegant generalisation to the two-fluid problem of the
standard expressions for the one-fluid polytrope ${\cal E} \propto \rho^2$. The
model equation of state includes entrainment and is general enough to allow for
realistic values for, say, mass and radius of the star. It also includes a term
that can be related to ``symmetry energy'' terms that appear in more realistic
equations of state. We use the analytic solution to explore how relative
rotation between the two fluids, the ``symmetry energy'' term, and entrainment
affect the neutron star's local distribution of particles, as well as global
quantities as the Kepler limit, ellipticity, and moments of inertia.
