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

HAL
open science

# Gas and dust hydrodynamical simulations of massive lopsided transition discs – I. Gas distribution

Zhaohuan Zhu, Clément Baruteau

► To cite this version:

Zhaohuan Zhu, Clément Baruteau. Gas and dust hydrodynamical simulations of massive lopsided transition discs – I. Gas distribution. Monthly Notices of the Royal Astronomical Society, 2016, 458 (4), pp.3918-3926. 10.1093/mnras/stw202. hal-02369693

HAL Id: hal-02369693

https://hal.science/hal-02369693

Submitted on 9 Sep 2021

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Gas and dust hydrodynamical simulations of massive lopsided transition discs – I. Gas distribution

Zhaohuan Zhu<fnref n="1" />*† and Clément Baruteau<fnref n="2" />*

Accepted 2016 January 21. Received 2016 January 21; in original form 2015 August 4

## ABSTRACT

Motivated by lopsided structures observed in some massive transition discs, we have carried out 2D numerical simulations to study vortex structure in massive discs, including the effects of disc self-gravity and the indirect force which is due to the displacement of the central star from the barycentre of the system by the lopsided structure. When only the indirect force is included, we confirm the finding by Mittal & Chiang that the vortex becomes stronger and can be more than two pressure scale heights wide, as long as the disc-to-star mass ratio is $\ge 1$ per cent. Such wide vortices can excite strong density waves in the disc and therefore migrate inwards rapidly. However, when disc self-gravity is also considered in simulations, self-gravity plays a more prominent role on the vortex structure. We confirm that when the disc Toomre $Q$ parameter is smaller than $\pi/(2h)$, where $h$ is the disc’s aspect ratio, the vortices are significantly weakened and their inward migration slows down dramatically. Most importantly, when the disc is massive enough (e.g. $Q \sim 3$), we find that the lopsided gas structure orbits around the star at a speed significantly slower than the local Keplerian speed. This sub-Keplerian pattern speed can lead to the concentration of dust particles at a radius beyond the lopsided gas structure (as shown in Paper II). Overall, disc self-gravity regulates the vortex structure in massive discs and the radial shift between the gas and dust distributions in vortices within massive discs may be probed by future observations.

**Key words:** hydrodynamics – instabilities – protoplanetary discs – stars: protostars.

## 1 INTRODUCTION

Transition discs are protoplanetary discs whose inner regions have undergone substantial clearing (see the review by Espaillat et al. 2014). Recent submm interferometric observations have suggested significant non-axisymmetric features in these discs (HD 142527 from Casassus et al. 2013; Oph IRS 48 from van der Marel et al. 2013; LkHα 330 from Isella et al. 2013). In the extreme case of Oph IRS 48, there is a highly asymmetric crescent-shaped dust structure between 45 and 80 au from the star. The peak emission from this dust structure is at least 130 times stronger than the upper limit of the opposite side of the disc. Observations at longer wavelengths (e.g. 8.8 mm with Australia Telescope Compact Array for HD 142527; Casassus et al. 2015) reveal an even more compact structure. These observations suggest that some dust trapping mechanism is operating in the azimuthal direction of the disc.

Theoretically, it has been known that anticyclonic vortices can be long lived (Godon & Livio 2000) and can efficiently trap

dust particles (Adams & Watkins 1995; Barge & Sommeria 1995; Tanga et al. 1996; Chavanis 2000; Johansen, Andersen & Bröndenburg 2004; Lyra et al. 2009a; Meheut et al. 2012; Zhu et al. 2014). Using the vortex gas structure derived by Kida (1981) and Goodman, Narayan & Goldreich (1987), particle distribution within the vortex has been derived and compared with observed asymmetric disc structures (Lyra & Lin 2013). By constructing realistic three-dimensional magnetohydrodynamics (MHD) simulations including dust particles, Zhu & Stone (2014) have found that dust trapping vortices can reproduce Atacama Large Millimeter/submillimeter Array (ALMA) observations reasonably well.

However, some transition discs are quite massive. For an example, the disc-to-star mass ratio is ~4 per cent in the HD 142527 system. Recently, Mittal & Chiang (2015) have suggested a different mechanism for azimuthal particle trapping in massive discs by allowing the star to move around the barycentre of the system due to the gravitational pull of the massive lopsided structure. They suggest that, if the azimuthal structure in these discs is massive enough, its gravitational force to the central star can displace the star from the barycentre of the system. Such displacement causes an $m = 1$ mode indirect force to the disc, which can lead to an asymmetric horseshoe flow pattern. This pattern can be self-sustaining as long as it can lead to enough gravitational force required to offset the

<footnote n="1">¹Department of Astrophysical Sciences, 4 Ivy Lane, Peyton Hall, Princeton University, Princeton, NJ 08544, USA</footnote>

<footnote n="2">²CNRS/Institut de Recherche en Astrophysique et Planétologie, 14 avenue Edouard Belin, F-31400 Toulouse, France</footnote>

<footnote>* E-mail: zhzhu@astro.princeton.edu (ZZ); clement.baruteau@irap.omp.eu (CB)
† Hubble fellow.</footnote>
star. However, Mittal & Chiang (2015) applied a prescribed indirect force to the disc, and it is unclear if this mode can sustain with the indirect force self-consistently calculated from the disc asymmetric structure itself.

Furthermore, the disc self-gravity has been ignored in Mittal & Chiang (2015). If the disc is so massive that the indirect force can become important, the disc self-gravity should play an even more important role in shaping the asymmetric structure. Suppose the lopsided pattern at a distance *r* from the star has a mass *m*<sub>LP</sub>, the acceleration of the central star is on the order of ~*Gm*<sub>LP</sub>*r*<sup>2</sup> and thus the indirect force felt by unit mass in the disc is ~*Gm*<sub>LP</sub>*r*<sup>2</sup>. On the other hand, the force due to the self-gravity of the asymmetric pattern is on the order of *Gm*<sub>LP</sub>*Δr*<sup>2</sup>, where *Δr* is the scale of the asymmetric pattern. Since *Δr* ~ ~*r*, the disc self-gravitating force should be larger than the indirect force. Previous simulations by Lin (2012) have shown that self-gravity can suppress large-scale vortices produced by Rossby wave instability (RWI; Lovelace et al. 1999). Linear instability analysis by Lovelace & Hohlfeld (2013) has also suggested that disc self-gravity can stabilize the *m* = 1 RWI mode when the disc's Toomre *Q* parameter *Q* < (π/2)(*r*/*H*), where *H* is the disc's scale height. On the other hand, the indirect force is a large-scale *m* = 1 driving force, which can facilitate *m* = 1 large-scale structure. Thus, a self-consistent model including both the indirect force and disc self-gravity is necessary.

