Ring Formation in Protoplanetary Disks Driven by an Eccentric Instability

We simulate a 2D self-gravitating gaseous disk with the hydrodynamics code LA-COMPASS (Li et al. 2005, 2009a). LA-COMPASS evolves the fluid equations with a Godunov, shock-capturing scheme. In addition to evolving the continuity equation and momentum equations for the gas, LA-COMPASS also evolves the gas energy equation with an additional cooling term,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {E}_{g}}{\partial t}+{\boldsymbol{\nabla }}\cdot \left(({E}_{g}+P){\boldsymbol{v}}\right)=-{\boldsymbol{v}}\cdot {\boldsymbol{\nabla }}{\rm{\Phi }}-\displaystyle \frac{{c}_{v}{\rm{\Sigma }}(T-{T}_{\mathrm{eq}})}{{t}_{c}},\end{eqnarray} \tag{ 1 }$$

where Σ, P, ${\boldsymbol{v}}=u\hat{{\boldsymbol{r}}}+v\hat{{\boldsymbol{\phi }}}$, and Φ are the gas surface density, pressure, 2D velocity in the r-ϕ midplane, and total gravitational potential, respectively. The specific heat capacity is c_v = k_b /(μ(γ − 1)), γ is the ratio of specific heats, k_b is the Boltzmann constant, and μ is the mean particle mass, which we fix to μ = 2.3m_p where m_p is the proton mass. The total gas energy is given by E_g = Σ∣ v ∣²/2 + P/(γ − 1). We assume that pressure is related to temperature, T, via the usual ideal gas equation of state, P = (k_b /μ)ΣT.

The gas cools to an equilibrium temperature T_eq on a timescale t_c that is a fixed multiple of the orbital timescale given by the parameter β = t_c Ω_K, where Ω_K is the Keplerian rotation rate (Gammie 2001). In the momentum equation, we include the gravitational potential of a central star with mass M_⋆, as well as the disk's self-gravity (Li et al. 2009b). We do not include the indirect potential associated with the acceleration of the star due to the pull of an off-center disk. We have checked that a simulation including the indirect potential finds nearly exactly the same results as those presented in Section 2.3.

The cooling model we adopt is meant to approximate the more accurate radiative cooling models thought to govern real protoplanetary disks. We parameterize the cooling timescale with the relatively simple β model that has been used extensively in studies of self-gravitating disks (e.g., Gammie 2001; Lodato & Rice 2004, 2005; Kratter & Lodato 2016). This model has been shown to produce similar results to more sophisticated cooling models and makes an analytic treatment of our results possible in Section 3. Additionally, we also specify an equilibrium temperature profile rather than self-consistently computing one by balancing cooling with stellar and internal heating. However, we do not expect this to impact our results significantly because more realistic equilibrium temperature profiles tend to be simple power laws (e.g., Chiang & Goldreich 1997).

While protoplanetary disks can be viscous in general, observations suggest many disks may have quite low viscosity (Flaherty et al. 2017, 2018, 2020). We focus on the inviscid limit by excluding an explicit disk viscosity term in the fluid equation, although there is implicit numerical viscosity in our simulations which is very small (Li et al. 2009b).

For the initial surface density, we choose a power-law profile with both an inner hole and an outer taper,

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Sigma }}(r)=2.03{{\rm{\Sigma }}}_{0}\mathop{\underbrace{\left(1-{e}^{-{\left(r/{R}_{0}\right)}^{6}}\right)}}\limits_{{\rm{inner}}\,{\rm{hole}}}\mathop{\overbrace{\left(\displaystyle \frac{{R}_{0}}{r}\right)}}\limits^{{\rm{power}}-{\rm{law}}}\\ \,\times \,\,\,\mathop{\underbrace{{e}^{-{\left(r/(2{R}_{0})\right)}^{2}}}}\limits_{\,{\rm{outer}}\,{\rm{taper}}}\,,\end{array}\end{eqnarray} \tag{ 2 }$$

We choose our reference radius to be R₀ and reference density to be Σ₀ = Σ(R₀).

The initial temperature is set via the gas sound speed, ${c}_{s}^{2}{(r)=\gamma ({k}_{b}/\mu )T={c}_{0}^{2}(r/{R}_{0})}^{-1/2}$. We set T_eq to this initial, power-law temperature. For all of our simulations, we set γ = 1.5 and ${c}_{0}=0.03\sqrt{{{GM}}_{\star }/{R}_{0}}$. Additionally, we set our unit of time to be the Keplerian frequency at R₀, i.e., ${{\rm{\Omega }}}_{{\rm{K}}}({R}_{0})=\sqrt{{{GM}}_{\star }/{R}_{0}^{3}}$.

