Revisiting gravitational instability in protostellar discs with improved radiative cooling models

Abstract: Young discs are expected to be significantly more massive than those observed at $&gt;1$ Myr and it is at this earliest stage that planet formation likely begins. Such massive discs may be susceptible to the gravitational instability (GI), therefore we need to determine the disc and stellar properties for which the GI is active to understand its role in early disc evolution and planet formation. Prior work has been limited by model assumptions and inaccuracies due to the complex nature of the thermodynamics of protostellar discs so we now revisit this question using an improved method to approximate radiative cooling within hydrodynamics simulations. We have explored a wide parameter space, representative of young protostellar discs of 0.1 to 1 M$<em>{\odot}$ and include irradiation from the host star. The parameters for which discs form spirals and fragment were found to differ to those obtained from earlier simulations. The outer regions of discs with radii of 50 au may be susceptible to fragmentation, meaning that GI-driven planet formation is not restricted to only the most extended discs. The additional thermal support due to stellar irradiation increases the disc mass that remains stable against GI: discs may reach up to $\gtrsim 0.4$ M$</em><em>$ without fragmenting, providing a considerable quantity of material for building planets. Large scale spiral arms only developed for $M_</em>\lesssim$ 0.3 M$_{\odot}$, except in the most compact discs. Furthermore, the long-lived spiral structures that form tend to be flocculent and compact, indicating that large-scale spiral arms should not be considered a typical outcome of GI.

• The paper shows that improved optical-depth estimates, especially the modified Lombardi method, substantially change GI outcomes by avoiding the excessive cooling produced by the Stamatellos approximation.
• Compact 50-au discs can fragment at 20–40 au when the disc-to-star mass ratio exceeds about 0.4, while irradiation stabilises many 100–200-au discs with masses of at least 0.4 times the stellar mass.
• GI-active discs usually produce faint, high-order flocculent spirals rather than prominent grand-design patterns, so simulations should use accurate cooling models and run for at least 10 outer rotation periods.

Young protostellar discs are expected to be substantially more massive than the discs observed at ages greater than 1 Myr, and it is precisely in this earliest phase that planet formation is thought to begin. Whether such massive discs undergo the gravitational instability (GI) — and if so, whether they fragment or settle into a self-regulated spiral state — depends sensitively on their thermodynamics. The paper "Revisiting gravitational instability in protostellar discs with improved radiative cooling models" (2608.13058) re-examines this question using a more accurate approximate radiative cooling treatment than has been used in most prior parameter studies, and finds that several established conclusions about where and when GI operates require revision.

Motivation and approach

The classical diagnostic for disc self-gravity is Toomre's $Q = c_s \kappa / (\pi G \Sigma)$, with instability setting in near $Q=1$ and observable effects appearing for $Q \lesssim 1.5$. A quasi-steady, self-regulated state can be maintained when cooling is slow relative to the orbital period ($$\beta_{\rm cool} = \Omega_{\rm K} t_{\rm cool} \gtrsim 10$$–20), whereas $\beta_{\rm cool} \lesssim 8$ typically leads to fragmentation. Analytical work had long suggested that fragmentation is confined to large radii — beyond roughly 70–100 au (Rafikov 2005; Clarke 2009) — implying that GI-driven planet formation yields few planetary-mass objects. However, these expectations rest on simplified thermodynamic treatments: imposed constant $\beta_{\rm cool}$, barotropic equations of state, or the widely used polytropic pseudo-cloud approximation of Stamatellos et al. (2007), which assumes each gas parcel sits within a spherical polytropic cloud and consequently overestimates column density in disc geometries.

The authors build on the improved column-density estimators introduced by Young et al. (2024), which tailor the pseudo-cloud approach specifically to self-gravitating discs. Four variants are compared: the original Stamatellos method (column density from the local gravitational potential); the Lombardi method (from local pressure and hydrodynamic acceleration); a combined method averaging the two scale-height estimates in inverse quadrature; and the "modified Lombardi" method, which replaces the Stamatellos estimate with an analytical scale height for a self-gravitating disc, $$H_0/H_* = \sqrt{\pi/2}/\sqrt{1 + 1/(Q_{3D}\sqrt{\pi/2})}$$, where $H_*=c_s/\Omega_{\rm K}$. The cooling prescription is coupled to flux-limited diffusion (FLD) in the hybrid scheme of Forgan et al. (2009), so that heat transport between optically thick neighbouring fluid elements is captured alongside radiative exchange at $\tau \sim 1$ surfaces. Stellar irradiation is included via a background temperature set by the stellar luminosity, attenuated by $\exp(-\Sigma_i \bar{\kappa}_i)$, with luminosities taken from MIST tracks at 0.5 Myr. Simulations use the SPH code {\sc phantom} with $Q=1$0 particles, sink-particle stars, and initial conditions constructed to be consistent with a self-gravitating, irradiated disc structure — an important safeguard against artificially unstable setups.