In this paper, we present a self-consistent model including both the indirect force and disc self-gravity, focusing on their effects on the gas structure of the vortex. In the subsequent paper (Baruteau & Zhu 2016), we will study how these effects impact particle concentration in vortices. In Section 2, we first follow Mittal & Chiang (2015) using test particle method to study asymmetric structures in pressureless non-self-gravitating discs. Then we introduce our hydrodynamical simulations including both disc pressure and self-gravity in Section 3. Results are presented in Section 4 and summarized in Section 5.

## 2 'FAST MODES' IN PRESSURELESS NON-SELF-GRAVITATING FLUIDS

Although the horseshoe solution led by the indirect force (Mittal & Chiang 2015) seems to be distinct from the traditional vortex solution (Goodman et al. 1987), both solutions are steady in a Keplerian rotating reference frame. Thus, they are both classified as 'fast modes', in contrast with the 'slow modes' (Tremaine 2001; Lin 2015) that is almost steady in an inertial reference frame. In this section, we present the traditional vortex solution and the new horseshoe solution together in the pressureless fluid.

In a frame which is centred on the star and rotates at the angular frequency $\Omega_f$, the Euler equation for the disc becomes

$$ \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla)\mathbf{v} = -\frac{1}{\rho} \nabla P - \nabla \Phi_* - \nabla \Phi_{\text{ind}} - \nabla \Phi_{\text{sg}} \\ - \Omega_f \times (\Omega_f \times \mathbf{r}) - 2\Omega_f \times \mathbf{v}, \quad (1) $$

where $\Phi_*$, $\Phi_{\text{ind}}$, and $\Phi_{\text{sg}}$ are the potential due to the direct gravitational force from the star, the indirect force due to the acceleration of the reference frame (or the star), and the disc self-gravity. Specifically $\Phi_{\text{ind}} = -a_c \cdot \mathbf{r}$, where $a_c$ is the acceleration of the central star. The last two terms in equation (1) are the centrifugal force and the Coriolis force. The pressure term, indirect force term, and self-gravity term are all related to the disc density distribution ($\rho$) which has to be solved with the continuity equation.

For a 'fast mode' which is steady in the rotating frame, we have $\partial\mathbf{v}/\partial t = 0$ so that the fluid streamlines and trajectories coincide. We can then use the trajectories of test particles to study the flow streamlines. Such test particle method has been widely used in studying discs in binary systems (Paczynski 1977) and galactic dynamics (Contopoulos 1979; Binney et al. 1991). The position and velocity of the test particle ($\mathbf{r}$, $\mathbf{u}$) are

$$ \frac{d\mathbf{r}}{dt} = \mathbf{u}, \quad (2) $$

$$ \frac{d\mathbf{u}}{dt} = -\frac{1}{\rho}\nabla P - \nabla\Phi_* - \nabla\Phi_{\text{ind}} - \nabla\Phi_{\text{sg}} \\ - \Omega_f \times (\Omega_f \times \mathbf{r}) - 2\Omega_f \times \mathbf{u}. \quad (3) $$

If the indirect force term is known beforehand and the density related terms (pressure and self-gravity) are small compared to other terms, we can ignore these density related terms to simplify the problem further and search all possible 'fast modes' as in Mittal & Chiang (2015). Equations (2) and (3) can be solved given initial position and velocity for a test particle. Although we can have infinite trajectories, for a given flow only trajectories that do not intercept each other or itself represent allowed streamlines. For a steady flow, the streamlines also have to be closed, unless there is external mass flowing in. Thus, in the rest of this section, we will look for closed non-intercepting trajectories through numerical integration of equations (2) and (3) using the fixed time step fifth-order Runge-Kutta method.

To search these possible 'fast mode' flow patterns, we choose $\Omega_p = \Omega_f = 1$, where $\Omega_p$ is the pattern speed. In this case, the pattern is steady in the rotating frame. The lopsided pattern is assumed to be symmetric to the x-axis so that the indirect force is along the x-axis. The indirect force is opposite to the direction of the acceleration of the central star, and it is a constant force anywhere in the disc. The prescribed indirect force is determined from hydrodynamical simulations in Section 3. In the unit that $GM_* = 1$, the indirect force is $-5.625 \times 10^{-3} \hat{x}$ (Section 3). Test particles are launched at the x-axis with $v_x = 0$. Since the force is symmetric to the x-axis and the particle is initialized with only $v_y$ at the x-axis, this symmetry implies that the particle trajectory is a non-intercepting closed orbit only if $v_x = 0$ when the particle crosses the x-axis for the first time after it has been launched. By varying $v_y$, we search these closed non-intercepting orbits.

We have searched closed non-intercepting orbits of particles with and without including indirect force due to the large-scale lopsided pattern. Without the indirect force, the most simple 'fast mode' is the traditional pressureless elliptical vortex having an aspect ratio($\chi$) of 2. This is shown in the left-hand panel of Fig. 1. Fluid streamlines in this mode are the trajectories of test particles undergoing epicyclic motion in the rotating frame. The guiding centre approximation suggests that, in the rotating frame whose angular frequency ($\Omega_f$) equals the test particles' mean motion, eccentric particles undergo anticyclonic epicyclic motions and form ellipses with aspect ratios of 2. All test particles in the left-hand panel of Fig. 1 have their mean motion equal to $\Omega_f$ but different eccentricities. This 'fast mode' is the only 'fast mode' when there is no indirect force.<fnref n="1" /> This mode