The simulations are performed in a 2D cylindrical domain from R_in = 0.2 to R_out = 6.0 on a uniformly spaced (N_r, N_ϕ ) = (4096, 4096) grid. The fluid field at the inner and outer boundaries is fixed at its initial value. A damping method for r ∈ [0.2, 0.3] ∪ [5.6, 6.0] is applied to prevent wave reflection (de Val-Borro et al. 2006). We set the initial azimuthal velocity v to provide the exact radial force balance, while a nonzero initial radial velocity u is sampled uniformly from [−10⁻³, +10⁻³]rΩ_K in each fluid cell as an initial perturbation. This perturbation could be caused by various processes, such as random fluid turbulence (Balbus & Hawley 1998; Flaherty et al. 2020). We have checked that the disk evolution is independent of the type and amplitude of the initial perturbation as long as there is power in the m = 1 azimuthal component.

Our Fiducial disk adopts ${{\rm{\Sigma }}}_{0}=1.6\times {10}^{-3}{M}_{\star }{R}_{0}^{-2}$ and β = 10⁻⁶. Hence, the disk has a total mass of 0.0186M_⋆ and is essentially locally isothermal. ³ With radiative cooling, typical protoplanetary disks can have β ≪ 1 (Zhang & Zhu 2020; see also Section 4.2). The effect of Σ₀ and β are explored in Section 3.1.

Figure 1 shows our initial disk configuration. The surface density peaks just outside of R₀ coinciding with a minimum in the Toomre Q = c_s Ω/(π GΣ) (Toomre 1964). We focus in this paper on disks that are gravitationally stable, Q > 2, to ensure that any instability we find is not due to a small Q.

Zoom In Zoom Out Reset image size

2.3.1. Modal Time Evolution

We let the disk evolve for 10⁴ orbits at R₀ to provide a general picture of the long-term evolution. This corresponds to roughly one million years around a solar-mass star if R₀ = 21.5 au.

Figure 2 shows three surface density snapshots. The disk evolves through three different stages: single-arm spiral (time = 1500 orbits), multiple closed eccentric rings (3500 orbits), and multiple nearly circular rings (9000 orbits). These features are all radially confined between r ∼ 1–2 R₀.

Zoom In Zoom Out Reset image size

To get a sense of the relative importance of the different azimuthal symmetries in the disk, we show the time evolution of the global Fourier amplitudes,

$$\begin{eqnarray}&&{C}_{m}(t)\equiv \displaystyle \frac{1}{{M}_{{\rm{d}}}}{\displaystyle \int }_{{R}_{\mathrm{in}}}^{{R}_{\mathrm{out}}}\left|{\displaystyle \int }_{0}^{2\pi }{\rm{\Sigma }}{e}^{-{im}\phi }d\phi \right|{rdr},\end{eqnarray} \tag{ 3 }$$

in the upper panel of Figure 3. The coefficients are normalized to the m = 0 coefficient, which corresponds to the total disk mass, M_d .

Zoom In Zoom Out Reset image size

At early times, the disk surface density variations are dominated by the m = 1 component. From the C₁ curve and the snapshots, this component appears as an exponentially growing, one-armed, trailing spiral with a growth rate of γ₁ ≈ 10⁻³. Upon measuring the time evolution of the m = 1 complex phase (shown in the bottom panel of Figure 3), we find that the spiral precesses coherently and in a prograde direction with a pattern speed of Ω_p ≈ 0.018. Because Ω_p ≪ Ω and the precession rate is coherent, what we are seeing is the growth of a slow eccentric mode (Tremaine 2001; Lin 2015; Lee et al. 2019a).

During this time, higher m modes are also exponentially growing with growth rates that are γ_m ≈ m γ₁, albeit at an amplitude much less than C₁. This relation suggests that these modes are harmonics of the m = 1 mode and may be driven by nonlinear, mode–mode interactions (Laughlin et al. 1997, 1998; Lee 2016).

By ∼1700 orbits, the growth stalls, and C₁ saturates to a strength of ∼10%. This phase of evolution coincides with a morphological transition of the one-armed spiral into a set of closed elliptical rings and a slowdown in the m = 1 pattern speed.

After ∼5000 orbits, the eccentricity transitions yet again to being predominantly peaked in the inner regions of the disk (r < R₀). All modes start to decay from their peak amplitude on a timescale of ∼10⁴/m orbits, and the rings transition to being mostly axisymmetric. While fading away, the m = 1 component maintains a coherent pattern speed of Ω_p ≈ 10⁻².

2.3.2. Ring Evolution

Figure 4 shows the azimuthally averaged Σ (denoted by $\left\langle {\rm{\Sigma }}\right\rangle$) in the Fiducial simulation. During the growth phase (upper panel), $\left\langle {\rm{\Sigma }}\right\rangle$ develops three well-defined overdensities that are at radii r ≈ 1.1, 1.5, 1.8R₀ and consistently grow over time. Once the instability saturates (middle panel), there is significant nonlinear evolution of 〈Σ〉 that changes the ring locations. In the final decay phase (bottom panel), $\left\langle {\rm{\Sigma }}\right\rangle$ settles into a quasi-steady-state configuration consisting of three prominent rings at r = 1.1, 1.2, 1.5R₀, and additional, weaker features at r =0.85R₀ and 1.8R₀. Even though $\left\langle {\rm{\Sigma }}\right\rangle$ increases in the rings, the disk maintains Q ≳ 2 throughout the simulation. We find no signatures of other instabilities, e.g., gravitational, Rossby wave, Rayleigh, or others.

