Simulating How Planets Are Born: New Guide for Radiation Hydrodynamics of Disks

Around every young star lies a swirling disk of gas and dust, the raw material from which planets are assembled. Understanding how these disks evolve, how turbulence stirs their contents, and how planets carve their first homes into the gas has long been one of the central challenges of modern astrophysics. Now, a comprehensive review published in Living Reviews in Computational Astrophysics offers researchers a detailed guidebook on how to build and test numerical simulations of protoplanetary disks, with a particular focus on radiation hydrodynamics, the coupled treatment of gas motion and the transport of light and heat.

The review, led by Hubert Klahr, Hans Baehr, Julio David Melon Fuksman, and Thomas Pfeil of the Max Planck Institute for Solar System Research, arrives at a moment when observational facilities such as the Atacama Large Millimeter/submillimeter Array are delivering breathtaking images of planet-forming disks. Rings, gaps, spiral arms, and hints of long-lived vortices appear routinely in the data, but interpreting these structures requires models that faithfully capture the underlying physics. The authors argue that the first step is deceptively simple: before adding magnetic fields, self-gravity, or dust feedback, a simulation must correctly reproduce the pure hydrodynamic behavior of an irradiated disk, and that behavior is far from trivial.

Disks around young stars are born as byproducts of star formation, a buffer for the excess angular momentum of a collapsing molecular cloud core. During the first tens of thousands of years, the star acquires most of its mass through the disk, a phase in which the disk’s own self-gravity supplies the main torque driving accretion. But after this main accretion phase, enough material remains to build planetary systems over the following ten million years, until winds and photoevaporation strip the gas away. Crucially, for planet formation, an overly vigorous accretion flow toward the star can wash away nascent planets, so understanding when and how disks settle into a calmer state matters enormously.

One of the review’s central themes is baroclinicity, the misalignment between surfaces of constant pressure and constant density in a disk. Because stellar irradiation heats the disk surface while viscous dissipation may warm the midplane, disks develop both radial and vertical temperature gradients. In such a configuration, pressure and density contours are inclined with respect to one another, a situation familiar from Earth’s atmosphere and oceans that also drives instabilities in disks. These effects were long missed in simulations because they are too weak to emerge at low resolution or with overly dissipative numerical schemes. Only recently, with increased computing power, have researchers been able to confirm the analytic predictions of thermal baroclinic instability theory with full numerical experiments.

The review catalogues the family of instabilities that radiation hydrodynamic simulations must capture. The vertical shear instability, or VSI, arises from the vertical gradient of rotation in a baroclinic disk and operates best when thermal relaxation is extremely fast, since rapid cooling prevents the stable vertical stratification from suppressing the unstable shear. The convective overstability, in contrast, is strongest when the thermal relaxation time is comparable to the orbital period; it grows from epicyclic oscillations of gas parcels displaced radially in a weakly convective environment. The Goldreich-Schubert-Fricke instability, a close cousin of the VSI inherited from the theory of rotating stars, operates alongside the convective overstability, and recent work has shown that any disk unstable to one is unstable to the other, because both depend equally on the disk’s baroclinicity. A third mechanism, the subcritical baroclinic instability, generates long-lived anticyclonic vortices from radial entropy gradients combined with thermal relaxation, though it lacks a linear growth-rate prediction and is therefore harder to use as a code benchmark.

To help researchers validate their codes, the authors lay out a systematic testing procedure. The method begins by constructing an equilibrium disk model, typically through a so-called 1+1-dimensional calculation that solves for vertical hydrostatic balance at each radius while enforcing energy conservation with flux-limited diffusion and physically motivated dust opacities. This equilibrium serves both as an initial condition for multidimensional simulations and as the basis for perturbation theory: analytic linear analysis yields predicted growth rates for unstable modes. A code that reproduces those growth rates in the linear regime can then be trusted to explore the nonlinear turbulence that follows. If a code fails to match the predicted growth, the review notes, this signals insufficient resolution or a numerically dissipative scheme rather than a failure of the underlying theory.