<footnote n="1">¹ In a frame rotating at $\Omega_f$, an eccentric particle must have its mean motion equal to $\Omega_f$ or in first-order mean motion resonance with $\Omega_f$ to form a non-intercepting closed orbit. However, streamlines cannot consist of different first-order mean motion resonant orbits (e.g. 2:1, 3:2, 4:3,...) since these orbits are discretized in space while streamlines should be continuous in</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Two closed non-crossing orbit configurations when the centre-of-mass is not at the star. Left-hand panel: the eccentric orbits forms a vortex pattern with an aspect ratio of 2. The solid curves have considered the indirect force while the dotted curves have not considered the indirect force. Clearly, the indirect force has small effect on these eccentric orbits. Right-hand panel: similar to Mittal & Chiang (2015) fig. 1, horseshoe shaped orbits form due to the indirect force. In both panels, the orbits are initiated at the same x positions. The lopsided pattern in the right-hand panel is much more elongated than the one in the left-hand panel.</figcaption></figure>

can also be derived from the traditional compressible vortex solution under the pressureless limit (Section 3).

When the indirect force due to the displacement of the star from the barycentre of the system is included, two classes of closed non-intercepting smooth orbits with dramatically different velocities have been found. The first class of orbits forms the $\chi \sim 2$ vortex solution above, with a slight modulation by the indirect force. In this mode, the test particles only need one orbital time $(2\pi/\Omega_f)$ to finish one circle around the vortex centre, and the acceleration due to the addition of the star's gravity and the inertial force is on the order of $\sim eGM/r^2$, much larger than the indirect force. Thus, the indirect force has very little effect on this class of orbits (red curves in Fig. 1). The aspect ratios of the ellipses are still $\sim 2$.

Another class of closed non-intercepting orbits in the rotating frame is from slowly moving test particles (the right-hand panel of Fig. 1),<fnref n="2" /> as originally shown by Mittal & Chiang (2015). Different from the vortex solution above where test particles undergo rapid epicyclic motion, test particles in this mode undergo slow horseshoe orbits induced by the indirect force. The libration frequency for the particle is around $\sqrt{3}\mu\Omega_f$ (Mittal & Chiang 2015), where $\mu$ is the ratio between the distance from the barycentre to the star and is equal to 0.0056 in our setup. Test particles in the right-hand panel of Fig. 1 thus need $\sim 10$ orbital time to finish one circle around the vortex centre. The motion of these test particles is so slow in the rotating frame that inertial forces almost balance the star's gravity, and the indirect force plays a crucial role in determining the orbits. During one horseshoe orbit, the test particle loses and gains angular momentum from the indirect force. In the limit that the indirect force becomes zero, these horseshoe orbits gradually become circular orbits around the central star.

space. Streamlines can neither consist of orbits in the same first-order mean motion resonance (e.g. 2:1) but having different eccentricities, since these orbits will intercept each other. Thus, the only possible fluid streamlines consist of trajectories of test particles having their mean motion equal to $\Omega_f$ but different eccentricities, as shown in the left-hand panel of Fig. 1.

## 3 HYDRODYNAMICAL SIMULATIONS

Although the two 'fast modes' in pressureless fluids have distinct aspect ratios, they become less distinct when the gas pressure is considered. The compressible vortex solution by Goodman et al. (1987, GNG solution) suggests that, even without either the indirect force or disc self-gravity, the gas pressure alone can lead to vortices with any aspect ratio larger than 2. In the special case that the vortex aspect ratio ($\chi$) equals 2, the GNG solution becomes the pressureless vortex solution with streamlines as those shown by the test particle method in the left-hand panel of Fig. 1. For a vortex with $\chi$ larger than 2, the pressure force starts to balance the tidal force, Coriolis force, and centrifugal force. When $\chi$ becomes larger, the net force which leads to vortex rotation becomes smaller, so that the vortex rotates slower and is more elongated. Eventually when $\chi$ keeps increasing, this net force can be so small that any additional force (such as the indirect force and self-gravity) starts to affect the vortex structure.

To study how the indirect force and disc self-gravity affect the vortex solution in compressible fluids, we have carried out two-dimensional hydrodynamical simulations including both of these forces.

### 3.1 Set-up

The hydrodynamical code we used is FARGO-ADSG (Baruteau & Masset 2008a,b) which is built on FARGO (Masset 2000) but with optional self-gravity and energy equation. Since the numerical grid is centred at the central star, the velocity equations in a rotating frame are

$$ \frac{\partial v_r}{\partial t} + v_r \frac{\partial v_r}{\partial r} + v_\phi \frac{\partial v_r}{r \partial \phi} - \frac{v_\phi^2}{r} = -\frac{\partial P}{\rho \partial r} - \frac{\partial \Phi_*}{\partial r} - \frac{\partial \Phi_{\text{ind}}}{\partial r} - \frac{\partial \Phi_{\text{sg}}}{\partial r} + \Omega_f^2 r + 2\Omega_f v_\phi, \quad (4) $$

$$ \frac{\partial v_\phi}{\partial t} + v_r \frac{\partial v_\phi}{\partial r} + v_\phi \frac{\partial v_\phi}{r \partial \phi} - \frac{v_\phi v_r}{r} = -\frac{1}{\rho r} \frac{\partial P}{\partial \phi} - \frac{1}{r} \frac{\partial \Phi_*}{\partial \phi} - \frac{1}{r} \frac{\partial \Phi_{\text{ind}}}{\partial \phi} - \frac{\partial \Phi_{\text{sg}}}{\partial \phi} - 2\Omega_f v_r, \quad (5) $$

where the indirect potential $\Phi_{\text{ind}} = -a_c \cdot r$ and $a_c$ is the acceleration of the central star due to the disc's gravity. It is calculated by integrating the gravitational acceleration from each grid cell on to the central star. Indirect forces are calculated at each time step and added to the equation of motion. The simulations are run in the inertial frame centred at the star, but analyses shown later are done in the rotating frame with $\Omega_f$ equal to the pattern speed of the asymmetric disc structure so that the time derivatives in these equations become almost zero.

To make each grid cell have the same length in both the radial and azimuthal directions, our grids are uniformly spaced in $\log r$ from $r_{\text{in}}$ to $r_{\text{out}}$, where $r_{\text{in}} = 0.1r_0$ and $r_{\text{out}} = 10r_0$. Our standard simulations have 752 grids in the radial direction, and 1024 grids in the azimuthal direction. When disc self-gravity is included, a softening length of 0.3 disc scale height is used in the self-gravity potential to mimic the effect of a finite disc thickness. Outflow boundary condition is used at both inner and outer boundaries.