Zoom In Zoom Out Reset image size

With our Fiducial simulation, we have shown that a disk with inner and outer tapers traps an unstable one-arm spiral. This growing spiral arm has been previously reported by Lin (2015) for locally isothermal disks (γ = 1 and β → 0) and in cooled disks with β ≲ 0.1, although the latter result was not discussed in detail and used preliminary simulations. We have shown, in detail, that a nonzero cooling time unambiguously results in a similar instability. But more importantly, we find that the nonlinear outcome of this spiral instability is the generation of concentric axisymmetric rings.

In reality, disk viscosity can affect the growth of the spiral and the lifetime of the rings. In our Fiducial simulation, viscosity is negligibly small. If we were to include a nonzero viscosity ν, it could suppress the growth of the spiral if the damping rate is faster than the growth rate. Moreover, a large ν may diffuse away the initial Σ profile before any growth could occur. Preliminary short-term simulations of the growth phase with constant viscosity find that ν ≲ 10⁻⁵, which corresponds to α ≲ 0.01 at R₀, gives the same growth rate and spiral profile as our Fiducial results. We expect that in a long-term viscous simulation, the lifetime of the axisymmetric rings will be longer than the lifetime of the spiral waves as a tightly wound spiral damps at a rate that is faster than the rings by a factor of (kw)² > 1, where k is the radial wavenumber of the spiral and w is the radial width of a typical ring.

In the next sections, we show from first principles why the disk is unstable and why that generically results in concentric rings.

Slow eccentric modes in self-gravitating, adiabatic disks have been thoroughly studied in Lee et al. (2019a). They found that there are three kinds of modes that are different in morphology and controlled by the profile of

$$\begin{eqnarray}&&g(r)=\displaystyle \frac{\pi G{\rm{\Sigma }}r}{{c}_{\mathrm{iso}}^{2}},\end{eqnarray} \tag{ 5 }$$

which measures the relative strength of pressure and self-gravity. ⁴

We show a collection of the different mode types in dispersion relation maps (DRM) in Figure 6. The DRM plots show contours of fixed ${\mathfrak{R}}(\omega )$ in r–kr space, where k is the radial wavenumber of the mode. The contours are calculated from a WKB-like approximation to the eccentricity equation (see Appendix B and Lee et al. 2019a for more details). Following the naming convention of Lee et al. (2019a), the three mode types are

• 1.  
  Pressure-dominated elliptical modes (e-p modes) occur when g < 1. These are long-wavelength retrograde (${\mathfrak{R}}(\omega )\lt 0$) modes that are trapped by two inner Lindblad resonances (locations were k = 0).
• 2.  
  Self-gravity-dominated elliptical modes (e-pg modes) occur when g > 1 and switch between long and short radial wavelength at Q-barriers (locations where kr = g).
• 3.  
  Spiral modes (s-modes) occur when g ≫ 1 and consist of leading (k < 0) and trailing (k > 0) spirals trapped by Q-barriers at short wavelengths. Each of these contours are separate, i.e., there are no Lindblad resonances connecting them. The mode seen in the Fiducial simulation is a trailing s-mode.

Zoom In Zoom Out Reset image size

To determine whether a given mode is unstable we look at the equation governing the modes angular momentum deficit, Σr³Ω∣E∣², which can be obtained from multiplying the eccentricity equation by the complex conjugate of E, E^*, and integrating in r over the whole disk. The growth rate can then be determined from the imaginary component of this equation. We show in Appendix C that doing so results in two terms that contribute to the growth rate ${\gamma }_{1}={\mathfrak{I}}(\omega )$:

$$\begin{eqnarray}\begin{array}{rcl}{\gamma }_{1}\displaystyle \int 2{r}^{3}{\rm{\Sigma }}{\rm{\Omega }}| E{| }^{2}{dr} & = & -\displaystyle \int {r}^{3}P\left[{\beta (\gamma -1)\left|\displaystyle \frac{{dE}}{{dr}}\right|}^{2}\right.\\ & & \left.-\displaystyle \frac{d\mathrm{ln}{c}_{\mathrm{iso}}^{2}}{{dr}}{\mathfrak{I}}\left\{(1-i\beta )E\displaystyle \frac{{{dE}}^{* }}{{dr}}\right\}\right]{dr},\end{array}\end{eqnarray} \tag{ 6 }$$

to first order in β. For an elliptical mode which has $E\propto \cos ({kr})$, the integrand on the right-hand side is proportional to $-{k}^{2}{\cos }^{2}({kr})$, so that there is no growth (γ₁ ≤ 0); hence, elliptical modes are stable.

For an s-mode, E ∝ e^ikr , so that ∣dE/dr∣² = k²∣E∣² and EdE^*/dr = − ik∣E∣². From the DRM, we see that s-modes are radially confined, so that we can approximate the integral at the radius where $g={g}_{\max }$. Doing so, we find that γ₁ is approximate,

