An Explicit M1 Radiation-hydrodynamics Scheme for 3D Protostellar Evolution

We present a radiation-hydrodynamics (RHD) scheme that enables 3D simulations resolving both protostellar interiors and their surrounding accretion flows within a single framework, to clarify how a protostar evolves while interacting with the accretion flow. The method builds on an explicit two-moment M1 closure scheme with a reduced speed of light approximation (RSLA) for massively parallel computation. Our scheme introduces a complementary non-RSLA radiation component that dominates in optically thick regions. This hybrid treatment restores physical energy conservation inside protostars, which would otherwise be violated under the RSLA, while retaining the advantage of large time steps. To overcome the limitation of the conventional M1 closure in solving radiative transfer in extremely optically thick regions inside protostars and across steep optical-depth gradients near their surfaces, we incorporate the optical-depth information of neighboring cells into the radiative transfer calculation. We further evolve photon-number densities in addition to radiation energy densities to reconstruct an effective local spectrum on the fly without resorting to costly multi-frequency transport. We implement this scheme in the adaptive mesh refinement code SFUMATO and verify its validity through a series of test calculations. As an application, we follow the early evolution of a massive protostar formed at high redshift, within a full cosmological context. The results reveal a continuous structure connecting the swollen protostar and its surrounding disk, which cannot be captured in conventional 1D models. This RHD scheme opens a path to studies of protostellar evolution and its interaction with the accretion flow in realistic 3D environments.

💡 Research Summary

This paper introduces a novel radiation‑hydrodynamics (RHD) scheme that enables fully three‑dimensional simulations capable of resolving both the interior of a protostar and its surrounding accretion flow within a single computational framework. The core of the method is an explicit two‑moment M1 closure combined with the reduced speed of light approximation (RSLA), which allows large time steps and excellent scalability on massively parallel architectures. However, the authors identify three critical shortcomings of a plain M1‑RSLA approach: (i) violation of total energy conservation in optically thick regions because the reduced light speed artificially suppresses gas‑radiation energy exchange; (ii) excessive numerical diffusion inside the protostar and inaccurate flux transmission across the steep optical‑depth gradient at the stellar surface; and (iii) the lack of spectral information when only a single radiation energy density is evolved.

To overcome these issues, the authors construct a hybrid scheme. First, they split the radiation field into two components: a conventional RSLA component that dominates in optically thin gas, and a “non‑RSLA” component that is evolved with the true speed of light in cells whose optical depth exceeds a prescribed threshold. This restores exact energy conservation inside the protostar while preserving the efficiency of RSLA elsewhere. Second, they augment the M1 flux calculation with a neighbor‑cell optical‑depth weighting, effectively incorporating information about the surrounding opacity distribution. This mitigates artificial leakage of radiation across sharp surface layers and reduces numerical diffusion in the deep interior. Third, they introduce a photon‑number density variable alongside the radiation energy density. By taking the ratio of the two, the local mean photon energy (and thus an approximate spectrum) can be reconstructed on the fly without the computational expense of multi‑frequency transport. The opacities (Planck, energy‑mean, and flux‑mean) are then updated dynamically using this locally estimated spectrum.

The scheme is implemented in the adaptive‑mesh‑refinement (AMR) code SFUMATO, which already supports self‑gravity, magnetohydrodynamics, and sink particles. The authors describe how refinement criteria are based on both the Jeans length and the local optical depth, ensuring that the protostellar interior is adequately resolved. Time integration follows the explicit RSLA CFL condition, while the non‑RSLA component is sub‑cycled where necessary. Parallelization uses MPI with OpenMP threading, and GPU kernels are provided for the most expensive radiation updates.

A suite of verification tests demonstrates the method’s accuracy and robustness. In a static diffusion test the numerical solution reproduces the analytic diffusion profile with sub‑percent errors even when the cell optical depth exceeds 10⁴. A shadow test shows that the hybrid scheme preserves sharp shadows behind an opaque obstacle, unlike a pure RSLA run that smears the shadow. Radiative shock tests confirm that both the shock structure and post‑shock temperature match analytical expectations, and that total energy is conserved to better than 10⁻⁸ relative error when the non‑RSLA component is active.

As a scientific application, the authors simulate the formation of a massive protostar at high redshift (z≈20) within a fully cosmological context. The initial conditions are extracted from a ΛCDM simulation and consist of a metal‑free minihalo collapsing under its own gravity. The simulation follows the protostar from its birth up to ∼10 M⊙. The results reveal a continuous, “puffy” protostellar envelope that smoothly connects to a rotating accretion disk. The hybrid radiation treatment captures the rapid transition from optically thick interior to optically thin exterior without artificial flux suppression, leading to a surface temperature evolution that differs markedly from 1‑D spherical models. Moreover, the photon‑number density method provides a time‑varying estimate of the emergent spectrum, showing an increase in UV and soft X‑ray output as the protostar grows, with direct implications for the ionization state of the surrounding primordial gas.

The discussion acknowledges remaining limitations: the choice of the optical‑depth threshold for the non‑RSLA component can affect results; the photon‑number approach yields only a mean photon energy and cannot resolve line features or non‑thermal spectra; and the M1 closure still struggles with intersecting radiation beams, suggesting future hybridization with Variable Eddington Tensor or Monte‑Carlo methods.

In summary, the paper delivers a practical, energy‑conserving, and spectrally aware RHD algorithm that bridges the gap between large‑scale star‑formation simulations and detailed protostellar evolution. By enabling fully resolved 3‑D protostellar interiors, the method opens new avenues for studying feedback, accretion geometry, and early‑Universe star formation with unprecedented fidelity.