To generate a vortex in the simulation, we initialize the disc surface density with a density bump in the radial direction which will later break into vortices through Papaloizou-Pringle instability

<footnote n="2">² Starting from some initial positions, we can find multiple smooth closed orbits while at some positions only roughly smooth closed orbits can be found.</footnote>
<table><caption>**Table 1.** Simulations.</caption><thead><tr><th>Run</th><th>Σ<sub>0</sub></th><th>Self-gravity</th><th>Indirect force</th><th>R<sub>c</sub><sup>a</sup></th><th>Σ<sub>max</sub>(R = R<sub>c</sub>)<sup>b</sup></th><th>Σ<sub>min</sub>(R = R<sub>c</sub>)<sup>b</sup></th><th>ΔR/R<sub>c</sub><sup>c</sup></th></tr></thead><tbody><tr><td>g5</td><td>5 × 10<sup>-3</sup></td><td>No</td><td>No</td><td>0.89</td><td>8.8 × 10<sup>-3</sup></td><td>3.8 × 10<sup>-3</sup></td><td>0.23</td></tr><tr><td>g5i</td><td>5 × 10<sup>-3</sup></td><td>No</td><td>Yes</td><td>0.78</td><td>1.3 × 10<sup>-2</sup></td><td>1.8 × 10<sup>-3</sup></td><td>0.30</td></tr><tr><td>g5g</td><td>5 × 10<sup>-3</sup></td><td>Yes</td><td>No</td><td>0.98</td><td>6.5 × 10<sup>-3</sup></td><td>5.0 × 10<sup>-3</sup></td><td>0.12</td></tr><tr><td>g5gi</td><td>5 × 10<sup>-3</sup></td><td>Yes</td><td>Yes</td><td>1.03</td><td>6.5 × 10<sup>-3</sup></td><td>4.8 × 10<sup>-3</sup></td><td>0.15</td></tr><tr><td>g5giH</td><td>5 × 10<sup>-3</sup></td><td>Yes</td><td>Yes</td><td>1.02</td><td>6.6 × 10<sup>-3</sup></td><td>4.8 × 10<sup>-3</sup></td><td>0.15</td></tr><tr><td>g0p2i</td><td>2 × 10<sup>-4</sup></td><td>No</td><td>Yes</td><td>0.88</td><td>3.7 × 10<sup>-4</sup></td><td>1.4 × 10<sup>-4</sup></td><td>0.24</td></tr><tr><td>g0p2g</td><td>2 × 10<sup>-4</sup></td><td>Yes</td><td>No</td><td>0.92</td><td>3.4 × 10<sup>-4</sup></td><td>1.6 × 10<sup>-4</sup></td><td>0.20</td></tr><tr><td>g0p2gi</td><td>2 × 10<sup>-4</sup></td><td>Yes</td><td>Yes</td><td>0.90</td><td>3.4 × 10<sup>-4</sup></td><td>1.6 × 10<sup>-4</sup></td><td>0.21</td></tr><tr><td>g0p5i</td><td>5 × 10<sup>-4</sup></td><td>No</td><td>Yes</td><td>0.88</td><td>9.1 × 10<sup>-4</sup></td><td>3.8 × 10<sup>-4</sup></td><td>0.22</td></tr><tr><td>g0p5g</td><td>5 × 10<sup>-4</sup></td><td>Yes</td><td>No</td><td>0.93</td><td>8.3 × 10<sup>-4</sup></td><td>4.3 × 10<sup>-4</sup></td><td>0.20</td></tr><tr><td>g0p5gi</td><td>5 × 10<sup>-4</sup></td><td>Yes</td><td>Yes</td><td>0.92</td><td>8.3 × 10<sup>-4</sup></td><td>4.1 × 10<sup>-4</sup></td><td>0.20</td></tr><tr><td>g1i</td><td>1 × 10<sup>-3</sup></td><td>No</td><td>Yes</td><td>0.87</td><td>1.9 × 10<sup>-3</sup></td><td>6.6 × 10<sup>-4</sup></td><td>0.25</td></tr><tr><td>g1g</td><td>1 × 10<sup>-3</sup></td><td>Yes</td><td>No</td><td>0.93</td><td>1.6 × 10<sup>-3</sup></td><td>8.9 × 10<sup>-4</sup></td><td>0.18</td></tr><tr><td>g1gi</td><td>1 × 10<sup>-3</sup></td><td>Yes</td><td>Yes</td><td>0.94</td><td>1.6 × 10<sup>-3</sup></td><td>8.6 × 10<sup>-4</sup></td><td>0.19</td></tr><tr><td>g1giH</td><td>1 × 10<sup>-3</sup></td><td>Yes</td><td>Yes</td><td>0.95</td><td>1.6 × 10<sup>-3</sup></td><td>8.6 × 10<sup>-4</sup></td><td>0.19</td></tr><tr><td>g2i</td><td>2 × 10<sup>-3</sup></td><td>No</td><td>Yes</td><td>0.86</td><td>3.8 × 10<sup>-3</sup></td><td>1.3 × 10<sup>-3</sup></td><td>0.25</td></tr><tr><td>g2g</td><td>2 × 10<sup>-3</sup></td><td>Yes</td><td>No</td><td>0.94</td><td>3.0 × 10<sup>-3</sup></td><td>1.9 × 10<sup>-3</sup></td><td>0.17</td></tr><tr><td>g2gi</td><td>2 × 10<sup>-3</sup></td><td>Yes</td><td>Yes</td><td>0.96</td><td>3.1 × 10<sup>-3</sup></td><td>1.8 × 10<sup>-3</sup></td><td>0.17</td></tr><tr><td>g10i</td><td>10 × 10<sup>-3</sup></td><td>No</td><td>Yes</td><td>0.51</td><td>5.8 × 10<sup>-2</sup></td><td>4.3 × 10<sup>-3</sup></td><td>0.26</td></tr><tr><td>g10g</td><td>10 × 10<sup>-3</sup></td><td>Yes</td><td>No</td><td>1.02</td><td>1.1 × 10<sup>-2</sup></td><td>9.4 × 10<sup>-3</sup></td><td>0.15</td></tr><tr><td>g10gi</td><td>10 × 10<sup>-3</sup></td><td>Yes</td><td>Yes</td><td>1.03</td><td>1.2 × 10<sup>-2</sup></td><td>9.5 × 10<sup>-3</sup></td><td>0.15</td></tr><tr><td>g10giH</td><td>10 × 10<sup>-3</sup></td><td>Yes</td><td>Yes</td><td>1.02</td><td>1.1 × 10<sup>-2</sup></td><td>9.6 × 10<sup>-3</sup></td><td>0.15</td></tr><tr><td colspan="8"><strong>Various disc α</strong></td></tr><tr><td>g5(α = 10<sup>-5</sup>)</td><td>5 × 10<sup>-3</sup></td><td>No</td><td>No</td><td>0.90</td><td>8.6 × 10<sup>-3</sup></td><td>3.8 × 10<sup>-3</sup></td><td>0.22</td></tr><tr><td>g5(α = 10<sup>-4</sup>)</td><td>5 × 10<sup>-3</sup></td><td>No</td><td>No</td><td>1.00</td><td>5.1 × 10<sup>-3</sup></td><td>4.8 × 10<sup>-3</sup></td><td>0.09</td></tr><tr><td>g5(α = 10<sup>-3</sup>)</td><td>5 × 10<sup>-3</sup></td><td>No</td><td>No</td><td>0.88</td><td>3.7 × 10<sup>-3</sup></td><td>3.7 × 10<sup>-3</sup></td><td>—</td></tr></tbody></table>