$$\begin{eqnarray}&&{\gamma }_{1}\approx -\displaystyle \frac{{g}_{\max }{h}^{2}{\rm{\Omega }}}{(2\gamma )}\left[(\gamma -1)\beta {g}_{\max }+\displaystyle \frac{d\mathrm{ln}{c}_{\mathrm{iso}}^{2}}{d\mathrm{ln}r}\mathrm{sign}(k)\right],\end{eqnarray} \tag{ 7 }$$

where h = c_s/(rΩ), and we approximate ${kr}\approx {g}_{\max }$. When β is small, γ₁ is approximately linear to ${g}_{\max }$.

One can see that only an s-mode can be unstable. In particular, only modes with kdc²/dr < 0 have the potential to be unstable, i.e., trailing (leading) spirals must be accompanied by decreasing (increasing) temperature profiles. For an instability to occur, β must satisfy

$$\begin{eqnarray}&&\beta \lesssim {\beta }_{{\rm{c}}}\equiv \displaystyle \frac{1}{1-\gamma }\displaystyle \frac{\mathrm{sign}(k)}{{g}_{\max }}\displaystyle \frac{d\mathrm{ln}{c}_{\mathrm{iso}}^{2}}{d\mathrm{ln}r}.\end{eqnarray} \tag{ 8 }$$

Because slower cooling acts to suppress the instability, the disk must cool fast enough for an initial spiral mode to grow.

Figure 7 compares the growth rates between our linear results (solid curve) and a set of hydrodynamical simulations (blue dots) for different values of ${g}_{\max }$ and β. Both the linear and the hydrodynamics results suggest that an unstable mode exists when ${g}_{\max }\,\gt \,9$ (top panel). The linear theory growth rates agree to 10% of hydrodynamics results for ${g}_{\max }\gt \,12$. Between ${g}_{\max }\,=\,9$ and 12, the results differ by 71% and 16%. This discrepancy is likely associated with the ${g}_{\max }$ being very close to the instability threshold. All of the larger ${g}_{\max }$ simulations remain gravitationally stable, Q ≳ 2, throughout their evolution.

Zoom In Zoom Out Reset image size

For a fast-enough cooling rate (bottom panel), the growth rate γ₁ is almost a constant. For β close to the critical value β_c, the growth rate drops sharply from nearly maximum to zero as β increases. The hydrodynamic results suggest that β_c is between 0.06 and 0.08. Equations (6) and (8) yield β_c = 0.08 and 0.07, respectively. The growth rates by Equation (6) agree within 13% of the hydrodynamics results for β ≤ 0.02.

Equation (7) (dashed line) provides a good approximation to γ₁ for ${g}_{\max }\gt 12$. The main error introduced when going from Equation (6) to (7), for ${g}_{\max }\gt 12$, comes from the approximation of the integral. Despite this, Equation (7), which does not require the E(r) profile, only overestimates γ₁ by ∼10%–20%. At a smaller ${g}_{\max }$, Equation (7) does not apply as the s-mode is either nonexistent or not prominent.

In this section, we have shown that an instability develops in a self-gravitating disk with cooling when (1) the ratio of self-gravity to pressure (i.e., g) is large enough to support a trapped spiral mode, (2) the spiral coincides with a nonzero temperature gradient, and (3) the cooling rate is faster than the critical rate given in Equation (8). Moreover, we have derived an analytical approximation for the resulting growth rate from the dispersion relation for the gas. In the limit of β → 0, our analytical growth rate agrees with the angular-momentum-based estimate of Lin (2015).

While we have focused on one particular choice of Σ and ${c}_{\mathrm{iso}}^{2}$ profiles, the spiral instability mechanism applies to a broad range of astrophysically interesting density and temperature profiles. This is because the appearance of the s-modes is controlled by the radial profile of g, not Σ or ${c}_{\mathrm{iso}}^{2}$ individually. As we discuss in detail in Section 4.1, this results in a wide range of density and temperature profiles that may be s-mode unstable.

From Figure 4, we know that when there is a growing one-armed spiral, the disk also produces axisymmetric rings. To understand this process, we examine how angular momentum is transferred from the spiral into the mean flow of the disk, and how this deposition process generically results in a nonconstant radial mass flux that drives the formation of multiple rings.

In general, angular momentum is stored in the mean flow of the disk as well as any waves present. Angular momentum exchange occurs between these two reservoirs during wave excitation and when there is nonzero dissipation (from e.g., viscosity or shocks) or baroclinicity (i.e., ∇ P × ∇Σ ≠ 0). These latter two processes are equivalent to the loss of conservation of the disk vortensity, $\xi =\hat{{\boldsymbol{z}}}\cdot ({\boldsymbol{\nabla }}\times {\boldsymbol{v}})/{\rm{\Sigma }}$.

In inviscid disks with cooling, baroclinic effects transfer angular momentum from the mean flow, whose mass and angular momentum evolve as

$$\begin{eqnarray}&&2\pi r\displaystyle \frac{\partial \left\langle {\rm{\Sigma }}\right\rangle }{\partial t}+\displaystyle \frac{\partial }{\partial r}(-\dot{M})=0,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&2\pi {r}^{2}\displaystyle \frac{\partial \left\langle {\rm{\Sigma }}\right\rangle \left\langle v\right\rangle }{\partial t}+\displaystyle \frac{\partial }{\partial r}\left(-\dot{M}r\left\langle v\right\rangle \right)=-{t}_{\mathrm{dep}}\end{eqnarray} \tag{ 10 }$$

