# Vortex survival in 3D self-gravitating accretion discs

Min-Kai Lin, A. Pierens

► To cite this version:

Min-Kai Lin, A. Pierens. Vortex survival in 3D self-gravitating accretion discs. Monthly Notices of the Royal Astronomical Society, 2018, 478 (1), pp.575-591. 10.1093/mnras/sty947. hal-01777647

HAL Id: hal-01777647

https://hal.science/hal-01777647v1

Submitted on 22 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Vortex survival in 3D self-gravitating accretion discs

Min-Kai Lin<fnref n="1" />* and Arnaud Pierens<fnref n="2" />,<fnref n="3" />

Accepted 2018 April 12. Received 2018 April 6; in original form 2018 February 22

## ABSTRACT

Large-scale, dust-trapping vortices may account for observations of asymmetric protoplanetary discs. Disc vortices are also potential sites for accelerated planetesimal formation by concentrating dust grains. However, in 3D discs vortices are subject to destructive ‘elliptic instabilities’, which reduces their viability as dust traps. The survival of vortices in 3D accretion discs is thus an important issue to address. In this work, we perform shearing box simulations to show that disc self-gravity enhances the survival of 3D vortices, even when self-gravity is weak in the classic sense (e.g. with a Toomre $Q \simeq 5$). We find a 3D self-gravitating vortex can grow on secular time-scales in spite of the elliptic instability. The vortex aspect ratio decreases as it strengthens, which feeds the elliptic instability. The result is a 3D vortex with a turbulent core that persists for $\sim 10^3$ orbits. We find when gravitational and hydrodynamic stresses become comparable, the vortex may undergo episodic bursts, which we interpret as an interaction between elliptic and gravitational instabilities. We estimate the distribution of dust particles in self-gravitating, turbulent vortices. Our results suggest large-scale vortices in protoplanetary discs are more easily observed at large radii.

**Key words:** accretion – hydrodynamics – instabilities – methods: numerical – protoplanetary discs – accretion discs.

## 1 INTRODUCTION

Vortex dynamics is playing an increasingly important role in protoplanetary discs from both observational and theoretical perspectives. Several observations of transition discs reveal lopsided asymmetries in dust (Casassus et al. 2013; Fukagawa et al. 2013; Isella et al. 2013; van der Marel et al. 2013, 2016; Pérez et al. 2014; Hashimoto et al. 2015; Marino et al. 2015; Ohta et al. 2016; Kraus et al. 2017). Disc vortices provide a possible explanation because they are associated with localized pressure bumps, which can act as dust traps (Barge & Sommeria 1995; Lyra & Lin 2013; Zhu et al. 2014). This interpretation suggests that large-scale vortices may be common in real protoplanetary discs.

The origin of protoplanetary disc vortices is usually attributed to hydrodynamic instabilities. To date, these include: the Rossby Wave Instability (RWI; Lovelace et al. 1999; Li et al. 2000; Meheut et al. 2010; Lin 2012b; Yellin-Bergovoy, Heifetz & Umurhan 2016; Ono et al. 2016); the Sub-critical Baroclinic Instability (SBI; Petersen, Julien & Stewart 2007; Lesur & Papaloizou 2010; Lyra & Klahr 2011; Raettig, Lyra & Klahr 2013; Barge, Richard & Le Dizès 2016); the Convective Overstability (ConO; Klahr & Hub-

bard 2014; Lyra 2014; Latter 2016); the Vertical Shear Instability (VSI; Nelson, Gressel & Umurhan 2013; Stoll & Kley 2014; Barker & Latter 2015; Lin & Youdin 2015; Richard, Nelson & Umurhan 2016); and the Zombie Vortex Instability (ZVI; Marcus et al. 2015; Lesur & Latter 2016; Umurhan, Shariff & Cuzzi 2016b). A common outcome of these instabilities is vortex formation and/or amplification.

Which of these hydrodynamic instabilities operate depends on the disc structure and thermal state, but taken together they pertain to a wide range of protoplanetary disc conditions (Umurhan, Estrada & Cuzzi 2016a). Since these mechanisms form or amplify vortices, we can expect vortices to be natural features of protoplanetary discs. These vortices may assist in planetesimal formation by concentrating dust particles (Inaba & Barge 2006; Meheut et al. 2012).

Given the relevance of protoplanetary disc vortices to planet formation, it is essential to understand how they evolve under realistic conditions. A well-known result is that in three-dimensional (3D) discs vortices can develop secondary instabilities that weaken or even destroy them (Lithwick 2009; Lesur & Papaloizou 2009; Chang & Oishi 2010; Railton & Papaloizou 2014). This would strongly threaten their survival in protoplanetary discs. It is therefore necessary to model disc vortices in 3D and account for such ‘elliptic instabilities’ to obtain realistic vortex lifetimes.

<footnote n="1">¹Institute of Astronomy and Astrophysics, Academia Sinica, Taipei 10617, Taiwan</footnote>

<footnote n="2">²Université de Bordeaux, Observatoire Aquitain des Sciences de l'Univers, BP89 F-33271 Floirac Cedex, France</footnote>

<footnote n="3">³CNRS, Laboratoire d'Astrophysique de Bordeaux, BP89 F-33271 Floirac Cedex, France</footnote>

<footnote>* E-mail: mklin@asiaa.sinica.edu.tw</footnote>
Another effect that is often neglected in modelling protoplanetary disc vortices is their self-gravity. Disc self-gravity is usually described by the Toomre parameter

$$Q = \frac{c_s \Omega}{\pi G \Sigma}, \qquad (1)$$

where $c_s$ is the sound speed, $\Omega$ is the disc's rotation frequency, $\Sigma$ is the surface density, and $G$ is the gravitational constant. Traditionally, self-gravity is considered important when $Q \lesssim 1$, which is the formal criterion for gravitational instability in laminar, axisymmetric thin discs (Toomre 1964). This typically translates to rather large disc-to-star mass ratios, $M_d/M_* \gtrsim 0.1$. However, the condition for gravitational instability, and hence the importance of self-gravity, can be relaxed when additional physics are considered (Lin & Kratter 2016).

The effect of self-gravity on the equilibrium structure of a disc vortex was first investigated by Adams & Watkins (1995) assuming circular vortices. However, the shear in protoplanetary discs likely only permits elliptical vortices (Lesur & Papaloizou 2009). Self-gravity has also been included in studies of vortex formation and evolution via the RWI (Lyra et al. 2009b; Lin & Papaloizou 2011; Yellin-Bergovoy et al. 2016; Regály & Vorobyov 2017) and in gravito-turbulent discs (Mamatsashvili & Rice 2009).

In the case of the RWI, recent studies show that disc self-gravity affects vortex evolution when $Q \lesssim \pi/2h$ (Lovelace & Hohlfeld 2013; Zhu & Baruteau 2016), where $h$ is the disc aspect ratio. For typical protoplanetary discs with $h \sim 0.05$, this means that self-gravity should be accounted when $Q \lesssim 20$; significantly larger than the classic condition. However, these calculations have been limited to razor-thin, two-dimensional (2D) discs.

The above discussion motivates us to study the evolution of elliptic vortices in 3D self-gravitating protoplanetary discs. This was first attempted by Lin (2012a) in global 3D numerical simulations of vortex formation via the RWI at planet gaps in self-gravitating discs. However, these simulations lacked the numerical resolution needed to capture the elliptic instability (EI). Thus, the survival of 3D vortices in self-gravitating discs remains unclear.

In this work, we study the evolution of 3D vortices in a local patch of a self-gravitating protoplanetary disc. We use the shearing box framework and insert a vortex to evolve as an initial value problem. While less realistic than Lin (2012a) in terms of vortex formation, this approach permits clean numerical experiments and higher resolution.

We find self-gravity enhances the survival of 3D vortices against the EI. The EI develops regardless of the strength of self-gravity, but with sufficient self-gravity the vortex subsequently undergoes a slow growth, with increasing levels of internal hydrodynamic turbulence. The vortex then undergoes bursts; alternating between coherent and turbulent cores. We interpret this as ‘gravito-elliptic’ feedback: a growing vortex spins up and develops hydrodynamic turbulence because of the EI, which contributes to further collapse by removing the large-scale internal rotation. The latter effect is similar to secular gravitational instabilities in viscous or dusty accretion discs (e.g. Youdin 2011; Takahashi & Inutsuka 2014; Lin & Kratter 2016).

This paper is organized as follows. We list the basic equations in Section 2. In Section 3 we give a qualitative discussion on the expected effect of self-gravity on disc vortices. We describe our simulation set-up and diagnostic measures in Section 4. We present simulation results in Section 5. We discuss implications of our simulations in Section 6 before summarizing in Section 7.

## 2 LOCAL DISC MODEL

We consider an inviscid, 3D self-gravitating accretion disc orbiting a central star. Cylindrical coordinates $(R, \phi, z)$ are centred on the star. The equilibrium rotation profile is described by $\Omega(R) \propto R^{-q}$. We are interested in Keplerian discs with $q = 3/2$ and a vertical oscillation frequency $\Omega_z = \Omega$. We adopt the shearing box framework to study a small patch of the disc (Goldreich & Lynden-Bell 1965). Local Cartesian coordinates $(x, y, z)$ are anchored on a fiducial point $(R_0, \phi_0, 0)$ in the global disc corotating with the background flow, so $\phi_0 = \Omega_0 t$ where $\Omega_0 \equiv \Omega(R_0)$. For simplicity, we drop the subscript zero hereafter.

The Cartesian shearing box fluid equations are

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0, \qquad (2)$$

$$\frac{\partial \mathbf{v}}{\partial t} + \mathbf{v} \cdot \nabla \mathbf{v} = -\frac{1}{\rho} \nabla p - \nabla \Phi - 2\Omega \hat{z} \times \mathbf{v} + \Omega^2(2qx, 0, -z), \qquad (3)$$

$$\nabla^2 \Phi = 4\pi G \rho, \qquad (4)$$

where $\rho$ is the mass density, $\mathbf{v}$ is the total fluid velocity in the box, $p$ is the pressure, and $\Phi$ is the gravitational potential of the gas. We adopt an isothermal equation of state so that

$$p = c_s^2 \rho, \qquad (5)$$

where $c_s = H\Omega$ is the constant isothermal sound speed, with $H$ being the disc scale height at the reference radius.

A simple, steady-state equilibrium solution to the above equations describes a stratified, axisymmetric shear flow with $\mathbf{v} = -q\Omega x \hat{y}$ and $\rho = \rho(z)$. In dimensionless variables $\hat{\rho} \equiv \rho/\rho_0$, where $\rho_0$ is the mid-plane density, and $Z \equiv z/H$, the density field satisfies

$$\frac{d^2 \ln \hat{\rho}}{dZ^2} + \frac{\hat{\rho}}{Q_{3D}} + 1 = 0, \qquad (6)$$

where

$$Q_{3D} \equiv \frac{\Omega^2}{4\pi G \rho_0} \qquad (7)$$

is a self-gravity parameter for 3D discs (Mamatsashvili & Rice 2010). The boundary conditions are $\hat{\rho}(0) = 1$ and $d\hat{\rho}/dz = 0$ at $z = 0$. In the non-self-gravitating limit, $Q_{3D} \to \infty$, we obtain the usual Gaussian solution, $\hat{\rho} = \exp(-Z^2/2)$.

However, for non-self-gravitating discs there exists another class of equilibrium flow solutions corresponding to a vortex, which we discuss below.

## 3 VORTEX SOLUTIONS AND THE QUALITATIVE EFFECT OF SELF-GRAVITY

