CASTRO: A New Compressible Astrophysical Solver. III. Multigroup Radiation Hydrodynamics

We present a formulation for multigroup radiation hydrodynamics that is correct to order $O(v/c)$ using the comoving-frame approach and the flux-limited diffusion approximation. We describe a numerical algorithm for solving the system, implemented in the compressible astrophysics code, CASTRO. CASTRO uses an Eulerian grid with block-structured adaptive mesh refinement based on a nested hierarchy of logically-rectangular variable-sized grids with simultaneous refinement in both space and time. In our multigroup radiation solver, the system is split into three parts, one part that couples the radiation and fluid in a hyperbolic subsystem, another part that advects the radiation in frequency space, and a parabolic part that evolves radiation diffusion and source-sink terms. The hyperbolic subsystem and the frequency space advection are solved explicitly with high-order Godunov schemes, whereas the parabolic part is solved implicitly with a first-order backward Euler method. Our multigroup radiation solver works for both neutrino and photon radiation.

💡 Research Summary

The paper presents a comprehensive formulation and implementation of multigroup radiation hydrodynamics (MGRHD) within the compressible astrophysics code CASTRO, achieving accuracy to order O(v/c) by employing a comoving‑frame approach combined with the flux‑limited diffusion (FLD) approximation. The authors begin by motivating the need for a fully relativistic treatment of radiation‑matter coupling in high‑energy astrophysical phenomena such as core‑collapse supernovae, neutron‑star mergers, and radiative accretion flows, where fluid velocities can approach a few percent of the speed of light and where neutrino transport plays a decisive role. Traditional single‑group or static‑frame methods either neglect O(v/c) terms or treat them only approximately, leading to systematic errors in energy and momentum exchange.

The theoretical development starts from the comoving‑frame radiation transfer equation, expanding all terms to first order in v/c. Radiation quantities—energy density (E_g), flux (\mathbf{F}_g), and pressure tensor (\mathbf{P}_g)—are defined for each energy group (g). The FLD closure replaces the full angular dependence with a diffusion‑like relation (\mathbf{F}_g = -D_g \nabla E_g), where the diffusion coefficient (D_g) incorporates a flux limiter that smoothly transitions between optically thick diffusion and optically thin free‑streaming regimes, thereby respecting the causal limit ( |\mathbf{F}_g| \le cE_g).

A key contribution of the work is the operator‑splitting strategy that isolates three physically distinct subsystems:

1. Hyperbolic radiation‑fluid coupling – The combined set of mass, momentum, and total‑energy conservation equations includes radiation pressure and energy‑exchange source terms. This subsystem is advanced explicitly using a high‑order Godunov scheme (piecewise‑parabolic reconstruction with an HLLC or Roe Riemann solver). The explicit treatment enables accurate capture of shocks, contact discontinuities, and rapid advection of radiation energy with the fluid.

2. Frequency‑space advection – The redistribution of radiation energy among groups, driven by absorption, emission, and inelastic scattering, is cast as an advection equation in frequency (or energy) space. The same Godunov framework is applied, ensuring that the group‑to‑group transfer conserves total radiation energy to machine precision.

3. Parabolic diffusion and source‑sink – The diffusion term and stiff emission/absorption processes form a parabolic subsystem. Because the diffusion time scale can be orders of magnitude smaller than the hydrodynamic Courant limit, this part is integrated implicitly with a first‑order backward‑Euler method. The resulting linear system couples all groups and is solved with a preconditioned Krylov solver (GMRES) accelerated by a multigrid preconditioner, achieving scalable performance on large parallel machines.

CASTRO’s block‑structured adaptive mesh refinement (AMR) infrastructure is fully leveraged. Grids are logically rectangular, and each refinement level advances with its own time step (sub‑cycling), preserving second‑order temporal accuracy across levels. This allows the code to concentrate resolution where radiation‑matter interaction is strongest (e.g., near the neutrinosphere in a supernova core) while keeping the overall computational cost manageable.

The implementation is deliberately agnostic to the nature of the radiation; the same multigroup framework can be used for photon transport (e.g., radiative stellar envelopes) or for neutrino transport (e.g., electron‑type, mu‑type, and tau‑type neutrinos). Group‑dependent opacities, emissivities, and scattering kernels are read from input tables, and the code automatically assembles the corresponding source terms.

A suite of verification tests validates each component. The authors first present a radiation‑hydrodynamic wave test that isolates the O(v/c) terms, demonstrating that inclusion of these terms yields phase speeds and damping rates that match analytic predictions, whereas a purely Newtonian treatment deviates noticeably. Next, a multigroup diffusion test confirms that energy is conserved across groups during pure diffusion, even in the presence of strong opacity gradients. Finally, a full core‑collapse supernova simulation showcases the practical impact: the multigroup neutrino transport reproduces the expected deleptonization burst, shock revival timing, and post‑bounce neutrino luminosities more faithfully than a single‑group approximation, while still maintaining acceptable wall‑clock times thanks to the AMR‑driven load balancing and the implicit diffusion solver.

In conclusion, the paper demonstrates that a carefully split, high‑order explicit–implicit algorithm can deliver O(v/c)‑accurate multigroup radiation hydrodynamics within an AMR framework, handling both photon and neutrino radiation with comparable efficiency. The authors outline future extensions, including the adoption of higher‑order moment closures (e.g., M1) for anisotropic radiation fields, incorporation of general relativistic corrections for compact objects, and coupling to nuclear reaction networks for fully self‑consistent astrophysical explosion modeling.