Massive Self-Gravitating Accretion Discs

• Massive self-gravitating accretion discs are defined by significant disc-to-star mass ratios (~0.1–0.5) where self-gravity competes with shear, pressure, and cooling to trigger gravito-turbulence and fragmentation.
• The topic employs diagnostics like the Toomre parameter, β-cooling times, and effective α-viscosity prescriptions to evaluate stability, thermal balance, and angular momentum transport.
• These discs are pivotal in varied astrophysical settings—from protostellar systems to AGN—shaping disc morphology, transport mechanisms, and the pathways for star or companion formation.

Searching arXiv for the cited SMD papers to ground the article in the literature. Massive self-gravitating accretion discs are accretion discs in which the disc’s own gravity is dynamically important, so that self-gravity competes with shear, pressure, cooling, and, in some settings, irradiation, magnetic transport, or binary torques. A common benchmark for “massive” is a disc-to-central-object mass ratio of order $q \equiv M_d/M_* \approx 0.1$, although several regimes discussed in the literature extend to $q \sim 0.5$ or higher [(Lodato et al., 2011); (Rice, 2016)]. In such discs, gravitational instability is commonly diagnosed with the Toomre parameter $Q = c_s \kappa / (\pi G \Sigma)$, with $\kappa \simeq \Omega$ in Keplerian flows; when $Q$ approaches unity, spiral structure, gravito-turbulence, non-local torques, and, under sufficiently rapid cooling, fragmentation into bound objects may occur [(Rice et al., 2011); (Rice, 2016)]. Across protoplanetary discs, discs around massive young stars, active galactic nuclei, and circumbinary black-hole discs, the central problem is how self-gravity modifies transport, thermal balance, morphology, and collapse thresholds [(Rice et al., 2011); (Roedig et al., 2012); (Illenseer et al., 2015); (Forgan et al., 2016)].

1. Dynamical definition and stability criteria

The standard local diagnostic of gravitational instability in an accretion disc is the Toomre parameter,

$Q = \frac{c_s \kappa}{\pi G \Sigma},$

where $c_s$ is the sound speed, $\kappa$ the epicyclic frequency, $G$ the gravitational constant, and $\Sigma$ the surface density. For razor-thin discs, axisymmetric linear instability occurs if $q \sim 0.5$0; in Keplerian discs, $q \sim 0.5$1 (Rice et al., 2011). The same criterion is used broadly across the SMD literature, with non-axisymmetric structure typically appearing when $q \sim 0.5$2–$q \sim 0.5$3 or, more generally, when the disc approaches a marginally self-regulated state with $q \sim 0.5$4 [(Lodato et al., 2011); (Rice, 2016)].

The distinction between local linear instability and nonlinear saturated states is important. In two-dimensional local shearing-sheet calculations with irradiation, quasi-steady self-gravitating states saturate at $q \sim 0.5$5–$q \sim 0.5$6, whereas comparable three-dimensional discs typically saturate near $q \sim 0.5$7 because finite thickness dilutes self-gravity by $q \sim 0.5$8 for the most unstable modes with $q \sim 0.5$9 (Rice et al., 2011). This implies that quoted threshold values are model-dependent, especially with respect to dimensionality and vertical structure.

Cooling introduces a second control parameter. In the standard $Q = c_s \kappa / (\pi G \Sigma)$0-cooling notation,

$Q = c_s \kappa / (\pi G \Sigma)$1

and fragmentation occurs if cooling is sufficiently rapid, conventionally written as

$Q = c_s \kappa / (\pi G \Sigma)$2

The review literature emphasizes that the fragmentation boundary is more physically interpreted as a maximum sustainable self-gravitating stress, often near $Q = c_s \kappa / (\pi G \Sigma)$3, rather than as a single universal $Q = c_s \kappa / (\pi G \Sigma)$4 (Rice, 2016). In two-dimensional irradiated local simulations with $Q = c_s \kappa / (\pi G \Sigma)$5, however, explicit thresholds are reported: $Q = c_s \kappa / (\pi G \Sigma)$6 without irradiation, declining to $Q = c_s \kappa / (\pi G \Sigma)$7 at the strongest irradiation studied (Rice et al., 2011).

