Instabilities at planetary gap edges in 3D self-gravitating disks

Numerical simulations are presented to study the stability of gaps opened by giant planets in 3D self-gravitating disks. In weakly self-gravitating disks, a few vortices develop at the gap edge and merge on orbital time-scales. The result is one large but weak vortex with Rossby number -0.01. In moderately self-gravitating disks, more vortices develop and their merging is resisted on dynamical time-scales. Self-gravity can sustain multi-vortex configurations, with Rossby number -0.2 to -0.1, over a time-scale of order 100 orbits. Self-gravity also enhances the vortex vertical density stratification, even in disks with initial Toomre parameter of order 10. However, vortex formation is suppressed in strongly self-gravitating disks and replaced by a global spiral instability associated with the gap edge which develops during gap formation.


Introduction
Gaps induced by planets in protoplanetary disks can become dynamically unstable if the disk viscosity is sufficiently small [1]. This is because planetary gap edges are associated with potential vorticity or vortensity extrema [2,3], the existence of which is necessary for instability [4]. Gap edges may undergo vortex formation in weakly self-gravitating disks (associated with vortensity minima) or a spiral instability (associated with vortensity maxima) in strongly self-gravitating yet Toomre-stable disks [5,6]. Development of such instabilities can significantly affect planetary migration [7] and dust evolution [8]. These studies have employed 2D disk models, but gap edges have characteristic widths of the disk local scale-height, so it is necessary to extend the study of gap stability to 3D.

Numerical simulations with ZEUS-MP
The system is an inviscid, non-magnetized 3D fluid disk embedded with a giant planet of mass M p , both rotating about a central star of mass M * . Spherical co-ordinates (r, θ, φ) centered about the star are adopted. Units are such that G = M * = 1 , where G is the gravitational constant.
The disk is governing by the Euler equations coupled with self-gravity through the Poisson equation. The equation of state is locally isothermal, so the sound-speed is c s = HΩ k , where H = hR is the isothermal scale-height with constant aspect-ratio h, Ω 2 k ≡ GM * /R 3 and R = r sin θ. Each disk model is labelled by its minimum Keplerian Toomre parameter Q 0 , located at the outer disk boundary. The planet is regarded as an external potential and held on a circular Keplerian orbit of radius r p at the a e-mail: mklin924@cita.utoronto.ca The Journal's name

Weakly self-gravitating disks (Q 0 = 8)
Two simulations with Q 0 = 8 were run, one with self-gravity and the other without. M p = 0.002M * and h = 0.07 are adopted. Both cases developed 2-3 vortices early on but the quasi-steady state is a single vortex with Ro ∼ −0.01, where the Rossby number is defined as Ro ≡ ω z / 2Ω , where ω z is the absolute vertical vorticity and Ω is the azimuthally-averaged angular speed.
The vortices differ noticeably in the (r, θ) plane. This is shown in Fig. 1. Self-gravity enhances the vertical density stratification of the vortex, with the midplane density enhancement being ≃ 50% larger than that in the non-self-gravitating run. The initial Keplerian Toomre parameter at the radius of vortex formation is about 10, but even this is sufficient to affect the vortex vertical structure. The creation of vortensity minima lowers the Toomre parameter, which further decreases with vortex formation since they are over-densities. Thus, self-gravity can become important in the perturbed state even if it is negligible initially. Fig. 2 shows the relative density perturbation and Rossby number at the end of the a simulation with Q 0 = 3. Self-gravity is included. The 5-vortex configuration is sustained from its initial development, unlike in the weakly self-gravitating case where merging occurred over the same time-scale. The preference for linear vortex modes with higher azimuthal wavenumber m with increasing strength of self-gravity was observed in 2D simulations [5,11], and persists in 3D. The smaller vortices here are stronger than the single vortex in previous case, with Rossby number Ro ∼ −0.2 and the relative density perturbation has significant vertical dependence.

Long term simulation
A smaller disk model, with r ∈ [2, 20], was simulated up to t ∼ 500P 0 . M p = 10 −3 M * and h = 0.05 were adopted for this run. Fig. 3 shows the relative density perturbation towards the end of the simulation. The multi-vortex configuration lasted ∼ 200 orbits at the vortex radius. Notice a vortex may reach comparable over-densities to the final post-merger vortex in the weakly self-gravitating disk. It was observed that |Ro| decreased from 0.2 at the onset of vortex formation, to 0.1 towards the end of the simulation, which may be due to limited numerical resolution.

Gap edge spiral instability (Q
The linear vortex instability can be suppressed by strong self-gravity. To demonstrate this, a disk model with Q 0 = 1.5, h = 0.05 and M p = 10 −3 M * was simulated. Fig. 4 shows the development of an m = 2 spiral mode associated with the outer gap edge. This instability occurs during gap formation and supplies positive disk-on-planet co-orbital torques, because the over-density protrudes the outer gap edge and approaches the planet from upstream. The disturbance is significantly stratified, with most of the perturbation confined near the midplane. The global spiral pattern appears transient, having decreased in amplitude by t = 50P 0 , but this is likely a radial boundary condition effect.

Discussion
Direct numerical simulations of 3D self-gravitating disk-planet systems confirm the stability properties of gap edges previously explored in 2D [5,6,11,12]. Vertical self-gravity enhances the density stratification of a vortex. Given the vortex instability is only expected to occur in low viscosity regions of protoplanetary disks -dead zones -which are overlaid by actively accreting layers [13], it may be advantageous to have self-gravity confining the over-density near the midplane, thereby mitigate upper disk boundary effects and make the instability a more robust mechanism for vortex formation.