<sup>a</sup>R<sub>c</sub> is the radial position where the maximum surface density lies at 150 orbits.

<sup>b</sup>Σ<sub>max</sub>(R = R<sub>c</sub>) and Σ<sub>min</sub>(R = R<sub>c</sub>) are the maximum and minimum density along the azimuthal direction at R = R<sub>c</sub> at 150 orbits.

<sup>c</sup>ΔR is the width of the vortex in the radial direction. It is the distance between two radial positions where Σ = (Σ<sub>max</sub> + Σ<sub>min</sub>)/2 at the azimuthal angle of the vortex centre.

(Papalozou & Pringle 1984, 1985) or RWI (Lovelace et al. 1999; Li et al. 2000, 2001). The initial density profile is

$$ \Sigma(r) = \Sigma_0 \left\{ \zeta + \exp \left[ - \frac{(r - r_0)^2}{2\sigma^2} \right] \right\}, \quad (6) $$

where ζ = 0.01 and r<sub>0</sub> = 1. Lyra et al. (2009b) and Regály et al. (2012) have showed that RWI can be excited in α discs only if the density jump is sharp enough with jump width less than about 2 pressure scale-height (H). H is defined as c<sub>s</sub>/Ω<sub>K</sub>, where c<sub>s</sub> is the disc sound speed, and Ω<sub>K</sub> is the orbital frequency. Thus, we choose σ = 2H(r<sub>0</sub>). To trigger m = 1 mode in the instability, we add a small perturbation to the gas surface density as in Heemskerk, Papalozou & Savonije (1992):

$$ \Sigma(r, \phi) = \Sigma(r) \left[ 1 + 10^{-3} \cos(\phi) \sin \left( \pi \frac{r - r_{\text{in}}}{r_{\text{out}} - r_{\text{in}}} \right) \right]. \quad (7) $$

We have run the simulations for 500 orbits, and in this paper we refer one orbit as the orbital time (2π/Ω) at r<sub>0</sub>.

We use the isothermal equation of state so that P = Σc<sub>s</sub><sup>2</sup>. The choice of an isothermal equation of state is for the sake of simplicity, but we note that the inclusion of an energy equation impacts the growth and saturation phases of the RWI (see e.g. Les & Lin 2015, where the RWI is induced by a gap-opening planet). Since the whole disc has the same temperature, the disc aspect ratio h ≡ H/r increases as √r. We take h = 0.1 at r = r<sub>0</sub>, which is typical for a protoplanetary disc at 10s of au. A small viscosity with α = 10<sup>-6</sup> is included in the simulation. We find that adding this small viscosity significantly improves the convergence of the simulation. To explore

the effect of α on the vortex structure, we have also increased α to 10<sup>-3</sup> for our fiducial run.

To study the effect of indirect force and self-gravity on the vortex structure, we choose discs which are relatively massive but still gravitationally stable. Our fiducial run assumes Σ<sub>0</sub> = 5 × 10<sup>-3</sup> (models g5 in Table 1) which has a total mass of 0.017 M<sub>☉</sub>, and the Toomre parameter, Q ≡ c<sub>s</sub>Ω/(πGΣ), is ~6 at r = r<sub>0</sub>. We also vary Σ<sub>0</sub> from 2 × 10<sup>-4</sup> to 0.01 (g0p2 to g10), and the corresponding Q varies from 160 to 3. The effects of the indirect force and self-gravity become more important in more massive discs. For runs which only include the indirect force, we add ‘i’ at the end of their model names. For runs which only include disc self-gravity, we add ‘g’ at the end of their names. For runs including both the indirect force and self-gravity, we add ‘gi’ at the end of their names. For three high-resolution runs having 1536 radial grids and 2048 azimuthal grids, the model names are ended with ‘H’. All the runs are summarized in Table 1.

## 4 RESULTS

After the simulations start, an m = 1 mode grows exponentially, and in most runs it saturates at ~40 orbits. For runs where self-gravity is important (e.g. g10g), it saturates at a slightly longer time (e.g. ~60 orbits). With our fiducial disc mass (g5), the disc surface density for runs having either the indirect force or disc self-gravity or both are shown in Fig. 2. In Fig. 2 and all figures below, we have rotate the images so that the maximum gas surface density is always at φ = π. Table 1 gives the radial position where the maximum surface density lies at 150 orbits, labelled as R<sub>c</sub>. Since R<sub>c</sub> = 1 initially, the
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Upper panels: the disc surface density contours at 150 orbits for case g5, g5i, g5g, and g5gi in the *r*–*φ* plane. Middle and bottom panels: the disc surface density contours for the same cases at 150 (middle panels) and 400 (bottom panels) orbits in the *x*–*y* plane.</figcaption></figure>

deviation of $R_c$ from 1 reflects the radial migration speed of the vortex. At 150 orbits, the maximum and minimum density along the azimuthal direction at $R = R_c$ are also given in Table 1, labelled as $\Sigma_{\max}$ and $\Sigma_{\min}$. The ratio between $\Sigma_{\max}$ and $\Sigma_{\min}$ represents the strength of the vortex. The vortex radial width based on the vortex density structure is also given in Table 1.

