Gravitational instability in protoplanetary disk with cooling: 2D global analysis

Gravitational instabilities (GI) play a significant role in both the evolution of protoplanetary disks (PPDs) and the planet formation process. (e.g., Cameron, 1978; Boss, 1997; Durisen et al., 2007; Boley & Durisen, 2010; Paardekooper, 2012; Xu & Kunz, 2021). They likely manifest during the early stages of protostellar disk evolution when the disk mass occupies more than 10% protostellar mass (Adams & Lin, 1993; Armitage, 2020). Disk mass estimates suggest that $50\%$ of Class 0 and $25\%$ of Class I disks exhibit GI instability (Kratter & Lodato, 2016).

The observational evidence suggests that GI may play a substantial role in PPDs. Recent advancements in observational facilities and instruments, such as ALMA, VLT/SPHERE, VLT/CRIRES, and GPI, have enabled spatially resolved observations of PPDs. Spiral structures have been detected in the sub-millimeter thermal emission reflected by the mm-sized dusts around the disk midplane (e.g., Pérez et al., 2016; Huang et al., 2018; Kurtovic et al., 2018). These observations are biased tracers of gas density profiles due to finite aerodynamic coupling between gas and dust. Moreover, optical and near-infrared (NIR) high contrast scattered-light images, arising from starlight scattered by $\mu$m-sized dusts suspended in the disk surface, have also revealed spiral structures (e.g., Garufi et al., 2013; Wagner et al., 2015; Avenhaus et al., 2017; Benisty et al., 2017; Uyama et al., 2018). The origin of these spirals remains a debate, with GI offering a plausible explanation (Pérez et al., 2016; Meru et al., 2017; Tomida et al., 2017; Hall et al., 2018; Huang et al., 2018; Speedie et al., 2024). Additionally, the presence of substellar companions at large radii (with semimajor axes $a>10$ AU) implies the importance of GI (e.g. Kratter et al., 2010; Forgan & Rice, 2013; Forgan et al., 2015).

The study of GI has a long history dated back to James Jeans. Regarding disk GI, in recent years, Goodman & Narayan (1988) explored the fundamental modes in self-gravitating incompressible tori, demonstrating the existence of I mode (with corotation outside the tori) and J mode (with corotation at the density maximum). Adams et al. (1989) showed that one-armed density wave can be unstable when the disk mass is sufficiently high ($M_{d}/M_{*}=1$). Noh et al. (1991) found nonaxisymmetric instability of high $m$’s when the disk is not so massive ($M_{d}=0.05\sim 0.5M_{\odot}$). Hadley et al. (2014) studied the global nonaxisymmetric instabilities in a thick disk to produce a large parameter space of star-disk system. Lodato & Rice (2005) revealed that sufficiently massive disks can exhibit short-lived arms. The series of studies (Laughlin & Rozyczka, 1996; Dong et al., 2015; Hall et al., 2019; Chen et al., 2021) demonstrated that GI-induced spirals can have more than two arms. In addition to classical GI, drag-mediated two-fluid GI (Coradini et al., 1981; Longarini et al., 2023) and secular GI (Ward, 2000; Youdin, 2005; Shariff & Cuzzi, 2011; Michikoshi et al., 2012; Takeuchi & Ida, 2012; Tominaga et al., 2020; Pierens, 2021) have been proposed for less massive disks.

The dynamical consequences of GI are strongly influenced by disk thermodynamics. When a disk becomes gravitationally unstable, GI-induced spirals spread across the disk, forming shocks that heat the disk to increase thermal pressure against GI. However, radiative cooling can help remove heat from the disk, making it more unstable. The final state is determined by the balance between heating and cooling. When the cooling timescale is shorter than the dynamical timescale, fragmentation may occur, leading to the formation of gravitationally bound objects (e.g., Boss, 2001; Gammie, 2001; Johnson & Gammie, 2003; Rice et al., 2005; Baehr et al., 2017). When the cooling timescale is larger than the dynamical timescale, the disk may settle into a quasi-steady, marginally unstable state (e.g., Goldreich & Lynden-Bell, 1965; Lin & Pringle, 1987; Gammie, 2001; Lodato & Rice, 2004; Boley et al., 2006; Michael et al., 2012; Steiman-Cameron et al., 2013), potentially sustaining long-lived substructures in the disk.

Motivated by the relationship of disk thermodynamics and GI, we perform a linear stability analysis to study the effects of cooling time for disk GI. In Section 2 we derive the linear perturbation equations and give the numerical methods. In Section 3 we show some representative modes in different GI regimes. In Section 4 we focus on the relationship of dominant modes and disk cooling timescales. In Section 5 we discuss some implications and caveats, especially the disk $\alpha$. In Section 6 and we summarize our results.