In this section we discuss the effect of disc self-gravity on equilibrium flows that correspond to a vortex. To do this, it is useful to work in coordinates that reflect such a flow pattern. We thus adopt the non-orthogonal ‘elliptico-polar’ coordinates $(s, \varphi, z)$ defined by

$$x = s \cos \varphi, \qquad (8)$$

$$y = \chi s \sin \varphi, \qquad (9)$$

where $\chi$ is a constant, and $z$ remains unchanged (see e.g. Kerswell 1994, for further discussion of this coordinate system). Here, we take $\chi > 1$ so that the coordinate $s \in [0, \infty)$ corresponds to the
semiminor axis of a series of similar ellipses, and $\varphi \in [0, 2\pi]$ corresponds to the azimuthal position along a given ellipse. We also define the fluid velocity components $v_s, v_\varphi$ such that

$$v_x \hat{x} + v_y \hat{y} = v_s s + v_\varphi \varphi, \quad (10)$$

where

$$s = \cos \varphi \hat{x} + \chi \sin \varphi \hat{y}, \quad (11)$$

$$\varphi = -\sin \varphi \hat{x} + \chi \cos \varphi \hat{y}. \quad (12)$$

The shearing box equations in elliptico-polar coordinates are listed in Appendix A.

## 3.1 The GNG vortex in the absence of self-gravity

Goodman, Narayan & Goldreich (1987) present exact vortex solutions to equations (2) and (3) in the absence of self-gravity (by setting $G = \Phi = 0$). Here, we review this ‘GNG’ vortex in elliptico-polar coordinates.

The GNG vortex corresponds to a steady axisymmetric flow in elliptico-polar coordinates in vertical hydrostatic equilibrium, so that $\partial_t = \partial_\varphi = v_s = v_z = 0$ and $v_\varphi = v_\varphi(s)$. The momentum equations then read

$$0 = \frac{v_{\varphi}^2}{s} - \frac{1}{2} (\xi_{-} \cos 2\varphi + \xi_{+}) \frac{\partial \eta}{\partial s} + \Omega v_{\varphi} \chi (\xi_{-} \cos 2\varphi + \xi_{+}) + q \Omega^2 s (\cos 2\varphi + 1), \quad (13)$$

$$0 = \frac{\xi_{-}}{2} \sin 2\varphi \frac{\partial \eta}{\partial s} - \Omega \chi \xi_{-} v_{\varphi} \sin 2\varphi - q \Omega^2 s \sin 2\varphi, \quad (14)$$

$$0 = z\Omega^2 + \frac{\partial \eta}{\partial z}, \quad (14)$$

where $\xi_{\pm} = 1 \pm 1/\chi^2$, and

$$\eta = c_s^2 \ln \rho \quad (15)$$

is the gas enthalpy.

Since $\partial_\varphi = 0$ by construction, the trigonometric coefficients in equations (13) and (14) must vanish. Then we obtain the flow solution

$$v_{\varphi}(s) = -\Omega_{v}s \text{ with } \Omega_{v} = \sqrt{\frac{2q}{\chi^2 - 1}} \Omega. \quad (16)$$

We have chosen the sign for $\Omega_v$ so that the flow is anticyclonic (Bodo et al. 2007). The corresponding density distribution is given via

$$\ln \frac{\rho}{\rho_c} = -\frac{s^2}{2H_{\text{eff}}^2} - \frac{z^2}{2H^2}, \quad (17)$$

where $\rho_c$ is the density at $s = z = 0$ and the characteristic horizontal size $H_{\text{eff}}$ is

$$H_{\text{eff}} \equiv \frac{H}{F}, \quad (18)$$

with

$$F^2(\chi; q) \equiv \frac{2}{\xi_{-}} \left[ \frac{\sqrt{2q(\chi^2 - 1)}}{\chi} - q \right]. \quad (19)$$

Note that, in a Keplerian disc, the requirement $F^2 > 0$ implies non-self-gravitating vortices only exist with aspect ratios $\chi > 2$.

<figure><img src="image_3.png" /><figcaption>Figure 1. Characteristic horizontal vortex size, given by equation (23), as a function of the aspect ratio, for a range of self-gravity parameters $Q_{3D}$.</figcaption></figure>

## 3.2 Effect of self-gravity on the GNG vortex density field

The GNG solution is not compatible with self-gravity because the trigonometric coefficients of the Poisson equation in elliptico-polar coordinates do not vanish (see equation A6). This is because the gravitational potential of an ellipsoidal density distribution does not have the same symmetry, except for $s \to 0$ (Chandrasekhar 1969).

In order to get a sense of the effect of self-gravity, let us consider the region about $s = 0$ (vortex centroid), and simplify further by neglecting vertical stratification, so $\partial_z = z = 0$. For this cylindrical problem, we assume the self-gravitational potential has the form

$$\Phi = \frac{1}{2}As^2 + \text{const.}, \quad (20)$$

for small $s$, where $A$ is a constant. Equation (20) is consistent with the Poisson equation near the origin provided that

$$A = \frac{\Omega^2}{Q_{3D}\xi_{+}}. \quad (21)$$

Including self-gravity in the momentum equations amounts to the replacement $\eta \to \eta + \Phi$. Then

$$\begin{aligned} \frac{\partial \ln \rho}{\partial s} &= -\frac{s}{H_{\text{eff}}^2} - \frac{1}{c_s^2} \frac{\partial \Phi}{\partial s} = -\frac{s}{H_{\text{eff}}^2} \left( 1 + \frac{1}{Q_{3D}\xi_{+}F^2} \right) \\ &= -\frac{s}{H_{\text{eff},sg}^2}, \end{aligned} \quad (22)$$

where we have inserted the above approximation for $\Phi$. Writing $H_{\text{eff,sg}} = H/F_{\text{sg}}$, we have

$$F_{\text{sg}}^2(\chi, Q_{3D}; q) \equiv \frac{H^2}{H_{\text{eff,sg}}^2} = F^2 + \frac{1}{Q_{3D}\xi_{+}}. \quad (23)$$

Thus, self-gravity reduces the characteristic vortex size by increasing $F^2$ to $F^2 + 1/Q_{3D}\xi_{+}$.

In Fig. 1 we plot the characteristic vortex size $H_{\text{eff,sg}}$ as a function of the aspect ratio and self-gravity. The characteristic vortex size is insensitive to its aspect ratio for $\chi \gtrsim 4$. At a fixed $\chi$ the vortex size is reduced by self-gravity, but this effect is small for elongated vortices. Notice the minimum allowed aspect ratio (defined by $F_{\text{sg}}^2(\chi) = 0$) is reduced by self-gravity. We show this in Fig. 2: as $Q_{3D} \to 0$ the minimum aspect ratio tends to unity.

Our estimate above differs from Adams & Watkins (1995), who found that self-gravity increases the size of a circular vortex. By contrast, Figs 1 and 2 indicate self-gravity reduces the size of elliptic vortices; and that circular vortices are only allowed in strongly self-gravitating discs.
<figure><img src="image_4.png" /><figcaption>Figure 2. Minimum allowed aspect ratio of the self-gravity-modified GNG vortex as a function of the self-gravity parameter $Q_{3D}$.</figcaption></figure>

### 3.3 Effect of self-gravity on the GNG vortex vorticity field

In Appendix B, we give a more careful discussion of the effect of self-gravity. There, we also consider how self-gravity modifies the distribution of vertical vorticity of the GNG vortex. This effect is not obvious since potentials do not directly source vorticity.

We find the vorticity perturbation due to self-gravity is

$$ \Delta\omega \propto (4\Omega_v - \Omega\chi\xi_+) \left(\frac{s}{H_{\text{eff}}}\right)^4 \quad (24) $$

near the vortex centre (see equation B37). For the GNG solution this implies that $\Delta\omega < 0$ for $\chi \gtrsim 2.5$. As we will be considering elongated vortices in practice, we expect that self-gravity enhances the GNG vortex by making its vorticity more negative. However, this is a small effect as the above estimate is restricted to $s \ll H_{\text{eff}}$. Adams & Watkins (1995) reached a similar conclusion for circular vortices, where they find self-gravity does not modify the vortex velocity field.

## 4 NUMERICAL SIMULATIONS

We simulate the evolution of a single vortex in the 3D, self-gravitating shearing box using the ATHENA code (Stone et al. 2008) to evolve equations (2)–(4). We adopt standard numerical configurations: third-order reconstruction, a Roe solver, and a Corner Transport Upwind (CTU) integrator. Disc self-gravity is solved via Fast Fourier Transform, but assuming a vacuum beyond the vertical boundaries. We also use orbital advection (Stone & Gardiner 2010).

### 4.1 Grid set-up and boundary conditions

The simulation domain is $x \in [-L_x, L_x]/2$, uniformly spaced with $N_x$ grid points; and similarly for the y and z directions. Our fiducial box size and resolution are $(L_x, L_y, L_z) = (16H, 32H, 6H)$ with $(N_x, N_y, N_z) = (512, 512, 192)$, or 32 cells per $H$ in the x, z directions and 16 cells per $H$ in y. We use a lower resolution in the azimuth since we consider elongated vortices. The box size and numerical resolution are a compromise between having a sufficiently large box to minimize the influence of neighbouring vortices due to periodic boundary conditions (see below), and having enough resolution to capture the EI. We find this set-up makes it feasible to simulate up to 2000 orbits. For selected cases we also double the x resolution.

We apply periodic boundary conditions in the y direction and shear-periodic boundary conditions in x. In the z direction we apply reflective boundary conditions for the momenta, and extrapolate the density field assuming vertical hydrostatic equilibrium. In addition,

we apply damping for $|x| \in [0.8L_x, L_x]$ in which the hydrodynamic variables are relaxed towards their initial values on a time-scale of 0.01 orbits. This is done to minimize spiral density waves launched by neighbouring vortices.

### 4.2 Units and self-gravity parameter

We adopt computational units such that $H = c_s = 1$, then $\Omega = 1$ as well. Time is quoted in units of the orbital period at the fiducial radius, $P_0 = 2\pi/\Omega$. The initial mid-plane density is $\rho_0 = 1$. In ATHENA, the strength of self-gravity is parametrized by the value of $4\pi G$. We specify it through $Q_{3D} = 1/4\pi G$ in computational units.

For the disc models considered in this work, the initial density distribution is approximately Gaussian, $\rho \simeq \rho_0\exp(-z^2/2H^2)$. Then our self-gravity parameter $Q_{3D}$ is related to the classic Toomre parameter $Q$ by

$$ Q_{3D} \approx \sqrt{\frac{\pi}{8}} Q. \quad (25) $$

Thus, using $Q$ or $Q_{3D}$ is immaterial as they only differ by a factor of the order of unity. We will label our simulations with $Q_{3D}$, however, since this is the natural choice for 3D discs.

### 4.3 Initial condition and vortex perturbation

The GNG vortex described above is useful for analytic discussion due to its simple form, but is inconvenient to set as initial conditions in shearing box simulations because it is unbounded in the horizontal plane, which does not respect the standard shearing box boundary conditions. In a global disc, we expect the flow far from the vortex to return to orbital motion about the star.

We thus initialize an axisymmetric, 3D shearing box as described in Section 2 and apply the GNG vortex flow as a perturbation:

$$ \Delta v_x = \frac{\Omega_v}{\chi} y \exp(-m^2), \quad (26) $$

$$ \Delta v_y = (q\Omega x - \Omega_v \chi x) \exp(-m^2), \quad (27) $$

with

