Cosmic Ray Magnetohydrodynamics: A New Two-Moment Framework with Numerical Implementation

Cosmic rays (CRs) play a pivotal role in various astrophysical systems, delivering feedback over a broad range of scales. However, modeling CR transport remains challenging due to its inherently multi-scale nature and complex microphysics. Recent advances in two-moment CR hydrodynamics have alleviated some of these challenges, improving understanding of CR feedback. Yet, current two-moment methods may not be able to directly incorporate all relevant CR transport processes, while the outcome of CR feedback sensitively depends on these underlying microphysics. Furthermore, numerical challenges persist, including instabilities from streaming terms and ambiguities in solver design for coupled CR-MHD systems. In this work, we develop a two-moment description for CR hydrodynamics from first principles. Beyond canonical CR streaming, our formulation accounts for CR pressure anisotropy and Alfvén waves propagating in both directions along the magnetic field, providing a general framework to incorporate more CR transport physics. We implement this framework as a new CR fluid module in the \textit{Athena}++ code, and validate it through a suite of benchmark tests. In particular, we derive the full dispersion relation of the two-moment CR-MHD system, identifying the CR-acoustic instability as well as other wave branches. These CR-MHD waves serve as rigorous benchmarks and also enable the use of realistic signal speeds in our Riemann solver. We propose a time step guideline to mitigate numerical instabilities arising from streaming source terms.

💡 Research Summary

This paper presents a comprehensive, first‑principles derivation and numerical implementation of a new two‑moment cosmic‑ray (CR) magnetohydrodynamics (MHD) framework that overcomes several limitations of existing CR transport models. Traditional one‑moment CR energy equations suffer from numerical instabilities because the streaming velocity, defined as v_st = −sgn(b·∇P_cr) v_A, is a discontinuous, sign‑changing term. Earlier two‑moment schemes (e.g., Jiang & Oh 2018; Thomas & Pfrommer 2019) improve stability by evolving a separate CR energy‑flux equation, but they typically assume Alfvén waves propagate only in one direction and do not treat CR pressure anisotropy. Consequently, they cannot capture the full range of microphysical processes such as the CR‑pressure‑anisotropy instability (CR‑PAI) that generates both forward and backward propagating Alfvén waves.

The authors start from the relativistic Boltzmann (Fokker‑Planck) equation for CRs and transform it into a frame moving with the component of the gas velocity perpendicular to the magnetic field (u_⊥). In this frame the CR distribution is gyrotropic, allowing a clean separation of pitch‑angle diffusion due to scattering off Alfvén waves. Crucially, they introduce separate scattering frequencies ν⁺ and ν⁻ for forward and backward waves, respectively, and retain the Doppler‑shifted wave velocities w_± = u ± v_A b̂. By expanding the collision term to leading order, they obtain a Fokker‑Planck operator that includes both pitch‑angle and momentum diffusion, preserving energy in the wave frame.

Taking the zeroth and first moments of the kinetic equation yields evolution equations for the CR energy density E_cr and the CR energy flux F_cr. The pressure tensor is split into parallel and perpendicular components (P_∥, P_⊥), giving rise to a pressure‑anisotropy term that appears as a non‑conservative source in the flux equation. The resulting system couples to the standard MHD equations through the total momentum and energy exchange terms, and it naturally incorporates the effects of both CR streaming (via v_st) and CR‑PAI (via pressure anisotropy).

A linear stability analysis is performed on the combined CR‑MHD system. By perturbing all variables with plane waves ∝ exp