A further distinction concerns whether self-gravity remains local. Global three-dimensional radiative SPH calculations show that low-mass, geometrically thin discs can be represented by a local $Q = c_s \kappa / (\pi G \Sigma)$8-parametrization when $Q = c_s \kappa / (\pi G \Sigma)$9, generally the case for $\kappa \simeq \Omega$0, but this approximation breaks down as $\kappa \simeq \Omega$1 approaches unity and low-$\kappa \simeq \Omega$2 global spirals dominate (Forgan et al., 2010). This suggests that “massive” is not a single regime but a continuum from local gravito-turbulence to strongly non-local, wave-mediated transport.

2. Thermal balance, effective stress, and transport laws

In non-fragmenting SMDs, transport is often described by an effective viscous stress whose magnitude is set by local thermal balance. Without irradiation, local equilibrium yields

$\kappa \simeq \Omega$3

With irradiation, the self-gravitating component of the stress is reduced by the fraction of pressure support provided by the irradiation floor,

$\kappa \simeq \Omega$4

with $\kappa \simeq \Omega$5 used in the irradiated shearing-sheet study (Rice et al., 2011).

That same study measures the transport coefficient directly from Reynolds and gravitational stresses as

$\kappa \simeq \Omega$6

where

$\kappa \simeq \Omega$7

The main result is that, in non-fragmenting runs, $\kappa \simeq \Omega$8 remains set by local thermal equilibrium even for very long cooling times, up to $\kappa \simeq \Omega$9 or $Q$0 (Rice et al., 2011). This is one of the clearest statements of the gravito-turbulent closure in irradiated local models.

Other works frame the same physics in alternative transport prescriptions. One-dimensional self-similar models for geometrically thin, viscous self-gravitating discs reduce the full vertically integrated system to a nonlinear advection–diffusion equation for $Q$1,

$Q$2

with the surface density coupled through

$Q$3

Three viscosity prescriptions are considered there: a self-regulated “LP” form, a DSB $Q$4-viscosity $Q$5, and an RZ form $Q$6 (Illenseer et al., 2015). In these similarity solutions, the outer rotation-law exponent $Q$7 is the key control parameter; flatter rotation laws at large radii yield higher accretion rates, and fully self-gravitating discs evolve faster than nearly Keplerian discs (Illenseer et al., 2015).

The transport locality question remains central. Global radiative SPH simulations find that for $Q$8, $Q$9 is consistent with the cooling-based $Q = \frac{c_s \kappa}{\pi G \Sigma},$0, whereas for $Q = \frac{c_s \kappa}{\pi G \Sigma},$1–$Q = \frac{c_s \kappa}{\pi G \Sigma},$2, global $Q = \frac{c_s \kappa}{\pi G \Sigma},$3 spirals produce non-local energy flux and systematic departures from local thermodynamic equilibrium (Forgan et al., 2010). A plausible implication is that the $Q = \frac{c_s \kappa}{\pi G \Sigma},$4-closure is robust in the thin-disc, modest-$Q = \frac{c_s \kappa}{\pi G \Sigma},$5 limit, but should be interpreted as a time-averaged phenomenology in the most massive discs.

3. Fragmentation, stochasticity, and numerical convergence

Fragmentation in SMDs is usually defined by the sustained formation of overdense, long-lived clumps. In the irradiated local simulations, fragmenting runs are identified by clumps with densities $Q = \frac{c_s \kappa}{\pi G \Sigma},$6 the mean that survive for many cooling times, while the disc-averaged $Q = \frac{c_s \kappa}{\pi G \Sigma},$7 rises rather than settling (Rice et al., 2011). In protostellar disc reviews, fragmentation is located predominantly in the outer disc, typically beyond $Q = \frac{c_s \kappa}{\pi G \Sigma},$8–$Q = \frac{c_s \kappa}{\pi G \Sigma},$9 au under standard opacities, and rarely inside $c_s$0–$c_s$1 au even when opacity is reduced (Rice, 2016).