$$ m^2(x, y) = \frac{x^2}{b^2} + \frac{y^2}{a^2} = \frac{1}{b^2} \left( x^2 + \frac{y^2}{\chi^2} \right), \quad (28) $$

where (a, b) are the semimajor and semiminor axes of the vortex perturbation, respectively. Near the origin, the total velocity field is approximately that of the GNG vortex, while far from it one recovers the Keplerian shear flow.

We introduce the above velocity perturbations through source terms in the horizontal momentum equations over a time-scale of $10P_0$. We specify the perturbation radial size $b$ and aspect ratio $\chi$, and obtain the appropriate value of $\Omega_v$ from equation (16). Figs 3 and 4 show that the initial vortex, formed by perturbing the disc with the GNG solution, is insensitive to $Q_{3D}$. The self-gravitating vortex is slightly more stratified (as found by Lin 2012a).

### 4.4 Diagnostics

Our key diagnostics are:

(i) The Rossby number

$$ \text{Ro} \equiv \frac{\hat{z} \cdot \nabla \times \mathbf{v}}{2\Omega} + \frac{q}{2}, \quad (29) $$
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Vortex formed by perturbing the shearing box with the GNG vortex solution. Left: self-gravitating disc with $Q_{3D} = 3$. Right: effectively non-self-gravitating disc with $Q_{3D} = 100$. The colour scale corresponds to the Rossby number $R_o = (\nabla \times \mathbf{v})/2\Omega + q/2$ at the mid-plane. Here, $\chi$ is the aspect ratio of the locus $R_o = 0.5 \min(R_o)$ marked by the white contour.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Vertical structure of vortices formed by perturbing the shearing box with the GNG vortex solution. Upper: self-gravitating disc with $Q_{3D} = 3$. Lower: non-self-gravitating disc with $Q_{3D} = 100$. The colour scale is the density perturbation relative to $t=0$ (thus removing the background disc stratification).</figcaption></figure>

as a dimensionless measure of the vertical vorticity relative to the background shear flow. The Rossby number can be used to measure the strength of the vortex, which can be defined as a coherent, extended region of $R_o < 0$ in the absence of turbulence (see Fig. 3).

(ii) The vertical Mach number

$$M_z \equiv \frac{|v_z|}{c_s} \qquad (30)$$

to characterize vertical motions and thus the EI, which is expected to be three-dimensional (Lesur & Papaloizou 2009).

(iii) The Reynolds and gravitational stresses

$$\alpha_R \equiv \frac{1}{c_s^2} v_x (v_y + q\Omega x), \qquad (31)$$

$$\alpha_G \equiv \frac{1}{c_s^2} \frac{\partial_x \Phi \partial_y \Phi}{4\pi G \rho}, \qquad (32)$$

to characterize turbulent activity. The total stress is $\alpha \equiv \alpha_R + \alpha_G$.

(iv) Vortex characteristics. We calculate the evolution of the vortex shape in post-processing, as follows. For this we use the vortex density field because it remains smooth throughout the simulation. (The vorticity field can involve small-scale structures so it becomes difficult to define a large-scale coherent pattern.) We first locate the horizontal coordinates $(x_v, y_v)$ of $\max(\rho)$ near the origin. We define the vortex boundary as the contour around $(x_v, y_v)$ on which $\rho = 0.5\max(\rho)$. This contour is usually close to an ellipse, except at late evolution. We then define the semiminor and major vortex axes, $b_\rho$ and $a_\rho$, as half the length and width of the contour, respectively. The corresponding aspect ratio is $\chi_\rho = a_\rho/b_\rho$.

We define the density-weighted average value of a quantity $f$ as

$$\langle f \rangle \equiv \frac{\int_V \rho f dV}{M_V}, \qquad (33)$$

where

$$M_V \equiv \int_V \rho dV, \qquad (34)$$

and the sample volume $V$ is $x \in [-b, b]$.

# 5 RESULTS

Table 1 summarizes the simulations we have carried out. Our fiducial set of runs, which are highlighted in bold and discussed in detail below, consists of Q3, Q4, Q5, and Q100. The initial vortex aspect ratio is $\chi = 10$ with a radial size $b = H/2$. Note that axisymmetric instability in 3D shearing boxes requires $Q_{3D} \lesssim 0.2$ (Mamatsashvili & Rice 2010). Thus, all cases are gravitationally stable in the absence of the vortex. (They correspond to Toomre parameters $Q \gtrsim 5$.)

## 5.1 Overview

Fig. 5 shows the evolution of average mass, Rossby number, vertical Mach number, and stresses in the sample volume $x \in [-b, b]$ and time-averaged over five orbits. Without self-gravity ($Q_{3D} = 100$) the vortex decays by $t \sim 700P_0$, after which vertical motions remain small ($|v_z| \sim 10^{-3}c_s$) and the average Rossby number is roughly constant, as is $M_V$. The box returns to a laminar state with $\alpha \sim 10^{-7}$.

By contrast, the self-gravitating cases with $3 \le Q_{3D} \le 5$ are characterized by an increase in $M_V$ for $t \gtrsim 300P_0$ with a corresponding increase in $|Ro|$, $|v_z|$, and $\alpha$. These indicate sustained 3D turbulence developed from the EI as the vortex spins up due to self-gravitational growth. The self-gravitating vortices survive longer than the non-self-gravitating case. However, the $Q_{3D} = 3$ vortex weakens again after becoming too strong ($t \gtrsim 1500P_0$, see Section 5.5).

Fig. 6 shows the evolution of the Toomre $Q$ parameter measured at the vortex centre, the density-based vortex semiminor axis $b_\rho$, and aspect ratio $\chi_\rho$. The Toomre parameter remains $Q > 2$, which would suggest self-gravity to be unimportant in the classic sense. Evidently, this is not the case: we observe vortex growth even with a Toomre $Q \sim 8$ initially.
<table><caption>**Table 1.** Summary of simulations. In all runs an initial phase ($t \lesssim 300P_0$) of EI weakens/lengthens the vortex. The subsequent evolution is described below. The fiducial runs are in bold and discussed in the text.</caption><thead><tr><th>Run</th><th>Q<sub>3D</sub></th><th>z<sub>max</sub></th><th>χ</th><th>N<sub>x</sub>/L<sub>x</sub></th><th>Vortex evolution</th></tr></thead><tbody><tr><td>Q100</td><td>100</td><td>3H</td><td>10</td><td>32H</td><td>Vortex destroyed by initial EI and decays into an axisymmetric density bump.</td></tr><tr><td>Q100z2HR</td><td>100</td><td>2H</td><td>10</td><td>64H</td><td>Vortex destroyed by the initial EI.</td></tr><tr><td>Q5</td><td>5</td><td>3H</td><td>10</td><td>32H</td><td>Vortex survives until the end of the simulation (t = 2000P<sub>0</sub>) with a turbulent core with 3D motions.</td></tr><tr><td>Q4</td><td>4</td><td>3H</td><td>10</td><td>32H</td><td>Similar to Q5, but the vortex undergoes an azimuthal collapse from t ~ 1500P<sub>0</sub> when α<sub>R</sub> ≳ α<sub>G</sub> in the vortex and oscillations in |v<sub>z</sub>| develop. For t ≳ 1700P<sub>0</sub> max|v<sub>z</sub>| appears outside the vortex core. An overdense vortex remains at the end of the simulation.</td></tr><tr><td>Q4HR</td><td>4</td><td>3H</td><td>10</td><td>64H</td><td>Vortex destroyed by the initial EI.</td></tr><tr><td>Q4z2HR</td><td>4</td><td>2H</td><td>10</td><td>64H</td><td>Vortex survives with a turbulent core.</td></tr><tr><td>Q4c5</td><td>4</td><td>3H</td><td>5</td><td>32H</td><td>Vortex destroyed by initial EI.</td></tr><tr><td>Q4c20</td><td>4</td><td>3H</td><td>20</td><td>32H</td><td>Similar to Q4, but azimuthal collapse occurs sooner at t ~ 1000P<sub>0</sub>.</td></tr><tr><td>Q3</td><td>3</td><td>3H</td><td>10</td><td>32H</td><td>Similar to Q3D = 4, but azimuthal collapse begins earlier (t ~ 1200P<sub>0</sub>). Vertical motions dominate outside the vortex core from t ~ 1400P<sub>0</sub>. Vortex decays after reaching a size of b ~ H at t ~ 1500P<sub>0</sub> and inducing spiral shock waves in the ambient disc.</td></tr><tr><td>Q3HR</td><td>3</td><td>3H</td><td>10</td><td>64H</td><td>Vortex survives with a turbulent core.</td></tr><tr><td>Q3z2HR</td><td>3</td><td>2H</td><td>10</td><td>64H</td><td>Vortex survives with a turbulent core. (Simulation terminated at t = 1200P<sub>0</sub>.)</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Top to bottom: total mass, Rossby number, vertical Mach number, and stresses measured in the sample volume $|x| \le H/2$ containing the vortex for different strengths of self-gravity. These correspond to runs labelled Q3, Q4, Q5, and Q100 in Table 1.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Time evolution of the Toomre Q parameter at the vortex centre (top), vortex semiminor axis (middle), and aspect ratio (bottom). Measurements are based on the vortex density field. For $Q_{3D} = 100$ the Toomre $Q \gtrsim 150$ and is therefore not shown. Dashed lines in the middle panel are linear fits to the growth in $b_\rho$.</figcaption></figure>

The vortex aspect ratio first increases as the EI weakens the vortex. Without self-gravity the vortex eventually decays into a ring. For $Q_{3D} = 5$ the aspect ratio is maintained at $\chi_\rho \sim 14–15$. For $Q_{3D} = 3, 4$ the aspect ratio decreases, first slowly then rapidly from $t \sim 1200P_0, 1500P_0$ for $Q_{3D} = 3, 4$, respectively. Notice this drop coincides with self-gravity becoming dominant ($\alpha_G > \alpha_R$).

Finally, in Fig. 7 we compare the stresses averaged over the simulation domain excluding buffer zones. These are qualitatively similar to stresses measured at the box centre, but the total stress
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Stresses averaged over the simulation domain for the fiducial runs.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Hydrodynamic turbulence due to the EI of a 3D vortex in a self-gravitating disc (upper panels), and in a non-self-gravitating disc (lower panels). Comparing the upper and lower panels shows that hydrodynamic turbulence can be sustained (and grows) in the presence of self-gravity. Colours show the Rossby number $Ro = (\nabla \times \mathbf{v})/2\Omega + q/2$ at $z=0$.</figcaption></figure>

is always dominated by the Reynolds contribution. This reflects the fact that all of the discs under consideration are weakly self-gravitating in the classic sense.

## 5.2 Elliptic instability

We find early vortex evolution and weakening ($t \lesssim 300P_0$) are driven by EI and independent of self-gravity.<fnref n="1" /> The EI is signified by the exponential increase in $\sqrt{\langle v_z^2 \rangle}$ in $t \lesssim 100P_0$ with a growth rate $\gamma \simeq 0.01\Omega$ and saturating at $\sqrt{\langle v_z^2 \rangle} \sim 0.01c_s$ in all cases. This is consistent with previous non-self-gravitating studies of the EI, and the weak dependence on self-gravity is expected since the EI is incompressible (Lesur & Papaloizou 2009). The left-hand panels in Fig. 8 show the result of the EI is small-scale hydrodynamic turbulence. The subsequent evolution then depends on self-gravity.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Secular vortex evolution, after the initial development of EI and corresponding hydrodynamic turbulence. Top: a self-gravitating disc. Bottom: a non-self-gravitating disc. Colours show the relative density perturbation at the mid-plane.</figcaption></figure>