Sensitivity to the cooling method

A central result is that disc evolution is acutely sensitive to the accuracy of the optical-depth estimate. In head-to-head comparisons, discs evolved with the Stamatellos method develop spirals or fragment under conditions where the Lombardi, combined, and modified Lombardi methods produce stable axisymmetric discs. The cause is traced directly to the overestimated optical depth of the Stamatellos method: its mid-plane optical depth remains above unity out to ~150 au versus within ~100 au for the other methods, producing cooler mid-planes, lower $Q=1$1, and steeper fitted cooling profiles ($Q=1$2 with $Q=1$3 for Stamatellos versus $Q=1$4–3 for the modified Lombardi). For marginally unstable systems, even the three improved methods differ in whether fragments form and when. This establishes plainly that parameter studies conducted with the older cooling approximations may misplace the fragmentation boundary, and that the modified Lombardi method should be preferred going forward.

The simulated mid-plane temperatures of ~20–30 K agree with observational estimates for Class 0/I discs, providing external validation of the thermodynamics.

Revised stability boundaries

The parameter study spans stellar masses of 0.1–1.0 $Q=1$5, disc-to-star mass ratios from 0.1 to beyond unity, and outer radii of 50, 100, and 200 au, with simulations now run to at least 10 outer rotation periods (ORPs) — longer than many earlier studies, which the authors show matters: one disc remained axisymmetric for nearly 9 ORPs before a spiral developed, suggesting Haworth et al. (2020) may have missed spirals by evolving for only 3 ORPs.

Three findings stand out:

• Fragmentation at small radii: 50 au discs fragment for $Q=1$6, producing clumps at 20–40 au — well inside the ~70 au threshold from analytical work. This contradicts the analytical expectation that compact discs cannot fragment unless accretion rates exceed those measured here by two orders of magnitude, and it widens the parameter space for GI-formed planets that could migrate inward rather than remaining on wide (>100 au) orbits.
• Enhanced stability of extended discs: stellar irradiation stabilises 100 and 200 au discs up to masses of $Q=1$7 without fragmentation, and stars below ~0.5 $Q=1$8 support discs as massive as themselves. A disc of 0.5 $Q=1$9 around a solar-mass star, which Haworth et al. (2020) found would fragment at 0.4 $Q \lesssim 1.5$0, remains stable here.
• Spiral morphology: large-scale, low-$Q \lesssim 1.5$1 grand-design spirals form only around low-luminosity stars ($Q \lesssim 1.5$2 $Q \lesssim 1.5$3) or in very high surface-density discs. More typically, GI-active discs ($Q \lesssim 1.5$4) display faint, flocculent, higher-$Q \lesssim 1.5$5 structures unlikely to be observationally detectable. This challenges the assumption that prominent spirals are the generic signature of GI and offers a possible explanation for the dearth of observed GI spirals.

Discs with sustained spirals exhibit the zig-zag pattern in $Q \lesssim 1.5$6 characteristic of thermal self-regulation, persisting over the full 10–20 ORPs followed (~7000 yr at 50 au). Accretion rates in spiral-bearing discs are elevated relative to axisymmetric cases, with quasi-steady self-gravitating discs expected to sustain $Q \lesssim 1.5$7 $Q \lesssim 1.5$8 yr$Q \lesssim 1.5$9.

Limitations and open questions

The authors are explicit about the boundaries of validity. The polytropic-cooling-plus-FLD hybrid cannot capture shadowing structures or complex geometries such as misaligned disc components, for which full ray-tracing radiative transfer would be required; nor does it treat local dust-density variations, since opacity is column-averaged (though comparable ray-tracing studies share this limitation). Temperatures are underestimated in inner disc regions, though judged not to affect the conclusions. Resolution tests at $$\beta_{\rm cool} = \Omega_{\rm K} t_{\rm cool} \gtrsim 10$$0 particles show broadly similar evolution but slightly earlier fragmentation, leaving some residual resolution dependence in marginal cases. The comparison with Haworth et al. (2020) shows that the choice of irradiation attenuation model shifts outcomes modestly, but the cooling-method choice dominates. Open questions include whether the fainter, higher-order spirals seen in irradiation-dominated discs require higher-resolution study, whether dust trapping enhances the observability of low-contrast spirals, and how infall-driven self-regulation interacts with the irradiated self-regulated states found here.

Conclusion

By replacing an overestimating optical-depth approximation with the modified Lombardi method across a broad parameter space, this work revises the map of gravitational instability in young protostellar discs: fragmentation extends inward to tens of au in compact massive discs, extended discs tolerate larger masses before destabilising, and conspicuous grand-design spirals are the exception rather than the rule. The practical recommendations — adopt the modified Lombardi cooling approximation and evolve simulations for at least 10 dynamical timescales — together with the prediction that GI activity is concentrated in discs younger than 1 Myr, position upcoming Class 0 observations to test directly how strongly GI shapes the earliest phases of planet formation (2608.13058).