  The faithful inclusion of the effects of bulk viscosity induced by the
presence of chemical reactions is an important issue for simulations of
core-collapse supernovae, binary neutron star mergers, and neutron star
oscillations, where particle abundances are locally pushed out of chemical
equilibrium by rarefaction and compression of the fluid elements. In this work,
we discuss three different approaches that can be used to implement bulk
viscosity in general relativistic hydrodynamic simulations of neutron stars:
the exact multi-component reacting fluid, and two M\"uller-Israel-Stewart
theories, namely the second order Hiscock-Lindblom model and its linear limit,
the Maxwell-Cattaneo model. After discussing the theory behind the three
approaches, we specialize their dynamics equations to spherical symmetry in the
radial gauge-polar slicing (i.e., Schwarzschild) coordinates. We also discuss a
particular choice for the equation of state of the fluid and the associated
neutrino emission rates, which are used in a companion paper for the numerical
comparison of the three frameworks, and we obtain the effective sound speed for
the Hiscock-Lindblom theory in the non-linear regime.