## 5.3 Secular vortex growth and sustained turbulence

Fig. 8 compares the secular evolution of the vortex after the initial EI. In both self-gravitating and non-self-gravitating cases, the initial EI results in an irregular vorticity field ($t = 250P_0$). However, it then reorganizes into a larger, coherent vortical structure by $t = 500P_0$, which serves as initial conditions for secular evolution. This initial vortex ‘revival’ may be due to the Kelvin–Helmholtz instability of vorticity strips that result from the shearing of vorticity seeds by the background shear flow left by the initial EI (Lithwick 2007). This second-generation vortex then undergoes EI again. However, in the non-self-gravitating case the vortex simply decays: we do not observe the reemergence of a large-scale coherent vortex due to Lithwick’s mechanism.

By contrast, we find that in self-gravitating discs the large-scale vortex persists with small-scale 3D hydrodynamic turbulence in its core (Fig. 8, upper panel). The vortex undergoes secular growth from $t \sim 300P_0$ to the end of the simulation for $Q_{3D} = 5$, but only to $t \sim 1500P_0$, $1200P_0$ for $Q_{3D} = 4$, 3, respectively. (In the latter cases the vortex then undergoes an azimuthal collapse, see Section 5.5.) Stresses $\langle \alpha \rangle$ and vertical velocity fluctuations $\sqrt{\langle v_z^2 \rangle}$ grow exponentially, but with a growth rate $\sim 2–9 \times 10^{-4} \Omega$, significantly slower than the initial EI because of the larger vortex aspect ratios than those imposed initially (Lesur & Papaloizou 2009). However, notice $\sqrt{\langle v_z^2 \rangle}$ grows faster with decreasing $Q_{3D}$. This indicates stronger EI, which is due to smaller vortex aspect ratios with increasing self-gravity (see Fig. 6).

Fig. 9 shows the evolution of the density field during secular growth. The self-gravitating vortex increases in density and size, but with decreasing aspect ratio. Notice the density field remains smooth even though the vorticity field is turbulent. This reflects the incom-

<footnote n="1">¹However, adjustment to the background Keplerian shear may also have minor contributions to the early vortex evolution.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Example of episodic EI in a self-gravitating 3D vortex. Top: Rossby number. Bottom: vertical Mach number.</figcaption></figure>

pressible nature of the EI. By contrast, the non-self-gravitating vortex is destroyed by the initial EI, leaving an axisymmetric density ring.

We note that during secular growth the monotonic increase in the vortex mass, rotation, and semiminor axis (Figs 5 and 6) in the self-gravitating cases proceed more rapidly than the decrease in its aspect ratio. (For $Q_{3D} = 5$ the vortex aspect ratio is nearly constant.) A possible reason is that vertical stratification due to self-gravity (see Fig. 4) makes the secular dynamics more similar to that of 2D vortices (Umurhan & Regev 2004; Johnson & Gammie 2005; Shen et al. 2006). These studies show that 2D vortices tend to merge into larger vortices (Godon & Livio 1999). Though merging cannot occur for our single, isolated vortex, secular evolution towards larger scales can be expected if the system behaves two-dimensionally due to the inverse cascade of energy (Kraichnan & Montgomery 1980). However, the additional stratification brought about by self-gravity is small; thus this effect is unlikely important for the present $Q_{3D}$ values. Here, the vortex growth in its radial size may be simply due to mass accretion.

## 5.4 Episodic bursts

We find secular vortex growth ends when $\alpha_R \sim \alpha_G \sim 10^{-5}$, at which point we observe *episodic bursts* of EI in self-gravitating cases. This is reflected in Fig. 5 as oscillations in $\sqrt{\langle v_z^2 \rangle}$, starting from $t \sim 1000P_0$, $1500P_0$ for $Q_{3D} = 3, 4$, respectively; and it only just begins to develop at the end of the $Q_{3D} = 5$ simulation. Fig. 10 shows snapshots of this episodic behaviour for $Q_{3D} = 4$ (where the simulation output frequency conveniently coincides with the burst). The vortex oscillates between having small-scale turbulence ($t = 1575P_0, 1625P_0$) and a coherent patch of negative vorticity ($t = 1600P_0, 1650P_0$). Notice in that $\max|Ro|$ also oscillates.

We comment that the oscillations are not reflected by the Rossby number in Fig. 5 because Ro has been averaged over $x \in [-b, b]$, which demonstrates that overall the vortex continues to strengthen. Nevertheless, we checked that $|Ro|$ does oscillate during the burst phase when the sampling volume is limited to near the vortex centre.

This episodic behaviour suggests competition between vortex growth due to self-gravity and vortex destruction by the EI. As the vortex grows, it spins up ($|Ro|$ increases) and its aspect ratio decreases. This feeds the EI, which converts some of the vortex's large-scale, planar rotation into small-scale 3D turbulence (signified by an increasing $|v_z|$). The EI activity decays, but is re-triggered

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Evolution of the $Q_{3D} = 3$ vortex during azimuthal collapse, when its aspect ratio drops rapidly. Left: mid-plane density perturbation. Middle: mid-plane Rossby number. Right: average vertical Mach number. Here, ellipses in Ro are drawn by inspection, and correspond to the loci $Ro = 0.65 \min(Ro)$ and $Ro = 0.3 \min(Ro)$ in the top and bottom snapshots, respectively.</figcaption></figure>

once the vortex reforms as it grows via self-gravity. This growth may be partly attributable to the removal of the vortex's internal rotation by the EI and enabling secular gravitational instabilities (see Section 6.2).

A similar bursty behaviour has been reported in 3D simulations of vortex amplification by the sub-critical baroclinic instability in competition with EI (Lesur & Papaloizou 2010, their fig. 16). We discuss this analogy further in Section 6.3.

## 5.5 Azimuthal collapse, gap opening, and self-destruction

As Reynolds stresses ($\alpha_R$) increase beyond the gravitational stresses ($\alpha_G$), the vortex in the $Q_{3D} = 3, 4$ discs enters a second, faster collapse phase where its aspect ratio rapidly drops, although its radial size continues to grow, see Figs 5 and 6. This azimuthal collapse phase corresponds to $t \in [1200, 1500]P_0$ and $t \in [1500, 2000]P_0$ for $Q_{3D} = 3, 4$, respectively. The total stress also increases rapidly (see Fig. 7).

Fig. 11 shows azimuthal collapse for the $Q_{3D} = 3$ simulation. The vortex exits the bursty phase with a coherent vorticity patch of an aspect ratio $\chi \sim 7$ (upper middle panel). Notice $\max|v_z|$ moves outside the vortex core as its aspect ratio decreases (lower right-hand panel). This is in fact consistent with Lesur & Papaloizou (2009), who find the EI occurs outside the vortex core for an intermediate range of aspect ratios in the absence of buoyancy, as is the case for isothermal gas considered here. In their specific case of Kida vortices (Kida 1981), this range is $4 \lesssim \chi \lesssim 6$, but the corresponding range may differ for the vortices in our simulations.

We find azimuthal collapse is associated with vortex-driven spiral density waves. Spiral density waves are naturally excited by vortex modes in compressible gas (Bodo et al. 2005; Heinemann & Papaloizou 2009; Mamatsashvili & Chagelishvili 2007). However, as the vortex half-width approaches $H$, these spiral density waves steepen into shocks (Paardekooper et al. 2010) and transfer angular momentum and energy away from the vortex (Bodo et al. 2007).
These spiral shocks can also be identified visually in Fig. 11 at $t = 1500P_0$ in both the density and vorticity maps.

Azimuthal collapse is terminated once the vortex half-width exceeds ~H. This is seen in Fig. 5 as the rapid drop in $|Ro|$ towards the end of the simulation with $Q_{3D} = 3$. The corresponding decrease in $M_V$ indicates mass-loss by the vortex to the ambient disc. This final vortex weakening is likely related to the strong spiral shocks induced by the vortex once it becomes too large in the radial direction. We observe a similar self-destruction behaviour as above in an extended run of the $Q_{3D} = 4$ case. However, a proper study of this final stage requires global simulations to minimize the influence of radial boundary conditions. We thus limit our discussion here to when $b_ρ \lesssim H$. (See Appendix C for examples of vortex evolution in 2D, global discs.)

Nevertheless, our observation that vortex radial size being limited by the associated density waves is consistent with previous studies of vortex evolution in 2D. The critical radial extent is expected to be ~H (Godon & Livio 1999), beyond which the vortex flow becomes sonic (Paardekooper et al. 2010). This is equivalent to the Jeans length $λ_J \sim O(QH)$ in a self-gravitating disc with $Q \sim 1$. Indeed, Mamatsashvili & Rice (2009) suggest 2D vortices are limited to $λ_J$ in their radial extent. They find in gravito-turbulent discs that upon reaching $λ_J$, 2D vortices weaken and are sheared away by the background flow (see also Lyra et al. 2009b). In our discs with $Q > 1$, $λ_J$ is somewhat larger than $H$. The vortex width limit in our 3D simulations is thus similar to that in corresponding 2D, non-self-gravitating discs.

Notice in Fig. 11 that two shallow gaps have developed on either side of the vortex. Similar double-gap structures appear in low-mass planet-disc simulations (e.g. Zhu et al. 2013; Dong et al. 2017). This suggests that the vortex can interact with the ambient disc in a similar way to disc-planet interaction. The strength of gravitational interaction between the vortex and the ambient disc can be estimated via the total mass $M_{\text{pert}}$ enclosed in the elliptical patch ($x^2 + y^2/χ_ρ^2)/b_ρ^2 \lesssim 1$:

$$q_v \equiv \frac{M_{\text{pert}}}{M_* h^3} = \frac{\chi_ρ}{Q} \left(\frac{b_ρ}{H}\right)^2, \quad (35)$$

where $M_*$ is the mass of the central star and $h \equiv H/R$ is the characteristic disc aspect ratio at the fiducial radius of the shearing box. Using Fig. 6 we find $q_v \sim 1.3$ for the $Q_{3D} = 3$ simulation at $t = 1500P_0$. This can be compared to the thermal condition for gap-opening by disc-planet interaction, $q_v \gtrsim 1$ (Crida et al. 2006). Thus, gap-opening by the vortex can indeed be expected. However, we expect gap-opening by an elliptical, extended vortex to be weaker than by a point-mass planet of the same mass.

The vortex-induced spiral shocks may also contributed to mass accretion on to the vortex, similar to shock-driven accretion in circumplanetary discs (Zhu et al. 2016). Both (1) the increasing vortex mass; and (2) the decreasing vortex aspect ratio are expected to enhance the vortex-driven spiral density waves. Thus, azimuthal collapse proceeds more rapidly than the secular growth, when the associated spiral waves are weak. Turbulence outside the vortex core may also contribute to angular momentum removal and hence mass accretion by the vortex.

## 5.6 Vortex aspect ratio, numerical resolution, and box height

We performed several variations of the fiducial runs described above. These are summarized in Table 1 and include different aspect ratios for the initial vortex perturbation ($χ$), the radial resolution ($L_x/N_x$), and the vertical domain size ($z_{\text{max}}$),

<figure><img src="image_14.png" /><figcaption>Figure 12. Comparison between vortex evolution in a $Q_{3D} = 3$ disc, but different radial resolutions: $N_x/L_x = 32/H$ (black) and $64/H$ (orange). The vortex persists in both simulations.</figcaption></figure>

For the $Q_{3D} = 4$ disc we ran additional cases with $χ = 20$ and $χ = 5$ for the initial vortex perturbation. The former behaves similarly to our fiducial run with $χ = 10$. However, for $χ = 5$ the vortex is destroyed by the initial phase of EI, similar to the non-self-gravitating case $Q_{100}$. This suggests that self-gravity cannot save vortices born with small aspect ratios.