Fig. 2 clearly shows that the indirect force strengthens the vortex. Compared with case g5 (the leftmost panels), the vortex in case g5i (the second panel from the left) is more roundish with a radial width larger than 2 disc scale height, and has a higher contrast between the vortex centre and the background (also shown in Table 1). This stronger vortex also excites stronger density waves. These density waves carry the angular momentum of the vortex away and enable vortex migration (Paardekooper, Lesur & Papaloizou 2010). With strong density waves excited, the vortex in g5i migrates inward fast. The vortex in case g5i migrates from $r=1$ to $0.5$ within 400 orbits.

On the other hand, disc self-gravity stabilizes the vortex. With only self-gravity included (case g5g, the third panel from the left), the saturated $m = 1$ mode breaks into two vortices at ~100 orbits and these two vortices are found to remain separated even at 400 orbits. This is consistent with Goodman & Narayan (1988), Lin (2012), and Yellin-Bergovoy, Heifetz & Umurhan (2015) that self-gravity inhibits RWI, especially for low $m$ modes. Quantitatively, linear instability analysis by Lovelace & Hohlfeld (2013) has suggested that disc self-gravity can suppress RWI modes with

$m < (\pi/2)(r/H)Q^{-1}$. The suppression of $m = 1$ mode in the case g5g whose $(\pi/2)(r/H)Q^{-1} = 2.5$ is consistent with this criterion.

However, when both disc self-gravity and the indirect force are included (case g5gi, the rightmost panel), the $m = 1$ mode persists in the simulation. The shift of power from the $m = 2$ mode to the $m = 1$ mode when the star is allowed to move has first been seen in numerical simulations by Christodoulou & Narayan (1992). Although disc self-gravity tries to stabilize the vortex and inhibits the $m = 1$ mode, the indirect force tries to maintain the $m = 1$ mode. Eventually, the vortex in case g5gi is weaker than the vortex in case g5, and it migrates considerably slower than the vortex in either g5 or g5i. Detailed analyses suggest that it migrates less than 0.5 per cent in the radial direction within 500 orbits.

Both the indirect force and disc self-gravity become more important for shaping the vortex structure when the disc becomes more massive. Fig. 3 shows the surface density of discs having different masses (increasing mass from left- to right-hand panels) if the indirect force or self-gravity or both are included (from upper to lower panels). For the smallest disc mass (g0p2), both disc self-gravity and the indirect force are weak, and they have negligible effects on the vortex structure. With the disc surface density increasing, the indirect force starts to strengthen the vortex and the vortex migrates faster in a more massive disc (upper panels). On the other hand, disc self-gravity starts to weaken the vortex (the middle row of Fig. 3) in massive discs. As mentioned above, RWI modes
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Upper panels: the disc surface density contours at 150 orbits for discs having different masses (from left- to right-hand panels) if the indirect force or self-gravity or both are included (from upper to lower panels). From left- to right-hand panels, $(\pi/2)(r/H)Q^{-1} = 0.1, 1, 2.5, 5$, and $m = 1$ mode is suppressed when $(\pi/2)(r/H)Q^{-1} > 1$.</figcaption></figure>

with $m < (\pi/2)(r/H)Q^{-1}$ (Lovelace & Hohlfeld 2013) will be affected by disc self-gravity. From left- to right-hand panels in Fig. 3, $(\pi/2)(r/H)Q^{-1} = 0.1, 1, 2.5, 5$, respectively. As expected, $m = 1$ mode is completely suppressed in the rightmost two cases, and only weak $m = 2$ perturbations have been observed. In the most massive case, even the $m = 2$ perturbation is significantly suppressed. When both the indirect force and self-gravity are included (lower panels), these two effects seem to be competing with each other, but eventually the effect of self-gravity seems to be dominant and the vortices become weaker than those in the top panels. The indirect force which has $m = 1$ symmetry also leaves its imprint by maintaining the $m = 1$ disc asymmetry. The middle and bottom panels suggest that, as long as disc self-gravity is included, the vortex migrates slowly in the disc.

To quantitatively compare the asymmetric structure in simulations with different masses, cuts of the disc’s surface density across the vortex centre along the *r* and *φ* directions are shown in Fig. 4. The peak density in the radial profile is at a larger position for a more massive disc, suggesting that the vortex migrates inward slower (or even migrates slightly outward) in a more massive disc. The density profile in the azimuthal direction clearly shows that the vortex becomes weaker in a more massive disc due to the stabilizing effect of disc self-gravity. The density ratio between the vortex centre and the background is ~2 for the least massive disc while it is ~1.2 for the most massive case. The similarity between low- and high-resolution runs suggests that the simulations are numerically converged. There is a slight difference for the azimuthal density

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Density cuts across the vortex centre along the *r* (the left-hand panel) and *φ* (the right-hand panel) directions for g0p2gi, g0p5gi, g1gi, g2gi, g5gi, and g10gi (from heavy to light curves). The dotted curves are from high-resolution runs g1giH, g5giH, and g10giH.</figcaption></figure>

profile for the most massive cases between g10gi and g10giH. In g10giH, some small vortices are present which do not merge with the big vortex at 150 orbits. However, the big vortex has similar structures as that in the lower resolution run.

Since the vortex structure sensitively depends on the disc viscosity (de Val-Borro et al. 2007; Fu et al. 2014; Zhu & Stone 2014), we have carried out three additional simulations with $\alpha = 10^{-5}$, $10^{-4}$, and $10^{-3}$ using the same set-up as our fiducial case g5 which
<figure><img src="image_7.png" /><figcaption>**Figure 5.** The upper left-hand panel: the streamlines in the rotating frame for horseshoe orbits in the pressureless fluid (the same as the right-hand panel of Fig. 1). Upper middle and right-hand panels: the density contours and streamlines in the frame corotating with the lopsided structure for g5, g5i, and g10gi at 150 orbits. Lower panels: the difference between the azimuthal velocity of the fluid element at $\phi = \pi$ and the local Keplerian velocity.</figcaption></figure>