Irradiation weakens instability but does not generally remove fragmentation. As irradiation is increased toward the quasi-linear stability threshold, the critical cooling time decreases only by about a factor of two, from $c_s$2 to $c_s$3 for $c_s$4, and the fragmentation boundary in the $c_s$5 plane does not follow contours of constant $c_s$6 (Rice et al., 2011). The interpretation advanced there is that stronger irradiation raises $c_s$7 and $c_s$8, but if other transport channels are weak, mass builds up until self-gravity turns on and fragmentation ensues anyway (Rice et al., 2011).

A long-standing controversy concerns convergence of the fragmentation threshold. SPH resolution studies argue that the measured threshold depends on $c_s$9, with finer resolution needed to fragment at larger $\kappa$0. Interpreting the Meru–Bate results as resolution-driven trends, Lodato and Clarke infer tentative convergence at $\kappa$1 a few $\kappa$2–$\kappa$3 particles for $\kappa$4, with a converged threshold around $\kappa$5–$\kappa$6, and a preferred estimate $\kappa$7 (Lodato et al., 2011). Two numerical origins are proposed: artificial-viscosity heating, quantified by

$\kappa$8

and smoothing of density peaks over finite $\kappa$9, leading to a saturating relation

$G$0

with best-fit $G$1 and $G$2 (Lodato et al., 2011).

The review literature treats this issue cautiously. It notes reports of non-convergence and proposed fragmentation at $G$3, but also cites high-resolution calculations showing no fragmentation for $G$4 and analyses indicating that simulations with $G$5 typically do not fragment (Rice, 2016). The robust consensus is narrower than the numerical debate: transport relations and the perturbation–stress connection are comparatively stable, whereas exact thermodynamic thresholds are sensitive to resolution and implementation [(Lodato et al., 2011); (Rice, 2016)].

Stochastic fragmentation has also been examined explicitly. Two-dimensional SPH studies of the waiting-time distribution between strong shocks find an exponential law,

$G$6

with $G$7, most probable waiting times of $G$8–$G$9, and negligible probability of shock-free intervals longer than $\Sigma$0 (Young et al., 2015). Combining this with $\Sigma$1, the survival probability of a contracting clump scales as $\Sigma$2 with $\Sigma$3, leading to the conclusion that stochastic fragmentation cannot move the fragmentation radius inward by more than $\Sigma$4 (Young et al., 2015). This suggests that stochasticity does not qualitatively alter the standard view that direct gravitational collapse is largely confined to the outer disc.

4. Morphology, spectra, and the universality of gravito-turbulent structure

In non-fragmenting SMDs, the turbulent structure is not arbitrary. The irradiated local simulations show that the power spectrum of surface-density perturbations is uniquely set by $\Sigma$5, not by $\Sigma$6 or irradiation level. All spectra share the same shape and peak at

$\Sigma$7

implying that most of the stress arises from wavelengths much smaller than the local box size, and simulations with the same $\Sigma$8 but different $\Sigma$9 have nearly identical spectra (Rice et al., 2011). The perturbation amplitude scales as

$q \sim 0.5$00

generalizing earlier results that $q \sim 0.5$01 (Rice et al., 2011).

Global radiative SPH calculations further distinguish between local, high-$q \sim 0.5$02 structure and global low-$q \sim 0.5$03 spirals. Thin low-$q \sim 0.5$04 discs distribute power over higher azimuthal mode number, with low variability in temperature and $q \sim 0.5$05, while massive thick discs are dominated by $q \sim 0.5$06 spirals, transient bursts, and outward wave-mediated energy transport (Forgan et al., 2010). In this sense, morphology is itself a diagnostic of whether the disc is in the local gravito-turbulent regime or the non-local global regime.

In massive protostellar systems, disc asymmetry feeds back dynamically on the central star. Radiation-hydrodynamic simulations including the indirect potential from stellar wobbling find that this backreaction makes discs smaller and rounder, delays and reduces fragmentation, changes angular momentum redistribution, and suppresses gaseous clump ejection (Meyer et al., 2024). Without wobbling, fragmentation begins between $q \sim 0.5$07–$q \sim 0.5$08 au by $q \sim 0.5$09 kyr and becomes violent by $q \sim 0.5$10 kyr, with many clumps out to $q \sim 0.5$11–$q \sim 0.5$12 au and later to $q \sim 0.5$13 au. With wobbling, fragmentation is delayed until $q \sim 0.5$14 kyr and remains milder, with only one migrating clump at $q \sim 0.5$15 kyr and no clump ejections reported (Meyer et al., 2024). This suggests that non-axisymmetric backreaction can act as an additional self-regulation channel in massive protostellar SMDs.