We also repeated some runs with twice the radial resolution. These are compared in Fig. 12 for the $Q_{3D} = 3$ disc. We find the initial EI growth rate is converged, but with better resolution the resulting turbulence leads to a weaker vortex to serve as initial conditions for secular growth. This is shown in Fig. 13 as the more elongated vortex at $t = 500P_0$ compared to that in Fig. 8. Since the EI turbulence is better resolved, it has a more destructive effect and the subsequent secular growth rate is slower and a 3D vortex with a turbulent core persists until the end of the simulation.

We find that reducing the vertical domain also enhances vortex survival, which is as expected as the system is forced to behave more two-dimensionally (Lithwick 2009). This is demonstrated in Fig. 14 with simulations using $Q_{3D} = 4$. Starting from the fiducial run with $z_{\text{max}} = 3H$, we double the radial resolution and find the vortex is destroyed by the initial EI. (Comparing with Fig. 13 confirms that self-gravity enhances vortex survival.) However, retaining the same resolution but setting $z_{\text{max}} = 2H$ results in a persistent 3D vortex with a turbulent core, also shown in Fig. 15.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Vortex evolution in the $Q_{3D} = 3$ disc but with double the radial resolution than the fiducial runs.</figcaption></figure>

These additional runs indicate that, in order to undergo self-gravitational growth, a 3D vortex needs to maintain its large-scale internal flow. Thus the initial EI cannot be too strong. This can be limited by diffusion (numerical in our case) and/or limited vertical domain. Notice in the simulations where a vortex persists that $\alpha_G > \alpha_R$ during secular growth, as observed in the standard run. We suggest this condition required for long-term vortex survival. If the initial EI turbulence destroys the large-scale vortex topology, then there is no vortex for self-gravity to maintain or amplify.

# 6 DISCUSSION

## 6.1 Self-gravitational vortex growth and hydrodynamic turbulence

Our main finding is that self-gravitating 3D vortices can survive the EI. Self-gravity enables vortex survival and growth even when the Toomre parameter $Q \sim 5$, which is significantly larger than that required for axisymmetric instability in 3D shearing boxes ($Q \lesssim 0.6$, Mamatsashvili & Rice 2010). Thus, the presence of a vortex reduces the gravitational stability of the flow.

We suggest that, provided the initial EI does not completely destroy the large-scale vortex flow topology (see Sections 5.6 and 6.6), the vortex can then slowly grow under its self-gravity. (Without self-gravity the vortex always decays.) The growth spins up the vortex and its aspect ratio decreases. This is expected for vortices in shearing discs, e.g. the GNG vortex spin increases with decreasing $\chi$ (see equation 16). Smaller aspect ratio vortices favour EI (Lesur & Palaoizou 2009), which generates further hydrodynamic turbulence.

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Vortex evolution in discs with $Q_{3D} = 4$ but different radial resolutions and vertical domain sizes.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Effect of vertical domain size on vortex evolution. The disc has $Q_{3D} = 4$ and the simulations were performed with twice the fiducial radial resolution. The vortex does not survive with $z_{max} = 3H$ but does with $z_{max} = 2H$.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 16.** The ‘gravito-elliptic’ interpretation of the bursty behaviour observed towards the end of secular vortex growth. A vortex spins up and its aspect ratio decreases, generating hydrodynamic turbulence via the EI. The associated turbulent viscosity favours gravitational instability by removing rotational support against self-gravity.</figcaption></figure>

The vortex persists during secular growth because the subsequent EI is weak, as the vortex has been elongated by the initial EI. Our numerical experiments indicate a necessary condition for vortex survival is having gravitational stresses exceed Reynolds stresses ($\alpha_G \gtrsim \alpha_R$). It signifies the evolution being dominated by self-gravity with hydrodynamic turbulence being a by-product.

To check that self-gravity can directly cause vortex growth, we ran several 2D simulations in Appendix C. We find for sufficiently strong self-gravity, the vortex can indeed grow without EI-induced turbulence. However, growth does not occur in razor-thin discs if $Q_{3D} \ge 4$, whereas we find growth in corresponding 3D simulations with EI. This suggests that EI-induced hydrodynamic turbulence may also contribute to gravitational instability, as discussed below.

## 6.2 Turbulence-enabled gravitational instability?

The EI converts the vortex’s internal large-scale, planar rotation into small-scale 3D turbulence. However, removing rotation in a self-gravitating flow allows gravitational collapse. This type of gravitational instability relies on dissipative processes to remove rotational stabilization (Lynden-Bell & Pringle 1974). Comparing our 3D simulations that have EI and show vortex growth and supplementary 2D simulations in Appendix C that show no EI or growth suggests that EI turbulence may provide a similar effect; thus enabling vortex growth at large $Q$. This picture is similar to secular gravitational instabilities in viscous or dusty discs, whereby friction allows gravitational instabilities to develop even when $Q > 1$ (e.g. Schmit & Tscharnuter 1995; Takahashi & Inutsuka 2014).

Gravitational instability catalysed by EI-induced turbulence leads to the possibility of ‘gravito-elliptic’ feedback, illustrated in Fig. 16. Vortex growth leads to hydrodynamic turbulence through the EI, which helps further collapse via secular gravitational instabilities, thus creating a cycle. The vortex attains a quasi-steady state if the destructive effect of EI is balanced by gravitational instability. This explains the ‘episodic bursts’ seen in our fiducial simulations when $\alpha_R \sim \alpha_G$ (Section 5.4).

## 6.3 Analogy with SBI in 3D

The episodic EI bursts observed in Section 5.4 bear some parallel to the SBI in 3D. The SBI is a mechanism for amplifying vortices in non-isothermal discs and requires an unstable radial entropy gradient and an appropriate thermal diffusion or cooling time-scale (Lesur & Papaloizou 2010). Neither are present in our standard, isothermal shearing boxes; Lesur & Papaloizou also found episodic EI in 3D simulations of the SBI.

In 3D, the SBI competes with the EI (Lesur & Papaloizou 2010; Lyra & Klahr 2011; Barge et al. 2016). These authors find the resulting vortex maintains its large-scale flow, but also contains small-scale, 3D hydrodynamic turbulence in its core. This is similar

to our self-gravitating vortices (e.g. Fig. 8). Self-gravity plays the role of SBI in our simulations by amplifying the vortex.

The important difference is that for SBI there exists an explicit baroclinic vorticity source, which comes from the global radial entropy gradient. Self-gravity cannot directly source vorticity. However, the vortex can spin-up as it collapses due to self-gravity. The result in both cases is that the vortex aspect ratio decreases, feeding the EI, which tries to remove the large-scale vortex rotation.

## 6.4 Application to Rossby vortices

We can relate our simulations to vortex formation by the RWI (Lovelace et al. 1999; Li et al. 2000, 2001). The RWI is a dynamical instability associated with radially structured discs. More specifically, it requires a minimum in the disc’s potential vorticity profile ($=\omega/\Sigma$). Candidate sites include gap edges opened by planets (Li et al. 2009; Lyra et al. 2009a; Lin 2012a) and dead zone boundaries (Lyra et al. 2008, 2009b). These isolated radial structures often involve length-scales of $O(H)$, and are thus consistent with our typical vortex half-widths ($\lesssim H$).

The expected initial aspect ratio of Rossby vortices is $\chi = \pi/hm_{max}$; here $m_{max} \sim 5-8$ is the most unstable azimuthal wavenumber of the RWI (Lin & Papaloizou 2011) and $h = H/R$ is the disc aspect ratio. For protoplanetary discs with $h = 0.05-0.1$ we expect $\chi \sim 4-12$. Our nominal choice $\chi = 10$ is consistent with the lower bound on lower $m_{max}$, which is applicable to thin and/or low-mass discs (Lin & Papaloizou 2011) as considered here. This is also the limit where Rossby vortices can merge. Thus our model could also correspond to a single, post-merger Rossby vortex.

## 6.5 Dust-trapping and observational implications

Our main result is that moderate disc self-gravity can help sustain 3D vortices in spite of the EI. This means that vortices are more likely to be observed in the outer parts of protoplanetary discs, where self-gravity becomes important. This could be relevant to vortex formation via the RWI due to gap-opening planets at large radii (Hammer, Kratter & Lin 2017) or the outer dead zone boundary (Lyra, Turner & McNally 2015).

For definiteness, let us consider the class II circumstellar disc models described in Kratter & Lodato (2016, their fig. 2a). For typical global accretion rates $\dot{M} \in [10^{-8}, 10^{-7}] M_{\odot} \text{yr}^{-1}$, the Toomre parameter ranges from $Q \simeq 10$ at 30 au to $Q \simeq 5$ at 100 au. These are similar to the $Q$ values in our simulations. Furthermore, at these radii, heating is dominated by stellar irradiation (Chiang & Goldreich 1997) and thus disc behaves almost isothermally; also consistent with our equation of state. Note also our simulation time-scales of $O(10^3)$ orbits correspond $10^6$ yr at 100 au, which is of the order of the disc lifetime.

Our fiducial simulation ($Q_{3D} = 3$ or $Q \simeq 5$) corresponds to 100 au in this physical disc model. Although $Q > 2$ is traditionally considered non-self-gravitating (Kratter & Lodato 2016), neglecting self-gravity altogether would result in a weak vortex with vanishing overdensity. However, properly including self-gravity gives a concentrated vortex (Fig. 9). This would suggest that the trapping of dust – often the observable – is favoured in self-gravitating discs.

On the other hand, self-gravitating 3D vortices have turbulent cores, which tends to expel dust particles due to turbulent dust diffusion (Youdin & Lithwick 2007). To examine this effect we use the dust-trapping model described in Lyra & Lin (2013) assuming an underlying GNG vortex. They find the distribution of small,
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Normalized (mid-plane) dust density distribution obtained from the dust-trap model of Lyra & Lin (2013) with vortex parameters taken from the fiducial simulations. We fix the Stokes number $St = 10^{-3}$.</figcaption></figure>

passive dust particles within the vortex core is given by

$$ \rho_d(s) = \rho_{d,max} \exp \left( -\frac{s^2}{2H_v^2} \right), \quad (36) $$

$$ H_v = \frac{H}{f(\chi)} \sqrt{\frac{\delta}{St + \delta}}, \quad (37) $$

where $St = \Omega\tau_s$ is the Stokes number and $\tau_s$ is the stopping time characterizing dust-gas friction and $\delta$ is a dimensionless parametrization of turbulent dust diffusion. The GNG vortex has $f(\chi) \sim 0.7$ for $\chi \gg 1$ (see Lyra & Lin 2013, for details). Recall $s$ is a coordinate that labels elliptical streamlines in the vortex by their semiminor axes (see Section 3).

We evaluate equation (37) using simulation results as input parameters. Assuming dust-gas interaction is dominated by drag forces, we set $\delta = \alpha_R$. We also use $\chi_\rho$ in place of $\chi$, but this is immaterial since the function $f$ is effectively constant for elongated vortices.

In Fig. 17 we plot the expected distribution of dust particles with $St = 10^{-3}$. We only consider the secular growth stage with EI turbulence occurring inside the vortex core, where Lyra & Lin's model is applicable. In all cases dust is initially concentrated towards the centre, but over time it diffuses outwards due to the increasing levels of hydrodynamic turbulence, which in turn is more vigorous with increasing self-gravity. Thus we can expect the dusty vortices to be wider in self-gravitating discs.