has $\alpha = 10^{-6}$. Table 1 shows that in the case with $\alpha = 10^{-5}$ the vortex has a similar strength as our fiducial case, but it becomes significantly weaker in cases with bigger $\alpha$. Intuitively, we would expect that the vortex cannot be generated when the growth time-scale of RWI is comparable with the viscous spreading time-scale of the density bump ($\Delta R^2/\nu$). In our simulations, the width of the density bump is $\Delta R = 0.2H$, and the vortex saturates at around 50 orbits. Thus, $\alpha = 0.01$ should suppress the vortex formation. In our simulations, a comparable but smaller $\alpha(\alpha = 10^{-3})$ has completely suppressed the vortex.

To understand how the gas pressure, the indirect force, and disc self-gravity affect the vortex structure, the density contours and the fluid streamlines in the frame corotating with the asymmetric structure are shown in the upper panels of Fig. 5. The indirect force alone in the leftmost panel can generate a lopsided structure which has a similar aspect ratio as the vortex in case g5 (the second panel from the left). However, they have very different velocity structures as shown in the bottom panels of Fig. 5, indicating the presence of the gas pressure has significantly altered the fluid velocity structure. When both the indirect force and the gas pressure are considered, the indirect force gives an additional force to spin the vortex, making the vortex stronger with a smaller aspect ratio (case g5i, the third panel from the left). The strong vortex in g5i has a large radial extend. The radial width measured from the vortex density structure is ~3.4 $H(R = R_c)$ (Table 1). The radial width measured from the largest elliptical vortex streamline in Fig. 5 is ~4 $H(R = R_c)$. These results suggest that the vortex can have an envelope extending beyond the

sonic point of the vortex by a factor of ~2 (also in Fu et al. 2014; Zhu et al. 2014).

However, when disc self-gravity is included and the disc is relatively massive (the rightmost panels in Fig. 5 for case g10gi), the streamlines remarkably go across the asymmetric structure without any rotary vortex motion around the density maximum. The vortex motion occurs at a much larger distance around $r \sim 1.2$, and the density contours do not coincide with the velocity streamlines. This is due to the fact that the lopsided structure orbits around the star at a speed significantly slower than the local Keplerian speed when the disc self-gravity starts to dominate the vortex dynamics. In these massive discs, the ‘fast mode’ appears to be transiting to the ‘slow mode’. As shown in Fig. 6, the red pluses suggest that the pattern speed for the lopsided structure is 80 per cent of the Keplerian speed in case g10gi. The corotation radius where the disc’s Keplerian angular speed<fnref n="3" /> matches the structure’s pattern speed is thus at $r \sim 1.2$ instead of $r \sim 1$. In the frame corotating with the disc asymmetric structure as shown in the rightmost panel of Fig. 5, the zero velocity point is thus at $r \sim 1.2$.

We would expect that such sub-Keplerian asymmetric disc structure cannot trap dust particles at the gas maximum. This is because dust particles always try to move at the local Keplerian speed and they will move in and out of the sub-Keplerian asymmetric structure

<footnote n="3">³ Strictly speaking, the disc does not rotate at the exact Keplerian speed due to disc self-gravity.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** The pattern speed with respect to the vortex position in case g5 (black pluses), g10i (green pluses), and g10gi (red pluses). The dotted curve represents the Keplerian speed. The squares are Keplerian speed at the same radial positions as the pluses. The dashed line labels the pattern speed for the g10gi case. The corotation radius is at $r \sim 1.2$, further than the vortex itself.</figcaption></figure>

quickly and cannot be trapped. But for particles at the corotation radius, they will remain at the same relative position with the asymmetric structure and can potentially be affected by the asymmetric disc structure. When we include dust particles in the simulations (Paper II), some particles indeed concentrate to the corotation radius at $r \sim 1.2$. Particles outside the vortex motion at $r \sim 1.2$ feel a non-zero gas drag and sink to the centre of this vortex motion. This interesting radial offset between the gas and dust lopsided structures will be shown in Paper II and may be observable with ALMA.

The detailed force balance for the vortex in various cases has been presented in Appendix A, where we can see the pressure gradient plays an important role in determining the vortex structure.

## 5 CONCLUSIONS

We have studied asymmetric flow patterns in massive discs. We have first searched possible 'fast modes' (lopsided structures orbiting around the central star at the Keplerian speed) in pressureless fluid using test particle methods. Without the indirect force that comes from the displacement of the central star by the lopsided structure, the only 'fast mode' is the traditional pressureless anticyclonic vortex with an aspect ratio ($\chi$) of 2. When the indirect force is considered, besides the pressureless vortex mode, another 'fast mode' whose streamlines consist of slowly liberating horseshoe orbits exist, as shown by Mittal & Chiang (2015). The streamlines in this mode are highly elongated with very large aspect ratios.

When the gas pressure is included, these two modes are less distinct. Together with the indirect force, the gas pressure significantly alters the lopsided structure. Vortices with any aspect ratio are now possible even due to the gas pressure alone (Goodman et al. 1987). Vortex with a larger aspect ratio rotates around the vortex centre slower, and can be more easily affected by other additional forces (such as the indirect force and self-gravity).

Using two dimensional global hydrodynamical simulations with both the indirect force and disc self-gravity self-consistently included, we have studied the effect of the indirect force on the vortex structure, which is initiated by the RWI. The simulations confirm that the indirect force alone can widen the vortex streamlines to be more than two disc scale heights wide in the radial direction when

the mass ratio between the disc and the central star is $\gtrsim 1$ per cent, and the vortex migrates faster when the indirect force becomes important. However, for such discs, disc self-gravity becomes equally (if not more) important than the indirect force. In massive discs where $(\pi/2)(r/H)Q^{-1} > 1$ (Lovelace & Hohlfeld 2013), disc self-gravity alone suppresses the $m=1$ mode. But the indirect force restores the $m=1$ lopsided structure. Still, overall, the vortex is significantly weakened by disc self-gravity. Vortices' inward migration slows down. Vortices in some cases can even migrate outwards.

One important observation is that, when the disc is massive enough (e.g. $Q \sim 3$), the lopsided gas structure orbits around the star at a speed significantly slower than the local Keplerian speed. In this case, there is a radial shift between the lopsided structure itself and the corotation radius where the disc's Keplerian rotation matches the structure's pattern speed. Since dust can be trapped at the corotation radius, it suggests that there could be a radial shift between the gas and dust distributions in vortices of massive discs. Dust distribution in vortices in massive discs is presented in Paper II. Overall, disc self-gravity is important to regulate the vortex structure in massive discs and it has observational signatures which may be probed by current and future observations.