In this section we test with the mass-averaged Toomre numbers $\overline{Q}=0.53,1.15,1.65$ across the entire range of $\beta$ and show some representative results of GI modes. We focus on $m=2$ modes because the two-armed spirals are commonly observed in both millimeter continuum (Andrews et al., 2018; Huang et al., 2018; Kurtovic et al., 2018) and near-infrared scattered-light observations (Garufi et al., 2013; Benisty et al., 2015, 2017; Canovas et al., 2018; Uyama et al., 2018). To characterize the modes, we compare the morphology of eigenfunctions, the oscillation rate (or pattern speed) $\omega_{p}$, the growth rate $\gamma_{m}$, and the maximum location of $|\delta\Sigma|$ ($r_{|\delta\Sigma|_{\rm{max}}}$) under different $\overline{Q}$ and $\beta$. Additionally, we provide the locations of the inner and outer Lindblad resonances within which the spiral density wave propagates.

We have known that the growth rate and the mode transition depend not only on $\overline{Q}$ but also on $\beta$, and therefore, we scan the parameter space ($\beta,\overline{Q}$) for $\beta$ from $10^{-4}$ to $10^{4}$ and $\overline{Q}$ from $0.5$ to $1.7$ with fixed disk-star mass ratio $M_{d}/M_{*}=0.4$. We focus on the $m=2,3,4$ modes since $m=1$ mode may be stable and higher modes ($m\geqslant 5$) may not be GI but shear instability (recall that GI prefers long wave).

We choose $r_{|\delta\Sigma|_{\rm{max}}}$ to distinguish the transitional regime ($r_{|\delta\Sigma|_{\rm{max}}}>0.5$) from the local and global regimes ($r_{|\delta\Sigma|_{\rm{max}}}\lesssim 0.5$). Fig.9 shows the contour of $r_{|\delta\Sigma|_{\rm{max}}}$ as a function of $\beta$ and $\overline{Q}$. The bold characters “$L,T,G$” denote respectively the local, transitional and global regimes. For $\overline{Q}\lesssim 0.9$, the dominant modes in the disk are local regardless of disk thermodynamics $\beta$. As $\overline{Q}$ increases, at $\beta\sim 1$ GI enters the transitional regime, and moreover, slower cooling ($\beta>1$) requires lower $\overline{Q}$ to enter the transitional regime than faster cooling ($\beta<1$). When $\overline{Q}$ continues to increase, GI with both large and small $\beta$ will eventually leave the transitional regime to enter the global regime. However, GI with $\beta\sim 1$ tends to retain in the transitional regime, and higher mode owns a wider transitional regime. Similar to growth rate, this behavior at $\beta\sim 1$ may also result from the resonance of orbital dynamics and thermodynamics.

Gravitational instability significantly influences the evolution of PPDs and the planet formation. Firstly, ALMA has detected various substructures, e.g., spirals, rings, and gaps (Andrews et al., 2018), and GI is a potential mechanism for these large-scale substructures (Lin & Shu, 1964; Binney & Tremaine, 2008; Bae et al., 2023; Meru et al., 2017; Speedie et al., 2024). Secondly, although core accretion is considered to be the classical model to understand the giant planet formation even at large distance $\gtrsim$ 10 AU (e.g., Nielsen et al., 2019; Vigan et al., 2021; Currie et al., 2023), GI is an alternative for the high-mass and long-period planet formation, especially to understand the metallicity non-correlation of brown dwarfs, debris disks, and their host stars (e.g., Nayakshin, 2017). Thirdly, besides turbulent viscosity, GI can be responsible for the angular momentum transport (Rice & Nayakshin, 2017).

The disk turbulent viscosity is related to accretion rate and angular momentum transport. Usually it is described by the disk $\alpha$ (Shakura & Sunyaev, 1973)

where $\nu$ is the turbulent viscosity, $\alpha$ is a dimensionless coefficient, $H$ is the disk scale height and $c_{s}$ is the local sound speed. The disk $\alpha$ is essentially the ratio of turbulent velocity of largest eddies to sound speed. The typical accretion rate $\sim 10^{-8}$ to $10^{-7}M_{\odot}$ in protoplanetary disks (e.g., Armitage, 2020) corresponds to $\alpha$ $\sim 10^{-3}$ to $10^{-2}$. Turbulent viscosity arises from Reynolds stress which is given by

In addition to Reynolds stress, GI can also play a role in transporting angular momentum. Gravitational stress is given by (e.g., Gammie, 2001; Lodato & Rice, 2004; Boley et al., 2006; Michael et al., 2012; Steiman-Cameron et al., 2023)

where $S$ is disk surface area. By Reynolds and gravitational stresses we can calculate the effective $\alpha$

where the stress $T_{r\phi}$ can be either $T^{\rm Rey}_{r\phi}$ or $T^{\rm grav}_{r\phi}$. By the thermal balance between viscous heating and radiative cooling, we can derive $\alpha_{\rm grav}\approx 4/[9\gamma(\gamma-1)\beta]$ (Armitage, 2020) such that $\alpha_{\rm grav}\approx 0.1$ at $\beta\approx 1$ and $\gamma\approx 2$. This argument is for a nonlinear self-sustained disk but our analysis is about linear stability (i.e., perturbation grows but not at the equilibrium balance) so that the disk $\alpha$ in our linear study is roughly proportional to disk mass (see Eq. (26)).