We comment that the above result for a single vortex may not be applicable to strongly self-gravitating discs with gravitoturbulence/fragmentation (Gammie 2001). However, in that limit it has also been shown that particle trapping is enhanced by self-gravity (Gibbons, Rice & Mamatsashvili 2012; Gibbons, Mamatsashvili & Rice 2014, 2015; Shi et al. 2016).

## 6.6 Caveats and future directions

We have limited our study to classically stable discs, i.e. $Q > 2$ (or $Q_{3D} > 1$). It would be interesting to study 3D vortex evolution in strongly self-gravitating discs with $Q \lesssim 2$. However, in this

limit the vortex must be initialized with care. Simply perturbing a strongly self-gravitating disc with the non-self-gravitating GNG vortex flow may result in transient effects. This problem may also be ill-posed because it is unclear if a single large-scale vortex can form in strongly self-gravitating discs in the first place. For example, it has been shown that the vortex-forming RWI can be suppressed by self-gravity (Lin & Papaloizou 2011; Lovelace & Hohlfeld 2013). Furthermore, classic gravitational instabilities can be expected in the ambient disc, which would complicate the evolution of an isolated vortex. We thus defer the strongly self-gravitating limit to a future study.

Returning to the $Q > 5$ discs considered in our study, the additional runs described in Section 5.6 indicate that 3D vortices undergo secular self-gravitational growth only if the initial EI-induced turbulence does not destroy the large-scale vortex topology. In our simulations the initial EI turbulence strength is limited by a reduced vertical domain (Lithwick 2009) and/or numerical diffusion (Lesur & Papaloizou 2009). Since our models are inviscid, true numerical convergence may be difficult to achieve.

However, these effects may not be unrealistic. Numerical diffusion may mimic other sources of low-level turbulent viscosity in a protoplanetary disc (Umurhan et al. 2016a). Vortices can also be limited in vertical extent if the background disc has a layered structure (Bai 2016). The presence of an active vortex-driving source may also help to retain the large-scale vortex flow against the initial EI. Nevertheless, future works should incorporate these physical effects explicitly.

Local shearing box simulations do not permit a proper study of vortex gap-opening and migration (Paardekooper et al. 2010), which stem from global spiral density waves launched by the vortex. Clearly, one must eventually perform global disc simulations of self-gravitating 3D vortices. This is computationally challenging because of the need to simultaneously resolve the small-scale EI. Though some recent non-self-gravitating, global simulations have done so (Barge et al. 2016; Richard et al. 2016), including 3D self-gravity can significantly increase the computational cost.

However, our simulations show that even with EI-induced hydrodynamic turbulence, the vortex maintains a smooth density field with moderate amplitude. Furthermore, elongated vortices do not require high resolution in azimuth. This means that an accurate global Poisson solver (with sub-$H$ resolution) is not essential for this problem.

We have purposefully neglected vortex-forming mechanisms to isolate the effect of self-gravity on vortex evolution. Nevertheless, future work should account for hydrodynamic instabilities that drive vortex formation (RWI, ConO/SBI, ZVI, or VSI), which require improved geometries and physics such as cooling.

Finally, we performed pure gas simulations and only estimated the extent of vortex dust trapping in post-processing. It will be desirable to conduct two-phase, gas-dust simulations to properly capture the evolution of dust particles and the back reaction of drag forces on the gas. The latter becomes important for high dust-to-gas ratios, which can be expected in vortex centres (Fu et al. 2014; Crnkovic-Rubsamen, Zhu & Stone 2015; Surville, Mayer & Lin 2016). A simple starting point to model these effects is the one-fluid, thermodynamic model of dusty gas recently developed by Lin & Youdin (2017).

## 7 SUMMARY

In this paper we study the evolution of 3D vortices in self-gravitating protoplanetary discs via customized numerical simulations in the
shearing box framework. We find self-gravity can help 3D vortices survive against elliptic instabilities that would otherwise destroy them in non-self-gravitating discs. With moderate self-gravity the result is a 3D vortex with a turbulent core. We emphasize that our disc models are gravitationally stable in the traditional sense that Toomre Q values safely exceed unity; yet we observe vortex growth, albeit on time-scales of $O(10^3)$ orbits. This indicates that the presence of a vortex reduces the long-term stability of the flow.

In addition to hydrodynamic turbulence as a by-product of self-gravitational vortex growth via elliptic instabilities, we suggest secular gravitational instabilities can be catalysed by the hydrodynamic turbulence that results from the EI. The EI removes rotational support against self-gravity, which favours gravitational instability; but the subsequent spin-up of the vortex feeds the EI. There is competition between the vortex destruction by the EI and vortex growth due to self-gravity (either directly or assisted by hydrodynamic turbulence); thus a quasi-steady state is possible. Our simulations suggest vortices survive more easily, and hence are more likely observed, in protoplanetary discs at large radii.

## ACKNOWLEDGEMENTS

We thank the anonymous referee for a detailed report that led to an improved physical interpretation and presentation of our results. This work is supported by the Theoretical Institute for Advanced Research in Astrophysics (TIARA) based in Academia Sinica Institute of Astronomy and Astrophysics (ASIAA). All 3D simulations were performed on the TIARA High Performance Computing cluster. The early stages of this project were carried out at the University of Arizona under the support of the Steward Theory Fellowship and made use of the El Gato and Ocelote clusters. This work is also part of the NASA Astrophysics Theory Programme NNX17AK59G. Computer time for the 2D runs presented in this paper was provided by HPC resources of Cines under allocation A0010406957 made by GENCI (Grand Equipement National de Calcul Intensif).

## REFERENCES

Adams F. C., Watkins R., 1995, ApJ, 451, 314

Bai X.-N., 2016, ApJ, 821, 80

Barge P., Sommeria J., 1995, A&A, 295, L1

Barge P., Richard S., Le Dizès S., 2016, A&A, 592, A136

Barker A. J., Latter H. N., 2015, MNRAS, 450, 21

Bodo G., Chagelishvili G., Murante G., Tevzadze A., Rossi P., Ferrari A., 2005, A&A, 437, 9

Bodo G., Tevzadze A., Chagelishvili G., Mignone A., Rossi P., Ferrari A., 2007, A&A, 475, 51

Casassus S. et al., 2013, Nature, 493, 191

Chandrasekhar S., 1969, Ellipsoidal Figures of Equilibrium, Yale University Press, New Haven

Chang P., Oishi J. S., 2010, ApJ, 721, 1593

Chiang E. I., Goldreich P., 1997, ApJ, 490, 368

Crida A., Morbidelli A., Masset F., 2006, Icarus, 181, 587

Crnkovic-Rubsamen I., Zhu Z., Stone J. M., 2015, MNRAS, 450, 4285

Dong R., Li S., Chiang E., Li H., 2017, ApJ, 843, 127

Fu W., Li H., Lubow S., Li S., 2014, ApJ, 788, L41

Fukagawa M. et al., 2013, PASJ, 65, L14

Gammie C. F., 2001, ApJ, 553, 174

Gibbons P. G., Rice W. K. M., Mamatsashvili G. R., 2012, MNRAS, 426, 1444

Gibbons P. G., Mamatsashvili G. R., Rice W. K. M., 2014, MNRAS, 442, 361

Gibbons P. G., Mamatsashvili G. R., Rice W. K. M., 2015, MNRAS, 453, 4232

Godon P., Livio M., 1999, ApJ, 523, 350

Goldreich P., Lynden-Bell D., 1965, MNRAS, 130, 125

Goodman J., Narayan R., Goldreich P., 1987, MNRAS, 225, 695

Hammer M., Kratter K. M., Lin M.-K., 2017, MNRAS, 466, 3533

Hashimoto J. et al., 2015, ApJ, 799, 43

Heinemann T., Papaloizou J. C. B., 2009, MNRAS, 397, 52

Inaba S., Barge P., 2006, ApJ, 649, 415

Isella A., Pérez L. M., Carpenter J. M., Ricci L., Andrews S., Rosenfeld K., 2013, ApJ, 775, 30

Johnson B. M., Gammie C. F., 2005, ApJ, 635, 149

Kerswell R. R., 1994, J. Fluid Mech., 274, 219

Kida S., 1981, J. Phys. Soc. Japan, 50, 3517

Klahr H., Hubbard A., 2014, ApJ, 788, 21

Kraichnan R. H., Montgomery D., 1980, Rep. Prog. Phys., 43, 547

Kratter K., Lodato G., 2016, ARA&A, 54, 271

Kraus S. et al., 2017, ApJ, 848, L11

Latter H. N., 2016, MNRAS, 455, 2608

Lesur G. R. J., Latter H., 2016, MNRAS, 462, 4549

Lesur G., Papaloizou J. C. B., 2009, A&A, 498, 1

Lesur G., Papaloizou J. C. B., 2010, A&A, 513, A60

Li H., Finn J. M., Lovelace R. V. E., Colgate S. A., 2000, ApJ, 533, 1023

Li H., Colgate S. A., Wendroff B., Liska R., 2001, ApJ, 551, 874

Li H., Lubow S. H., Li S., Lin D. N. C., 2009, ApJ, 690, L52

Lin M.-K., 2012a, MNRAS, 426, 3211

Lin M.-K., 2012b, ApJ, 754, 21

Lin M.-K., Kratter K. M., 2016, ApJ, 824, 91

Lin M.-K., Papaloizou J. C. B., 2011, MNRAS, 415, 1426

Lin M.-K., Youdin A. N., 2015, ApJ, 811, 17

Lin M.-K., Youdin A. N., 2017, ApJ, 849, 129

Lithwick Y., 2007, ApJ, 670, 789

Lithwick Y., 2009, ApJ, 693, 85

Lovelace R. V. E., Hohlfeld R. G., 2013, MNRAS, 429, 529

Lovelace R. V. E., Li H., Colgate S. A., Nelson A. F., 1999, ApJ, 513, 805

Lynden-Bell D., Pringle J. E., 1974, MNRAS, 168, 603

Lyra W., 2014, ApJ, 789, 77

Lyra W., Klahr H., 2011, A&A, 527, A138

Lyra W., Lin M.-K., 2013, ApJ, 775, 17

Lyra W., Johansen A., Klahr H., Piskunov N., 2008, A&A, 491, L41

Lyra W., Johansen A., Klahr H., Piskunov N., 2009a, A&A, 493, 1125

Lyra W., Johansen A., Zsom A., Klahr H., Piskunov N., 2009b, A&A, 497, 869

Lyra W., Turner N. J., McNally C. P., 2015, A&A, 574, A10

Mamatsashvili G. R., Chagelishvili G. D., 2007, MNRAS, 381, 809

Mamatsashvili G. R., Rice W. K. M., 2009, MNRAS, 394, 2153

Mamatsashvili G. R., Rice W. K. M., 2010, MNRAS, 406, 2050

Marcus P. S., Pei S., Jiang C.-H., Barranco J. A., Hassanzadeh P., Lecoanet D., 2015, ApJ, 808, 87

Marino S., Casassus S., Perez S., Lyra W., Roman P. E., Avenhaus H., Wright C. M., Maddison S. T., 2015, ApJ, 813, 76

Meheut H., Casse F., Varniere P., Tagger M., 2010, A&A, 516, A31

Meheut H., Meliani Z., Varniere P., Benz W., 2012, A&A, 545, A134

Nelson R. P., Gressel O., Umurhan O. M., 2013, MNRAS, 435, 2610

Ohta Y. et al., 2016, PASJ, 68, 53

Ono T., Muto T., Takeuchi T., Nomura H., 2016, ApJ, 823, 84