Our simulations have several limitations. First, we have only explored a small parameter space of disc structures. Besides $H/r=0.1$, we have also explored discs with $H/r = 0.05, 0.07, 0.08$, and $0.09$ at $r_0$ while keeping $\sigma = 2H(r_0)$ in equation (6). Significant growth of the RWI has been observed only when $H/r > \sim 0.09$ despite the fact that the density bump is always 2 scale height wide in each case, implying the vortex production by the RWI is sensitive to the disc and density bump structure. Second, our 2D simulations cannot capture vortex instabilities occurring in 3D, such as the elliptical instability (Lesur & Papaloizou 2009). The effects of indirect force and disc self-gravity on these instabilities need further study. Moreover, in massive discs where the gas and dust lopsided structures have a radial offset, how such offset affects the parametric instability from dust to gas feedback also needs to be revisited (Railton & Papaloizou 2014).

## ACKNOWLEDGEMENTS

All hydrodynamical simulations are carried out using computer supported by the Princeton Institute of Computational Science and Engineering, and the Texas Advanced Computing Center (TACC) at The University of Texas at Austin through XSEDE grant TG-AST130002. ZZ acknowledges support by NASA through Hubble Fellowship grant HST-HF-51333.01-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555.

## REFERENCES

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

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

* Baruteau C., Masset F., 2008a, ApJ, 672, 1054

* Baruteau C., Masset F., 2008b, ApJ, 678, 483

* Baruteau C., Zhu Z., 2016, MNRAS, 458, 3927 (Paper II)

* Binney J., Gerhard O. E., Stark A. A., Bally J., Uchida K. I., 1991, MNRAS, 252, 210

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

* Casassus S. et al., 2015, ApJ, 812, 126

* Chavanis P. H., 2000, A&A, 356, 1089

* Christodoulou D. M., Narayan R., 1992, ApJ, 388, 451
Contopoulos G., 1979, in Evans D. S., ed., Photometry, Kinematics and Dynamics of Galaxies. University of Texas, Austin, p. 425

de Val-Borro M., Artymowicz P., D'Angelo G., Peplinski A., 2007, A&A, 471, 1043

Espaillat C. et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. University of Arizona Press, Tucson, p. 497

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

Godon P., Livio M., 2000, ApJ, 537, 396

Goodman J., Narayan R., 1988, MNRAS, 231, 97

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

Heemskerk M. H. M., Papaloizou J. C., Savonije G. J., 1992, A&A, 260, 161

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

Johansen A., Andersen A. C., Brandenburg A., 2004, A&A, 417, 361

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

Les R., Lin M.-K., 2015, MNRAS, 450, 1503

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

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

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

Lin M.-K., 2015, MNRAS, 448, 3806

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

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

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

Masset F., 2000, A&AS, 141, 165

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

Mittal T., Chiang E., 2015, ApJ, 798, L25

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

Paczynski B., 1977, ApJ, 216, 822

Papaloizou J. C. B., Pringle J. E., 1984, MNRAS, 208, 721

Papaloizou J. C. B., Pringle J. E., 1985, MNRAS, 213, 799

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

Regály Z., Juhász A., Sándor Z., Dullemond C. P., 2012, MNRAS, 419, 1701

Tanga P., Babiano A., Dubrulle B., Provenzale A., 1996, Icarus, 121, 158

Tremaine S., 2001, AJ, 121, 1776

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

Yellin-Bergovoy R., Heifetz E., Umurhan O. M., 2015, preprint (arXiv:1503.08470)

Zhu Z., Stone J. M., 2014, ApJ, 795, 53

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

## APPENDIX A: FORCE BALANCE WITHIN THE VORTEX

The force balance at $\theta = \pi$ for various cases is shown in Fig. A1. For a circular orbit, the tidal force $-\partial\Phi_*/\partial r + \Omega_t^2 r$, Coriolis force ($2\Omega_f v_\phi$), and centrifugal force ($v_\phi^2/r$) balance each other. When only the indirect force is included (the leftmost panel), the velocity at $r=1$ becomes slightly super-Keplerian since the Coriolis force and the centrifugal force need to balance the indirect force. Thus, the velocity in the rotating frame at $r=1$ becomes $v_\phi \sim v_K + f_{ind}/2\Omega$ where $v_K$ is the Keplerian velocity. The leftmost panel of Fig. A1 also suggests that the force which turns the vortex ($-v_r \frac{\partial v_r}{\partial r} - v_\phi \frac{\partial v_r}{\partial \phi}$) is at a similar amplitude as the indirect force. With the gas pressure included (right-hand panels), the pressure gradient (the cyan dotted curve) is quite large, and its balance with other forces leads to the vortex rotation. When disc self-gravity has been included (the rightmost panels), the self-gravitational force becomes zero at $r \sim 0.95$. The pressure gradient becomes zero at $r \sim 1$. To balance the self-gravitational force at $r \sim 1$, the addition of tidal force, Coriolis force, and centrifugal force has to be larger than zero at $r \sim 1$. This is possible if the corotation radius is at $r > 1$.

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

**Figure A1.** The disc surface density and force balance with respect to $r$ at $\phi = \pi$ for horseshoe orbits in pressureless fluid, g5, g5i, and g10gi at 150 orbits (from left to right). Bottom panels: various force components. The red curve is $-\frac{\partial\Phi_+}{\partial r} + \Omega_t^2 r + 2\Omega_f v_\phi$, the blue dotted curve is $\frac{v_\phi^2}{r}$, the cyan dotted curve is $-\frac{\partial P}{\partial r}$, the green dotted curve is $-v_r \frac{\partial v_r}{\partial r} - v_\phi \frac{\partial v_r}{\partial \phi}$, the black dotted curve is $-\frac{\partial\Phi_{ind}}{\partial r}$, and the orange dotted curve is $-\frac{\partial\Phi_{sg}}{\partial r}$. The solid curves are the total forces which add one additional force represented by the dotted curve in the same colour.

This paper has been typeset from a TeX/LA TeX file prepared by the author.