Abstract 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 multifrequency 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.

Abstract Direct collapse of supermassive stars is a possible pathway to form supermassive black hole seeds at high redshifts. Whereas previous three-dimensional (3D) simulations demonstrate that supermassive stars form via rapid mass accretion, those resolving the stellar interior have been limited. Here, we report 3D radiation-hydrodynamic (RHD) simulations following the evolution of rapidly accreting protostars resolving the stellar interior. We use an adaptive mesh refinement code with our newly developed RHD solver employing an explicit M1 closure method. We follow the early evolution until the stellar mass reaches ∼10 M_⊙ from two different initial configurations of spherical and turbulent clouds. We demonstrate that, in both cases, a swollen protostar whose radius is 100–1000 R_⊙ appears, as predicted by the stellar evolution calculations. Its effective temperature remains a few thousand Kelvin, and the radiative feedback by ionizing photons is too weak to disturb the accretion flow up to the epoch examined in this work. In the turbulent case, the protostar rotates rapidly at more than 0.4 times the Keplerian velocity owing to the angular momentum provided by the initial turbulence. The protostar approximates an oblate spheroid, and its equatorial radius is more than twice the polar radius. Our results suggest that we need to consider the rapid stellar rotation to elucidate the realistic 3D protostellar evolution in the supermassive star formation.