  Cosmic rays (CRs) play an important role in many astrophysical systems.
Acting on plasma scales to galactic environments, CRs are usually modeled as a
fluid, using the CR energy density as the evolving quantity. This method comes
with the flaw that the corresponding CR evolution equation is not in
conservative form as it contains an adiabatic source term that couples CRs to
the thermal gas. In the absence of non-adiabatic changes, instead evolving the
CR entropy density is a physically equivalent option that avoids this potential
numerical inconsistency. In this work, we study both approaches for evolving
CRs in the context of magneto-hydrodynamic (MHD) simulations using the
massively parallel moving-mesh code AREPO. We investigate the performance of
both methods in a sequence of shock-tube tests with various resolutions and
shock Mach numbers. We find that the entropy-conserving scheme performs best
for the idealized case of purely adiabatic CRs across the shock while both
approaches yield similar results at lower resolution. In this setup, both
schemes operate well and almost independently of the shock Mach number. Taking
active CR acceleration at the shock into account, the energy-based method
proves to be numerically much more stable and significantly more accurate in
determining the shock velocity, in particular at low resolution, which is more
typical for astrophysical large-scale simulations. For a more realistic
application, we simulate the formation of several isolated galaxies at
different halo masses and find that both numerical methods yield almost
identical results with differences far below common astrophysical
uncertainties.
