Ring Formation in Protoplanetary Disks Driven by an Eccentric Instability
DDraft version February 5, 2021
Typeset using L A TEX twocolumn style in AASTeX63
Ring Formation in Protoplanetary Disks Driven by an Eccentric Instability
Jiaru Li ,
1, 2
Adam M. Dempsey , Hui Li , and Shengtai Li Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA Center for Astrophysics and Planetary Science, Department of Astronomy, Cornell University, Ithaca, NY 14853, USA (Accepted January 29, 2021)
Submitted to ApJABSTRACTWe find that, under certain conditions, protoplanetary disks may spontaneously generate multiple,concentric gas rings without an embedded planet through an eccentric cooling instability. Using bothlinear theory and non-linear hydrodynamics simulations, we show that a variety of background statesmay trap a slowly precessing, one-armed spiral mode that becomes unstable when a gravitationally-stable disk rapidly cools. The angular momentum required to excite this spiral comes at the expenseof non-uniform mass transport that generically results in multiple rings. For example, one long-termhydrodynamics simulation exhibits four long-lived, axisymmetric gas rings. We verify the instabilityevolution and ring formation mechanism from first principles with our linear theory, which showsremarkable agreement with the simulation results. Dust trapped in these rings may produce observablefeatures consistent with observed disks. Additionally, direct detection of the eccentric gas motions maybe possible when the instability saturates, and any residual eccentricity left over in the rings at latertimes may also provide direct observational evidence of this mechanism.
Keywords: accretion, accretion disks — hydrodynamics — instabilities — methods: analytical —methods: numerical — protoplanetary disks INTRODUCTIONObservations of protoplanetary disks have revealed adiversity of substructures, ranging from annular rings(e.g., ALMA Partnership et al. 2015; Andrews et al.2016; Isella et al. 2016; Huang et al. 2018a), to lop-sided vortices (e.g., Fukagawa et al. 2013; van der Marelet al. 2013) and spiral arms (e.g., Huang et al. 2018b,c).Among these, the most common substructures are con-centric rings and gaps (Andrews et al. 2018; Garufi et al.2018), which may provide crucial insights to our un-derstanding of protoplanetary disk evolution and planetformation. Yet, their origins remain open questions.Unseen planets are commonly invoked to explain thesefeatures (e.g., Dipierro et al. 2015; Dong et al. 2015,2017; Jin et al. 2016; Zhang et al. 2018; Wafflard-Fernandez & Baruteau 2020), as they may excite spi-ral waves (Goldreich & Tremaine 1979) and carve an-
Corresponding author: Jiaru Lijiaru [email protected] nular gaps in the disk under certain conditions (Lin& Papaloizou 1986; Li et al. 2005). To produce ringsat the observed location in their host disks, many pro-posed gap-opening planets must be at a few tens of AU(or even a hundred) from the star (e.g. Dipierro et al.2015; Zhang et al. 2018). However, without pre-exsitingsubstructures, dust and solid material can quickly driftfrom large radii to the inner parts of the disk (Weiden-schilling 1977). Given a low solid density environmentand a long orbital timescale at large disk radii, the cur-rent core-accretion model of planet formation takes morethan a million years to form planets (Helled & Boden-heimer 2014; Morbidelli 2020). Given this, some ob-served ringed disks may be too young to assemble theproposed planets (e.g., ALMA Partnership et al. 2015;Segura-Cox et al. 2020).On the other hand, there are several ring-formingmechanisms that do not require planets. These in-clude a dust-gas viscous gravitational instability (Tomi-naga et al. 2019), and a magneto-hydrodynamical wind-driven instability (Riols & Lesur 2019; Riols et al. 2020). a r X i v : . [ a s t r o - ph . E P ] F e b Li et al.
Additionally, dust rings may occur at the edges of deadzones (Flock et al. 2015) and snow lines (Okuzumi et al.2016).Another interesting possibility is that rings may beformed via angular momentum exchange with an unsta-ble non-axisymmetric mode. One candidate for this isan eccentric mode (Lubow 1991; Ogilvie 2001; Lee et al.2019a,b). In this context, an eccentric mode refers tothe coherent precession of gas on eccentric orbits acrossdifferent radii. The precession can be caused by bothpressure and self-gravity.It has been shown recently that almost any disk witha realistic density profile can sustain long-lived eccen-tric modes (Lee et al. 2019b). In non-self-gravitatingdisks, the linearized equation of motion for the complexeccentricity can be mapped to a Schr¨onger-like equa-tion, which shows that free eccentric eigenmodes canbe trapped in the disk like a particle in a potentialwell (Ogilvie 2008; Saini et al. 2009; Lee et al. 2019b;Mu˜noz & Lithwick 2020). These modes have discretespectra and appear as apse-align elliptical rings. Withdisk self-gravity, Tremaine (2001) and Lee et al. (2019a)found that disks can also support degenerate eccentricmodes which appear as single-arm spirals. Forced ec-centric modes have also been found to exist in diskscoupled with both low-mass (e.g., Papaloizou 2002;Teyssandier & Ogilvie 2016) and high-mass secondaryobjects (Mu˜noz & Lithwick 2020).A global eccentric mode can grow its amplitudethrough several mechanisms, such as: the SLING mech-anism (Adams et al. 1989; Shu et al. 1990) that ampli-fies an eccentric perturbation through the wobble of thecentral star, and instantaneous cooling (Lin 2015).In this paper, we present a theory of spontaneous(planet-free) gas ring formation through a spiral eccen-tric instability due to thermodynamic cooling. We findthat this instability can produce a variety of gas featuresin just a few thousand orbits: from an initial spiral, totransient ellipses, and finally to a series of concentricrings.The rest of this paper is organized as follows. InSection 2, we present a hydrodynamical simulationthat demonstrates the three substructures in one diskthrough its evolution. In Section 3, we formulate a lin-ear theory for eccentric modes in disks with a finite cool-ing rate, and we use this theory to explain the growthand saturation of the spiral in the simulation. In partic-ular, Section 3.2 discusses how an unstable one-armedspiral generically results in multiple, concentric rings.Section 4 discusses the astrophysical applicability anddetectability of our mechanism. In Section 5, we give asummary of our findings. FIDUCIAL EVOLUTION2.1.
Disk model
We simulate a 2D self-gravitating gaseous disk withthe hydrodynamics code
LA-COMPASS (Li et al. 2005,2009a).
LA-COMPASS evolves the fluid equations witha Godunov, shock-capturing scheme. In addition toevolving the continuity equation and momentum equa-tions for the gas,
LA-COMPASS also evolves the gas energyequation with an additional cooling term, ∂E g ∂t + ∇ · (( E g + P ) v ) = − v · ∇ Φ − c v Σ( T − T eq ) t c , (1)where Σ, P , v = u ˆ r + v ˆ φ , and Φ are the gas surfacedensity, pressure, 2D velocity in the r - φ midplane, andtotal gravitational potential, respectively. The specificheat capacity is c v = k b / ( µ ( γ − γ is the ratio of spe-cific heats, k b is the Boltzmann constant, and µ is themean particle mass which we fix to µ = 2 . m p where m p is the proton mass. The total gas energy is givenby E g = Σ | v | / P/ ( γ − T , via the usual ideal gasequation of state, P = ( k b /µ )Σ T .Gas cools to an equilibrium temperature T eq on atimescale t c that is a fixed multiple of the orbitaltimescale given by the parameter β = t c Ω K , where Ω K isthe Keplerian rotation rate (Gammie 2001). In the mo-mentum equation, we include the gravitational potentialof a central star with mass M (cid:63) , as well as the disk’s selfgravity (Li et al. 2009b). We do not include the indirectpotential associated with the acceleration of the star dueto the pull of an off-center disk. We have checked that asimulation including the indirect potential finds nearlyexactly the same results as those presented in Section2.3.The cooling model we adopt is meant to approximatethe more accurate radiative cooling models thought togovern real protoplanetary disks. We parameterize thecooling timescale with the relatively simple β model thathas been used extensively in studies of self-gravitatingdisks (e.g., Gammie 2001; Lodato & Rice 2004, 2005;Kratter & Lodato 2016). This model has been shownto produce similar results to more sophisticated cool-ing models and makes an analytic treatment of our re-sults possible in Section 3. Additionally, we also spec-ify an equilibrium temperature profile rather than self-consistently computing one by balancing cooling withstellar and internal heating. However, we do not expectthis to impact our results significantly because more re-alistic equilibrium temperature profiles tend to be simplepower-laws (e.g., Chiang & Goldreich 1997).While protoplanetary disks can be viscous in general,observations suggest many disks may have quite low vis- ings formed by an eccentric instability Initial conditions
For the initial surface density, we choose a power-lawprofile with both an inner hole and an outer taper,Σ( r ) = 2 . (cid:16) − e − ( r/R ) (cid:17)(cid:124) (cid:123)(cid:122) (cid:125) inner hole power − law (cid:122) (cid:125)(cid:124) (cid:123)(cid:18) R r (cid:19) e − ( r/ (2 R )) (cid:124) (cid:123)(cid:122) (cid:125) outer taper , (2)We choose our reference radius to be R and referencedensity to be Σ = Σ( R ).The initial temperature is set via the gas sound speed, c s ( r ) = γ ( k b /µ ) T = c ( r/R ) − / . We set T eq to thisinitial, power-law temperature. For all of our simula-tions we set γ = 1 . c = 0 . (cid:112) GM (cid:63) /R . Ad-ditionally, we set our unit of time to be the Keplerianfrequency at R , i.e, Ω K ( R ) = (cid:112) GM (cid:63) /R .The simulations are performed in a two-dimensionalcylindrical domain from R in = 0 . R out = 6 . N r , N φ ) = (4096 , r ∈ [0 . , . ∪ [5 . , .
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 non-zero initial radial velocity u is sampled uniformly from U nn o r m a li z e d , c s c r / R T oo m r e Q Q = 2 = 1.6 × 10 M R Figure 1.
Initial Σ and c (top) and radial Q profiles (bot-tom) used in our study. Note that we only consider gravita-tionally stable disks with Q > Figure 2.
Snapshots of the
Fiducial hydrodynamics simu-lation at time = 1500 orbits (top), 3500 orbits (middle), and9000 orbits (bottom). The left column shows disk surfacedensity. The right column shows the m (cid:54) = 0 component ofthe density deviation in the r - φ coordinates. [ − − , +10 − ] r Ω K in each fluid cell as an initial per-turbation. This perturbation could be caused by vari-ous processes, such as random fluid turbulence (Balbus& Hawley 1998; Flaherty et al. 2020). We have checkedthat the disk evolution is independent of the type andamplitude of the initial perturbation so long as there ispower in the m = 1 azimuthal component.2.3. Fiducial evolution
Our
Fiducial disk adopts Σ = 1 . × − M (cid:63) R − and β = 10 − . Hence, the disk has a total mass 0 . M (cid:63) and is essentially locally isothermal . With radiativecooling, typical protoplanetary disks can have β (cid:28) and β are explored in Section 3.1.Figure 1 shows our initial disk configuration. The sur-face density peaks just outside of R coinciding witha minimum in the Toomre Q = c s Ω / ( πG Σ) (Toomre Despite the name, essentially locally isothermal is fundamentallydifferent from strictly locally isothermal, which does not evolvean energy equation, as the latter does not conserve wave angularmomentum (Lee 2016; Miranda & Rafikov 2019).
Li et al.
Q >
2, to ensure that any instability wefind is not due to a small Q .2.3.1. Modal time evolution
We let the disk evolve for 10 orbits at R to providea general picture of the long-term evolution. This corre-sponds to roughly one million years around a solar massstar if R = 21 . r ∼ − R .To get a sense of the relative importance of the differ-ent azimuthal symmetries in the disk, we show the timeevolution of the global Fourier amplitudes, C m ( t ) ≡ M d (cid:90) R out R in (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) π Σ e − imφ dφ (cid:12)(cid:12)(cid:12)(cid:12) rdr, (3)in the upper panel of Figure 3. The coefficients arenormalized to the m = 0 coefficient which correspondsto the total disk mass, M d .At early times, the disk surface density variations aredominated by the m = 1 component. From the C curve and the snapshots, this component appears as anexponentially growing, one-armed, trailing spiral withgrowth rate γ ≈ − . Upon measuring the space-timeevolution of the m = 1 complex phase (shown in the bot-tom panel of Figure 3), we find that the spiral precesses coherently and in a prograde direction with a patternspeed of Ω p ≈ . p (cid:28) Ω and the preces-sion rate is coherent, what we are seeing is the growthof a slow eccentric mode (Tremaine 2001; Lin 2015; Leeet al. 2019a).During this time, higher m modes are also exponen-tially growing with growth rates that are γ m ≈ mγ ,albeit at an amplitude much less than C . This relationsuggests that these modes are harmonics of the m = 1mode and may be driven by non-linear, mode-mode in-teractions (Laughlin et al. 1997, 1998; Lee 2016).By ∼ C saturatesto a strength of ∼ m = 1 pattern speed.After ∼ r < R ). All modes start to decay fromtheir peak amplitude on a timescale of ∼ /m orbits,and the rings transition to being mostly axisymmetric. While fading away, the m = 1 component maintains acoherent pattern speed of Ω p ≈ − .2.3.2. Ring evolution
Figure 4 shows the azimuthally-averaged Σ (denotedby (cid:104) Σ (cid:105) ) in the Fiducial simulation. During the growthphase (upper panel), (cid:104) Σ (cid:105) develops three well-definedover-densities that are at radii r ≈ . , . , . R andconsistently grow over time. Once the instability satu-rates (middle panel), there is significant non-linear evo-lution of (cid:104) Σ (cid:105) that changes the ring locations. In the finaldecay phase (bottom panel), (cid:104) Σ (cid:105) settles into a quasi-steady-state configuration consisting of three prominentrings at r = 1 . , . , . R , and additional, weaker fea-tures at r = 0 . R and 1 . R . Even though (cid:104) Σ (cid:105) in-creases in the rings, the disk maintains Q (cid:38) Fiducial simulation, we have shown thata disk with inner and outer tapers traps an unstableone-arm spiral. This growing spiral arm has been previ-ously reported by Lin (2015) for locally isothermal disks( γ = 1 and β →
0) and in cooled disks with β (cid:46) .
1, al-though the latter result was not discussed in detail andused preliminary simulations. We have shown, in detail,that a non-zero cooling time unambiguously results in asimilar instability. But more importantly, we find thatthe non-linear outcome of this spiral instability is thegeneration of concentric axisymmetric rings.In reality, disk viscosity can affect the growth of thespiral and the lifetime of the rings. In our
Fiducial simulation, viscosity is negligibly small. If we were to in-clude a non-zero viscosity ν , it could suppress the growthof the spiral if the damping rate is faster than the growthrate. Moreover, a large ν may diffuse away the initialΣ profile before any growth could occur. Preliminaryshort-term simulations of the growth phase with con-stant viscosity find that ν (cid:46) − , which correspondsto α (cid:46) .
01 at R , gives the same growth rate and spi-ral profile as our Fiducial results. We expect that ina long-term viscous simulation, the lifetime of the ax-isymmetric rings will be longer than the lifetime of thespiral waves as a tightly-wound spiral damps at a ratethat is faster than the rings by a factor of ( kw ) > k is the radial wave number of the spiral and w is the radial width of a typical ring.In the next sections, we show from first principles whythe disk is unstable and why that generically results inconcentric rings. THEORY OF SPIRAL INSTABILITY ANDCOMPARISON TO THE FIDUCIAL RESULTS ings formed by an eccentric instability C m m=1m=2m=3m=4m=5 Figure 3. Top:
Time evolution of the global Fourier mode amplitudes for m = 1 − Fiducial simulation. The redtriangles marks the time of the snapshots in Figure 2. The dashed black vertical lines divide the whole evolution into (a) theinitial growth of the spiral, (b) the saturation and morphology transition to closed eccentric rings, and (c) final multi-ring stage.
Bottom:
Space-time diagrams of the local m = 1 Fourier mode phase, φ , from 1000 to 6000 orbits. The opaqueness of thecolormap scales with the eccentricity at ( r, t ). The m = 1 pattern is precessing at a coherent rate of Ω p = 0 .
018 when its is aspiral and a different rate of Ω p = 0 .
010 after saturation.
Our theory is based on the linearized equation of mo-tion for eccentricity in an adiabatic disk with a finitecooling rate (Equation 1). In Appendix A we show that normal modes of the disk, e.g. u ∝ e − iωt + imφ , that areboth slow ( (cid:60) ( ω ) (cid:28) Ω K ) and eccentric ( m = 1) satisfythe eigenvalue equation,2 r Ω K Σ ω = 11 + β ddr (cid:26) r P (cid:20) c ddr (cid:18) Ec (cid:19) + γβ dEdr (cid:21) + iβr P (cid:20) ( γ − dEdr + d ln c dr E (cid:21)(cid:27) + r dPdr E (cid:124) (cid:123)(cid:122) (cid:125) pressure with β cooling − disk self-gravity (cid:122) (cid:125)(cid:124) (cid:123) Σ r ddr (cid:0) r d Φ dr (cid:1) E − Σ ddr (cid:0) r Φ (cid:1) (4)where ω is the complex mode frequency (the “eigen-value”), E is the complex eccentricity (the “eigenvec- tor”), and c = P/ Σ is the isothermal sound speed
Li et al. / / r / R / T i m e ( o r b i t ) Figure 4.
Time evolution of the azimuthally averaged den-sity profile (cid:104) Σ (cid:105) in the Fiducial simulation. The three panelsroughly correspond to the three phases of evolution shownin Figure 2: one-armed spiral (top), elliptical rings (middle),circular rings (bottom). squared. The eccentricity equation includes the contri-bution from the m = [0 ,
1] self gravity potentials, Φ [0 , ,and the effects of cooling through β .Many previous authors have derived the governingequation of the disk eccentricity by adopting either theadiabatic (e.g., Goodchild & Ogilvie 2006; Lee et al.2019a,b) or the locally isothermal approximation (e.g.,Teyssandier & Ogilvie 2016; Mu˜noz & Lithwick 2020).Their results can be recovered from our Equation (4) bytaking the limits β → + ∞ and β →
0, respectively.We numerically solve the boundary-eigenvalue prob-lem given by Equation (4) for the complex eccentric-ity and mode frequency using a finite difference ma-trix method (see e.g., Adams et al. 1989; Lee et al.2019a). We use a uniformly log-spaced grid that cov-ers r ∈ (0 . ,
50) with N = 2048 cells. At the bound-aries we set dE/dr = 0, but note that since our densitygoes to zero at the boundaries, our choice of boundarycondition has little effect on the modes. Solving the ec-centricity equation gives a set of N modes. We focusonly on the mode which has the fastest growth rate asthis one should be present in the hydrodynamical simu-lation.Figure 5 compares the eigenfunction E that we ob-tain for the Fiducial disk model to the eccentricity inthe simulation E sim (at 1500 orbits). The linear theorymatches the hydrodynamics result almost perfectly. Theargument of periapse, arg( E ), of the eccentricity varies r / R E cc e n t r i c i t y Im( E )Re( E ) Eigenmode E ( r )Hydro result E sim Figure 5.
Radial profile of the (unstable) fundamental ec-centric mode E ( r ) in the Fiducial disk from the linear the-ory (Equation 4; solid green) and the eccentricity E sim fromhydrodynamics simulation at t = 1500 orbits (dashed blue).The real and imaginary parts are shown as opaque and trans-parent lines, respectively. The eigenfunction is multiplied bya real amplitude and a complex phase factor to fit the sim-ulation result. with r , which produces the single-arm spiral that we seein the 2D gas density.3.1. Characteristics and growth of slow eccentric modes
Slow eccentric modes in self-gravitating, adiabaticdisks have been thoroughly studied in Lee et al. (2019a).They found that there are three kinds of modes that aredifferent in morphology and controlled by the profile of g ( r ) = πG Σ rc , (5)which measures the relative strength of pressure andself-gravity .We show a collection of the different mode types indispersion relation maps (DRM) in Figure 6. The DRMplots show contours of fixed (cid:60) ( ω ) in r - kr space, where k is the radial wave number of the mode. The contours arecalculated from a WKB-like approximation to the eccen-tricity equation (see Appendix B and Lee et al. 2019afor more details). Following the naming convention ofLee et al. (2019a) the three mode types are: • Pressure dominated elliptical modes (e-p modes)occur when g <
1. These are long wavelengthretrograde ( (cid:60) ( ω ) <
0) modes that are trapped bytwo inner Lindblad resonances (locations were k =0). In disks with cooling, the sound-speed used in the definition of g is a mix of the isothermal and adiabatic values controlled by β (see e.g., Equation (B12) in Appendix B). Since we mainly focuson rapidly cooled disks, we use the isothermal sound speed in g . ings formed by an eccentric instability K k r g max = 0.1 K k r g max = 1.2 Lindblad resonancesQ-barriers10 r / R K k r g max = 13 trailing ( k > 0) onlyleading ( k < 0) only Figure 6.
Dispersion relation map (DRM) for the eccen-tric modes in Equation (4) based on the WKB dispersionrelation derived in Appendix B. The closed curves show thecontours of constant ω . Each panel adopts a different disk g max (Equation 5). The darkest contour represents the fun-damental mode, while the higher harmonics are displayed asmore transparent contours. The elliptical modes are shownin black. The leading and trailing spiral modes are shownin blue and red, respectively. We also mark the Lindbladresonances (blue dots) and the Q-barriers (orange dots) forthe fundamental mode defined by Equation (B12). • Self-gravity dominated elliptical modes (e-pgmodes) occur when g > kr = g ). • Spiral modes (s-modes) occur when g (cid:29) k <
0) and trailing ( k >
0) spi-rals trapped by Q-barriers at short wavelengths.Each of these contours are separate, i.e, there areno Lindblad resonances connecting them. Themode seen in the
Fiducial simulation is a trailings-mode.To determine whether a given mode is unstable welook at the equation governing the modes angular mo-mentum deficit, Σ r Ω | E | , which can be obtained frommultiplying the eccentricity equation by complex con- jugate of E , E ∗ , and integrating in r over the wholedisk. The growth rate can then be determined fromthe imaginary component of this equation. We show inAppendix C that doing so results in two terms whichcontribute to the growth rate γ = (cid:61) ( ω ): γ (cid:90) r ΣΩ | E | dr = − (cid:90) r P (cid:2) β ( γ − | dEdr | (6) − d ln c dr (cid:61) (cid:26) (1 − iβ ) E dE ∗ dr (cid:27) (cid:3) dr, to first order in β . For an elliptical mode which has E ∝ cos( kr ), the integrand on the right-hand side isproportional to − k cos ( kr ), so that there is no growth( γ ≤ 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 approximatethe integral at the radius where g = g max . Doing so, wefind that γ is approximately, γ ≈ − g max h Ω(2 γ ) [( γ − βg max + d ln c d ln r sign( k )] , (7)where h = c s / ( r Ω) and we approximate kr ≈ g max .When β is small, γ is approximately linear to g max .One can see that only an s-mode can be unstable . Inparticular, only modes with kdc /dr < β must satisfy, β (cid:46) β c ≡ − γ sign( k ) g max d ln c d ln r . (8)Because slower cooling acts to suppress the instability,the disk must cool fast enough for an initial spiral modeto grow.Figure 7 compares the growth rates between our linearresults (solid curve) and a set of hydrodynamical sim-ulations (blue dots) for different values of g max and β .Both the linear and the hydrodynamics results suggestthat an unstable mode exists when g max > g max >
12. Between g max = 9 and12, the results differ by 71% and 16%. This discrepancyis likely associated with the g max being very close to theinstability threshold. All of the larger g max simulationsremain gravitationally stable, Q (cid:38)
2, throughout theirevolution.For fast enough cooling rate (bottom panel), thegrowth rate γ is almost a constant. For β close tothe critical value β c , the growth rate drops sharply fromnearly maximum to zero as β increases. The hydrody-namic results suggest that β c is between 0 .
06 and 0 . Li et al. g max Fiducial
Eq. 6, same as Im( ) of Eq.4Approximation (Equation 7)Hydro Simulation Result0.0008 0.0012 0.0016 0.002 0.0024 [ M R ] m a x ( , ) Eq. 6, same as Im( ) of Eq.4Approximation (Equation 7)Hydro Simulation Result
Figure 7.
Growth rate of the fundamental eccentric modeas a function of g max (top panel; all disks have β = 10 − )and β (bottom panel; all disks have g max = 12 .
6. The greenline shows the rate calculated from the linear theory (Equa-tion 4 and 6). The blue star and dots represent the ratesmeasured from non-linear hydrodynamics simulations. Thedashed green line shows the approximated prediction of g max and β -dependence (Equation 7). We plot max(0, γ ) when β > β c because our simulations do not measure the suppres-sion rate. Equation (6) and Equation (8) yields β c = 0 .
08 and 0 . β ≤ . γ for g max >
12. The main error introducedwhen going from Equation (6) to (7), for g max > E ( r ) pro-file, only overestimates γ by ∼ − g max , Equation (7) does not apply as the s-mode is eithernon-existent or not prominent.In this section we have shown that an instability devel-ops in a self-gravitating disk with cooling when (1) theratio of self-gravity to pressure (i.e., g ) is large enoughto support a trapped spiral mode, (2) the spiral co-incides with a non-zero temperature gradient, and (3)the cooling rate is faster than the critical rate given inEquation (8). Moreover, we have derived an analytical approximation for the resulting growth rate from thedispersion relation for the gas. In the limit of β → c profiles, the spiral instability mechanism appliesto a broad range of astrophysically interesting densityand temperature profiles. This is because the appear-ance of the s-modes is controlled by the radial profile of g , not Σ or c individually. As we discuss in detail inSection 4.1, this results in a wide range of density andtemperature profiles that may be s-mode unstable.3.2. Ring formation in response to a growing spiral
From Figure 4, we know that when there is a growingone-armed spiral, the disk also produces axisymmetricrings. To understand this process, we examine how an-gular momentum is transferred from the spiral into themean flow of the disk, and how this deposition processgenerically results in a non-constant radial mass fluxthat drives the formation of multiple rings.In general, angular momentum is stored in the meanflow of the disk as well as any waves present. Angularmomentum exchange occurs between these two reser-voirs during wave excitation and when there is non-zero dissipation (from e.g., viscosity or shocks) or baro-clincity (i.e., ∇ P × ∇ Σ (cid:54) = 0). These latter two processesare equivalent to the loss of conservation of the diskvortensity, ξ = ˆ z · ( ∇ × v ) / Σ.In inviscid disks with cooling, baroclinic effects trans-fer angular momentum from the mean flow, whose massand angular momentum evolve as,2 πr ∂ (cid:104) Σ (cid:105) ∂t + ∂∂r ( − ˙ M ) = 0 (9)2 πr ∂ (cid:104) Σ (cid:105) (cid:104) v (cid:105) ∂t + ∂∂r (cid:16) − ˙ M r (cid:104) v (cid:105) (cid:17) = − t dep (10)into a wave, whose angular momentum follows πr ∂ (cid:104) Σ (cid:48) v (cid:48) (cid:105) ∂t + ∂∂r (cid:0) πr (cid:104) Σ uv (cid:48) (cid:105) + F G (cid:1) = t dep (11)where ˙ M = − πr (cid:104) Σ u (cid:105) is the radial mass flux and F G isthe flux of angular momentum due to self-gravity (Gol-dreich & Tremaine 1979). The exchange torque den-sity between the wave and mean flow is (Dempsey et al.2020), t dep πr = − (cid:104) u (cid:48) Σ (cid:48) (cid:105) ∂ ( r (cid:104) v (cid:105) ) ∂r + (cid:104) Σ (cid:105) (cid:28) u (cid:48) ∂ ( rv (cid:48) ) ∂r (cid:29) − (cid:28) Σ (cid:48) Σ ∂P∂φ (cid:48) (cid:29) . (12) Note that 2 πr (cid:104) Σ (cid:48) v (cid:48) (cid:105) is the angular momentum in the wave inthe Eulerian frame, i.e., at fixed radius. This is different than thewave angular momentum, − πr ΣΩ | E | , which can be calculatedusing the Lagrangian perturbations (Lin 2015). ings formed by an eccentric instability t dep for an eccentric wave can be shown tobe (Lubow 1990, see also Appendix D), t dep πr = − γ r Σ ∂ (cid:104) ξ (cid:105) ∂r | E | + r ∂c ∂r Σ (cid:61) (cid:18) E ∗ r ∂E∂r (cid:19) + β ( γ − r c Σ (cid:12)(cid:12)(cid:12)(cid:12) ∂E∂r (cid:12)(cid:12)(cid:12)(cid:12) , (13)in disks with β (cid:28)
1. The first term depends on thebackground vortensity gradient and is zero for stablemodes, while the remaining terms are exactly the sameas the terms on the right-hand-side of Equation (6) thatset the growth rate. We see now that they correspond toa non-zero, persistent background torque proportionalto the disk’s temperature gradient (Lin & Papaloizou2011; Lin 2015), and a stabilizing torque proportionalto the cooling rate.The loss of angular momentum from the mean flowinduces a non-zero radial mass flux in the location ofthe wave. From Equations (9)-(10) this ˙ M is,˙ M = 2 r Ω K (cid:20) πr (cid:104) Σ (cid:105) ∂ (cid:104) v (cid:105) ∂t + t dep (cid:21) . (14)The ∂ t (cid:104) v (cid:105) term is mostly determined by the disk eccen-tricity, which is ultimately sculpted by t dep . We esti-mate ∂ t (cid:104) v (cid:105) by assuming the disk maintains an approxi-mate radial force balance by modifying the pressure andself-gravity corrected Keplerian rotation, Ω , by a smallamount, Ω = Ω + Ω . From the azimuthally averagedradial velocity equation the Ω required to maintain bal-ance is (Lubow 1990),2 r Ω Ω = 1 (cid:104) Σ (cid:105) ∂ (cid:104) P (cid:105) ∂r + ∂ (cid:104) Φ (cid:105) ∂r − r Ω − (cid:28) v (cid:48) r (cid:29) + (cid:28) v (cid:48) r ∂u (cid:48) ∂φ (cid:29) + (cid:42)(cid:18) (cid:19) (cid:48) ∂P (cid:48) ∂r (cid:43) + (cid:28) u (cid:48) ∂u (cid:48) ∂r (cid:29) . (15)The first three terms on the right-hand side sum to zeroas they constitute the dominate radial force balance. Forslow eccentric modes the pressure term is ∼ h smallerthan the velocity terms and so we can neglect it. Be-cause we are interested in unstable modes, the remain-ing RHS terms grow as e γ t , so that the ∂ t (cid:104) v (cid:105) term inEquation (14) is approximately, ∂ (cid:104) v (cid:105) ∂t ≈ γ r Ω (cid:20) − (cid:10) v (cid:48) (cid:11) + (cid:28) v (cid:48) ∂u (cid:48) ∂φ (cid:29) + (cid:28) u (cid:48) r ∂u (cid:48) ∂r (cid:29)(cid:21) ≈ − γ r Ω K | E | + γ r Ω K ∂∂r | E | (16) r / R M [ R K ( R )] Ring forming spot (convergent flow) M ( theory; -- hydro) v / t term (Eq.17a) / r effect (Eq.17b)baroclinic effect (Eq.17c) Figure 8.
Radial mass flux, ˙ M , at 1500 orbits in our Fidu-cial simulation (dashed black) and theoretical predictions(solid lines) from Equation 17a-17c). The colored curvesshow the contributions to ˙ M from the ∂ t (cid:104) v (cid:105) component(Equation 17a), the ∂ r (cid:104) ξ (cid:105) component (Equation 17b), andthe baroclinic component (Equation 17c). The E ( r ) profileis calculated as in Section 3 with the disk initial profiles andscaled to | E | in the snapshot. In Equations (17a) - (17b),we use the measured γ from the Fiducial simulation. Thevertical lines represent the ring forming spots in the simula-tion. where in the last line we replaced u (cid:48) and v (cid:48) with thecomplex eccentricity for an s-mode. The final induced˙ M for rapid cooling is,˙ M ≈ + 2 πr Σ γ (cid:20) −| E | + 2 r ∂∂r | E | (cid:21) (17a) − πr Σ γ Ω K ∂ (cid:104) ξ (cid:105) ∂r | E | (17b)+ 8 πr Ω K (cid:34) r Σ dc dr (cid:61) ( E ∗ ∂∂r E ) + βc ( γ − (cid:12)(cid:12)(cid:12)(cid:12) ∂E∂r (cid:12)(cid:12)(cid:12)(cid:12) (cid:35) . (17c)The terms proportional to γ were previously derived byLubow (1990) for general m , while the last two terms areunique to a β -cooled disk with an eccentric spiral. Tocomplete the calculation of ˙ M we can use the linear so-lution for E . One can then use Equation (9) to computethe mean density evolution.We show in Figure 8 that ˙ M as given by Equa-tions (17a)-(17c) agrees remarkably well with the ˙ M inthe exponential growth phase of the Fiducial simula-tion. In particular, we recover the approximate locationsof the rings, which coincide with convergent points inthe axisymmetric flow, i.e., at rising inflection points of ˙ M . The primary ring locations are set by the ∂ t (cid:104) v (cid:105) terms in Equation (17a). The vortensity (Equation 17b)and baroclinic (Equation 17c) terms mostly shift the ˙ M profile and slightly modify the ring locations. The smalldisagreement between the linear and hydrodynamical re-0 Li et al. sults is primarily due to the m (cid:54) = 1 perturbations andthe terms of O ( γ h , γ β ) that we have dropped.Overall, this excellent agreement demonstrates thatthe initial mass redistribution leading to the formationof multiple rings occurs because of the exchange of angu-lar momentum from the background disk and the spiralwave. 3.3. Post-saturation evolution
While the initial phase is linear and fully describableby our linear theory, non-linear evolution takes over atlater times (cf. Figure 4). In our
Fiducial simulation, (cid:104) Σ (cid:105) eventually evolves to the point where the instabilitysaturates at t ∼ Fiducial simulation, because an ellip-tical eccentric mode dominated by self-gravity, unlike ans-mode, tends to have larger eccentricity near the inneredge of the disk (Lee et al. 2019a). The lifetime of sucha mode will thus be short as any dissipation – whetherphysical or numerical – will operate on a faster timescalethan the dissipation of the rings at larger radii. DISCUSSION4.1.
Disk shape and temperature profile
In Section 3.1, we showed that disks are unstable whenthey can both cool faster than a critical rate and trapa slow spiral mode. Spiral modes appear when the dis-tance (in phase-space) between two Q-barriers is largeenough 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 someradius and wide enough to provide sufficient separationof the Q-barriers. This suggests that the larger scaleproperties of the Σ and c profiles are unimportant,and that only the local properties in the vicinity of the g maximum are relevant.Consider temperature and density profiles which pro-duce a peak in g at the radius r max . Taylor expanding g ( r ) around this peak produces a simple Gaussian g-profile, g ( r ) = g max e − ( r − r max ) /σ g , (18)where σ g = (cid:112) d ln g/dr ) − is the Gaussian width.For a given r max , g max , and σ g the density and tempera-ture profiles need only satisfy the following requirements g max g / R = . = . = . un s t a b l e f o r = . < . n o s - m o d e s ( g a n d g m a x t oo s m a ll ) m o r e i n s t a b li t y a s d e c r e a s e s s t a b l e Figure 9.
Stability of the eccentric mode for different com-binations of g max , σ g , and β . Bluer regions correspond tofaster growth rates, while white regions are stable. Spiralmodes may only exist above the solid line g max σ g /R ∼ γ = 0 contour) for the cooling rate β as specified.The unstable region expands as β decreases. The curves arecalculated from the linear equation (Equation 4) for g asEquation (18). The red star represents the g profile of the Fiducial disk. at r max : g max = πG Σ rc (cid:12)(cid:12)(cid:12)(cid:12) r max , (19) dd ln r ln (cid:18) r Σ c (cid:19) r max = 0 , σ − g = 12 d dr ln (cid:18) c Σ (cid:19) r max . (20)A wide range of Σ and c profiles can be constructedto fulfill these requirements.Using this g ( r ) profile, we solve for the linear γ in the( g max , σ g , β ) parameter space. Our results are shownin Figure 9 for c ∝ r − / and Σ given by Equa-tions (5) and (18). We find that there is a region of σ g - g max parameter space that is unstable to one-armedspirals. This region is bounded at large g max by Equa-tion (8) where the cooling is too slow, and below by g max σ g /R ∼
3, where the g -profile is too narrow to trapspiral modes. Lowering β extends the right boundary ofthe unstable region to higher g max values. Using Figure9 and Equations (19)-(20) we can thus estimate whetherany given Σ and c profiles will become unstable.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 transi-tion disks, which are characterized by large inner cavities(e.g., Muzerolle et al. 2010; Espaillat et al. 2014; Ansdellet al. 2016, 2017). However, transition disks make uponly (cid:46)
20% of observed protoplanetary disks. In Figure10, we show that monotonic density profiles , i.e. diskswithout inner holes, can also produce peaked g ( r ) pro-files capable of trapping unstable s-modes. Hence, the ings formed by an eccentric instability r / R U nn o r m a li z e d c g | E | Figure 10.
Example of a disk with a monotonic density,power-law sound speed, and spiral-trapping g profile with g max = 9, σ g = 0 . R , r max = R . The black curve showsthe unnormalized linear | E | profile. s-mode instability may be a generic instability capableof producing multiple rings across a variety of observedprotoplanetary disks.Our results are independent of what physical lengthunits we give to R . Thus, we can take the Σ profilesshown in Figures 1 and 10 and apply them to realisticdisks. Recall that the s-mode instability requires botha large g max and small β . For a given disk mass we canplace a constraint on what H/r is required to producea g max (cid:38)
10. For example, a transition disk (similar toFigure 1) with a cavity inside of 50 AU and disk mass ∼ . M (cid:12) would need H/r (cid:46) .
06 at 50 AU to attain g max (cid:38)
10. Similarly, a non-transition disk (similar toFigure 10) with an exponential cutoff outside 150 AU atthe same mass would require
H/r (cid:46) .
05 at 150 AU toattain g max . To see whether such disks will actually beunstable requires an estimate of the cooling rate, whichwe focus on in the next subsection.4.2. Cooling rate
So far, we have focused exclusively on coolingtimescales with constant β . In contrast, more realistic,radiative cooling that takes into account various opacitysources in the disk produces a non-constant β profile. Amore general β function affects both the ability of thedisk to trap spiral modes, through, the g-profile, as wellas the criterion for instability (Equation 8).From Zhang & Zhu (2020), we can assign an effective β to a more realistic radiative cooling model, β = 0 .
015 ( f .
01 ) − ( κ R / g ) − ( L (cid:63) L (cid:12) ) / ( φ .
02 ) − / × ( M (cid:63) M (cid:12) ) / (1 + τ ) , (21)which depends on the central star properties ( M (cid:63) , L (cid:63) ),the dust-to-gas ratio f , the flaring angle, φ , the disk opacity, κ , and the optical depth, τ . For typical pro-toplanetary disk parameters, β can easily be (cid:28)
1, i.e,realistic cooling may destabilize a disk with a trappedspiral mode. Even smaller β can be achieved by rais-ing the dust-to-gas ratio or the opacity, or by looking atdisks around less massive or less luminous stars.The above estimation does not take into account theeffects of non-constant β on the disk stability. In Ap-pendix A we derive the correction to the disk normalmode equation (Equation A9) for a radius-dependent β . To the lowest order in dβ/dr , there is an ex-tra torque on a slow eccentric mode proportional to( dβ/dr )( dc /dr )Σ | E | . This torque modifies the crite-rion for instability to, β < β c = 1 g ( γ − (cid:18) − d ln c d ln r (cid:19) (cid:20) sign( k ) − g r dβdr (cid:21) . (22)Depending on the sign of dβ/dr , this extra term can ei-ther help stabilize or further destabilize a spiral-mode.Note that this additional torque is independent of themode type, i.e., elliptical modes also experience a torqueproporitional to ( dβ/dr )( dc /dr ). This raises the in-triguing possibility that elliptical modes may also beunstable if ( dβ/dr )( dc /dr ) <
3D effects
In this work, we have only considered 2D disks. How-ever, the 3D simulations in Lin (2015) mostly confirmedthe 2D results. He found that the s-mode instabilitycan still occur in 3D disks, but with a smaller growthrate. The spiral pattern varies with height, while themid-plane pattern is similar to the 2D results.For the linear theory, an additional pressure-like termon the order of ∼ h E should be added to Equation(4) in 3D (Ogilvie 2008). Hence, the g max and β c crite-ria for the spiral instability may change. But, our ringgeneration mechanism is derived for an arbitrary m = 1perturbation until we plug in the eccentricity formalismin the last step, thus we expect the ring formation tooccur in 3D disks as well. Nevertheless, long-term 3Dsimulations still need to be done in the future to ex-plore the important non-linear outcomes and final ringconfiguration.4.4. Implications for observations
Our model manifests a variety of disk substructuresthat may be observable in dust continuum images. Toget a sense of the dust response, we have performed a lo-cally 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 distri-bution as the gas with a dust-to-gas ratio of 0 .
01 and is2
Li et al. g a s , d u s t × gasdust01234 g a s , d u s t × r [20AU]01234 g a s , d u s t × Figure 11.
2D Snapshots of the hydrodynamics simulationwith dust (left) and (cid:104) Σ (cid:105) (right) at time = 1500 orbits (top),3500 orbits (middle), and 9000 orbits (bottom). The simula-tion includes 0.2mm dust initialized to the same distributionas the gas with a 0.01 dust-to-gas ratio. The Stokes numberat R is ∼ − initially. The unit of the colorbars in theleft panel and the vertical axes in the right panel is Σ . not allowed to feedback onto the gas. For simplicity, weonly consider dust with a diameter of 0.2mm and inter-nal density of 1.25g/cm . By assuming a solar-type starand a reference radius R = 20AU, the 0.2mm grainshave a Stokes number of ∼ − at R .Figure 11 shows that the dust largely resembles thegas pattern. The (cid:104) Σ (cid:105) plots demonstrate that the pri-mary dust ring is formed at the primary gas ring ( r ∼ . R ). The two comparable dust rings coincide withthe secondary and tertiary gas ring ( r ∼ . , . R ). Weexpect the dust rings to be wider in viscous disks. Thisresult suggests that our substructure formation mech-anism may create significant observational features inthe dust continuum. Further work including a range ofdust sizes and dust-feedback is required to make morerealistic observational predictions.Kinematical signatures of disk eccentricity may alsobe observable. During the spiral instability stage in ourmodel, the eccentric motions of the gas reach ∼
10% ofthe Keplerian motion, and in the final stage, the ringsmay still possess some residual eccentricity. Gas kine-matics traced by any molecular line emission can verifywhether or not the mechanism we present in this pa-per is the true origin of some observed disk substruc-tures. A caveat here is that we have considered only two-dimensional disks. As we discussed in Section 4.3,while the starlight is supposedly scattered at the topand bottom surfaces of the disk, it is not clear if thedisk eccentricity should be independent of height.Observations of a spiral instability in disks may alsoprovide new constraints to the mass, opacity, and otherbasic disk parameters because the instability depends onthe disk properties such as the self-gravity-to-pressureratio and cooling efficiency. SUMMARYIn this paper, we have explored the properties andlong-term evolution of a spiral instability in 2D self-gravitating (but gravitationally stable) disks with cool-ing. Our main results can be summarized as follows:1. With our
Fiducial long-term hydrodynamicssimulation, we demonstrated that there is atrapped one-armed spiral instability that evolvesinto a set of axisymmetric rings over longtimescales when the disk cools rapidly.2. We presented a linear equation that governs theslow global eccentric modes. We showed that thedisk must have large enough g max and must coolfaster than β c ∼ g − for the instability to occur.The eigensolutions we obtain from our linear the-ory can accurately describe our non-linear hydro-dynamics simulations in terms of the mode mor-phology, growth rate, and initial ring locations.3. We showed that the ring-forming mass flux relieson a persistent angular momentum transfer be-tween the axisymmetric background and the non-axisymmetric wave. Furthermore, we presenteda first-principle method for estimating the loca-tions and amplitudes of the rings using only theinitial equilibrium of the disk and the eccentricmode from our linear theory.4. We showed that the s-mode instability does notrely on a specific model of density or temperatureprofile. As long as the disk profile of g exhibits s-mode trapping regions, the spiral and subsequentring forming process can be expected to be generic.To conclude, our results demonstrate a new mecha-nism to form disk substructures without planets . Themorphology and kinematics of these substructures canbe quantitatively predicted by our analytical theory.Our theory may also provide new interpretation to theprotoplanetary disk observations. ings formed by an eccentric instability Software:
LA-COMPASS (Li et al. 2005, 2009a),Matplotlib (Hunter 2007), NumPy (van der Walt et al.2011), SciPy (Virtanen et al. 2020) APPENDIX A. ECCENTRIC NORMAL MODE EQUATIONWe start with re-writing the energy equation (Equation 1) in terms of the gas entropy, S ≡ ln( P/ Σ γ ),Σ T ( ∂S∂t + u ∂S∂r + vr ∂S∂φ ) = − Σ T − T eq t c , (A1)Now we assume that the background state is steady and axisymmetric, while the perturbation is small and has m = 1 symmetry. After removing the equilibrium background from Equation (A1) and linearizing, we get − i ( ω − Ω) (cid:104) P (cid:48) P − γ Σ (cid:48) Σ (cid:105) + u (cid:48) (cid:104) P dPdr − γ d Σ dr (cid:105) = Σ (cid:48) T − P (cid:48) t c P , (A2)where ω is the mode frequency and we have approximated T eq = T . Unprimed quantities denote the axisymmetricbackground, while primed quantities are the m = 1 perturbations.Using the slow-mode approximation Ω = Ω K = (cid:112) GM (cid:63) /r (cid:29) | ω | and defining β ≡ Ω K t c , the pressure perturbation P (cid:48) can be expressed as P (cid:48) = iβiβ + 1 (cid:18) − γrP dEdr − r dPdr E (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) adiabatic pressure + 1 iβ + 1 (cid:18) − rP dEdr − rP Σ d Σ dr E (cid:19)(cid:124) (cid:123)(cid:122) (cid:125) isothermal pressure , (A3)where we have expressed the perturbations using the complex eccentricity, E , defined as (Lee et al. 2019a) u (cid:48) = i Ω rE, v (cid:48) = −
12 Ω rE and Σ (cid:48) = − r ddr (Σ E ) . (A4)Note that P (cid:48) is a linear combination of the pressure perturbations in the adiabatic limit and the locally isothermallimit.4 Li et al.
Following Lee et al. (2019a), using P (cid:48) and the linearized equations of motion for u (cid:48) , v (cid:48) , and Σ (cid:48) , we arrive at thenormal mode eccentricity equation with cooling,2 r Ω K Σ ωE = [ iβiβ + 1 M adi + 1 iβ + 1 M iso + M dsg + M β ] E. (A5)where the four operators are, M adi E = ddr (cid:18) γr P dEdr (cid:19) + r dPdr E, (A6) M iso E = ddr (cid:18) r P dEdr (cid:19) + r dPdr E − ddr (cid:18) Σ dc dr r E (cid:19) , (A7) M dsg E = − Σ r ddr (cid:18) r d Φ dr (cid:19) E − Σ ddr ( r Φ ) , (A8) M β E = P r (cid:20) ddr (cid:18) β + iβ β (cid:19) (cid:18) γr dEdr + rP dPdr E (cid:19) + ddr (cid:18) − iβ β (cid:19) (cid:18) r dEdr + r Σ d Σ dr E (cid:19)(cid:21) . (A9)These quantify the adiabatic, isothermal, self-gravity, and non-constant β effects. Expressions for Φ and Φ can befound in Appendix A.3 of Lee et al. (2019a). Setting M β = 0 and simplifying results in Equation (4). B. WKB DISPERSION RELATIONFollowing Lee et al. (2019a,b), the dispersion relation governing an eccentric wave can be derived using a WKB-likeapproximation. By letting y = ( r P ) / E , one can transform the pressure operators into M adi E = γ ( r P ) / d dr y + (cid:20)(cid:114) rP dPdr − γ d dr ( r P ) / (cid:21) y, (B10) M iso E = ( r P c ) / d dr (cid:18) yc iso (cid:19) + (cid:20)(cid:114) rP dPdr − c iso d dr ( r P c ) / (cid:21) y. (B11)For the self-gravity operator, we use the same second-order expression as Lee et al. (2019a), and set M β = 0.Transforming E → y removes all first derivative terms from the pressure operators. Taking d /dr → − k , we have ω = (cid:18) γβ β (cid:19) (cid:18) − k c (cid:19) + ω p (cid:124) (cid:123)(cid:122) (cid:125) pressure with β cooling + πG ΣΩ | k | + ω g (cid:124) (cid:123)(cid:122) (cid:125) disk self-gravity , (B12)where ω p = c (cid:20) rP dPdr −
11 + β (cid:18) β γ ( r P ) − / d dr ( r P ) / + ( r P c ) − / d dr ( r P c ) / (cid:19)(cid:21) (B13) ω g = − r ddr (cid:18) r d Φ dr (cid:19) − (cid:18) πG ΣΩ r (cid:19) (cid:112) k r + 9 / (cid:20) dd ln r ln Σ( r P ) / − (cid:21) (B14)are the “potential” profiles with corrections of order | kr | − . We have dropped all the imaginary terms in the finalexpressions. Our pressure part is different from Lee et al. (2019a,b) because of cooling, but the self-gravity part is thesame.The DRMs (e.g., Figure 6) show the contours of ω in r - kr space. Each contour that satisfies the quantum conditionfor standing waves as in Lee et al. (2019a) represents an eccentric mode allowed in the WKB theory. Analogous to theenergy level diagram of the solutions to wave functions in a potential well, Figure 12 is a frequency level diagram ofthe first few modes from the DRM of our Fiducial disk. The two modes at the top are the WKB ground state s-modeand its first harmonic. The ground state frequency shows excellent agreement with the simulated pattern speed. Therest are e-p/pg modes, which are characterized by their encounters with the Lindblad resonances. ings formed by an eccentric instability Figure 12.
Frequency level diagram of the
Fiducial disk. The black curves represent the first few standing modes fromour WKB theory (fundamental mode at the top). Overplotted are the Lindblad resonances ( ω | k =0 ; blue curve) and Q-barriers( ω | kr = g ; orange curve). Intersection points between these curves and a given mode’s ω mark the radial locations of that mode’sLindblad resonances or Q-barriers.Modes can only propagate in the shaded regions and are exponentially decaying outside theseregions. The pattern speed from the simulation is shown in red to compare with the fundamental frequency.C. GROWTH RATEThe equation governing the modes angular momentum deficit can be obtained, to lowest order in β and dβ/dr , bymultiplying Equation (4) by E ∗ and integrating over the disk, (cid:90) rLω | E | dr = (cid:90) (cid:18) E ∗ β ddr (cid:26) r P (cid:20) c ddr (cid:18) Ec (cid:19) + γβ dEdr (cid:21) + iβr P (cid:20) ( γ − dEdr + d ln c dr E (cid:21)(cid:27) + r dPdr | E | (cid:124) (cid:123)(cid:122) (cid:125) pressure with β cooling − disk self-gravity (cid:122) (cid:125)(cid:124) (cid:123) Σ r ddr (cid:0) r d Φ dr (cid:1) | E | − E ∗ Σ ddr (cid:0) r Φ (cid:1) + non-constant β (cid:122) (cid:125)(cid:124) (cid:123) i dβdr r P (cid:20) ( γ − rE ∗ dEdr + d ln c d ln r | E | (cid:21) (cid:19) dr, (C15)where L = r Ω K Σ.We first ignore the ‘non-constant β ’ contribution to derive the growth rate in Section 3. In the adiabatic limit, ithas been shown that our form of self-gravity only contributes to the real part of ω , i.e. the mode pattern speed (Leeet al. 2019a). Hence, the imaginary component of this equation is γ (cid:90) rL | E | dr = (cid:90) (cid:61) (cid:18) E ∗ β ddr (cid:26) r P (cid:20) c ddr (cid:18) Ec (cid:19) + γβ dEdr (cid:21) + iβr P (cid:20) ( γ − dEdr + d ln c dr E (cid:21)(cid:27) (cid:19) dr. (C16)The right-hand side can be further simplified through integration by parts and assuming that the boundary contributionis zero: γ (cid:90) rL | E | dr = −
11 + β (cid:90) (cid:61) (cid:18) dE ∗ dr (cid:26) r P c ddr (cid:18) Ec (cid:19) + iβr P (cid:20) ( γ − dEdr + d ln c dr E (cid:21)(cid:27) (cid:19) dr = −
11 + β (cid:90) r P (cid:20) β ( γ − | dEdr | − (cid:61) ((1 − iβ ) d ln c dr E dE ∗ dr ) (cid:21) dr, (C17)as Equation (6) shows.The correction due to the ‘non-constant β ’ part in Eqution (C15) is (cid:61) ( ω β ) (cid:90) rL | E | dr = (cid:90) r P dβdr (cid:18) ( γ − r (cid:60) ( E ∗ dEdr ) + d ln c d ln r | E | (cid:19) dr. (C18)Hence, with dE/dr = ikE and k = g max /r , the growth rate is γ ≈ − g max h Ω(2 γ ) [( γ − βg max + d ln c d ln r sign( k ) − rg max d ln c d ln r dβdr ] , (C19)6 Li et al. which has the criterion for instability: β < β c = 1 g ( γ − (cid:18) − d ln c d ln r (cid:19) (cid:20) sign( k ) − rg dβdr (cid:21) (C20)as mentioned in Section 4.2. D. MODE ANGULAR MOMENTUM DEPOSITIONThe torque from the mean flow to the wave is (Dempsey et al. 2020) t dep πr = − (cid:104) u (cid:48) Σ (cid:48) (cid:105) ∂r (cid:104) v (cid:105) ∂r + (cid:104) Σ (cid:105) (cid:28) u (cid:48) ∂rv (cid:48) ∂r (cid:29) − (cid:28) Σ (cid:48) Σ ∂P∂φ (cid:48) (cid:29) . (D21)The first two terms are related to the disk vortensity, while the third term is β dependent and vanishes in the β → β → + ∞ limits.To evaluate the torque in term of the disk eccentricity, we first rewrite the first two terms as t dep ,ξ πr ≡ − (cid:104) u (cid:48) Σ (cid:48) (cid:105) ∂r (cid:104) v (cid:105) ∂r + (cid:104) Σ (cid:105) (cid:28) u (cid:48) ∂rv (cid:48) ∂r (cid:29) = 2Σ (cid:60) ( u (cid:48)∗ ∂ ( rv (cid:48) ) ∂r − u (cid:48)∗ Σ (cid:48) rξ ) , (D22)where ξ = ˆ z · ∇ × v / Σ = ∂ r ( rv ) / ( r Σ) is the vortensity to the zeroth order.The disk vortensity evolution is governed by ∂ξ∂t + u ∂ξ∂r + vr ∂ξ∂φ = 1Σ ∇ Σ × ∇ P. (D23)For a mode ξ (cid:48) ∝ e i ( φ − ωt ) with | ω | (cid:28) Ω and (cid:61) ( ω ) = γ , the governing equation can be simplified to( i Ω + γ ) r Σ ξ (cid:48) + u (cid:48) r Σ ∂ξ∂r = iT f (cid:48) , (D24)where r Σ ξ (cid:48) = ∂(cid:96) (cid:48) ∂r − iu (cid:48) − Σ (cid:48) rξ and f (cid:48) = µ Σ T (cid:48) T − η Σ (cid:48) , (D25)and where η = d ln T /dr and µ = d ln Σ /dr for T ≡ c . Multiplying Equation (D24) by ( − i Ω + γ ) u (cid:48)∗ , we get(Ω + γ ) r Σ u (cid:48)∗ ξ (cid:48) + ( − i Ω + γ ) | u (cid:48) | r Σ ∂ξ∂r = (Ω + iγ ) T u (cid:48)∗ f (cid:48) . (D26)The real part of the above equation says(Ω + γ ) t dep ,ξ πr + γ r Σ ∂ξ∂r | u (cid:48) | = Ω T (cid:60) ( u (cid:48)∗ f (cid:48) ) − γ T (cid:61) ( u (cid:48)∗ f (cid:48) ) . (D27)Hence, a non-zero t dep ,ξ can be caused by the growth/suppression of an eccentric mode (second LHS term) and thebaroclinic effect (the RHS terms).To proceed, we replace all primed quantities as u (cid:48) = i Ω rE, Σ (cid:48) = − ( Ω + i Ω γ Ω + γ ) r ∂∂r (Σ E ) , and T (cid:48) T = iβ ( γ −
1) Σ (cid:48) Σ − β Ω ∂S∂r u (cid:48) . (D28)The velocity u (cid:48) is from the eccentric formalism. Σ (cid:48) here is different from Equation (A4) because we need to keepthe correction of order γ . The expression for T (cid:48) is from Equation (A1) in the limit of β (cid:28)