The same work links morphology to observability through synthetic millimetre imaging. Post-processing with RADMC-3D and CASA shows that wobbling models yield images in better agreement with ALMA observations of AFGL 4176 mm1, G17.64+0.16, and G353.273 than fixed-star models (Meyer et al., 2024). The significance is not that self-gravity simply produces spirals, but that the detailed nonlinear morphology of SMDs depends on feedbacks internal to the star–disc system.

5. Astrophysical settings

SMDs occur in several astrophysical environments, with the same basic instability criteria but different thermal and dynamical consequences.

In protoplanetary discs, self-gravity can dominate angular momentum transport at large radii. Under self-luminous conditions, fragmentation is expected beyond $q \sim 0.5$16–$q \sim 0.5$17 au; in more general protostellar-disc discussions, a practical fragmentation zone of $q \sim 0.5$18–$q \sim 0.5$19 au is emphasized under realistic opacities [(Rice et al., 2011); (Rice, 2016)]. Early, massive protostellar discs with $q \sim 0.5$20 around a solar-mass star and $q \sim 0.5$21 are globally stable but still self-gravity-modified. Radiation-hydrodynamic SPH calculations of $q \sim 0.5$22 seeds inserted into such discs show an initial rapid inward migration phase with $q \sim 0.5$23–$q \sim 0.5$24 yr, followed by gap opening and either slower inward Type II–like migration or outward migration if the gap edges become gravitationally unstable (Stamatellos et al., 2018). Without radiative feedback from the protoplanet, gap edges reach $q \sim 0.5$25–2 and drive outward migration; with feedback, edges remain at $q \sim 0.5$26 and migration remains inward (Stamatellos et al., 2018). This is not fragmentation of the disc itself, but it shows how proximity to self-gravity changes the migration and mass-growth pathways of embedded objects.

The solids budget is also altered in self-gravitating protostellar discs. Quasi-steady $q \sim 0.5$27 models with radiative cooling and no irradiation indicate that mm-sized pebbles are fragmentation-limited to $q \sim 0.5$28, that midplane pebble-to-gas density ratios are typically $q \sim 0.5$29, and that the streaming instability is therefore generally suppressed (Forgan, 2019). By contrast, GI fragments with initial masses $q \sim 0.5$30 can accrete pebbles efficiently while continuing to migrate, since they open gaps in the pebble component but generally fail to open gas gaps (Forgan, 2019). A plausible implication is that early self-gravitating phases redistribute solids in a way that imprints later core-accretion conditions.

Around massive young stars, SMDs are both a transport engine and an observational challenge. Semi-analytic $q \sim 0.5$31 models of candidate discs show that continuum-based masses can underpredict true masses by factors of $q \sim 0.5$32–$q \sim 0.5$33 because high optical depth suppresses mm emission (Forgan et al., 2016). In G11.92−0.61 MM1 and NGC 6334 I(N) SMA1b, self-gravitating models match observed continuum masses within a factor $q \sim 0.5$34 and predict outer-disc fragmentation into low-mass stellar companions, with fragment masses typically $q \sim 0.5$35 (Forgan et al., 2016). By contrast, AFGL 4176 mm1 and IRAS 16547−4247 are interpreted there as gravitationally stable because irradiation keeps $q \sim 0.5$36 above the instability regime (Forgan et al., 2016). The same theme appears in radiation-hydrodynamic collapse calculations, where self-gravity-generated torques in massive circumstellar discs sustain mean stellar accretion rates of $q \sim 0.5$37, comparable to explicit $q \sim 0.5$38-disc models, even though the three-dimensional accretion is episodic (Kuiper et al., 2011).