Paardekooper S., Lesur G., Papaloizou J. C. B., 2010, ApJ, 725, 146

Pérez L. M., Isella A., Carpenter J. M., Chandler C. J., 2014, ApJ, 783, L13

Petersen M. R., Julien K., Stewart G. R., 2007, ApJ, 658, 1236

Raettig N., Lyra W., Klahr H., 2013, ApJ, 765, 115

Railton A. D., Papaloizou J. C. B., 2014, MNRAS, 445, 4409

Regály Z., Vorobyov E., 2017, MNRAS, 471, 2204

Richard S., Nelson R. P., Umurhan O. M., 2016, MNRAS, 456, 3571

Schmit U., Tscharnuter W. M., 1995, Icarus, 115, 304

Shen Y., Stone J. M., Gardiner T. A., 2006, ApJ, 653, 513

Shi J.-M., Zhu Z., Stone J. M., Chiang E., 2016, MNRAS, 459, 982

Stoll M. H. R., Kley W., 2014, A&A, 572, A77

Stone J. M., Gardiner T. A., 2010, ApJS, 189, 142

Stone J. M., Gardiner T. A., Teuben P., Hawley J. F., Simon J. B., 2008, ApJS, 178, 137
Surville C., Mayer L., Lin D. N. C., 2016, ApJ, 831, 82

Takahashi S. Z., Inutsuka S.-i., 2014, ApJ, 794, 55

Toomre A., 1964, ApJ, 139, 1217

Umurhan O. M., Regev O., 2004, A&A, 427, 855

Umurhan O. M., Estrada R. R., Cuzzi J. N., 2016a, in Lunar and Planetary
Science Conference, 47. p. 2887

Umurhan O. M., Shariff K., Cuzzi J. N., 2016b, ApJ, 830, 95

van der Marel N. et al., 2013, Science, 340, 1199

van der Marel N., Cazzoletti P., Pinilla P., Garufi A., 2016, ApJ, 832, 178

Yellin-Bergovoy R., Heifetz E., Umurhan O. M., 2016, Geophys. Astrophys.
Fluid Dyn., 110, 274

Youdin A. N., 2011, ApJ, 731, 99

Youdin A. N., Lithwick Y., 2007, Icarus, 192, 588

Zhu Z., Baruteau C., 2016, MNRAS, 458, 3918

Zhu Z., Stone J. M., Rafikov R. R., 2013, ApJ, 768, 143

Zhu Z., Stone J. M., Rafikov R. R., Bai X.-n., 2014, ApJ, 785, 122

Zhu Z., Ju W., Stone J. M., 2016, ApJ, 832, 193

APPENDIX A: SHEARING BOX EQUATIONS IN
ELLIPTICO-POLAR COORDINATES

In elliptico-polar coordinates (s, φ, z) the divergence and advective derivative operators have identical form to that in the usual cylindrical coordinates. That is,

$$
\nabla \cdot \mathbf{W} = \frac{1}{s} \frac{\partial}{\partial s} (s W_s) + \frac{1}{s} \frac{\partial W_{\varphi}}{\partial \varphi} + \frac{\partial W_z}{\partial z}, \quad (A1)
$$

$$
\mathbf{W} \cdot \nabla = W_s \frac{\partial}{\partial s} + \frac{W_{\varphi}}{s} \frac{\partial}{\partial \varphi} + W_z \frac{\partial}{\partial z}, \quad (A2)
$$

for any vector field W. The continuity equation thus retains the
same form as in the Cartesian box, as does the vertical momen-
tum equation. The horizontal momentum equations, however, are
modified to read

$$
\begin{align}
\frac{\partial v_s}{\partial t} + \mathbf{v} \cdot \nabla v_s - \frac{v_\psi^2}{s} & \nonumber \\
&= \frac{1}{2} (\xi_+ + \xi_- \cos 2\varphi) \left( 2\Omega \chi v_\psi - \frac{\partial \tilde{\eta}}{\partial s} \right) \nonumber \\
&\quad + \frac{\xi_-}{2} \sin 2\varphi \left( 2\Omega \chi v_s + \frac{1}{s} \frac{\partial \tilde{\eta}}{\partial \varphi} \right) + q\Omega^2 s (\cos 2\varphi + 1) \tag{A3}
\end{align}
$$

$$
\begin{align}
\frac{\partial v_{\varphi}}{\partial t} + \mathbf{v} \cdot \nabla v_{\varphi} + \frac{v_{\varphi} v_s}{s} \nonumber \\
&= \frac{\xi_{-}}{2} \sin 2\varphi \left( \frac{\partial \tilde{\eta}}{\partial s} - 2\Omega \chi v_{\varphi} \right) \nonumber \\
&\quad - \frac{1}{2} (\xi_{+} - \xi_{-} \cos 2\varphi) \left( \frac{1}{s} \frac{\partial \tilde{\eta}}{\partial \varphi} + 2\Omega \chi v_s \right) - q\Omega^2 s \sin 2\varphi, \tag{A4}
\end{align}
$$

where $\tilde{\eta} = \eta + \Phi$ and

$$
\xi_{\pm} = 1 \pm \frac{1}{\chi^2}. \tag{A5}
$$

Finally, the Poisson equation becomes

$$
\begin{bmatrix} \frac{1}{2}(\xi_{-} \cos 2\varphi + \xi_{+}) \frac{\partial^2}{\partial s^2} + \frac{1}{2s^2}(\xi_{+} - \xi_{-} \cos 2\varphi) \frac{\partial^2}{\partial \varphi^2} \\ -\frac{\xi_{-}}{s} \sin 2\varphi \frac{\partial^2}{\partial \varphi \partial s} + \frac{1}{2s}(\xi_{+} - \xi_{-} \cos 2\varphi) \frac{\partial}{\partial s} + \frac{\xi_{-}}{s^2} \sin 2\varphi \frac{\partial}{\partial \varphi} \\ + \frac{\partial^2}{\partial z^2} \end{bmatrix} \Phi = 4\pi G \rho. \quad (A6)
$$

APPENDIX B: LINEAR RESPONSE OF THE
GNG VORTEX TO ITS GRAVITATIONAL
POTENTIAL

We explore the perturbative effect of introducing self-gravity on the
GNG vortex solution. We simplify the problem by considering an
unstratified disc with no vertical velocities. Then $v_z = \partial_z = 0$. This
should be an adequate approximation for regions close to the vortex
centre ($s, z \to 0$). We set

$$
G \to 0 + \delta G. \tag{B1}
$$

This induces perturbations

$$
\rho \rightarrow \rho(s) + \delta\rho, \tag{B2}
$$

$$
v_s \to 0 + \delta v_s, \quad (\text{B3})
$$

$$
v_{\varphi} \rightarrow -\Omega_{\nu}s + \delta v_{\varphi}, \quad (B4)
$$

$$
\Phi \to 0 + \delta\Phi. \tag{B5}
$$

The background flow, $v_\varphi(s)$ and $\rho(s)$, is given by the non-self-
gravitating GNG vortex solution (equations 16 and 17) with the
vertical dimension suppressed. Since the basic state is symmetric in
$\varphi$, we Fourier analyse in azimuth and write

$$
\begin{bmatrix} \delta\rho \\ \delta v_s \\ \delta v_\varphi \\ \delta\Phi \end{bmatrix} = \text{Re} \sum_{l=0}^{\infty} \begin{bmatrix} \rho \frac{Q_l}{c_s^2} \\ iU_l \\ V_l \\ \Phi_l \end{bmatrix} \exp il\varphi. \quad \text{(B6)}
$$

We insert the above perturbations into the steady-state continuity,
horizontal momentum, and Poisson equations and linearize. We then
multiply by $e^{-in\varphi}$, and integrate in $\varphi$ to obtain

$$
\frac{1}{s} (s \rho U_n)' + \frac{n}{s} \rho V_n - n \Omega_v \frac{\rho Q_n}{c_s^2} = 0 \quad (B7)
$$

from the continuity equation, where ′ denotes d/ds; and the horizon-
tal momentum equations give

$$
\begin{equation}
\begin{split}
2\Omega_v V_n + n\Omega_v U_n ={}& -\frac{1}{2} \left[ \xi_+ \tilde{\eta}'_n + \frac{\xi_-}{2} (\tilde{\eta}'_{n-2} + \tilde{\eta}'_{n+2}) \right] \\
& + \frac{\xi_-}{4s} [(n-2)\tilde{\eta}_{n-2} - (n+2)\tilde{\eta}_{n+2}] \\
& + \frac{1}{2} \Omega \xi_- \chi (U_{n-2} - U_{n+2}) \\
& + \Omega \chi \left[ \xi_+ V_n + \frac{\xi_-}{2} (V_{n-2} + V_{n+2}) \right],
\end{split}
\tag{B8}
\end{equation}
$$

$$
\begin{equation}
\begin{split}
2\Omega_v U_n + n\Omega_v V_n ={}& \frac{\xi_-}{4} (\tilde{\eta}'_{n-2} - \tilde{\eta}'_{n+2}) + \frac{1}{2s} \left\{ n\xi_+ \tilde{\eta}_n \right. \\
& \left. - \frac{\xi_-}{2} [(n-2)\tilde{\eta}_{n-2} + (n+2)\tilde{\eta}_{n+2}] \right\} \\
& - \frac{1}{2} \Omega \xi_- \chi (V_{n-2} - V_{n+2}) \\
& + \Omega \chi \left[ \xi_+ U_n - \frac{\xi_-}{2} (U_{n-2} + U_{n+2}) \right],
\end{split}
\tag{B9}
\end{equation}
$$

where

$$
\tilde{\eta}_n \equiv Q_n + \Phi_n. \tag{B10}
$$
Adding and subtracting the horizontal momentum equations give the more compact form

$$ \frac{\xi_+}{2} \left( \tilde{\eta}'_n + \frac{n}{s} \tilde{\eta}_n \right) + \frac{\xi_-}{2} \left( \tilde{\eta}'_{n-2} - \frac{n-2}{s} \tilde{\eta}_{n-2} \right) = \left[ (2-n) \Omega_v - \Omega \chi \xi_+ \right] Y_n + \Omega \chi \xi_- X_{n-2}, \quad (B11) $$

$$ \frac{\xi_+}{2} \left( \tilde{\eta}'_n - \frac{n}{s} \tilde{\eta}_n \right) + \frac{\xi_-}{2} \left( \tilde{\eta}'_{n+2} + \frac{n+2}{s} \tilde{\eta}_{n+2} \right) = \left[ \Omega \chi \xi_+ - (n+2) \Omega_v \right] X_n - \Omega \chi \xi_- Y_{n+2}, \quad (B12) $$

where

$$ X_n \equiv U_n + V_n, \quad (B13) $$

$$ Y_n \equiv U_n - V_n. \quad (B14) $$

The linearized Poisson equation is

$$ \begin{aligned} & \frac{\xi_-}{4} \left[ \frac{d^2}{ds^2} - \frac{(2n-3)}{s} \frac{d}{ds} + \frac{n(n-2)}{s^2} \right] \Phi_{n-2} \\ & \qquad + \frac{\xi_+}{2} \left( \frac{d^2}{ds^2} + \frac{1}{s} \frac{d}{ds} - \frac{n^2}{s^2} \right) \Phi_n \\ & \qquad + \frac{\xi_-}{4} \left[ \frac{d^2}{ds^2} + \frac{(2n+3)}{s} \frac{d}{ds} + \frac{n(n+2)}{s^2} \right] \Phi_{n+2} \\ & = 4\pi \delta G \rho(s) \delta_{n0}. \end{aligned} \quad (B15) $$

Here $\delta_{ij}$ is the Kronecker delta.

## B1 The $n=2$ problem

To make analytic progress we truncate the series expansion (equation B6) at $n=2$. Then for $n=0$ we have

$$ \frac{\xi_+}{2} \tilde{\eta}'_0 = (\Omega \chi \xi_+ - 2 \Omega_v) V_0, \quad (B16) $$

$$ \begin{aligned} & \frac{\xi_+}{2} \tilde{\eta}'_0 + \frac{\xi_-}{2} \left( \tilde{\eta}'_2 + \frac{2}{s} \tilde{\eta}_2 \right) \\ & = (\Omega \chi \xi_+ - 2 \Omega_v) V_0 - \Omega \chi \xi_- Y_2, \end{aligned} \quad (B17) $$

