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

TL;DR

This work develops a general, first‑principles two‑moment framework for cosmic ray hydrodynamics that evolves CR energy density, flux, and pressure with a finite CR pressure anisotropy and bidirectional Alfvén waves, integrated with MHD in Athena++. By deriving the full CR–MHD dispersion relation and validating the code against diverse benchmarks, the authors demonstrate robust wave behavior, including CR‑modified acoustic modes, and provide a practical stability criterion for time stepping. The framework unifies previously separate two‑moment approaches, extends them to arbitrary wave configurations, and offers subgrid prescriptions for wave–driven scattering and anisotropy. The implementation enables flexible, physically motivated CR transport modeling across astrophysical environments, with potential extensions to multi‑energy CRs and dynamic wave evolution. Overall, the work advances predictive CR feedback modeling by coupling a comprehensive transport physics module to high‑fidelity MHD simulations, supported by rigorous numerical benchmarks.

Abstract

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.

Figures (12)

• Figure 1: Numerical results of 1D streaming tests from \ref{['subsec:1Dstreaming']}. Top left:$\mathcal{E}_{\rm cr}$ and $F_{\rm cr}$ in a static background gas $(u=0)$. Colored lines denote snapshots at increasing simulation times, while dashed curves show analytical solutions, validating agreement between numerical and theoretical results. Top right: Identical setup with a uniformly advected background gas $(u=1)$, introducing additional $F_{\rm cr}$ from gas motion. The $\mathcal{E}_{\rm cr}$ profile drifts at $4u/3$, consistent with expectations. Bottom: Spatial convergence of $L1$ errors $(\sum_{i=1}^{N_x}|\mathcal{E}_{\rm cr,i}-\mathcal{E}_{\rm cr,theory}|/N_x)$ at $t=0.3$. Errors roughly scale as $\propto N_x^{-1.34}$.
• Figure 2: Numerical results of the 1D streaming test using the same setup as in \ref{['fig:1Dstream']}, except for elevated scattering coefficients $\widetilde{\sigma}_{0,\rm st}$ to illustrate the performance of our stability criterion. Time steps $\Delta t$ progressively decrease from left to right, as indicated by the titles. We see numerical instability and deviation from analytical solutions persist until the stability condition $\Delta t<\sqrt{3}\Delta x/(1+4\widetilde{\sigma}_{0,\rm st})V_m$ is satisfied in the rightmost column.
• Figure 3: Numerical results of 1D streaming tests from \ref{['subsec:1Dnlld']}. This is similar to \ref{['fig:1Dstream']} but the damping mechanism follows the form expected for non-linear Landau damping, with the scattering coefficients defined in \ref{['eq:sigma_nll']}. Dashed lines denote our theoretical predictions, which align well with the simulation data at the propagating fronts. Deviations occur in the central plateau, where $\nabla P_{\rm cr}\rightarrow 0$ violates our theoretical assumption of a non-negligible $\nabla P_{\rm cr}$. As a result $\mathcal{E}_{\rm cr}$ and $F_{\rm cr}$ smoothly transition to the propagating fronts.
• Figure 4: Numerical results of 1D diffusion tests described in \ref{['sec:diff_test']}. Top panels: Evolution of $\mathcal{E}_{\rm cr}$ profiles in a static MHD background. Bottom panels: Evolution of $\mathcal{E}_{\rm cr}$ profiles in a uniformly moving MHD background (MHD velocity $u=-1$).
• Figure 5: 2D tests for the anisotropic transport of $\mathcal{E}_{\rm cr}$ along magnetic field lines on Cartesian grids, as described in \ref{['sec:aniso_flux']}. An initial patch with high CR energy density streams along circular background magnetic fields. The bottom right panel shows the analytical solution while other three panels demonstrate numerical results with different spatial resolutions at the same time.