Comparing energy and entropy formulations for cosmic ray hydrodynamics

Cosmic rays (CRs) play an important role in many astrophysical systems. Acting on plasma scales to galactic environments, CRs are usually modelled 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 magnetohydrodynamic (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 set-up, 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.