In active galactic nuclei and related massive black-hole environments, SMDs can be both transport layers and star-forming reservoirs. One-dimensional similarity solutions applied to AGN conclude that flatter outer rotation laws yield higher accretion rates and that fully self-gravitating discs evolve faster than nearly Keplerian discs, a point linked there to supermassive black-hole growth and quasar evolution (Illenseer et al., 2015). The review literature likewise notes that AGN discs may self-regulate or fragment to form stars depending on cooling conditions (Rice, 2016).

Circumbinary massive black-hole discs provide a distinct dynamical realization. Three-dimensional SPH simulations of binaries in self-gravitating circumbinary discs with $q \sim 0.5$39, $q \sim 0.5$40, and $q \sim 0.5$41–2 show a leaky cavity, edge overdensities, and spiral arms at $q \sim 0.5$42–$q \sim 0.5$43, with the net gravitational torque changing sign across corotation at $q \sim 0.5$44 (Roedig et al., 2012). The dominant net torque is described there as kinematic and non-resonant, arising from cavity streams on eccentric orbits rather than from a simple resonant tidal picture (Roedig et al., 2012). In retrograde circumbinary discs, where outer Lindblad resonances are absent, the cavity is smaller, $q \sim 0.5$45, yet circular binaries shrink at essentially the same rate as in prograde discs; for eccentric retrograde binaries, both $q \sim 0.5$46 decay and $q \sim 0.5$47 growth are approximately exponential while the binary remains coplanar (Roedig et al., 2013). These results situate SMDs within binary hardening theory and show that disc self-gravity can influence not only transport but also the orbital architecture of compact-object binaries.

6. Caveats, controversies, and theoretical limits

Several caveats recur across the SMD literature. Dimensionality matters: two-dimensional razor-thin calculations can overestimate the effective strength of self-gravity, and finite thickness lowers the saturated $q \sim 0.5$48 in three dimensions (Rice et al., 2011). Local shearing-sheet models omit global modes, non-local torques, and infall, whereas global models with realistic radiation often have limited resolution in the inner disc or in collapsing clumps [(Rice et al., 2011); (Kuiper et al., 2011); (Meyer et al., 2024)].

Cooling prescriptions are another major uncertainty. Idealized $q \sim 0.5$49-cooling is analytically and numerically convenient, but realistic discs experience opacity transitions, irradiation, and radiative diffusion. This matters acutely for fragmentation thresholds. Resolution studies explicitly argue that thermodynamic thresholds are less robust than transport diagnostics (Lodato et al., 2011), and the review literature continues to regard convergence, multiple fragmentation modes, and stochastic fragmentation as open issues even while judging the overall physical picture robust (Rice, 2016).

The validity of local transport closures has a similarly conditional status. For $q \sim 0.5$50, local $q \sim 0.5$51-models can be a good approximation; by $q \sim 0.5$52, transient low-$q \sim 0.5$53 spirals become important; and by $q \sim 0.5$54, only time-averaged local descriptions remain plausible because global wave transport is significant (Rice, 2016). Global radiative simulations sharpen this by associating the transition with $q \sim 0.5$55 and $q \sim 0.5$56 (Forgan et al., 2010).

Observational interpretation is also nontrivial. Massive discs may appear light in the continuum because they are optically thick (Forgan et al., 2016); spirals may indicate self-gravity only if the outer disc midplane is not overly stabilized by irradiation (Haworth et al., 2020); and compact structure in massive protostellar discs may depend sensitively on whether stellar wobbling is included (Meyer et al., 2024). The literature therefore repeatedly warns against direct inference of disc mass or GI state from continuum data alone.

Taken together, these caveats do not negate the core framework. Rather, they delimit it. Self-gravitating discs are consistently described as systems that hover near $q \sim 0.5$57, transport angular momentum through gravitationally driven stresses or global torques, fragment when cooling overwhelms self-regulation, and shift from local to non-local behaviour as mass ratio and thickness increase [(Rice et al., 2011); (Rice, 2016); (Forgan et al., 2010)]. What remains under active refinement is the precise placement of boundaries between these regimes, and the degree to which irradiation, resolution, inflow, magnetic transport, or stellar backreaction modify them.