Comparing energy and entropy formulations for cosmic ray hydrodynamics

TL;DR

This paper compares energy and entropy formulations for cosmic ray hydrodynamics in MHD simulations, finding that entropy-based methods are more stable and accurate for shock velocity, while both methods agree in galaxy formation scenarios.

Contribution

It introduces and evaluates the performance of energy versus entropy formulations for cosmic ray evolution in MHD simulations, highlighting the advantages of entropy-based schemes.

Findings

Entropy-conserving scheme performs best for adiabatic CRs across shocks.

Both methods yield similar results at low resolution in idealized tests.

Energy-based method is more stable and accurate for shock velocity at low resolution.

Abstract

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.