into a wave, whose angular momentum follows ⁵

$$\begin{eqnarray}&&2\pi {r}^{2}\displaystyle \frac{\partial \left\langle {\rm{\Sigma }}^{\prime} v^{\prime} \right\rangle }{\partial t}+\displaystyle \frac{\partial }{\partial r}\left(2\pi {r}^{2}\left\langle {\rm{\Sigma }}{uv}^{\prime} \right\rangle +{F}_{G}\right)={t}_{\mathrm{dep}},\end{eqnarray} \tag{ 11 }$$

where $\dot{M}=-2\pi r\left\langle {\rm{\Sigma }}u\right\rangle$ is the radial mass flux and F_G is the flux of angular momentum due to self-gravity (Goldreich & Tremaine 1979). The exchange torque density between the wave and mean flow is (Dempsey et al. 2020)

$$\begin{eqnarray}&&\displaystyle \frac{{t}_{\mathrm{dep}}}{2\pi r}=-\left\langle u^{\prime} {\rm{\Sigma }}^{\prime} \right\rangle \displaystyle \frac{\partial (r\left\langle v\right\rangle )}{\partial r}+\left\langle {\rm{\Sigma }}\right\rangle \left\langle u^{\prime} \displaystyle \frac{\partial ({rv}^{\prime} )}{\partial r}\right\rangle -\left\langle \displaystyle \frac{{\rm{\Sigma }}^{\prime} }{{\rm{\Sigma }}}\displaystyle \frac{\partial P^{\prime} }{\partial \phi }\right\rangle .\end{eqnarray} \tag{ 12 }$$

From the equation governing the evolution of the disk vortensity, t_dep for an eccentric wave can be shown to be (Lubow 1990, see also Appendix D)

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{t}_{\mathrm{dep}}}{4\pi r} & = & -{\gamma }_{1}{r}^{3}{{\rm{\Sigma }}}^{2}\displaystyle \frac{\partial \left\langle \xi \right\rangle }{\partial r}| E{| }^{2}+r\displaystyle \frac{\partial {c}_{\mathrm{iso}}^{2}}{\partial r}{\rm{\Sigma }}{\mathfrak{I}}\left({E}^{* }r\displaystyle \frac{\partial E}{\partial r}\right)\\ & & +\beta (\gamma -1){r}^{2}{c}_{\mathrm{iso}}^{2}{\rm{\Sigma }}{\left|\displaystyle \frac{\partial E}{\partial r}\right|}^{2},\end{array}\end{eqnarray} \tag{ 13 }$$

in disks with β ≪ 1. The first term depends on the background vortensity gradient and is zero for stable modes, while the remaining terms are exactly the same as the terms on the right-hand side of Equation (6) that set the growth rate. We see now that they correspond to a nonzero persistent background torque proportional to the disk's temperature gradient (Lin & Papaloizou 2011; Lin 2015) and a stabilizing torque proportional to the cooling rate.

The loss of angular momentum from the mean flow induces a nonzero radial mass flux in the location of the wave. From Equations (9)–(10), this $\dot{M}$ is

$$\begin{eqnarray}&&\dot{M}=\displaystyle \frac{2}{r{{\rm{\Omega }}}_{K}}\left[2\pi {r}^{2}\left\langle {\rm{\Sigma }}\right\rangle \displaystyle \frac{\partial \left\langle v\right\rangle }{\partial t}+{t}_{\mathrm{dep}}\right].\end{eqnarray} \tag{ 14 }$$

The ${\partial }_{t}\left\langle v\right\rangle$ term is mostly determined by the disk eccentricity, which is ultimately sculpted by t_dep. We estimate ${\partial }_{t}\left\langle v\right\rangle$ by assuming the disk maintains an approximate radial force balance by modifying the pressure and self-gravity-corrected Keplerian rotation, Ω₀, by a small amount, Ω = Ω₀ + Ω₁. From the azimuthally averaged radial velocity equation, the Ω₁ required to maintain balance is (Lubow 1990)

$$\begin{eqnarray}\begin{array}{rcl}2r{{\rm{\Omega }}}_{0}{{\rm{\Omega }}}_{1} & = & \displaystyle \frac{1}{\left\langle {\rm{\Sigma }}\right\rangle }\displaystyle \frac{\partial \left\langle P\right\rangle }{\partial r}+\displaystyle \frac{\partial \left\langle {\rm{\Phi }}\right\rangle }{\partial r}-r{{\rm{\Omega }}}_{0}^{2}\\ & & -\left\langle \displaystyle \frac{v{{\prime} }^{2}}{r}\right\rangle +\left\langle \displaystyle \frac{v^{\prime} }{r}\displaystyle \frac{\partial u^{\prime} }{\partial \phi }\right\rangle \\ & & +\left\langle {\left(\displaystyle \frac{1}{{\rm{\Sigma }}}\right)}^{{\prime} }\displaystyle \frac{\partial P^{\prime} }{\partial r}\right\rangle +\left\langle u^{\prime} \displaystyle \frac{\partial u^{\prime} }{\partial r}\right\rangle .\end{array}\end{eqnarray} \tag{ 15 }$$