Fig.10 shows the distribution of $\alpha_{\rm grav}$ and $\alpha_{\rm Rey}$ at different $\overline{Q}$ and $\beta$ in the three regimes. The deep valleys are caused by the sign reversals of torque. In our 2D study, $\alpha_{\rm grav}$ and $\alpha_{\rm Rey}$ are comparable in most regions of disk, but it is reported that $\alpha_{\rm Rey}$ is less significant than $\alpha_{\rm grav}$ in 3D studies (e.g., Michael et al., 2012; Steiman-Cameron et al., 2013; Bae et al., 2016; Béthune & Latter, 2022). In the local and global regimes $\alpha_{\rm grav}$ at $r_{|\delta\Sigma|_{\rm max}}$ reaches $\sim 10^{-2}$, but in the transitional regime it reaches $\sim 1$. Although density perturbation becomes highest at $r_{|\delta\Sigma|_{\rm max}}$ ($|\delta\Sigma|/\Sigma\gtrsim 50\%$ in the local regime, 20-40% in the transitional regime and 15% in the global regime), because of the other quantities in Eq.(26) (e.g., $\delta\Psi$, $r$, etc.), $\alpha_{\rm grav}$ at $r_{|\delta\Sigma|_{\rm max}}$ in the transitional regime is much stronger than in the other two regimes.

Our linear studies in both local and global regimes show strong similarities with 3D simulations of disks governed by GI (e.g., Boss, 2017). The rapid growth of instability in the local regime with low $\overline{Q}$ tends to create fragmentations and gravitationally bound objects in the inner disk, such as high-mass planets and brown dwarfs, regardless of cooling timescales. However, since GI saturates in the subsequent nonlinear phase and cooling is inefficient in the inner disk (e.g., $<1$ AU), fragmentation can be suppressed. For disks in the global regime, the lower growth rate leads to long-lived non-axisymmetric spirals. Notably, even for high $\overline{Q}$ ($\gtrsim 1.5$) and high $\beta$ ($\gtrsim 10$), the effective $\alpha$ can still reach $\sim 10^{-3}$, consistent with results from marginally gravitationally unstable disk studies (Boss, 2012), which suggest that GI could be the driving mechanism for generic accretion disk evolution.

The transitional regime is particularly intriguing, as robust instability develops in the outer disk, where $\beta\lesssim 1$ is typically satisfied (Lin & Youdin, 2015; Pfeil & Klahr, 2019). With the efficient cooling, substructures, fragmentations and substellar companions can be created by GI. Since the transitional regime dominates at $\beta\approx 1$ (see Fig.9), we expect that spirals with pitch angles ranging from a few degrees to over 10 degrees can be observed in the outer disk. Moreover, the large effective $\alpha_{\rm{grav}}$ in the transitional regime implies that GI is responsible for the angular momentum transport in the outer disk rather than viscosity.

In this paper we explore the characteristics of GI in a protoplanetary disk at various $\overline{Q}$ and $\beta$, and our key results are summarized as follows. Firstly, faster/slower cooling timescale (smaller/larger $\beta$) leads to faster/slower growth rate, which is consistent with numerical results. Moreover, growth rate changes sharply at $\beta\approx 1$. Secondly, different $\overline{Q}$ corresponds to different modes: low $\overline{Q}$ to local modes, intermediate $\overline{Q}\approx 1$ to transitional modes, and high $\overline{Q}$ to global modes. Density perturbation maximum $r_{|{\delta\Sigma}|_{\rm max}}$ of the local and global modes is located in the inner disk and close to equilibrium density maximum $r_{\Sigma_{\rm max}}$, but $r_{|{\delta\Sigma}|_{\rm max}}$ of the transitional modes in the outer disk and far away from $r_{\Sigma_{\rm max}}$. Moreover, the transitional modes lead to much stronger $\alpha_{\rm grav}$ than the other two modes, suggesting much more efficient transport of angular momentum. Thirdly, the parameter scan reveals that the transitional modes dominate at $\beta\approx 1$. In summary, we infer that GI works in the outer disk at $\overline{Q}\approx 1$ and $\beta\approx 1$ to trigger the transitional modes for efficient transport of angular momentum, disk substructures (e.g., spirals and fragmentations), and planet/brown dwarf formation.

We assume a 2D thin disk. Although our calculations involve a relatively large aspect ratio, it should be noted that an early-age disk can be fairly thick even with an envelope, which will induce the additional mechanisms (Kratter & Lodato, 2016). Moreover, our linear analysis cannot capture the nonlinear effects such as shocks, wave amplitude saturation, etc. In addition, some other factors such as accretion and magnetic fields can be taken into account.

We thank Hongping Deng and Cong Yu for the help on equilibrium, Qiang Hou for his initial assistance, Haiyang Wang for self-gravity, Yuru Xu for boundary conditions, and Cathie Clarke for disk $\alpha$.