The practical details of such tests are demanding. In local axisymmetric simulations carved from a global disk model, the authors demonstrate that resolving roughly 256 cells per pressure scale height is required to reproduce growth rates down to ten thousandths of the orbital frequency. Numerical schemes matter as well: high-order reconstruction methods combined with accurate Riemann solvers recover the predicted linear growth, whereas more diffusive approximate solvers can suppress instability entirely. Even the handling of cooling, implemented through thermal relaxation of the pressure toward an equilibrium value, must be treated carefully, with operator splitting modified so that relaxation times shorter than the dynamical step remain stable. These benchmarks, the authors argue, should become standard practice for any group embarking on radiation hydrodynamic disk simulations.

Beyond linear tests, the review surveys how different radiative transfer approximations shape the outcomes of full three-dimensional simulations. Flux-limited diffusion, the workhorse of earlier decades, is computationally efficient and accurate in optically thick regions but smears out shadows and introduces unphysical diffusion where radiation streams freely. The M1 two-moment method, implemented in codes such as PLUTO, preserves the direction of radiative fluxes and captures shadowing, but it artificially merges crossing beams of light, which can overestimate midplane temperatures by more than 40 percent in single-group calculations. A newer half-moment scheme reduces that error to a few percent, while discrete ordinates and Monte Carlo methods offer the greatest accuracy at the highest cost. The review’s message is that no single method suits every problem, and the choice must weigh computational expense against the physics one needs to capture, whether self-shadowing, scattering, or frequency-dependent heating.

The consequences of getting radiation transport right extend deep into planet formation theory. Global simulations of the vertical shear instability show that the resulting turbulence generates stresses whose strength depends on the square of the radial temperature gradient and on the local cooling time, challenging the classical assumption that turbulent viscosity scales simply with gas pressure. VSI turbulence can spawn long-lived anticyclonic vortices, provided simulations span the full azimuthal extent of the disk, and these storm systems are expected to act as efficient traps for pebbles and dust, potentially accelerating planetesimal formation. In irradiated disks, radiation hydrodynamical studies reveal that the instability can be localized to the surface layers when dust depletion lengthens midplane cooling times, producing a quiescent midplane beneath a vigorously turbulent atmosphere, a stratification with direct consequences for where dust can settle and planets can grow.

The review closes with a forward-looking agenda. The authors recommend abandoning fixed-temperature disk models in favor of self-consistent thermal evolution, since every realistic disk structure is subject to thermal baroclinic instabilities that fixed-temperature setups either suppress artificially or misrepresent. Future work must couple evolving dust populations to the opacity and cooling calculations, treat the separate temperatures of gas, dust, and radiation, and develop well-balanced or low-Mach-number schemes that resolve subsonic fluctuations in a supersonically rotating medium. The payoff is substantial: radiation hydrodynamic simulations calibrated against linear theory will underpin the interpretation of molecular line kinematics and scattered-light images from current and next-generation telescopes, transforming stunning pictures of planet-forming disks into quantitative tests of how worlds are born.

Subject of Research: Numerical radiation hydrodynamics methods and stability tests for simulating protoplanetary disks around young stars

Article Title: Numerical radiation hydrodynamics for circumstellar disks

Article References: Klahr, H., Baehr, H., Melon Fuksman, J. D., & Pfeil, T. (2026). Numerical radiation hydrodynamics for circumstellar disks. Living Reviews in Computational Astrophysics, 12(1), Article 3. https://doi.org/10.1007/s41115-026-00026-6

Image Credits: AI Generated

DOI: 10.1007/s41115-026-00026-6

Keywords: protoplanetary disks, radiation hydrodynamics, planet formation, vertical shear instability, convective overstability, baroclinic instability, radiative transfer, turbulence, numerical simulations, stellar irradiation, computational astrophysics, disk stability