The first three terms on the right-hand side sum to zero as they constitute the dominant radial force balance. For slow eccentric modes, the pressure term is ∼h² smaller than the velocity terms, and so we can ignore it. Because we are interested in unstable modes, the remaining right-hand-side terms grow as ${e}^{2{\gamma }_{1}t}$, so that the ${\partial }_{t}\left\langle v\right\rangle$ term in Equation (14) is approximately

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial \left\langle v\right\rangle }{\partial t} & \approx & \displaystyle \frac{{\gamma }_{1}}{r{{\rm{\Omega }}}_{0}}\left[-\left\langle v{{\prime} }^{2}\right\rangle +\left\langle v^{\prime} \displaystyle \frac{\partial u^{\prime} }{\partial \phi }\right\rangle +\left\langle u^{\prime} r\displaystyle \frac{\partial u^{\prime} }{\partial r}\right\rangle \right]\\ & & \approx -\displaystyle \frac{1}{2}{\gamma }_{1}r{{\rm{\Omega }}}_{{\rm{K}}}| E{| }^{2}+{\gamma }_{1}{r}^{2}{{\rm{\Omega }}}_{{\rm{K}}}\displaystyle \frac{\partial }{\partial r}| E{| }^{2},\end{array}\end{eqnarray} \tag{ 16 }$$

where in the last line we replaced $u^{\prime}$ and $v^{\prime}$ with the complex eccentricity for an s-mode. The final induced $\dot{M}$ for rapid cooling is

$$\begin{eqnarray}&&\dot{M}\approx +2\pi {r}^{2}{\rm{\Sigma }}{\gamma }_{1}\left[-| E{| }^{2}+2r\displaystyle \frac{\partial }{\partial r}| E{| }^{2}\right]\end{eqnarray} \tag{ 17a }$$

$$\begin{eqnarray}&&-8\pi {r}^{3}{{\rm{\Sigma }}}^{2}\displaystyle \frac{{\gamma }_{1}}{{{\rm{\Omega }}}_{{\rm{K}}}}\displaystyle \frac{\partial \left\langle \xi \right\rangle }{\partial r}| E{| }^{2}\end{eqnarray} \tag{ 17b }$$

$$\begin{eqnarray}&&+\displaystyle \frac{8\pi r}{{{\rm{\Omega }}}_{{\rm{K}}}}\left[r{\rm{\Sigma }}\displaystyle \frac{{{dc}}_{\mathrm{iso}}^{2}}{{dr}}{\mathfrak{I}}\left({E}^{* }\displaystyle \frac{\partial }{\partial r}E\right)+\beta {c}_{\mathrm{iso}}^{2}(\gamma -1){\left|\displaystyle \frac{\partial E}{\partial r}\right|}^{2}\right].\end{eqnarray} \tag{ 17c }$$

The terms proportional to γ₁ were previously derived by Lubow (1990) for general m, while the last two terms are unique to a β-cooled disk with an eccentric spiral. To complete the calculation of $\dot{M}$, we can use the linear solution for E. One can then use Equation (9) to compute the mean density evolution.

We show in Figure 8 that $\dot{M}$ as given by Equations 17(a)–(c) agrees remarkably well with the $\dot{M}$ in the exponential growth phase of the Fiducial simulation. In particular, we recover the approximate locations of the rings, which coincide with convergent points in the axisymmetric flow, i.e., at rising inflection points of $\dot{M}$. The primary ring locations are set by the ${\partial }_{t}\left\langle v\right\rangle$ terms in Equation 17(a). The vortensity (Equation 17(b)) and baroclinic (Equation 17(c)) terms mostly shift the $\dot{M}$ profile and slightly modify the ring locations. The small disagreement between the linear and hydrodynamical results is primarily due to the m ≠ 1 perturbations and the terms of ${ \mathcal O }({\gamma }_{1}{h}^{2},{\gamma }_{1}\beta )$ that we have dropped.

Zoom In Zoom Out Reset image size

Overall, this excellent agreement demonstrates that the initial mass redistribution leading to the formation of multiple rings occurs because of the exchange of angular momentum from the background disk and the spiral wave.

While the initial phase is linear and fully describable by our linear theory, nonlinear evolution takes over at later times (see Figure 4). In our Fiducial simulation, $\left\langle {\rm{\Sigma }}\right\rangle$ eventually evolves to the point where the instability saturates at t ∼ 1700 orbits. By this time, the disk can no longer trap the s-mode, and instead, the fundamental mode becomes elliptical.

The ring formation driven by the s-mode also stops with the saturation of the instability. Follow-up nonlinear evolution can still modify the sharpness and locations of the rings until the elliptical mode also fades away. We suspect this is caused by the inner boundary condition in the Fiducial simulation, because an elliptical eccentric mode dominated by self-gravity, unlike an s-mode, tends to have a larger eccentricity near the inner edge of the disk (Lee et al. 2019a). The lifetime of such a mode will thus be short as any dissipation—whether physical or numerical—will operate on a faster timescale than the dissipation of the rings at larger radii.