where we have utilized the fact that $U_0 = 0$ as implied by the continuity equation (equation B7) for $n=0$. Next, for $n=2$ we have

$$ \frac{\xi_+}{2} \left( \tilde{\eta}'_2 + \frac{2}{s} \tilde{\eta}_2 \right) + \frac{\xi_-}{2} \tilde{\eta}'_0 = \Omega \chi \xi_- V_0 - \Omega \chi \xi_+ Y_2, \quad (B18) $$

$$ \frac{\xi_\pm}{2} (\tilde{\eta}'_2 - \frac{2}{s} \tilde{\eta}_2) = (\Omega \chi \xi_+ - 4 \Omega_v) X_2. \quad (B19) $$

Equations (B16)–(B19) imply

$$ \tilde{\eta}'_0 = V_0 = 0, \quad (B20) $$

$$ 2\Omega_v V_2 = \frac{\xi_+}{s} \tilde{\eta}_2 + (\Omega \chi \xi_+ - 2\Omega_v) U_2. \quad (B21) $$

Hence, an ODE for $\tilde{\eta}_2$ is

$$ \tilde{\eta}'_2 + \frac{1}{s} \left( 2 - \frac{\Omega \chi \xi_+}{\Omega_v} \right) \tilde{\eta}_2 + \Omega \chi \left( 4 - \frac{\Omega \chi \xi_+}{\Omega_v} \right) U_2 = 0. \quad (B22) $$

Next, the $n=2$ continuity equation, after using equations (B11) and (B21), becomes

$$ \begin{aligned} U_2' + \left[ \frac{\rho'}{\rho} + \frac{1}{s} \left( \frac{\Omega \chi \xi_+}{\Omega_v} - 1 \right) \right] U_2 + \left( \frac{\xi_+}{s^2 \Omega_v} - \frac{2 \Omega_v}{c_s^2} \right) \tilde{\eta}_2 + \frac{2 \Omega_v}{c_s^2} \Phi_2 = 0. \end{aligned} \quad (B23) $$

Further eliminating $U_2$ gives a second-order ODE for $\tilde{\eta}_2$

$$ \begin{aligned} \tilde{\eta}_2'' + \left( \frac{\rho'}{\rho} + \frac{1}{s} \right) \tilde{\eta}_2' + \left[ \frac{\lambda \rho'}{s \rho} + 2 \Omega \chi (\lambda + 2) \frac{\Omega_v}{c_s^2} - \frac{4}{s^2} \right] \tilde{\eta}_2 \\ = 2 \Omega \chi (\lambda + 2) \frac{\Omega_v}{c_s^2} \Phi_2, \end{aligned} \quad (B24) $$

where $\lambda \equiv 2 - \Omega \chi \xi_+ / \Omega_v$. We require $\tilde{\eta}_2'(0) = \tilde{\eta}_2(0) = 0$.

Equation (B24) describes the non-axisymmetric response of the GNG vortex to the potential $\Phi_2$. We are interested in the case where the potential arises from the vortex mass itself. Thus, $\Phi_2$ is given by a solution to the linearized Poisson solution.

## B2 Approximate solution as $s \to 0$

Here, we consider the approximate solution to equation (B24) for small $s$. We first obtain the potential solution $\Phi_2$. To do so in a tractable manner, we apply the ‘source term approximation’ to solve the Poisson equation as follows.

We first solve the $n=0$ cylindrical Poisson equation neglecting the $n=2$ term,

$$ \Phi_0'' + \frac{1}{s} \Phi_0' = \frac{2\Omega^2}{\xi_+ Q_{3D}} \exp(-s^2/2H_{\text{eff}}^2). \quad (B25) $$

The solution is

$$ \Phi_0' = \frac{2\Omega^2 H_{\text{eff}}^2}{\xi_+ Q_{3D}} \frac{1}{s} \left[ 1 - \exp\left(-\frac{s^2}{2H_{\text{eff}}^2}\right) \right]. \quad (B26) $$

We then use the above solution for $\Phi_0'$ as a source term in the $n=2$ Poisson equation

$$ \begin{aligned} \Phi_2'' + \frac{1}{s} \Phi_2 - \frac{4}{s^2} \Phi_2 = & -\frac{\xi_-}{2\xi_+} \left( \Phi_0'' - \frac{1}{s} \Phi_0' \right) \\ & \approx \frac{\xi_- \Omega^2}{4\xi_+^2 Q_{3D} H_{\text{eff}}^2} s^2 \equiv \mathcal{A} s^2, \end{aligned} \quad (B27) $$

where for the approximation we have expanded the solution for $\Phi_0'$ to first order in $s^2/2H_{\text{eff}}^2$. The approximate solution for $\Phi_2$ is then

$$ \Phi_2 \approx \frac{\mathcal{A}}{12} s^4. \quad (B28) $$

For small $s$ equation (B24), with the above solution for $\Phi_2$, is approximately

$$ \tilde{\eta}''_2 + \frac{1}{s} \tilde{\eta}'_2 - \frac{4}{s^2} \tilde{\eta}_2 = \mathcal{B} s^4, \quad (B29) $$

where

$$ \mathcal{B} \equiv \mathcal{A} \Omega \chi (\lambda + 2) \Omega_v / 6c_s^2. \quad (B30) $$

The solution is

$$ \tilde{\eta}_2 = \frac{\mathcal{B} s^6}{32}. \quad (B31) $$

The corresponding solution for $U_2$ and $V_2$ can be obtained from equations (B23) and (B21), respectively. We find

$$ U_2 = - \frac{\mathcal{A} \Omega_v}{192 c_s^2} (6 + \lambda) s^5, \quad (B32) $$
and

$$V_2 - U_2 = \frac{A\Omega_v}{48c_s^2}(2 + \lambda)s^5. \qquad (B33)$$

## B3 Vorticity perturbation

The perturbation to the vertical component of vorticity, $\delta\omega$, is defined in Cartesian coordinates as

$$\delta\omega \equiv \frac{\partial}{\partial x}\delta v_y - \frac{\partial}{\partial y}\delta v_x. \qquad (B34)$$

Expanding $\delta\omega$ similarly to equation (B6) and writing in terms of the perturbed elliptico-polar velocity components $U_n$, $V_n$ we have

$$\omega_n = \frac{\chi \xi_+}{2} V'_n + \frac{\chi \xi_+}{2s} (V_n + nU_n) \\ + \frac{\chi \xi_-}{4} \left( \frac{d}{ds} - \frac{n-1}{s} \right) (U_{n-2} + V_{n-2}) \\ - \frac{\chi \xi_-}{4} \left( \frac{d}{ds} + \frac{n+1}{s} \right) (U_{n+2} - V_{n+2}). \qquad (B35)$$

Then

$$\omega_0 = \frac{\chi \xi_-}{4} \left( \frac{d}{ds} + \frac{1}{s} \right) (V_2 - U_2). \qquad (B36)$$

Inserting the above solution for small $s$, equation (B33), we find

$$\omega_0 = \frac{\xi_+^2 \chi}{128 \xi_+^2 Q_{3D} F^2} (4\Omega_v - \Omega \chi \xi_+) \left(\frac{s}{H_{\text{eff}}}\right)^4. \qquad (B37)$$

Because of the second term in brackets, the sign of the vorticity perturbation depends on $\chi$. For $\chi \gtrsim 2.56$ the vorticity is made more negative by self-gravity and vice versa.

## APPENDIX C: VORTEX EVOLUTION IN 2D DISCS

We also simulated vortex evolution in 2D, razor-thin self-gravitating discs. These simulations were carried out in global cylindrical geometry ($R, \phi$) using the GENESIS code in a set-up similar to Pierens & Lin (in preparation). However, we adopt the same self-gravity parameter and vortex initialization as in our 3D simulations.

Fig. C1 compares the vortex evolution in discs with $Q_{3D} = 3, 4, 5$ and 100. Being 2D, the elliptic instability is suppressed and thus plays no role in these simulations. The $Q_{3D} = 100$ vortex decays due to numerical diffusion. This also occurs for the $Q_{3D} = 4, 5$, unlike in the corresponding 3D fiducial simulations where we find secular growth. The $Q_{3D} = 3$ case does show vortex growth, similar to the corresponding 3D simulation.

In Fig. C2 we plot the time evolution of vortex properties in the $Q_{3D} = 3, 2D$ disc. Similar to the 3D case, vortex growth is limited to $\lesssim H$ in its half-width. The global set-up here allows us to simulate beyond this stage without interference from boundary conditions. We find the vortex induces strong spiral shocks and is weakened, signified by the drop in $[Ro]$. The vortex is eventually destroyed by $t = 500P_0$, i.e. it does not readjust into a new configuration (Bodo et al. 2007). This self-limited growth is qualitatively similar to Mamatsashvili & Rice (2009). However, Mamatsashvili & Rice (2009) considered non-isothermal discs with $Q$ of the order of unity, which results in a gravito-turbulent state with multiple vortices, each lasting only two orbits; whereas our isolated vortex persists for $\sim 300$ orbits.

For the $Q_{3D} = 3, 2D$ disc we also experimented with different parameters, including the disc aspect ratio and flaring index, surface

<figure><img src="image_20.png" /><figcaption>Figure C1. Razor-thin, 2D global self-gravitating disc simulations of vortex evolution. Top: Rossby number; bottom: surface density.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure C2. Vortex evolution in the razor-thin, 2D disc with $Q_{3D} = 3$ shown in Fig. C1. Top to bottom: Toomre $Q$ at the vortex centre, Rossby number, aspect ratio, and radial size.</figcaption></figure>

density profile, and the initial vortex perturbation aspect ratio. The results are shown in Fig. C3. We find vortex growth is robust, unless the vortex is initially weak (left-hand panel). This supports the idea that a sufficiently strong vortex reduces the local gravitational stability of the disc.

Comparing these 2D simulations to our 3D cases indicates that for small $Q_{3D}$ (e.g. 3) EI turbulence is not necessary for the vortex to undergo self-gravitational growth. However, at larger $Q_{3D}$, self-gravity in itself is insufficient for vortex growth. Instead, we suggest
<figure><img src="image_22.png" /></figure>

**Figure C3.** Razor-thin, 2D global self-gravitating disc simulations of vortex evolution. We fix $Q_{3D} = 3$ but use different aspect ratios $h$ and flaring indices $f \equiv \partial \ln h / \partial \ln R$, surface density profiles $\sigma \equiv -\partial \ln \Sigma / \partial \ln R$, and the initial vortex aspect ratios $\chi$.

for $Q_{3D} \gtrsim 3$, EI turbulence can mediate vortex growth via secular gravitational instabilities, since for $Q_{3D} = 4, 5$ we only observe growth in the 3D simulations with EI, and not in the laminar, 2D discs.

This paper has been typeset from a TEX/LATEX file prepared by the author.