In Section 3.1, we showed that disks are unstable when they can both cool faster than a critical rate and trap a slow spiral mode. Spiral modes appear when the distance (in phase space) between two Q-barriers is large enough to satisfy a quantum condition (Lee et al. 2019a; see also Figure 6). Because Q-barriers occur where kr = g(r), the g profile needs to be both peaked at some radius and wide enough to provide sufficient separation of the Q-barriers. This suggests that the larger-scale properties of the Σ and ${c}_{\mathrm{iso}}^{2}$ profiles are unimportant and that only the local properties in the vicinity of the g maximum are relevant.

Consider temperature and density profiles that produce a peak in g at the radius ${r}_{\max }$. Taylor-expanding g(r) around this peak produces a simple Gaussian g profile,

$$\begin{eqnarray}&&g(r)={g}_{\max }{e}^{-{\left(r-{r}_{\max }\right)}^{2}/{\sigma }_{g}^{2}},\end{eqnarray} \tag{ 18 }$$

where ${\sigma }_{g}=\sqrt{2{({d}^{2}\mathrm{ln}g/{{dr}}^{2})}^{-1}}$ is the Gaussian width. For a given ${r}_{\max }$, ${g}_{\max }$, and σ_g , the density and temperature profiles need only satisfy the following requirements at ${r}_{\max }$:

$$\begin{eqnarray}&&{g}_{\max }\,=\,{\left.\displaystyle \frac{\pi G{\rm{\Sigma }}r}{{c}_{\mathrm{iso}}^{2}}\right|}_{{r}_{\max }},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\displaystyle \frac{d}{d\mathrm{ln}r}\mathrm{ln}{\left(\displaystyle \frac{r{\rm{\Sigma }}}{{c}_{\mathrm{iso}}^{2}}\right)}_{{r}_{\max }}=0,\qquad {\sigma }_{g}^{-2}=\displaystyle \frac{1}{2}\displaystyle \frac{{d}^{2}}{{{dr}}^{2}}\mathrm{ln}{\left(\displaystyle \frac{{c}_{\mathrm{iso}}^{2}}{{\rm{\Sigma }}}\right)}_{{r}_{\max }}.\end{eqnarray} \tag{ 20 }$$

A wide range of Σ and ${c}_{\mathrm{iso}}^{2}$ profiles can be constructed to fulfill these requirements.

Using this g(r) profile, we solve for the linear γ₁ in the $({g}_{\max },{\sigma }_{g},\beta )$ parameter space. Our results are shown in Figure 9 for ${c}_{\mathrm{iso}}^{2}\propto {r}^{-1/2}$ and Σ given by Equations (5) and (18). We find that there is a region of σ_g –${g}_{\max }$ parameter space that is unstable to one-armed spirals. This region is bounded at large ${g}_{\max }$ by Equation (8) where the cooling is too slow, and below by ${g}_{\max }{\sigma }_{g}\sim 3$, where the g profile is too narrow to trap spiral modes. Lowering β extends the right boundary of the unstable region to higher ${g}_{\max }$ values. Using Figure 9 and Equations (19)–(20), we can thus estimate whether any given Σ and ${c}_{\mathrm{iso}}^{2}$ profiles will become unstable.

Zoom In Zoom Out Reset image size

One limitation of our Fiducial simulation is that our Σ profile has both an inner hole and an outer taper. Such a model may be representative of observed transition disks, which are characterized by large inner cavities (e.g., Muzerolle et al. 2010; Espaillat et al. 2014; Ansdell et al. 2016, 2017). However, transition disks make up only ≲20% of observed protoplanetary disks. In Figure 10, we show that monotonic density profiles, i.e., disks without inner holes, can also produce peaked g(r) profiles capable of trapping unstable s-modes. Hence, the s-mode instability may be a generic instability capable of producing multiple rings across a variety of observed protoplanetary disks.

Zoom In Zoom Out Reset image size

Our results are independent of what physical length units we give to R₀. Thus, we can take the Σ profiles shown in Figures 1 and 10 and apply them to realistic disks. Recall that the s-mode instability requires both a large ${g}_{\max }$ and small β. For a given disk mass, we can place a constraint on what H/r is required to produce a ${g}_{\max }\gtrsim \,10$. For example, a transition disk (similar to Figure 1) with a cavity inside of 50 au and disk mass ∼0.05 M_⊙ would need H/r ≲ 0.06 at 50 au to attain ${g}_{\max }$ ≳ 10. Similarly, a nontransition disk (similar to Figure 10) with an exponential cutoff outside 150 au at the same mass would require H/r ≲ 0.05 at 150 au to attain ${g}_{\max }$. To see whether such disks will actually be unstable requires an estimate of the cooling rate, which we focus on in the next subsection.

So far, we have focused exclusively on cooling timescales with constant β. In contrast, more realistic, radiative cooling that takes into account various opacity sources in the disk produces a nonconstant β profile. A more general β function affects both the ability of the disk to trap spiral modes, through the g profile, as well as the criterion for instability (Equation 8).

From Zhang & Zhu (2020), we can assign an effective β to a more realistic radiative cooling model,

$$\begin{eqnarray}\begin{array}{rcl}\beta & = & 0.015{\left(\displaystyle \frac{f}{0.01}\right)}^{-1}{\left(\displaystyle \frac{{\kappa }_{R}}{1{\mathrm{cm}}^{2}/{\rm{g}}}\right)}^{-1}{\left(\displaystyle \frac{{L}_{\star }}{{L}_{\odot }}\right)}^{1/2}{\left(\displaystyle \frac{\phi }{0.02}\right)}^{-3/4}\\ & & \times {\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{1/2}(1+{\tau }^{2}),\end{array}\end{eqnarray} \tag{ 21 }$$

which depends on the central star properties (M_⋆, L_⋆), the dust-to-gas ratio f, the flaring angle, ϕ, the disk opacity, κ, and the optical depth, τ. For typical protoplanetary disk parameters, β can easily be ≪1, i.e., realistic cooling may destabilize a disk with a trapped spiral mode. Even smaller β can be achieved by raising the dust-to-gas ratio or the opacity, or by looking at disks around less massive or less luminous stars.

The above estimation does not take into account the effects of nonconstant β on the disk stability. In Appendix A, we derive the correction to the disk normal mode equation (Equation A9) for a radius-dependent β. To the lowest order in d β/dr, there is an extra torque on a slow eccentric mode proportional to (d β/dr)(dc²/dr)Σ∣E∣². This torque modifies the criterion for instability to

$$\begin{eqnarray}&&\beta \lt {\beta }_{c}=\displaystyle \frac{1}{g(\gamma -1)}\left(-\displaystyle \frac{d\mathrm{ln}{c}^{2}}{d\mathrm{ln}r}\right)\left[\mathrm{sign}(k)-\displaystyle \frac{1}{g}r\displaystyle \frac{d\beta }{{dr}}\right].\end{eqnarray} \tag{ 22 }$$

Depending on the sign of d β/dr, this extra term can either help stabilize or further destabilize a spiral mode. Note that this additional torque is independent of the mode type, i.e., elliptical modes also experience a torque proportional to (d β/dr)(dc²/dr). This raises the intriguing possibility that elliptical modes may also be unstable if (d β/dr)(dc²/dr) < 0.

In this work, we have only considered 2D disks. However, the 3D simulations in Lin (2015) mostly confirmed the 2D results. He found that the s-mode instability can still occur in 3D disks, but with a smaller growth rate. The spiral pattern varies with height, while the midplane pattern is similar to the 2D results.

For the linear theory, an additional pressure-like term on the order of ∼h² E should be added to Equation (4) in 3D (Ogilvie 2008). Hence, the ${g}_{\max }$ and β_c criteria for the spiral instability may change. But, our ring-generation mechanism is derived for an arbitrary m = 1 perturbation until we plug in the eccentricity formalism in the last step, thus we expect the ring formation to occur in 3D disks as well. Nevertheless, long-term 3D simulations still need to be done in the future to explore the important nonlinear outcomes and final ring configuration.

Our model manifests a variety of disk substructures that may be observable in dust continuum images. To get a sense of the dust response, we have performed a locally isothermal simulation with dust (see, e.g., Fu et al. 2014 for a description of the dust implementation in LA-COMPASS). The dust is initialized to the same distribution as the gas with a dust-to-gas ratio of 0.01 and is not allowed to feedback onto the gas. For simplicity, we only consider dust with a diameter of 0.2 mm and internal density of 1.25 g cm⁻³. By assuming a solar-type star and a reference radius R₀ = 20 au, the 0.2 mm grains have a Stokes number of ∼10⁻³ at R₀.

Figure 11 shows that the dust largely resembles the gas pattern. The $\left\langle {\rm{\Sigma }}\right\rangle$ plots demonstrate that the primary dust ring is formed at the primary gas ring (r ∼ 1.2R₀). The two comparable dust rings coincide with the secondary and tertiary gas ring (r ∼ 1.5, 2.3R₀). We expect the dust rings to be wider in viscous disks. This result suggests that our substructure formation mechanism may create significant observational features in the dust continuum. Further work including a range of dust sizes and dust feedback is required to make more realistic observational predictions.

Zoom In Zoom Out Reset image size

Kinematical signatures of disk eccentricity may also be observable. During the spiral instability stage in our model, the eccentric motions of the gas reach ∼10% of the Keplerian motion, and in the final stage, the rings may still possess some residual eccentricity. Gas kinematics traced by any molecular line emission can verify whether or not the mechanism we present in this paper is the true origin of some observed disk substructures. A caveat here is that we have considered only 2D disks. As we discussed in Section 4.3, while the starlight is supposedly scattered at the top and bottom surfaces of the disk, it is not clear if the disk eccentricity should be independent of height.

Observations of a spiral instability in disks may also provide new constraints on the mass, opacity, and other basic disk parameters because the instability depends on the disk properties such as the self-gravity to pressure ratio and cooling efficiency.