# Morphology and dynamical stability of self-gravitating vortices. Numerical simulations

S. Rendon Restrepo, P. Barge

► To cite this version:

S. Rendon Restrepo, P. Barge. Morphology and dynamical stability of self-gravitating vortices. Numerical simulations. Astronomy and Astrophysics - A&A, 2022, 666, 10.1051/0004-6361/202243518. hal-03948190

HAL Id: hal-03948190

https://hal.science/hal-03948190

Submitted on 26 Jan 2023

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
# Morphology and dynamical stability of self-gravitating vortices

## Numerical simulations*

S. Rendon Restrepo<fnref n="1" />,<fnref n="2" />* and P. Barge<fnref n="1" />

Received 10 March 2022 / Accepted 20 June 2022

## ABSTRACT

**Context.** Theoretical and numerical studies have shown that large-scale vortices in protoplanetary discs can result from various hydrodynamical instabilities. Once produced, such vortices can survive nearly unchanged over a large number of rotation periods, slowly migrating towards the star. Lopsided asymmetries recently observed at sub-millimetre and millimetre wavelengths in a number of transition discs could be explained by the emission of the solid particles trapped by vortices in the outer disc. However, at such a distance from the star, disc self-gravity (SG) may affect the vortex evolution and must be included in models.

**Aims.** Our first goal is to identify how vortex morphology is affected by its own gravity. Next, we look for conditions that a self-gravitating disc must satisfy in order to permit vortex survival at long timescales. Finally, we characterise as well as possible the persistent self-gravitating vortices we have found in isothermal and non-isothermal discs.

**Methods.** We performed 2D hydrodynamic simulations using the RoSSBi 3.0 code. The outline of our computations was limited to Euler's equations assuming a non-homotropic and non-adiabatic flow for an ideal gas. A series of 45 runs were carried out starting from a Gaussian vortex-model; the evolution of vortices was followed during 300 orbits for various values of the vortex parameters and the Toomre parameter. Two simulations, with the highest resolution thus far for studies of vortices, were also run to better characterise the internal structure of the vortices and for the purpose of comparison with an isothermal case.

**Results.** We find that SG tends to destabilise the injected vortices, but compact small-scale vortices seem to be more robust than large-scale oblong vortices. Vortex survival critically depends on the value of the disc's Toomre parameter, but may also depend on the disc temperature at equilibrium. Disc SG must be small enough to avoid destruction in successive splitting and an approximate 'stability' criterion is deduced for vortices. The self-gravitating vortices that we found persist during hundreds of rotation periods and look like the quasi-steady vortices obtained in the non-self-gravitating case. A number of these self-gravitating vortices are eventually accompanied by a secondary vortex with a horseshoe motion. These vortices reach a new rotational equilibrium in their core, tend to contract in the radial direction, and spin faster.

**Conclusions.** We propose an approximate 'robustness criterion', which states that, for a given morphology, a vortex appears stable provided that the disc's Toomre parameter overcomes a fixed threshold. Global simulations with a high enough numerical resolution are required to avoid inappropriate decay and to follow the evolution of self-gravitating vortices in protoplanetary discs. Vortices reach a nearly steady-state more easily in non-isothermal discs than in isothermal discs.

**Key words.** protoplanetary disks – hydrodynamics – methods: numerical – planets and satellites: formation

## 1. Introduction

The formation of planetesimals from the dust-grains embedded in the gas of protoplanetary discs (hereafter PPDs) remains one of the key problems of planet formation. Standard scenarios start from the sub-layer that the dust forms in the mid-plane of the disc after settling under the vertical component of the star's gravity. They are based on the collisional growth of the solid particles through coagulation and/or sticking mechanisms, generally followed by the formation of self-gravitating dust-clumps that are thought to compact into kilometre-sized planetesimals (Safronov 1972; Goldreich & Ward 1973). However, these models are known to stumble over a number of 'barriers' (metre-size, fragmentation, bouncing; Whipple 1972; Weidenschilling 1977; Dullemond & Dominik 2005; Blum & Wurm 2008; Zsom et al. 2010) that prevent them from giving a complete picture of the

problem. More recent models are based on resonant drag instabilities such as the streaming instability (Youdin & Goodman 2005; Jacquet et al. 2011), in which pressure bumps can marry with dust-clumps to deliver self-gravitating dust-clouds that are dense enough to form big planetesimals (up to Ceres-sized planetesimals; Johansen et al. 2007). Once planet cores have grown due to the accumulation of planetesimals or the sweeping of residual pebbles (Lambrechts & Johansen 2012), another important requirement for models is to reach the critical mass for run-away gas accretion before the solid material is accreted by the star and the gas is dissipated by winds and photoevaporation.

Another scenario proposes that anti-cyclonic vortices, formed during the early history of PPDs, could play an important role in the disc evolution, with a redistribution of angular momentum and a clumping of the solid material that (possibly) results in the formation of planetesimals or planetary cores (Barge & Sommeria 1995). Indeed, such vortices are known to very efficiently trap the dust particles, rapidly increasing the dust-to-gas ratio in the core of the vortices (this mass-ratio ~1

<footnote n="1">¹ Aix-Marseille Univ., CNRS, CNES, LAM, Marseille, France
e-mail: steven.rendon-restrepo@lam.fr; pierre.barge@lam.fr</footnote>

<footnote n="2">² Aix-Marseille Univ., CNRS, Centrale Marseille, IRPHE, Marseille, France</footnote>

<footnote>* Movies are only available at https://www.aanda.org</footnote>
in a few orbital periods for an optimal Stokes parameter, e.g. Raettig et al. 2015).

Theoretical works and numerical simulations have shown that gaseous vortices in PPDs can be the outcomes of hydro-dynamical instabilities, such as the Rossby wave instability (hereafter RWI; Lovelace et al. 1999; Li et al. 2000). This instability can grow when the generalised potential vortensity distribution reaches a local extrema<fnref n="1" />, which has also been studied in 3D (Meheut et al. 2010; Lin 2012; Richard et al. 2013). Such a situation is favoured at the boundaries of the magnetically inactive region (Varnière & Tagger 2006; Terquem 2008; Regály et al. 2012) or at the edges of the gap carved by a giant planet (de Val-Borro et al. 2006, 2007; Lyra et al. 2009; Lin & Papaloizou 2011; Baruteau et al. 2019). Besides RWI, many other instabilities were found to generate vortices, such as the baroclinic instability (Klahr & Bodenheimer 2003; Klahr 2004; Petersen et al. 2007a,b; Lesur & Papaloizou 2010; Barge et al. 2016), the convective overstability (Teed & Latter 2021), the vertical shear instability (Richard et al. 2016), or the zombie vortex instability (Marcus et al. 2015). Klahr et al. (2018) give a detailed review of instabilities occurring on PPDs, and further details about the origin and evolution of vortices can be found in Meheut et al. (2012, 2013), Regály et al. (2012, 2013), Fu et al. (2014a), Hammer et al. (2017), Baruteau et al. (2019), Robert et al. (2020) and Rometsch et al. (2021).

On the other hand, observations with high-resolution instruments, such as the Atacama Large Millimeter/submillimeter Array (ALMA; Andrews et al. 2018), have revealed strong azimuthal asymmetries in the dust thermal emission of transition discs around young stars such as Oph IRS 48 (van der Marel et al. 2013), MWC 758 (Dong et al. 2018), HD 135344B (Cazzoletti et al. 2018), HD 143006 (Pérez et al. 2018), TW Hya (Tsukagoshi et al. 2019), and HD 142527 (Soon et al. 2019). It is striking that such asymmetries can be interpreted as the thermal emission of dust confined in an undetected gaseous vortex. The possibility that such asymmetries are actually tracks of vortices in PPDs has been explored by various authors, but more solid arguments will have to wait for observations from instruments like the James Webb space telescope (JWST; with the Mid-InfraRed Instrument (MIRI) instrument in the thermal mid-IR region) or the future Extremely Large Telescope (ELT; using High-REsolution Spectrograph (HIRES) in synergy with ALMA).

From a theoretical point of view, strong dust concentrations in PPDs are known to be the location of instabilities due to the dust feedback. Similar situations also occur in the strong dust concentrations produced by a vortex, as observed in numerical 2D simulations (Fu et al. 2014b; Surville et al. 2016). In this case various dust and gas instabilities are at work, such as the Kelvin-Helmoltz instability or specific resonant-drag instabilities (Squire & Hopkins 2018). This evolution is complex and seems to end in a dissipation of the vortex. Of course, a more realistic description of the problem should include the vortex generation mechanism; for example, Miranda et al. (2017) claim that dust-laden vortices could last for thousands of orbits, which would be sufficient for dust growth or gravitational collapse. In any case, in 2D models, the dissipation timescales are difficult to estimate since the back-reaction of the dust particle becomes problematic at a high dust-to-gas ratio. 3D simulations are scarce

and seem to show that vortices remain stable and behave as Taylor columns, even if the mid-plane region is perturbed by a thin dust-layer (Raettig et al. 2021). These authors also find that dust densities in these vortices can be higher than Roche's density, possibly leading to a gravitational collapse. Another issue of the vortex scenario is the possible growth of the elliptical instability, which is known to destroy 3D-vortices (Lesur & Papaloizou 2009). However, vortices with a large enough aspect ratio could survive due to a weaker growth rate of this instability (Richard et al. 2013), and furthermore, the exact role of compressibility on the growth rate also appears unclear.

The impact of self-gravity (SG) on the formation and evolution of giant vortices in PPDs is key in the vortex scenario, particularly in order to explain the observed disc asymmetries. This was first addressed in a number of two-dimensional studies: Lovelace & Hohlfeld (2013) have shown that SG can affect RWI only if the Toomre parameter is less than the inverse of the non-dimensional geometric aspect ratio *h*; other authors have claimed that 2D vortices are either stretched by gravitational torques (Regály & Vorobyov 2017b), strengthened by the indirect potential induced by the vortex itself (Mittal & Chiang 2015; Zhu & Baruteau 2016), or weakened by thermal diffusion (Tarczay-Nehéz et al. 2020). On the other hand, three-dimensional studies (Lin & Pierens 2018) seem to show that 3D vortices with a turbulent core could survive the growth of the elliptic instability. Thus, the exact role of SG in the vortex evolution is still an open question.

In this paper we revisit the problem using the same strategy as Surville & Barge (2015), looking for approximate numerical vortex solutions of the compressible Euler equation when it is coupled to the Poisson equation. Thus, our simulations consisted of following up a Gaussian vortex initially superimposed on a disc at equilibrium. The vortex persistence was checked all along the evolution, particularly against splitting into secondary vortices. Following this procedure, a series of simulations were carried out and our first task was to organise the results to build up a 'stability' map of the vortices, which we used to identify the most suitable conditions for vortex survival. Finally, high-resolution (HR) simulations were carried out to characterise as well as possible the internal structure of a self-gravitating vortex hosted in an isothermal and a non-isothermal disc.

This work is organised as follows. Section 2 is devoted to setting up the disc model and presenting the methodology used in the paper. Then, in Sect. 3, we compare the evolution of a fiducial quasi-stationary vortex for different values of the Toomre parameters, which permits us to distinguish between three different cases according to the importance of SG. In Sect. 4, a large parametric study is performed in order to evaluate the stability of vortices with respect to the axisymmetric disc's Toomre parameter and the vortex internal structure. This enabled us to choose, in Sect. 5, the most promising structures and compare the HR simulations for vortices evolving in isothermal and non-isothermal discs. Finally, in Sect. 6 we suggest a discussion on observational constraints, numerical resolution, and spiral waves, followed by a conclusion.

## 2. Theoretical context

This section is devoted to introducing the main theoretical background of this paper. We begin with the description of the equations governing the disc evolution.

<footnote n="1">¹ This quantity is $\mathcal{L}(r) = (\sigma\Omega/\kappa^2) (\frac{p}{\sigma\gamma})^{2/\gamma}$ and the inverse of the vortensity, $\sigma/(\nabla \times \mathbf{v}) \cdot e_z$, for non-homentropic and for homeontropic 2D discs, respectively.</footnote>
## 2.1. Disc model and evolution equations

We investigate the evolution of an ideal gas (molecular dihydrogen, $\mu \sim 2.34$) in a thin (2D), non-homotropic disc orbiting a solar-type star. Compressibility must be taken into account since the vortex size often exceeds the scale height of the disc. On the other hand, viscosity is neglected since we are working in the magnetically dead zone of the protoplanetary discs. The equations governing gas dynamics in a 2D protoplanetary disc are thus the continuity equation, Euler's equation, and the conservation of energy equation, which are:

$$ \partial_t \sigma + \nabla \cdot (\sigma \mathbf{v}) = 0, \quad (1) $$

$$ \partial_t \mathbf{v} + (\mathbf{v} \cdot \nabla) \mathbf{v} = -\frac{\nabla P}{\sigma} - \nabla (\Phi_{\odot} + \Phi_{\text{SG}} + \Phi_{\text{ind}}), \quad (2) $$

$$ \sigma [\partial_t e + (\mathbf{v} \cdot \nabla) e] = -P \nabla \cdot \mathbf{v}, \quad (3) $$

where $P$, $\sigma$, and $\mathbf{v}$ are the vertically integrated pressure, density, and the gas velocity, respectively, $e = \frac{1}{\gamma-1} \frac{P}{\sigma}$ is the internal energy of an ideal gas, and $\gamma = 1.4$ is the adiabatic index for a diatomic gas. The gravitational potential has three different contributions: $\Phi_{\odot}$, $\Phi_{\text{SG}}$, and $\Phi_{\text{ind}}$ due to the central object, the disc SG and the vortex itself (indirect term), respectively. The last term is a virtual potential that accounts for the offset between the centre of the frame of reference, located at the star's position, and the {star+vortex} centre of mass. No viscosity is expected in the non-active region of the PPDs, although Stoll & Kley (2014) noted that, even in this region, vertical-shear instability can generate turbulence associated with an effective viscosity. According to Manger et al. (2020, 2021), Formula (14), our simulations correspond to an infinite cooling time and to $\alpha/(10^{-4}) \ll 1$, so that the viscosity is very small and compatible with our inviscid disc assumption. In the case of a finite cooling time, we should account for a non-negligible viscosity in the computations.

At equilibrium the background flow is assumed to be axisymmetric with radial profiles for the temperature $T_0(r)$, the surface density $\sigma_0(r)$, and the pressure $P_0(r)$ chosen in the form of simple power laws. The model only differ from the minimum mass solar nebulae (MMSN) hypothesis (Hayashi 1981) by the choice of the local values of the surface density index $\beta_{\sigma}$ and the temperature index $\beta_T$ to account for the observational constraints. With these assumptions, the unperturbed axisymmetric disc's Toomre parameter, $Q_0(r) = \Omega_k c_{s,0}/(\pi G \sigma_0)$, also obeys power laws that depend on the surface density in the centre of the simulation box. Thus, the main variables defining the stationary, axisymmetric flow are:

$$ \sigma_0(r) = \sigma_0 (r/r_0)^{\beta_{\sigma}}, \quad \text{with } \beta_{\sigma} = -1.1, $$

$$ T_0(r) = T_0 (r/r_0)^{\beta_T}, $$

$$ H/r = 0.05 (r/r_0)^{\beta_H - 1}, \quad \text{with } \beta_H = \frac{\beta_T + 3}{2}, $$

$$ Q_0(r) = 0.0279 \left( \frac{1700 \text{ g cm}^{-2}}{\sigma_0} \right) (r/r_0)^{\beta_Q}, \quad \text{with } \beta_Q = \beta_H - 3 - \beta_{\sigma}, \quad (4) $$

where $r$ is the normalised radius at 1 AU, $r_0 = 50$, $\sigma_0(r)$ and $\sigma_0$ are the surface density (in g cm$^{-2}$) at radial positions $r$ and $r_0$, respectively, and $H$ is the pressure scale height. $T_0(r)$ and $T_0 = 6.2$ K refer to the temperature at radial positions $r$ and $r_0$, respectively. In particular, the temperature was chosen such that the disc flaring, $h$, in the centre of the simulation box is equal to 0.05. For an ideal gas, the axisymmetric pressure also obeys a

<table><caption>**Table 1.** SG parameter and surface density at $r_0 = 50$ AU.</caption><tr><td>$q_0$</td><td>0.25</td><td>0.5</td><td>1</td><td>2</td><td>4</td></tr><tr><td>$\sigma_0$ (g cm<sup>-2</sup>)</td><td>9.38</td><td>4.69</td><td>2.35</td><td>1.17</td><td>0.59</td></tr></table>

power law with index $\beta_p = \beta_T + \beta_{\sigma}$. In Sects. 3 and 4, and for the non-isothermal disc studied in Sect. 5, we have set $\beta_T = -0.5$. On the other hand, we have chosen $\beta_T = 0$ for the isothermal disc considered in Sect. 5.

## 2.2. Gravity terms

In polar coordinates the SG potential reads:

$$ \Phi_{\text{SG}}(\mathbf{r}) = -G \int_{\text{disc}} \frac{\sigma_0(r')}{|\mathbf{r} - \mathbf{r}'|} d^2 \mathbf{r}' - G \int_{\text{disc}} \frac{\sigma(r', \theta') - \sigma_0(r')}{|\mathbf{r} - \mathbf{r}'|} d^2 \mathbf{r}', \quad (5) $$

where the first term is the axisymmetric contribution of the disc, while the second is the contribution of all other symmetry deviations. For example, the last term can correspond to the presence of an anti-cyclonic vortex or to a density wave. In order to account for the disc thickness in SG computations, we considered a softening length equal to 0.6 H (Müller et al. 2012). The unperturbed disc potential can be expressed analytically as an infinite sequence (Huré et al. 2008), which, truncated at second order, allows the axisymmetric disc potential to be estimated at an accuracy better than 99% (Surville 2013, in French). The presence of a massive vortex in the outer disc can lead to an offset of the mass barycentre, which may strengthen the vortex itself (Regály & Vorobyov 2017a; Zhu & Baruteau 2016). To account for this possible effect, the gradient of the indirect potential has been taken into account; this additional term can be written as:

$$ \nabla \Phi_{\text{ind}} = \left( \frac{A \cos \theta + B \sin \theta}{-A \sin \theta + B \cos \theta} \right), \quad (6) $$

$$ \text{where } A = \int_{\text{disc}} \frac{G \cos \theta' dm'}{r'^2} \quad \text{and} \quad B = \int_{\text{disc}} \frac{G \sin \theta' dm'}{r'^2}. $$

For convenience, all of our simulations are referenced and labelled with the Toomre parameter of the unperturbed disc. Following Lovelace & Hohlfeld (2013), it is also useful to introduce in the discussion the critical parameter $Q_c = 1/h$. Indeed, RWI is significantly changed by SG when $q_0 = Q_0/Q_c \le 1$. The normalised value of the Toomre parameter, $q_0$, will be used through the rest of the paper and will be called the SG parameter. Our investigation was restricted to $q_0 = [0.25, 0.5, 1, 2, 4]$ because below the lower limit, any vortex is unstable, whereas beyond $q_0 = 4$ (the upper limit), SG has no observable effect on the vortex. The relationship between $q_0$ and $\sigma_0$, the surface density at 50 AU, is presented in Table 1.

## 2.3. Numerical method and physical parameters

The system of Eqs. (1)–(3) was solved numerically thanks to the latest version of the RoSSBi (Rotating Systems Simulation for Bi-fluids) code, which was specifically developed to study the evolution of PPDs (Inaba et al. 2005; Inaba & Barge 2006; Surville 2013). This new version has been extended to 3D and structured for high-performance parallelism; it is presented in detail elsewhere (Rendon Restrepo et al. 2022) with many tests
<table><caption>**Table 2.** Vortex parameters at *t* = 50 *t*₀ with respect to the initial Gaussian model parameters.</caption><thead><tr><th>δ<br>χ</th><th>1.0</th><th>1.5</th><th>2.5</th></tr></thead><tbody><tr><td>5</td><td>(7.50, 0.79)</td><td>(6.21, 1.43)</td><td>(5.65, 2.78)</td></tr><tr><td>8</td><td>(8.78, 1.02)</td><td>(9.07, 1.68)</td><td>(9.00, 3.06)</td></tr><tr><td>14</td><td>(10.73, 1.12)</td><td>(13.6, 1.65)</td><td>(12.62, 2.69)</td></tr></tbody></table>

*Notes.* Column and row header: Gaussian model parameters at *t* = 0.

and validations. The code includes SG, which is computed using the Fastest Fourier Transform in the West 3 library (FFTW3; Frigo & Johnson 2005) with a logarithmic and a linear grid in the radial and azimuthal direction, respectively.

In the absence of both SG and viscosity, a number of authors (Godon & Livio 1999; Barranco & Marcus 2005) started numerical simulations from crude vortex models. Later, other studies started from approximate vortex models of the steady state equations (Lesur & Papaloiou 2009; Surville & Barge 2015). A Gaussian vortex model is particularly suitable in the case of compressible simulations; the injected vortices were found to rapidly relax to long-lived vortices, which can then be studied in detail (Surville & Barge 2015). Here, the same strategy is used to address the problem of vortices in a disc that has SG. Thus, instead of generating vortices through RWI or baroclinic instability, we preferred to inject a Gaussian solution as modelled in Surville & Barge (2015) since this allows a large variety of vortices to be explored by tuning their control parameters (radial width, Rossby number, and aspect ratio). This method allows the most appropriate vortex parameters to be selected, but also to save computational time since the timescale necessary to generate a steady vortex from a hydrodynamic instability generally exceeds hundreds of orbits. Since the geometrical parameters of the vortex can change during the relaxation phase, we gathered in Table 2 all of the parameters used to initialise our numerical simulations, as well as their value at *t* = 50 *t*₀, the time at which SG is plugged in. It is interesting to note that in presence of viscosity, *α* = 10⁻⁴, and simulations with or without SG have shown that vortices decay at a rate proportional to the disc mass (Regály & Vorobyov 2017b).

### 2.3.1. A smooth activation of SG

In the presence of SG, numerical tests (at low-resolution, hereafter LR) have shown that persistent vortices can also survive, but after a transition period that is longer and more complex than in the non-self-gravitating case. This transient evolution is due to the difference between the injected and the exact vortex-solutions; of course it reduces to the relaxation already discussed in the previous section for vanishing SG. For these reasons, and to keep the same strategy as previously used, it was decided to gradually introduce SG during the computations, in order to avoid spurious vortex bursts, which occur when SG is introduced with a Heaviside time step function. To this end we have implemented in the code the time function:

$$g(t) = \begin{cases} 0 & \text{if } t \le 60 t_0 \\ 1 - \exp\left(-\left(\frac{t-60 t_0}{T_{\text{act}}}\right)^2\right) & \text{otherwise} \end{cases}. \quad (7)$$

The time constant was chosen to be equal to the vortex relaxation time in the non-SG case; that is $T_{\text{act}} \approx 15 t_0$, where $t_0 = 353$ years is the orbital period at $r_0 = 50$ AU. With this procedure, SG reaches 98% of its final value in ≈30 orbital periods.

### 2.3.2. Choice of computational box

Besides the smooth activation of SG, another requirement is to correctly account for the vortex migration. In a self-gravitating disc, this migration not only results from the asymmetry of the vortex wake between the inner and outer regions (e.g. Paardekooper et al. 2013) but also from all the Lindblad torques produced by the density waves that are excited by the vortex.

Since in this work we are mainly interested by variations in vortices' structure and not in their type I migration, the initial vortex should keep its radial position during the whole run. This is possible provided that the computational box includes the same number of inner and outer Lindblad resonances. Therefore, prior to the main study, for each Gaussian vortex structure, we ran a preliminary study at LR, ($N_r, N_{\theta}$) = (360, 720), during 150 orbits and for a large computational window ($r \in [15, 100]$ AU) in order to catch all Lindblad resonances and forecast their location for the main study. This allowed us to choose the suited simulation window, $r \in [25, 72]$ AU, such that there was the same number of inner and outer Lindblad resonances during the whole simulation, to keep the vortex in the centre of the simulation box and avoid spurious vortex migration. Figure A.1 illustrates the extent of the full computational box and shows an example of compressibility waves emitted by a quasi-steady vortex.

### 2.3.3. Suitable physical parameters

A number of physical parameters have been chosen to characterise as well as possible the coherent gaseous vortices and to make easier the discussion. It is, for example, convenient to distinguish between morphological and global parameters such as the radial position, the bulk-mass, and the spin. Since elliptical streamline models (e.g. Kida 1981a; Goodman et al. 1987; Chavanis 2000) approximately describe the central region of 2D vortices, the most natural morphological parameters are the ellipse parameters such as the vortex aspect ratio $\chi$ and its radial extent (scaled to the pressure scale height), $\delta = \Delta r/H$. It is useful to recall that when the radial extent lies below unity, the flow is subsonic (incompressible), while when it lies above unity, the flow is supersonic (compressible). These structure parameters are determined by fitting ellipses on the iso-contours of the Rossby number at the $R_0 = -0.08$ level, where $R_0 = \nabla \times v'/2\Omega_k$ and $v'$ is the relative gas velocity with respect to the local pressure supported rotation rate (see Table 3 and Fig. C.2). Appendix C.4 gives the details of the procedure and explains why the Rossby number is chosen instead of density or pressure.

Global parameters are well suited to describe the bulk-vortex evolution once a nearly steady-state is reached, such as a long-term dynamical evolution. In particular, we numerically computed the distance from the vortex to the star, tracking the pressure maximum as a function of time. A debated question is whether vortices can be long-lived and robust against instabilities, and for this reason we computed the mean Rossby number in the vortex core or spin, $R_0$, namely the vortex strength. According to Kida (1981a) and GNG (Goodman et al. 1987) models, the vortex strength is inversely proportional to the aspect ratio. Numerical simulations by Surville & Barge (2015) have also shown that, even if the dependence of the Rossby number as a function of the vortex aspect ratio differs from the standard GNG relation ($1/\chi$), the general trend is the same with a saturation for elongated vortices: $\lim_{\chi \to \infty} R_0 \approx -0.1$. This general trend motivated some authors to use the inverse of the aspect ratio as a measure of vortex strength. Yet, comforted by
<table><caption>**Table 3.** Parameter definitions.</caption><thead><tr><th>Symbol</th><th>Definition/value</th><th>Description</th></tr></thead><tbody><tr><td>r</td><td></td><td>Normalised radius at 1 AU</td></tr><tr><td>σ</td><td></td><td>Gas surface-density</td></tr><tr><td>P</td><td></td><td>Vertically integrated pressure</td></tr><tr><td>T</td><td>μm<sub>u</sub> P / k<sub>B</sub> σ</td><td>Gas temperature</td></tr><tr><td>v</td><td></td><td>Gas velocity-field</td></tr><tr><td>v'</td><td>v - Ω<sub>K</sub> re<sub>θ</sub></td><td>Relative gas-velocity (with respect to the angular gas velocity)</td></tr><tr><td>Φ<sub>☉</sub></td><td></td><td>Central object potential</td></tr><tr><td>Φ<sub>SG</sub></td><td></td><td>Gas SG potential</td></tr><tr><td>Φ<sub>ind</sub></td><td></td><td>Indirect potential (offset of Φ<sub>☉</sub>)</td></tr><tr><td>Ω<sub>K</sub></td><td>√(GM<sub>☉</sub> / r + r<sup>2</sup> P<sub>0</sub> / σ<sub>0</sub><sup>2</sup>)</td><td>Pressure-supported rotation-rate</td></tr><tr><td>γ</td><td>1.4</td><td>Adiabatic index</td></tr><tr><td>c<sub>s</sub></td><td>√(γ<sup>2</sup> / σ<sup>2</sup>)</td><td>Gas sound-speed</td></tr><tr><td>H</td><td>√(2 / γ) c<sub>s</sub> / Ω<sub>k</sub></td><td>Pressure scale-height</td></tr><tr><td>h(r)</td><td>H / r = 0.05 (r/r<sub>0</sub>)<sup>βH-1</sup></td><td>Disc aspect ratio</td></tr><tr><td>Q</td><td>Ω<sub>k</sub> c<sub>s</sub> / (πGσ)</td><td>Local Toomre's parameter</td></tr><tr><td>Q<sub>0</sub></td><td>Ω<sub>k</sub> c<sub>s,0</sub> / (πGσ<sub>0</sub>)</td><td>Toomre's parameter at 50 AU for the Axis-symmetric disc</td></tr><tr><td>Q<sub>c</sub></td><td>1 / h = 20</td><td>Critical Toomre parameter</td></tr><tr><td>q<sub>0</sub></td><td>Q<sub>0</sub> / Q<sub>c</sub></td><td>SG parameter</td></tr><tr><td>Ro</td><td>∇ × v' / 2Ω<sub>k</sub></td><td>Rossby number</td></tr><tr><td>Ro</td><td>⟨Ro⟩<sub>core</sub></td><td>Vortex spin</td></tr><tr><td>χ</td><td></td><td>Vortex aspect-ratio</td></tr><tr><td>Δr</td><td></td><td>Vortex radial-width</td></tr><tr><td>δ</td><td>Δr / H</td><td>Vortex radial-extent</td></tr><tr><td>t<sub>0</sub></td><td>353 yr</td><td>Orbital period at 50 AU</td></tr><tr><td>C<sub>σ</sub></td><td>σ / σ<sub>0</sub> - 1</td><td>Density contrast</td></tr><tr><td>C<sub>P</sub></td><td>P / P<sub>0</sub> - 1</td><td>Pressure contrast</td></tr></tbody></table>

**Notes.** Quantities defined for the unperturbed axis-symmetric gas disc are denoted with a 0 subscript. The values provided for $t_0$, h, and $Q_c$ correspond to the values at 50 AU.

Surville & Barge (2015) and by results shown in Sect. 3.3.1, we decided to use Ro as the vortex strength through the rest of the paper. Finally, since we are mainly interested in SG effects, it is important to evaluate the mass that can be attributed to the vortex. This requires some caution since a vortex is a highly distributed object in contrast to planets, which are usually described as material points. Furthermore, a vortex is also accompanied by spiral waves, an annular over-density (see Sect. 5), and, eventually, some secondary eddies (see Sects. 3.4 and 5), which leads to a prudent calculation. We encourage the interested reader to refer to Appendix C, where we explain in detail the way we define the mass of a gaseous vortex and also its mean Rossby number. In Appendix D we estimated uncertainties. The definition of all the parameters and key quantities used in this paper are gathered in Table 3.

## 3. Vortex evolution

Six numerical simulations were run to address the problem, starting from a single reference vortex that was followed over 300 orbital periods. All runs performed in this section started

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Gaussian vortex model injected at $t = 0$. Top: Rossby number. Bottom: normalised density $(\sigma(r, \theta)/\sigma_0(r))$.</figcaption></figure>

from a Gaussian vortex with aspect ratio $\chi = 14$ and radial extent $\delta = 1.5$ which was injected at $r_0 = 50$ AU and $\theta_0 = \pi$. We exhibit in Fig. 1 the initial Gaussian vortex. Disc SG was quantified thanks to the SG parameter, $q_0$, whose values were selected among six different values: 0.25, 0.5, 1, 2, 4 and $\infty$ (non-SG case). The results of the simulations are presented in Figs. 2, 3, 4, and 5, and the temporal evolution can be found in online movies 1 and 2. The end-of-run snapshots of the surface density and the Rossby number, obtained for the six values of $q_0$, are presented in Figs. 2 and 3. Figures 4 and 5 focus on the evolution of physical quantities such as migration, mass, and spin. In the whole section, the case $q_0 = \infty$ and values at $t = 50 t_0$ are considered as the references on which comparisons are based. It is interesting to note that the slow increase in the vortex spin over the whole run is likely an artefact related to the numerical resolution as testified in Fig. E.1.

### 3.1. Numerical procedure for the run series

Computations started at $t = 0$ with a LR defined by ($N_r$, $N_\theta$) = (500, 720). During the first 50 orbits, the vortex morphological and global parameters did not stay constant but
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Comparison of the end-of-run surface density for increasing SG ($t = 300 t_0$). Top row and from left to right: $q_0 = \infty, 4, 2$; bottom row and from left to right: $q_0 = 1, 0.5, 0.25$. The initial vortex parameters are: $\delta = 1.5$ and $\chi = 14$. The simulation window was centred in $r \in [42.5, 57.5]$ AU in order to better appreciate the inner vortex structure.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Comparison of the end-of-run Rossby number for increasing SG ($t = 300 t_0$). Top row and from left to right: $q_0 = \infty, 4, 2$; bottom row and from left to right: $q_0 = 1, 0.5, 0.25$. The initial vortex parameters are: $\delta = 1.5$ and $\chi = 14$. The simulation window was centred in $r \in [42.5, 57.5]$ AU in order to better appreciate the inner vortex structure.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Morphological evolution for different values of the SG parameter. Top: semi-minor axis with respect to the pressure scale height at the vortex core for $q_0 = 0.25, 0.5, 1, 2, 4, \infty$. Bottom: aspect ratio at the $Ro = -0.08$ level. Initial vortex parameters are: $\delta = 1.5$ and $\chi = 14$. For readability, uncertainties are plotted at every fourth time steps (see Appendix D).</figcaption></figure>

slowly changed as a function of time. Indeed, the injected vortex was not an exact solution of the fluid equations but less than 50 orbits were sufficient for the injected solution to relax and reach a quasi-stationary regime (see Table 2). Then, at $t = 50 t_0$, the resolution of the simulations was increased using a 2D cubic interpolation<fnref n="2" /> into an intermediate grid resolution of ($N_r, N_\theta$) = (1500, 3000), which is equivalent to ~71 cells/H and ~24 cells/H in the radial<fnref n="3" /> and azimuth directions, respectively. At $t = 60 t_0$, SG was activated gradually thanks to the function described in Eq. (7). Finally, all the simulations were stopped at the 300th orbit. The successive steps of our procedure are recapped in a few lines:

$t = 0$ : Gaussian vortex at LR, with : ($N_r, N_\theta$) = (500, 720),

$0 < t < 50 t_0$ : Vortex relaxes to a quasi-steady state,

$t = 50 t_0$ : Resolution is increased to (1500, 3000),

$t = 60 t_0$ : Plug-in of SG,

$60 t_0 < t < 300 t_0$ : Simulation at intermediate resolution.

The rapid relaxation of the Gaussian solution was eased by the LR step, which introduced a significant numerical

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Bulk-vortex evolution for different values of the SG-parameter. Top: radial distance to the star. Middle: bulk of the vortex mass; $M_i$ is the vortex mass at $t = 0$. Bottom: absolute value of the mean Rossby number $Ro$. The initial vortex parameters are: $\delta = 1.5$ and $\chi = 14$. For readability, uncertainties are plotted at every fourth time steps (see Appendix D).</figcaption></figure>

viscosity. On the other side, during the last step, the choice of a higher resolution allowed a numerical convergence to be reached. Appendix E.1, stresses that LR can significantly speed up the decay of vortices. More quantitatively, during a complete simulation with $t \in [0, 300 t_0]$, we found that the vortex mass-loss is ~15% for every 250 orbits, if the resolution is ($N_r, N_\theta$) = (500, 720), while mass-loss is negligible if resolution is increased to ($N_r, N_\theta$) = (1500, 3000).

## 3.2. Weak SG: $q_0 = 2, 4$

There is no significant difference compared to the non-SG case where the vortex is slightly stretched by the shear. Only a

<footnote n="2">² With this interpolation, residuals are highly damped in less than two orbits.</footnote>

<footnote n="3">³ For a logarithmic mesh, the resolution at $r_0$ is $\delta_* = r_0(\alpha^{1/2} - \alpha^{-1/2})$ where $\alpha = (r_{\text{out}}/r_{\text{in}})^{1/N_r}$.</footnote>
marginal contraction and flattening of the vortex is observed; its
radial extent and aspect ratio change by less than 15 and 5%,
respectively. On the other hand, spin is not unchanged and vor-
tex mass is weakly affected (less than 5% with respect to initial
mass). Finally, as desired, migration is very slow, only 0.1 AU
in 300 orbits, which is nearly negligible with regards to the
estimated numerical precision.

3.3. Intermediate SG

This case corresponds to the transition between stretching and
splitting under gravitational forces. In the first occurrence, the
initial vortex is deformed but keeps a coherent structure, whereas
in the second one, coherence is lost and the vortex tends to split
into independent vortices.

3.3.1. Vortex shaped by the SG: $q_0 = 1$

If $q_0 = 1$, the vortex is deeply modified by the SG. It is stronger
than the initial vortex with a 1% increase in the mean Rossby
number, but it is also: (i) stretched in the azimuthal direction with
a 50% increase in aspect ratio; and (ii) contracted in the radial
direction with a 15% decrease in $\delta$. In particular, the aspect ratio
rapidly reaches a maximum at $t \sim 100 t_0$, followed by a slow
decrease until the end of the simulation. It is important to stress
that this increase in vorticity contrasts with the evolution of a
vortex only stretched by the shear. Migration is slightly accel-
erated to 0.15 AU in 300 orbits but, similarly to previous cases,
it remains very weak. Vortex mass declines by 5%. It is impor-
tant to mention that, in this case, both vortex strength and aspect
ratio increase simultaneously. This observation, in disagreement
with the statement that vortex strength is inversely proportional
to $\chi$, confirms the choice we made in Sect. 2.3.3. Finally, the
migration of the vortex is 50% faster than in the previous case
(which remains negligible over 300 orbits) and its mass declines
by $\sim 5\%$.

3.3.2. Vortex splitting: $q_0=0.5$

At $t = 60 t_0$, a sudden increase of 92% occurs in the aspect ratio, and a sudden decrease of 30% occurs in the radial extent. These sharp variations are followed, at the 90th orbit, by a splitting. Just before the splitting, $|Ro|$ sharply varies as a function of time. This likely indicates that the release of a secondary vortex results from the changes introduced by SG in the fluid-force equilibrium (pressure gradient, star attraction, and fictive forces); in the present case, ejection of secondaries could result from the predominance of the centrifugal forces. In particular, this secondary vortex picks up mass to the primary which, at $t = 90 t_0$, can lose more than 25% of its initial mass.

After the first release, the secondary vortex remains in the co-
orbital region and crosses the main vortex at ~117 $t_0$ and ~145 $t_0$.
During each encounter, the main vortex can pick up part of the
secondary mass and, at the end of the run, reaches 60% of its
initial value. After the 220th orbit the flow stabilises in a no-
encounter state. The secondary, in gravitational interaction with
the main vortex, has horseshoe oscillations of period $T_{hs} \sim 65t_0$.
The exchange of angular momentum between the two vortices
can easily explain the oscillatory behaviour observed in the
migration rate with the same period (see the top panel of Fig. 5).

3.4. Strong SG: $q_0 = 0.25$

In this case the injected vortex is stiffly modified with a drastic
increase in its spin and a strong radial contraction; it produces

two secondary vortices, which tend to cascade into smaller
vortices or transient eddies.

Similarly to the previous case, after $t = 75 t_0$, the vortex
is strongly deformed and loses more than 60% of its mass
during the ejection of two secondaries. The evolution of the
secondaries is also quite complex; it is governed by three-body
interactions with the primary-vortex, but also by successions
of splitting, merging, or dissipation processes. After the 190th
orbit, the vortex strongly decays and finally disintegrates at the
300th orbit in an annular overdensity and two residual vortices
that are located at [20°, 180°]. At this stage the density maxi-
mum in the annulus reaches half the peak density of the initial
vortex. At the end of the run, the main-vortex has lost 75% of
its mass.

3.5. Summary of vortex evolution

Three different regimes have been identified in our simulations.
The first evolution is dominated by (A) shear ($1 \le q_0$): there
is no significant departure from the non-SG case. The outer
regions of the vortex are slightly stretched by shear and SG.
The vortex reaches a quasi-steady regime with nearly constant
mass and migration. The second regime is the domain of (B)
self-gravitating vortices ($0.5 \le q_0 \le 1$): SG is able to change the
morphology and spin of the vortex: vortices spin slightly faster,
contract radially, and elongate. This is a new regime that vortices
can reach for intermediate values of $q_0$. Our simulations show
that such vortices have a stable and long-term evolution, similar
to the non-SG case. This suggests the existence of quasi-steady
vortex solutions of the Euler and Poisson equations, but further
investigations are needed for firm conclusions on this point. We
must stress that regime (B) includes the case of vortices releas-
ing a single secondary vortex, such as observed when $q_0 = 0.5$;
it corresponds, in fact, to a marginal case. Indeed, no splitting of
the main vortex is observed when the same run (at $q_0 = 0.5$) is
performed using the methodology presented in Sect. 5. The last
occurrence refers to (C) splitting by SG ($q_0 < 0.5$): SG strongly
affects the vortex shape, up to destruction, but also the global
flow. The initial vortex releases more than one secondary vortex,
which remain trapped in the co-orbital region of the primary. In
this regime the global parameters of the primary significantly
change as a function of time. For smaller values of $q_0$, either
the vortex directly releases multiple secondaries, or the process
repeats on the primary with successive releases of weaker sec-
ondaries and/or eddies. The result is a progressive fading of the
main vortex that correlates with an increasing number of eddies
and the growth of an annular overdensity.

In this last case $Q_0 = 5$ whereas, in the vortex core, the mean
value of Toomre parameter is $Q \sim 10$. These values of $Q_0$ and
$Q$ indicate that vortex splitting, although due to SG, cannot be
the result of a standard gravitational instability with a criterion
similar to that derived in the case of rotating discs. In the follow-
ing we suggest a possible way to extend this standard criterion to
address vortex stability.

4. Vortex stability

This section is devoted to better characterising the stability of
vortices submitted to their own gravity. This was achieved thanks
to two different methods: a study based on the compilation of
an extended series of numerical simulations, and a theoretical
estimate of the vortex parameters when constrained by Toomre's
criterion analogue.
<figure><img src="image_7.png" /><figcaption>Fig. 6. Stability domain of the self-gravitating vortices as a function of $q_0$. Left: when estimated from numerical simulations and right, when deduced from Eq. (9). Vortices are stable above the surface and unstable below (blue and orange, respectively). Both surfaces were projected on the left-hand side and the right-hand side of the box. $\chi$ and $\delta$ are, respectively, the vortex aspect ratio and the radial extent at $t = 50 t_0$.</figcaption></figure>

## 4.1. Numerical study

Here we used the same method as in Sect. 3, with the same numerical resolution but with initial parameters selected among the range of values gathered in Table 2. Our definition of vortex stability is a simple but empirical one: vortices are said to be ‘stable’ if (1) they survive until the end of the run ($300 t_0$), (2) there is a possible ejection of small scale eddies, and (3) there is a relative density contrast decrease no higher than 20%: $|C_\sigma(t_f) - C_\sigma(t_i)|/C_\sigma(t_i) \le 0.2$, where $t_i = 50 t_0$ and $t_f = 300 t_0$. A total of 45 different runs were carried out and provide a set of dots characterising the end-of-run state of the vortex in the $(\chi, \delta, q_0)$ space. Raw results are reported in Fig. B.1. They are presented another way in Fig. 6, in the form of a stability map in the left panel, and are compared to our theoretical estimates in the right panel (see Sect. 4.2). When the aspect ratio is constant: (i) if $\chi \sim 5$, vortex stability is independent of $\delta$; the vortex is stable if $q_0 > 0.375$ and unstable otherwise and (ii) if $\chi \sim 14$, vortex stability depends on $\delta$; the critical value of $q_0$ increases with $\delta$ (from 0.375 to 0.65), and therefore large vortices are more easily destabilised than small vortices (even at lower surface-densities). When the radial-extent is constant: (i) if $\delta \sim 1$, vortex stability is independent of the aspect ratio; the vortex is stable if $q_0 > 0.375$ and unstable otherwise and (ii) if $\delta \sim 2.5$, vortex stability depends on the aspect ratio; the critical value of $q_0$ increases with $\chi$ (from 0.375 to 0.65); in this case, elongated vortices are more easily destabilised than more compact vortices (even at lower surface-densities).

In summary, SG is able to significantly change the stability of gaseous vortices. Assuming standard disc configurations with $h(r_0) = 0.05$ and $\beta_\sigma = -1.1$, we find that SG destabilises all vortices if $q_0 < 0.375$ (or again $Q_0 < 7.5$) and all vortices remain stable for moderate SG when $q_0 > 0.65$. In a general way, small and compact vortices are more stable against SG splitting than large and elongated ones. Finally, a destructive splitting of the vortices can occur if $Q_0$ and $Q$ are approximately of the order of $\sim 8$.

## 4.2. Theoretical estimate

Gravitational stability of circumstellar discs is constrained by the famous Toomre criterion $c_s\kappa/(\pi G\sigma) > 1$, where $\kappa$ is the epicyclic

frequency (Toomre 1964); in contrast, little is known on the gravitational stability of massive gaseous vortices. Here, in an initial attempt to resolve the problem, we extrapolate this criterion by replacing the local epicyclic frequency with the vortex vorticity $\omega = 2|\overline{Ro}|\Omega_k$; assuming that $|\overline{Ro}| \sim |\overline{Ro}|_{\text{max}}$, we get:

$$2 \frac{c_s |\overline{Ro}|_{\text{max}} \Omega_k}{\pi G \sigma} \ge 1. \qquad (8)$$

If this condition is satisfied, we say that the vortex is stable against SG (hereafter ‘SG stable’). Of course such a definition is somewhat inappropriate since, stricto sensu, the criterion refers to the stability of a Keplerian disc in which the physical conditions would be the same as in the vortex core. Nonetheless, we use it below to make the discussions easier. Condition (8) can be reformulated using the fact that vortices are also pressure and density bumps, namely $P = P_0(1 + C_P)$ and $\sigma = \sigma_0(1 + C_\sigma)^4$. Accounting for a scaling factor of $\sim 1.65$ to better fit the numerical data, we propose an empirical condition for SG-stable vortices:

$$q_0 \ge q_{0,c} = \frac{0.83 h (1 + C_\sigma)^{3/2}}{|\overline{Ro}|_{\text{max}} (1 + C_P)^{1/2}} \quad (\text{hosting disc criterion}), \quad (9)$$

or again,

$$1.2 |\overline{Ro}|_{\text{max}} Q \ge 1 \qquad (\text{vortex core criterion}), \quad (10)$$

where $q_{0,c}$ is the value above which vortices are SG stable and $Q$ is the Toomre parameter estimated in the vortex core. Consistency of this condition with our numerical simulations was tested by computing $q_{0,c}$ at the end of the first relaxation phase (at $t = 50 t_0$). The values of $|\overline{Ro}|_{\text{max}}$, $C_\sigma$, and $C_P$ necessary for the computations are gathered in Table 4, and the values of $q_{0,c}$ are presented in a map plotted in Fig. 6 (right panel). It is striking that the two maps (right and left panels) are nearly parallel to one another in the $(\chi, \delta, q_0)$ space. Applying the above criteria to the vortex of Sect. 3, we find that such a vortex is SG stable

<footnote n="4">⁴ Using definitions exposed in Table 3 we obtain the sound speed in the vortex core, $c_s = c_{s,0} \sqrt{(1 + C_P)/(1 + C_\sigma)}$, where $c_{s,0}$ is the unperturbed sound speed.</footnote>
<table><caption>**Table 4.** Triplets ($|R_o|_{\text{max}}$, $C_\sigma$, $C_p$) for simulated vortices at $t = 50 t_0$.</caption><thead><tr><th>δ<br>χ</th><th>1.0</th><th>1.5</th><th>2.5</th></tr></thead><tbody><tr><td>5</td><td>(0.14, 0.20, 0.31)</td><td>(0.17, 0.50, 0.78)</td><td>(0.21, 1.73, 3.27)</td></tr><tr><td>8</td><td>(0.13, 0.32, 0.48)</td><td>(0.15, 0.72, 1.19)</td><td>(0.19, 2.49, 5.12)</td></tr><tr><td>14</td><td>(0.12, 0.38, 0.57)</td><td>(0.13, 0.82, 1.37)</td><td>(0.16, 2.72, 5.72)</td></tr></tbody></table>

if $Q_0 \ge 10.9$, or $Q \ge 6.9$, a condition that is consistent with our conclusion in the aforementioned section, stating that vortices begin to split even if the hosting disc and the vortex core are stable following the standard Toomre criterion.

### 4.3. Summary

The simple theoretical estimate we used to study the stability of the vortex core against SG allowed us to define a simple criterion based on a modified form of Toomre's criterion. The new criterion was formulated in two different ways (see Eqs. (9) and (10)) to make the comparison easier. The widened numerical experiment we carried out also allowed us to build up a mapping of the SG-parameter, $q_0$, in the $(\chi, \delta)$ space. Vortex stability is thus presented in the form of a simple map of the threshold value $q_{0,c}$ as a function of $\delta$ and $\chi$, selected in the range $0.8 \le \delta \le 3$ and $5.7 \le \chi \le 13.6$, respectively. In light of our numerical data, the stability criterion appears satisfactory but remains a somewhat empirical relation that requires further analysis and comparison. Our simulations showed that, when vortices are unstable with respect to SG, their destruction is preceded by a strong stretching that is likely explained by the gravitational torque exerted by the vortex itself (Regály & Vorobyov 2017b). Finally, this stability criterion could enable to identify upper density limits in hosting discs where a vortex presence is confirmed. This will be the subject of a discussion in Sect. 6.5.1.

Thanks to the threshold parameter $q_{0,c}$ obtained in this section, we were able to choose the suited parameters for studying self-gravitating vortices in greater detail at HR.

## 5. HR simulations

We performed HR simulations to better characterise the internal structure of the self-gravitating vortices identified in Sect. 4. We also studied two different cases, one in which the disc at equilibrium was isothermal and one in which it was not. In Figs. 7 and 8, the Rossby number, density, and compressibility at $t = 140 t_0$ are plotted. For the density plots (left panel of Fig. 8), the main vortex was masked to highlight the weak contribution of other substructures associated with a vortex. Besides the vortex internal structure, other features of interest are the spiral waves and the annular overdensity, which are intrinsically associated with vortices. Indeed, such features are possible clues for theoretical models and observational studies as well. For reference comparison, the case of a vortex in a non-SG and non-isothermal case has been added to these figures. In Fig. 9 the radial and azimuthal profiles of the density and the Rossby number are plotted for a non-isothermal disc in order to show the differences between quasi-steady vortices where SG is present and where it is not. The analysis of local spin inside the vortex is completed by Fig. 10, where the iso-contours of the Rossby number have been plotted for a non-isothermal disc. Finally, for completeness, the SG potential of the vortex and the induced forces in the radial

and azimuthal directions are reported in Fig. F.1 in the case of the non-isothermal disc.

### 5.1. Numerical resolution in the vortex core

Numerical resolution is important for relevant descriptions of SG and gas evolution in the inner region of vortices. In a circular vortex patch around a constant-density core, the enclosed mass is $m \sim \pi (\Delta r)^2 \sigma$. Thus, the gravitational influence of this region of the vortex extends to a mean distance given by Hill's radius, $R_H = r_0 (m/3 M_\odot)^{1/3}$. Then, if the radius of the circular patch reaches $\sim R_H$, the enclosed mass is $m \sim \pi R_H^2 \sigma$ and we get:

$$R_H = \frac{H}{3Q} \qquad (11)$$

For instance, if $Q = 10$, we need to resolve (at least) 30 times the pressure scale height to correctly describe the vortex core. In the foregoing results (cf. Sect. 3), this condition is better satisfied in the radial than in the azimuthal direction. This could contribute to the stronger contraction of the vortex observed in the radial than in the azimuthal direction. Motivated by Eq. (11), we undertook HR simulations, the steps for which are described in next subsection.

### 5.2. Activating SG for HR simulations

From a number of preliminary tests, we learned that the procedure to smoothly activate SG must be adapted to the numerical resolution. Indeed, when activating SG, the vortex evolution is two fold: higher resolution may lead to undesired splitting, while lower resolution, at long timescales, leads to enhanced decay. In fact, an appropriate choice of the time steps and the resolution, at which numerical convergence is reached (see Appendices E.1 and E.2), can avoid spurious evolutions. Thus, a new sequence for SG activation at HR has been defined with the following steps:

$t = 0$ : Gaussian vortex at LR with $(N_r, N_\theta) = (500, 640)$,

$t = 60 t_0$ : SG is smoothly activated,

$t = 120 t_0$ : Resolution increased to (1800, 16 000),

$120 t_0 < t < 140 t_0$ : HR simulation.

The final resolution is equivalent to 146 cells/H and 127 cells/H in the radial and azimuthal directions, respectively. Such a resolution was chosen so that both directions are almost equally resolved, which is important to correctly account for the SG. For obvious reasons of computational costs, the last HR-stage was stopped after only 20 orbits; despite this drawback, results are sufficient to distinguish substructures and to guess a possible steady state in the evolution of the main vortex. The outline for our HR simulations are the following: $r \in [34, 63]$ AU, $h_0 = 0.05$, $r_0 = 50$ AU, and $(\chi, \delta) = (14, 0.92)$ at $t = 120 t_0$. Furthermore, in order to choose the Toomre parameter most suited to obtaining self-gravitating vortices, several values of $q_0$ were tested close to the critical limit provided by our stability Criterion 9. This relation leads to values of $q_{0,c}$ equal to 0.74 and 0.44 for the isothermal and the non-isothermal cases, respectively. Finally, we found that the best critical values of the parameters are $(q_{0,c} = 0.7, \sigma_0 = 3.35 \text{ g cm}^{-2})$ and $(q_{0,c} = 0.425, \sigma_0 = 5.52 \text{ g cm}^{-2})$ in the isothermal and the non-isothermal cases, respectively. For both simulations we get $Q \sim 6$, which means that, in our simulations, the Hill radius of the vortex is resolved at least eight times.
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Rossby number at $t = 140 t_0$ for an isothermal and a non-isothermal disc with SG. Radial zoom with $r \in [46, 54]$ AU. **Top:** self-gravitating vortex in an isothermal disc. **Middle:** self-gravitating vortex in a non-isothermal disc. **Bottom:** reference vortex, without SG and in a non-isothermal disc.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Density and velocity-divergence at *t* = 140 *t*₀ for an isothermal and a non-isothermal disc (with and without SG). **Left column:** relative density (σ/σ₀). The main vortex is masked to distinguish the annular over-density and the spiral waves. **Right column:** normalised relative velocity field divergence (∇·v'/(2Ωₖ)). From top to bottom: self-gravitating vortex in an isothermal disc, self-gravitating vortex in a non-isothermal disc, and the reference vortex (without SG and in a non-isothermal disc), respectively.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Radial and azimuthal profiles of the Rossby number and the relative density at $t = 140t_0$ for a non-isothermal disc. **Top:** radial profile of the Rossby number and the relative density ($\sigma/\sigma_0$), from left to right, respectively. **Bottom:** azimuth profile of the Rossby number and the relative density ($\sigma/\sigma_0$), from left to right, respectively. The dotted green line corresponds to the level $R_0 = -0.08$ used to compute the vortex geometrical parameters.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Iso-contours of the Rossby number in a non-isothermal disc at $t = 140t_0$. **Left:** without SG. **Right:** with SG.</figcaption></figure>

We will now use our HR results to discriminate between vortices evolving in an isothermal or a non-isothermal disc. We will also use the results to better describe both the internal structure of a SG vortex and its co-orbital environment.

### 5.3. Isothermal versus non-isothermal

In both these two cases, SG is present. Comparisons are based on Figs. 7 (the two upper panels) and 8.

#### 5.3.1. Vortex core

In the vortex core, the most striking difference between the isothermal and the non-isothermal cases is the presence of a

rotating spiral in the Rossby number and the density distributions, which both evoke Taylor-Green vortices. This spiral pattern, which barely appears during the LR-stage, tends to strengthen during the HR-stage. It is also associated with the emission of weak compression waves from the core that are spiralling outwards. On the other hand, the outer region that directly surrounds the vortex-core seems to keep a stable elliptical shape. In contrast, in the non-isothermal case, the vortex keeps a smooth and steady elliptical shape with a nearly constant vorticity (or $\overline{R_0} \approx -0.12$). The same spiral structure is present in this case but is very weak and barely visible, likely due to the near incompressibility of the flow in the vortex core. These results suggest that stationary SG vortices have to be found in non-isothermal discs.
### 5.3.2. Annular bump and waves

Gas flowing in the horseshoe region of the main vortex is one of the peculiarities of self-gravitating vortices. This is, indeed, obvious in the two upper panels of Fig. 7. The flow appears rather laminar in the non-isothermal case, but a more complex evolution is possible in the isothermal case, such as observed in the top panel of Fig. 7. This snapshot at $t = 140 t_0$ shows six small vortices in the horseshoe region, located at [50°, 90°, 95°, 135°, 250°, 270°]. These vortices are emitting small amplitude spiral waves with density contrasts lower than 5%, which radially damp in less than 10 AU. The corresponding evolution sequence, available in the online movie 3, shows that these small vortices are not ejected by the main vortex but are instead formed in the shear of the horseshoe region (this effect amplifies with increasing resolution). In the non-isothermal case, the flow structure at $r_0 = 50$ and near $\theta = [95°, 265°]$ could be misinterpreted as eddies (see Fig. 8). Indeed, because of the presence of the massive vortex, the trajectories of small scale eddies are expected to be horseshoe. However, a detailed examination of our data shows that these weak overdensities are stationary in the frame centred on the vortex and don't correspond to vortical structures (see Fig. 7). The above arguments suggest that these two structures correspond instead to transition regions between the annular bump and the vortex itself.

## 5.4. With and without SG

This section focuses on vortices in non-isothermal discs and stresses the difference between self-gravitating and gravity-free vortices. Comparisons are based on Figs. 7 and 8 (the two lower panels), and on Figs. 9 and 10.

### 5.4.1. Vortex

Vortex shape can be estimated from the distribution in space of the vorticity or the density. The mapping and the iso-contours of the Rossby number in Figs. 7 and 10, respectively, clearly show that the vortex shape significantly differs when SG is or isn't taken into account. Under SG, the vortex is found to contract in the radial direction (by a factor of ~1.5) and to remain nearly unchanged in the azimuthal direction. Self-gravitating vortices appear flatter than the gravity-free vortices, which are only constrained by the shear. All these findings are consistent with previous works and particularly with the stretching mechanism revealed by Regály & Vorobyov (2017b). When comparing the radial and azimuthal profiles of the Rossby number and the density in Fig. 9, we find, in addition to the abovementioned vortex flattening, a much steeper increase in the vorticity and a lower density peak correlated with a higher annular bump. When SG is at work, the uniform distribution of vorticity and the strong tightening in the Rossby iso-contours for the external vortex regions suggest that gas compression is more important for self-gravitating vortices, as will be discussed in the next subsection.

Finally, it must be noted that a non-stationary spiral pattern is always present in the density distribution, near the vortex centre. However, this is hardly perceptible.

### 5.4.2. Annular bump and waves

In the gravity-free case, the vortex is commonly associated with an annular bump in the density and the vorticity. On the other hand, when SG is taken into account, the vorticity bump splits

into the two branches of a horseshoe. Thus, in this annular region, the flow is either a simple shear with systematic interactions with the vortex, or a complete horseshoe oscillation between a U-turn on each side of the vortex. This is a situation already observed in the LR simulations of Sect. 3, where secondary vortices coexist with the massive primary, which is an obvious consequence of three-body interactions. Finally, both vortices emit spiral waves with low density contrast (≤7%). Interestingly, when SG is absent, the spiral pattern is made of a rarefaction region located between two compression zones, while when SG is present, it is dominated by compression (Li et al. 2001). The compression of the gas flow in the spiral waves (right panel of Fig. 8) is five times stronger with SG than without SG, likely due to a stiffer connection between the vortex and the background.

## 5.5. Summary on HR results

Our HR simulations enabled detailed comparisons between SG vortices in the isothermal and the non-isothermal cases, but also between gravity-free and self-gravitating vortices in the non-isothermal case. Important conclusions are reached concerning the shape, the evolution, and the co-orbital environment of the vortices. When SG is taken into account: (i) vortices in an isothermal-disc tend to distort and may be easily destabilised in a chain of eddies moving along a horseshoe (a lot of waves are spiralling outwards), and the vorticity distribution of the initial vortex presents a rotating spiral pattern in its core; (ii) vortices in a non-isothermal disc keep a stable elliptical shape until the end of run, and outside the vortex the flow of gas in the co-orbital region is horseshoe shaped. In the case of non-isothermal discs: (i) self-gravitating vortices are flatter than gravity-free vortices (with only a radial contraction), gas compression in spiral waves is stronger (in relation with the vortex distortions), and the vorticity bump splits into the two branches of a horseshoe; (ii) if gravity is neglected, the vortex remains in a steady state and is associated with an annular bump in the density and vorticity distributions. We want to stress that the flow of gas in this annular region can be part of the vortex solution at steady state (cf. Sect. 6.4). The gas compression in the spiral waves is weak.

HR simulations have shown that large-scale vortices reach a quasi-steady state in non-isothermal discs, which is not the case in isothermal discs. The disc temperature seems to play a central role in the dynamical evolution of the vortices, particularly when SG is taken into account. Finally, this study shows that higher resolutions are needed for correctly and accurately studying SG in vortices, but also to describe small-scale structures, such as the ones observed in the isothermal run of Fig. 7, or in future dust trapping simulations.

## 6. Discussion

This section comes back to the strategy used in this paper to simulate the evolution of self-gravitating vortices. We discuss the interest and drawbacks of the method, but also its limitations and the possible improvements. Our results are compared to those obtained by Zhu & Baruteau (2016), Tarczay-Nehéz et al. (2022), and Lin & Pierens (2018), particularly where they concern vortex stability. Finally, we discuss possible observational implications.

### 6.1. Single-vortex strategy

In the absence of SG, previous studies showed the possibility of finding approximate vortex solutions of the stationary Euler
equations. When injected in the numerical simulations, these models relax into long-lived and quasi-steady vortices. This strategy was successful in the case of non-SG discs (Surville & Barge 2015) and is also appealing when SG is taken into account. However, this generalisation requires analytical vortex solutions of equations that include SG. We didn't find analytical solutions in this case, even in an approximate form, but nonetheless decided to keep the same strategy, with the idea of reaching the solution numerically following a succession of steady states. This is the reason for the construction of numerical procedures in Sects. 3 and 5 to get a gradual activation of SG starting from a model of a gravity-free vortex.

### 6.1.1. Learning from HR runs

In fact, the numerical procedures used in Sects. 3 and 5 manage two different effects: a slow variation of SG as a function of time (which does not correspond to a true physical situation), and an increase in the numerical resolution, which necessarily introduces noise in the computations. Of course an increase in the numerical resolution also corresponds to a decrease in the numerical viscosity, which at the small sizes scales as ~$\delta_*^2\Omega$, where $\delta_*$ is the mesh size. Therefore, in the HR case, the relaxation process is longer and often more complex than in the LR case. For example, in a number of HR runs, vortices are observed splitting into smaller vortices. This clearly stresses that numerical resolution is key in our SG-activation procedure.

### 6.1.2. Required time steps for SG activation

Moreover, following our strategy, HR should not be increased before a quasi-steady state is reached. Thus, so that our procedure does not lead to spurious vortex evolutions, a nearly steady state must be reached before any increase in the numerical resolution. In other words, the procedure should satisfy the following sequence: (i) activate smoothly SG at LR, (ii) proceed until a quasi-steady state is reached and (iii) increase resolution.

## 6.2. Instability-generated vortices

Another way to study the evolution of self-gravitating vortices in PPDs is to assume that they are first formed following one of the numerous instabilities known to produce large-scale vortices, such as RWI (Lovelace & Hohlfeld 2013; Lin & Pierens 2018) or the baroclinic instabilities (Klahr & Bodenheimer 2003; Klahr 2004). The constraint, in this case, is that the outline of the computations is imposed by the specific conditions necessary for the growth of the instability. Thus, the vortex evolution is model-dependant, which is not the case in the single-vortex approach developed above.

## 6.3. Comparisons with other works

We now compare our results with those obtained in three different papers. These are Tarczay-Nehéz et al. (2022), Pierens & Lin (2018) and Zhu & Baruteau (2016).

### 6.3.1. Comparison with Tarczay-Nehéz et al. (2022)

The authors performed a similar study to our own using a wider sample (1980 different runs). In contrast to our work, their outline is different since vortices are generated by RWI at the edge of a magnetically dead zone in a viscous $\alpha$-disc model with $\alpha = (10^{-4}, 10^{-5})$, and their simulations had

<table><caption>**Table 5.** Simulation parameters for isothermal runs in Lin & Pierens (2018) at the 500th orbit.</caption><thead><tr><th>f</th><th>Q</th><th>|Ro|</th><th>1.2 |Ro| Q</th></tr></thead><tbody><tr><td>1</td><td>8</td><td>0.12</td><td>1.15</td></tr><tr><td>2</td><td>4</td><td>0.10</td><td>0.48</td></tr></tbody></table>

**Notes.** Here isothermal refers to the convention used in the aforementioned article.

LR with ($N_r, N_\theta$) = (256, 512). In spite of the differences, mainly due to a stronger viscosity (turbulent and numerical, both), the results they found are quite similar to ours with the discrimination of various classes: vortex survival, vortex splitting at various scales, and ring formation.

Following Fig. 5 of their paper, we notice that vortex splitting occurs for $Q \approx 5$; however, to check the relevance or not of the stability criterion reported in our Eq. (10), the corresponding value of Ro is also required. This value is not provided by the authors, but it can be estimated in a simple way using the approximate dependence between Ro and the aspect ratio (Surville 2013). For elongations between five and 20, we have $0.1 \le |Ro| \le 0.2$ and we get $0.6 \le 1.2|Ro|Q \le 1.2$, which is consistent with our criterion in Eq. (9).

### 6.3.2. Comparison with Pierens & Lin (2018)

As in the previous paper, the authors study Rossby vortices in a disc with SG. These vortices are generated with or without $\beta$-cooling at the frontier of the dead-zone, where the alpha-viscosity jumps from $10^{-4}$ inside the dead zone to $10^{-2}$ in the magnetically active region. Their simulations were performed with a radial and an azimuthal resolution at the vortex position equal to 16 cells/H and 10 cells/H, respectively. The obtained vortices are either similar to our SG vortices, $-0.13 \le Ro \le -0.10$ and $5 \le \chi \le 7$, or very different with $Ro = -0.05$ and $\chi \ge 25$. These very weak and elongated structures are likely due to the artificial $\alpha$-viscosity and the low numerical resolution. Care should be paid to the fact that the authors named theirs studies without and with $\beta$-cooling and thermal radiation as isothermal and non-isothermal, respectively. For meaningful comparisons and possible testing of Criterion 10, we focus on their isothermal case and excluded the cases $f = [4, 8]$ because in our simulations, we never encountered vortices with such low Rossby numbers and such high elongations ($Ro = -0.05, \chi \ge 25$).

The estimated vortex parameters and the values of $|Ro|$ and $Q$ are gathered in Table 5 at the 500th orbit, which is the time at which RWI saturated and formed a unique vortex. In this table $f$ is the density scaling factor used by the authors, and the last column shows that our criterion is also comforted by this work, since the two vortices satisfying $1.2|Ro|Q \le 1$ are decaying.

### 6.3.3. Comparison with Zhu & Baruteau (2016)

The authors studied the effect of SG on Rossby vortices but also considered the effect of the indirect gravitational potential of the {star+disc} system. While SG alone tends to delay RWI and inhibit low $m$ modes, they found that the indirect potential tends to strengthen vortices. Such conclusions were also found by Regály & Vorobyov (2017b).

In their Fig. 3 (bottom panels), they report the results of various runs (g5gi, g0p2gi, g2gi and g10gi), which can be fruitfully
<table><caption>**Table 6.** Our SG parameter, $q_0$, corresponding to the density factor, $\Sigma_0$, of Zhu & Baruteau (2016).</caption><thead><tr><th>Σ<sub>0</sub></th><th>q<sub>0</sub></th></tr></thead><tbody><tr><td>0.0002</td><td>15.7</td></tr><tr><td>0.002</td><td>1.6</td></tr><tr><td>0.005</td><td>0.63</td></tr><tr><td>0.01</td><td>0.31</td></tr></tbody></table>

compared to ours. In all these runs, no splitting or vortex destruction is observed; this suggests, following our criterion, that $q_0$ lies above the threshold value $q_{0,c}$. The morphological and SG parameters necessary to test the criterion can be found in Fig. 3 of their paper and in our Table 6; thus, we get $(\chi, \delta) = (6, 0.85)$ and with our Fig. 6, we estimated that $q_{0,c} \le 0.375^5$. The results gathered in our Table 6 are found to be consistent with our criterion ($q_0 \ge 0.375 \ge q_{0,c}$), except in the last line where the inequality is indefinite.

### 6.3.4. Comparison summary

Equations (9) and (10) were compared to the three independent works and exhibit satisfactory results. It is worth noticing that when more physical parameters are taken into account, there should exist a more general stability criterion under the form: $f(Q, |Ro|, h, Re, \tau_c) \ge 1$, where Re is the Reynolds number and $\tau_c$ is the beta-cooling time normalised with respect to the Keplerian period. Finding the general function $f$ could be a very difficult task, if not impossible, but as a first step, it could be really interesting if researchers possessing a large simulation sample and more physical ingredients, such as Tarczay-Nehéz et al. (2022), checked if a simple stability criterion involving power laws, $Q^a|Ro|^bRe^c\tau_c^d\eta^k \ge 1$, can succeed.

## 6.4. Theoretical models for steady vortices in SG-free discs

Simple vortex solutions of the incompressible Euler equations at steady state have been studied by a number of authors (Goodman et al. 1987, hereafter GNG; Kida 1981b; Chavanis 2000). GNG's model provides a rough solution of the Euler and continuity equations for a polytropic<fnref n="6" /> sheared flow. It also indicates that, for a vortex patch of constant vorticity, the pressure and density contours should be ellipses with equal aspect ratios. More generally, a stationary, incompressible, barotropic<fnref n="7" /> vortex solution in a sheared flow should satisfy the requirement that the surface density, pressure, and vorticity are functions of the stream function, which indirectly implies that these three quantities must have the same aspect ratios.

Nonetheless, the nearly steady state reached after the relaxation of the Gaussian model in the present work (non-isothermal and gravity-free run studied in Sect. 5) and in Surville & Barge (2015) suggests the existence of an analytical solution of the steady, compressible Euler equations in the case of non-barotropic fluids. Such a solution should satisfy the following requirements: (i) distributions of pressure, density, and vorticity, which can differ from one another. For instance, we found different aspect ratios for these quantities in the gravity-free runs (i.e. $\chi_{Ro} = 9$, $\chi_{\sigma} = 7.7$, and $\chi_{P} = 7.3$, respectively), (ii) an analytical

way to connect the vortex and the annular bump, (iii) a steady spiral pattern that can account for the vortex wake and (iv) compressibility in the vortex core. This is a difficult task and the necessary developments to improve GNG's model are out of the scope of the present paper.

## 6.5. Possible observational consequences

Although the capture of dust and its concentration in vortices has not been addressed in the present paper, we think that some of our results could already provide interesting clues to test the presence of vortices in PPDs or infer disc properties from the survival of SG vortices. More details on observational signatures of vortices can be found in Regály et al. (2012), Robert et al. (2020) and the references therein.

### 6.5.1. An estimation for HD 163296 surface density

Around this Herbig Ae star of mass $2.3 M_\odot$, a dust clump orbiting at a distance of ~0.3 AU (period ~40 d) was recently detected by the VLTI/MATISSE and GRAVITY instruments (Lopez et al. 2014; GRAVITY Collaboration 2017; Varga et al. 2021). The authors tested different possible observational scenarios that could account for this feature, such as a stellar or a sub-stellar companion, and an inner disc misalignment with the line of sight. Following their conclusions, the most likely mechanism to explain such a dynamic asymmetry is growing dust grains being captured by a gaseous vortex. Indeed, the proximity of the star could produce axisymmetrical step jumps in the density and the pressure that induce a local extrema in the generalised potential vortensity distribution (Li et al. 2000); this is known to be favourable to RWI and the generation of a large-scale vortex. The presence of a vortex seems to be supported by the shadow variations revealed by the Hubble Space Telescope (Rich et al. 2020), since vortices are also pressure bumps whose vertical extent easily exceeds the mean scale height.

The results we get in this paper have potential implications on the disc density that observational predictions could infer. If a large-scale vortex is definitely present in the HD 163296 disc and at $r_0 = 0.3$ AU from the star, it should also be stable against SG forces; thus, following our criterion, we should have $q_0 \ge q_{0,c}$. Of course, we are aware that, in the case of HD 163296, the physical context is at the limit of validity of the method used in the present paper, but we have found it interesting to test this idea nonetheless. For this test the numerical outline is the following: an isothermal gas at $T \sim 1400$ K with a similar molecular weight as H₂ and a vortex located at $r_0 = 0.3$ AU from the star with $\chi \sim 7$ and $\delta \sim 1$. Thanks to previous assumptions, we were able to estimate the disc flaring equal to $h = 0.04$ and $Q_c = 25$. Then, Fig. 6 provides the threshold value $q_{0,c} = 0.375$, which fixes an upper limit to the density:

$$ \sigma_{0,\text{max}}(0.3\text{AU}) = \frac{c_s \Omega_k}{0.375 \cdot 4\pi G Q_c} = 62 \times 10^3 \text{ g cm}^{-2}. $$

Finally, if a vortex is really present, the surface density of the disc should not be larger than $62 \times 10^3 \text{ g cm}^{-2}$ at 0.3 AU to avoid being destroyed by SG splitting. This is, of course, a very crude approach to the problem, but our goal was only to stress that new constraints can emerge when SG is taken into account.

### 6.5.2. Vortex-induced spiral waves

In the absence of SG, vortices are known to excite spiral waves that cause them to migrate inwards towards the star. These waves,

<footnote n="5">⁵ The map of Fig. 6 is valid only in the range $0.8 \le \delta \le 3$.</footnote>

<footnote n="6">⁶ $P = A\sigma^{1+1/n}$.</footnote>

<footnote n="7">⁷ $P = P(\sigma)$.</footnote>
in contrast to the spiral density waves produced in disc or planet
interactions, originate from the compression and expansion of
the gas layers sandwiched between the vortex core and the back-
ground flow. Following Huang et al. (2019), the density contrast
of these waves (lower than 20%) is much too weak to be observed
in scattered light.

On the other hand, when SG is taken into account, the mass
of the vortex also contributes to the wave emission with the
excitation of spiral density waves. Of course, the more massive
the vortex, the higher the density contrast of the waves. In fact,
the maximum mass a vortex can reach is limited by its stabil-
ity against SG, which is constrained according to Criteria (9)
or (10). Anyway, it can be concluded that spiral waves are always
much weaker if excited by a gaseous vortex than if excited by a
planet. The density contrast of the waves associated with the vor-
tex remains very weak (see Figs. 8 and A.1), and the only way
to obtain stronger waves is to have large amount of solid mate-
rial trapped in the core of the vortex, mimicking the effect of a
planet.

6.5.3. Annular substructures

Vortices in PPDs are known to easily capture solid particles with
optimal Stokes number (St ~ 1) and to confine them in the core
(Barge & Sommeria 1994, 1995; Tanga et al. 1996; Bracco et al.
1999). As the vortices are always associated with an annular
bump (see Fig. 8), we have found it interesting to explore the
possible evolution of the solid particles in this annular region
during the trapping-in-vortex stage. Indeed, dust concentrations
in this region could be observational markers of the vortices. As
noted in Sect. 5, this annular region could be a simple annular
pressure and density bump in a nearly Keplerian shear, for weak
SG, or a more complex region with a horseshoe flow and pos-
sible substructures, for intermediate SG. Since pressure bumps
are possible sinks and vortical substructures are potential traps,
we expect that solid particles with St ≤ 1 should form a ring
strewed with various clumps. These clumps may have horseshoe
motions. Finally, these annular bumps, which can trap a non-
negligible amount of solid particles, can be considered as a tank
of intermediate-sized particles, which possibly ‘feed’ the vortex
during its evolution.

More precisely, since the density amplitudes between the
vortex core and the ring differ, we expect that a segregation
with respect to dust size should occur, which implies that dust
asymmetries and rings located at the same radial position should
be observable at different wavelengths in scattered light. Here,
we conjecture that this mechanism of gap cleaning and the
ring distribution of dust induced by vortices are due to the
small amplitude overdensity ring, which always goes along a
vortex. This is out of the scope of current paper, but would
be an interesting work to explore and to be compared with
observations.

7. Conclusions and perspectives

In this paper we reported the evolution, in an inviscid disc, of vor-
tices under the effect of SG thanks to HR numerical simulations.
As an initial state, vortices were injected considering a Gaussian
model (Surville & Barge 2015) and relaxed by numerical viscos-
ity. Our study shows that, in the presence of SG, vortices evolve
following three possible regimes:

**Shear dominated.** vortices are modelled by the shear but stay nearly unchanged with respect to the non-SG case.

**Self-gravitating vortices.** vortices stretch out, contract radi-
ally, and have a core spinning slightly faster. The main vortex
eventually releases small secondary vortices.

**Splitting cascade.** the initial vortex is gradually destroyed and, after a few hundred orbits, spreads into a residual annular bump that eventually contains small-scale vortex remnants; this can occur even if Toomre's parameter $Q_0 \sim 10$.

From these simulations it seems that for self-gravitating
vortices, vortex strength is not necessarily proportional to the
aspect ratio inverse. This should be confirmed by other inde-
pendent works. The other conclusions we get concern vortex
stability, particularly in regard to numerical resolution and disc
temperature. The most interesting points are the following:

**Approximate stability criterion.** from an extended numerical
study of vortices submitted to their own gravity, and after testing
with simple analytical arguments, we propose a new criterion,
which states that in inviscid discs, vortices resist SG destabilis-
ing effects if $1.2|Ro|_{max}Q \ge 1$, or again if $q_0 \ge q_{0c}$. This criterion
is supported by comparison with the results of three indepen-
dent studies (Tarczay-Nehéz et al. 2022; Zhu & Baruteau 2016;
Pierens & Lin 2018). In addition, the aforementioned criteria
were also applied to a possible vortex candidate in HD 163296
disc (Varga et al. 2021) which allowed an estimation to be made
of the upper limit for the gas density.

**Numerical resolution.** In the absence of viscosity, we found
that simulations with at least (71,24) cells/$H$ in the radial
and azimuthal directions, respectively, are needed to accurately
describe vortices. If such a requirement is not fulfilled, vortices
lose 15% of their mass during 250 orbits (see Appendix E.1)
and this trend is accentuated in the presence of SG. On the
other hand, simulations with ($N_r, N_\theta$) = (1800, 16 000), have
shown that numerical resolution is key to correctly and finely
describe the shape and evolution of SG vortices. Nearly uni-
form resolution, in $r$ and $\theta$, is necessary and simple estimates
show that simulations should resolve distances of the order of
~$H/(3Q)$.

**Disc temperature.** Comparisons between the HR simula-
tions of Sect. 5 clearly show that SG vortices in non-isothermal
discs: (i) evolve more slowly than vortices in isothermal discs,
and (ii) keep a smooth elliptical shape throughout the run. On
the other hand, vortices in isothermal discs are distorted by
the SG forces and vortical substructures are generated in the
horseshoe region; they also excite more compressive waves than
in the case of non-isothermal discs. The key point here is that
a temperature gradient of the disc could be necessary for SG
vortices to reach a nearly steady-state.

Zhu & Baruteau (2016), Regály & Vorobyov (2017b) and
Tarczay-Nehéz et al. (2022) found that vortices cannot form in
massive discs since lower mode numbers are inhibited by SG.
Although the strategy in this paper was different, our results are
consistent with their findings: not only can vortices not form in
massive discs, but even if they are well established (as the ini-
tial state in our simulations), they cannot be sustained. Thus, if
observed crescent-shaped asymmetries are indeed vortex pres-
ence tracks, these vortices should be hosted in non-massive
discs. However, in light of self-gravitating vortices stability
(as defined in this work), the definition of non-massive disc
should be clarified with respect to the other physical ingre-
dients governing the dynamics, such as viscosity or thermal
cooling. In a roundabout way, we provided an attempt of this
definition for inviscid discs through the unique Eq. (9) (or Eq. (10)).

Finally, the results of this paper also provide useful input to prepare dust-trapping simulations in permanent self-gravitating vortices. For example, these simulations should be carried out in non-isothermal discs in which $q_0 \ge 1.2 q_{0,c}$ minimises the destabilising effects of SG, in anticipation of their possible enhancement by the dust back-reaction; but also a $q_0$ value relatively close to $q_{0,c}$ would imply non-negligible effects from SG in the gas and dusty phases.

*Acknowledgements.* We thank the referee, Zsolt Regály, for his helpful comments which improved the quality of this work. We would like to warmly acknowledge Stéphane Le Dizès for fruitfull discussions during the preparation of the paper and also for his help to find funding of part of the project. This research has made use of computing facilities operated by CeSAM data center at LAM, Marseille, France. The Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high performance computing resources. This work was granted access to the HPC resources of CINES/IDRIS under the allocation 2020-A0090411982 made by GENCI. This work was granted access to the HPC/AI resources of CINES under the allocation 2018-A0050407407 made by GENCI.

## References

Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41

Barge, P., & Sommeria, J. 1994, in *Circumstellar Dust Disks and Planet Formation*, eds. R. Ferlet, & A. Vidal-Madjar (Berlin: Springer), 10, 295

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

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

Barranco, J. A., & Marcus, P. S. 2005, ApJ, 623, 1157

Baruteau, C., Barraza, M., Pérez, S., et al. 2019, MNRAS, 486, 304

Blum, J., & Wurm, G. 2008, ARA&A, 46, 21

Bracco, A., Chavanis, P. H., Provenzale, A., & Spiegel, E. A. 1999, *Phys. Fluids*, 11, 2280

Cazzoletti, P., van Dishoeck, E. F., Pinilla, P., et al. 2018, A&A, 619, A161

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

de Val-Borro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529

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

Dong, R., Liu, S.-y., Eisner, J., et al. 2018, ApJ, 860, 124

Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971

Frigo, M., & Johnson, S. G. 2005. Special issue on “Program Generation, Optimization, and Platform Adaptation” *Proceedings of the IEEE*, 93, 216

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

Fu, W., Li, H., Lubow, S., Li, S., & Liang, E. 2014b, ApJ, 795, L39

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

Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051

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

GRAVITY Collaboration (Abuter, R., et al.) 2017, A&A, 602, A94

Hammel, B., & Sullivan-Molina, N. 2020, https://doi.org/10.5281/zenodo.3723294

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

Hayashi, C. 1981, *Prog. Theor. Phys. Suppl.*, 70, 35

Huang, P., Dong, R., Li, H., Li, S., & Ji, J. 2019, ApJ, 883, L39

Huré, J. M., Hersant, F., Carreau, C., & Busset, J. P. 2008, A&A, 490, 477

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

Inaba, S., Barge, P., Daniel, E., & Guillard, H. 2005, A&A, 431, 365

Jacquet, E., Balbus, S., & Latter, H. 2011, MNRAS, 415, 3591

Johansen, A., & Youdin, A. 2007, ApJ, 662, 627

Johansen, A., Henning, T., & Klahr, H. 2006, ApJ, 643, 1219

Johansen, A., Oishi, J. S., Mac Low, M.-M., et al. 2007, *Nature*, 448, 1022

Kida, S. 1981a, J. Fluid Mech., 112, 397

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

Klahr, H. 2004, ApJ, 606, 1070

Klahr, H. H., & Bodenheimer, P. 2003, ApJ, 582, 869

Klahr, H., Pfeil, T., & Schreiber, A. 2018, *Handbook of Exoplanets, Instabilities and Flow Structures in Protoplanetary Disks: Setting the Stage for Planetary Formation*, eds. D., Hans J., & B., Juan Antonio (Berlin: Springer), 138

Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32

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

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

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

Lin, M.-K., & Pierens, A. 2018, MNRAS, 478, 575

Lopez, B., Lagarde, S., Jaffe, W., et al. 2014, *The Messenger*, 157, 5

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., Johansen, A., Klahr, H., & Piskunov, N. 2009, A&A, 493, 1125

Manger, N., Klahr, H., Kley, W., & Flock, M. 2020, MNRAS, 499, 1841

Manger, N., Pfeil, T., & Klahr, H. 2021, MNRAS, 508, 5402

Marcus, P. S., Pei, S., Jiang, C.-H., et al. 2015, ApJ, 808, 87

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

Meheut, H., Lovelace, R., & Lai, D. 2013, MNRAS, 430, 3

Miranda, R., Li, H. Li, S., & Jin, S. 2017, ApJ, 835, 118

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

Müller, T. W. A., Kley, W., & Meru, F. 2012, A&A, 541, A123

Paardekooper, S.-J., Lesur, G., & Papaloizou, J. C. B. 2013, *Eur. Phys. J. Web Conf.*, 46, 05001

Pérez, L. M., Benisty, M., Andrews, S. M., et al. 2018, ApJ, 869, L50

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

Petersen, M. R., Stewart, G. R., & Julien, K. 2007b, ApJ, 658, 1252

Pierens, A., & Lin, M.-K. 2018, MNRAS, 479, 4878

Raettig, N., Klahr, H., & Lyra, W. 2015, ApJ, 804, 35

Raettig, N., Lyra, W., & Klahr, H. 2021, ApJ, 913, 92

Regály, Z., & Vorobyov, E. 2017a, A&A, 601, A24

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

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

Regály, Z., Sándor, Z., Csomós, P., & Ataiee, S. 2013, MNRAS, 433, 2626

Rendon Restrepo, S., Barge, P., & Vavrik, R. 2022, A&A, in press, https://doi.org/10.1051/0004-6361/202243518

Rich, E. A., Wisniewski, J. P., Sitko, M. L., et al. 2020, ApJ, 902, 4

Richard, S., Barge, P., & Le Dizès, S. 2013, A&A, 559, A30

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

Robert, C. M. T., Méheut, H., & Ménard, F. 2020, A&A, 641, A128

Rometsch, T., Ziampras, A., Kley, W., & Béthune, W. 2021, A&A, 656, A130

Safronov, V. 1972, Evolution of the protoplanetary cloud and formation of the earth and the planets, NASA technical translation (Israel Program for Scientific Translations)

Soon, K.-L., Momose, M., Muto, T., et al. 2019, PASJ, 71, 124

Squire, J., & Hopkins, P. F. 2018, MNRAS, 477, 5011

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

Surville, C. 2013, Ph.D. Thesis, Aix-Marseille University, ed. P. Barge, Marseille, France

Surville, C., & Barge, P. 2015, A&A, 579, A100

Surville, C., Mayer, L., & Lin, D. N. C. 2016, ApJ, 831, 82

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

Tarczay-Nehéz, D., Regály, Z., & Vorobyov, E. 2020, MNRAS, 493, 3014

Tarczay-Nehéz, D., Rozgonyi, K., & Regály, Z. 2022, MNRAS, 511, 6055

Teed, R. J., & Latter, H. N. 2021, MNRAS, 507, 5523

Terquem, C. E. J. M. L. J. 2008, ApJ, 689, 532

Toomre, A. 1964, ApJ, 139, 1217

Tsukagoshi, T., Muto, T., Nomura, H., et al. 2019, ApJ, 878, L8

van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2013, *Science*, 340, 1199

Varga, J., Hogerheijde, M., van Boekel, R., et al. 2021, A&A, 647, A56

Varnière, P., & Tagger, M. 2006, A&A, 446, L13

Weidenschilling, S. J. 1977, MNRAS, 180, 57

Weidenschilling, S. J., & Cuzzi, J. N. 1993, in *Protostars and Planets III*, ed. E. H. Levy, & J. I. Lunine (Tucson: University of Arizona Press), 1031

Whipple, F. L. 1972, in *From Plasma to Planet*, ed. A. Elvius (Hoboken: John Wiley & Sons Inc), 211

Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459

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

Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57
## Appendix A: Full simulation box and emitted spiral

Figure A.1 presents the divergence of the velocity field normalised with respect to the gas rotation rate in order to show the full extent of the simulation box used in Section 3.

<figure><img src="image_12.png" /><figcaption>Fig. A.1: Compressibility waves, $\nabla \cdot v' / (2\Omega_k)$, emitted by a quasi-steady vortex in its full simulation box, $r \in [25, 72]$ AU, as used in Section 3.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. B.1: Vortex stability in the $(\chi, \delta, q_0)$ space. Each point is an intermediate resolution simulation: $(N_r, N_\theta) = (1500, 3000)$. $\chi$ and $\delta$ are, respectively, the vortex aspect ratio and the radial extent at $t = 50 t_0$.</figcaption></figure>

## Appendix B: Stable and unstable vortex structures with respect to Toomre's parameter

This section refers to the definition of stability provided in Sect. 4.1. Figure B.1 illustrates for each triplet $(\chi, \delta, q_0)$ the vortex stability with respect to SG. In order to improve results, three additional simulations were performed for the following triplets: $(\chi, \delta, q_0) = [(12.6, 2.7, 0.8), (12.6, 2.7, 0.65), (9.00, 3.06, 0.6)]$. The blue surface of Figure 6 (left panel) was obtained estimating the transition zone between stable and unstable vortices. Below are the prescriptions used to estimate, for a given couple $(\chi, \delta)$, the transition, in the $q_0$ direction, from the stable to the unstable region.

* No blue cross in the $q_0$ direction indicates that the transition between stable and unstable vortices is located in the mid-value between the lowest green triangle and the upper red dot.

* A blue cross in the $q_0$ direction indicates that the transition between stable and unstable vortices is the lowest blue cross.

## Appendix C: Computed quantities

This section is devoted to introducing in detail the way in which we computed different quantities, such as the distances, the vortex mass, the mean Rossby number, and the vortex geometrical parameters. We begin describing the way in which we estimated the vortex position.

### Appendix C.1. Main vortex position and distance to the star

We found that even in presence of secondary vortices the position of the main vortex is given by the pressure maximum. For

the above reason, the distance to the star and the angular position of the vortex are computed by tracking the position of the maximum of pressure contrast $C_p = \frac{P}{P_0} - 1$. We denote the vortex angular position by $\theta_{max}$.

### Appendix C.2. Mass of the vortices

For comparison of the vortices studied in Section 3, it was decided to define the core ‘mass’ of a vortex as the region where the density contrast is higher than a given threshold $\mu$. In order to choose appropriately this threshold we plotted in Figure C.1 contours for different values of $\mu$ and provided the associated vortex masses in the caption. Of course, the higher the threshold, the smaller the captured mass that leads to the lowest (preferred) value for $\mu$. However, a lower threshold, such as $\mu = [0.05, 0.15]$, is not convenient since the numerical filter captures the annular overdensity and a small share of the spiral waves emitted by the vortex (see top panels of Figure C.3). In order to overcome this issue, we decided to set the threshold parameter at a value of $\mu = 1.25 - 1 = 0.25$. Yet even with this clip, another problem arises because of the presence of secondary vortices (see cases $q_0 = [0.25, 0.5]$ in Section 3 or the vortex hosted in an isothermal disc in Section 5). We get rid of them by applying a second filter that extracts only the region located at a maximal angular position of $45^\circ$ from the pressure maximum. The bottom panels of Figure C.3 show how the density is filtered without (left) and with (right) this second filter. Despite our best efforts, the reader should keep in mind that this second filter is not optimal when secondary vortices enter the region $|\theta - \theta_{max}| \le 45^\circ$. In summary, the mass of a vortex core is computed numerically in our code as:

$$m_{\text{core}} = \int_{\tau} \sigma(r, \theta) r dr d\theta$$

where $\tau = \{(r, \theta) \in \mathbb{R}^+ \times [0, 2\pi] : C_\sigma \ge 0.25 \text{ and } |\theta - \theta_{max}| \le 45^\circ\}$.
<figure><img src="image_14.png" /><figcaption>Fig. C.1: **Mass of the vortex core with respect to the μ threshold at t = 158 t₀.**</figcaption></figure>

μ + 1 = [1.25, 1.40, 1.50, 1.60] (from the exterior to the vortex core). The respective vortex-core masses are 0.18, 0.1, 0.06, and 0.02 (in M<sub>J</sub>). The initial vortex parameters are: (χ, δ) = (14, 1.5) (t = 0) and q<sub>0</sub> = 0.5.

<figure><img src="image_15.png" /><figcaption>Fig. C.2: **Definition of aspect ratio, χ, and radial extent, Δr, fitting the Ro = -0.8 iso-contours.**</figcaption></figure>

In the (r, θ) coordinate system, Rossby iso-contours are ellipses, while in the Cartesian coordinate system, ellipses are curved.

### Appendix C.3. Mean Rossby number

Similarly to the previous subsection, we applied a clip to the Rossby number in order to extract the 'vortical' core region: $Ro \le -0.8$. The mean Rossby number is computed numerically in the extracted region:

$$ \overline{Ro} = \left( \int_{\Gamma} R o r dr d\theta \right) / S, $$

where $S = \int_{\Gamma} r dr d\theta$ is the total surface of the extracted region and $\Gamma = \{(r, \theta) \in \mathbb{R}^{+} \times [0, 2\pi] : Ro \le -0.8 \text{ and } |\theta - \theta_{max}| \le 45^{\circ}\}$.

### Appendix C.4. Aspect ratios and radial width

Aspect ratios were computed thanks to the python routine `Least Squares` fitting of ellipses (Hammel & Sullivan-Molina 2020), which provides the ellipse parameters (semi-major and semi-minor axes, inclination) which better fits a given set of points.

<table><caption>Table E.1. Resolution for numerical convergency tests</caption><thead><tr><th rowspan="2">Name</th><th colspan="2">Sections 3 and 4</th><th colspan="2">Section 5</th></tr><tr><th>(N<sub>r</sub>, N<sub>θ</sub>)</th><th>cells/H</th><th>(N<sub>r</sub>, N<sub>θ</sub>)</th><th>cells/H</th></tr></thead><tbody><tr><td>Res1</td><td>(240, 320)</td><td>(11, 2.6)</td><td>(1500, 3000)</td><td>(122, 24)</td></tr><tr><td>Res2</td><td>(500, 720)</td><td>(24, 5)</td><td>(1400, 11520)</td><td>(113, 92)</td></tr><tr><td>Res3</td><td>(750, 1536)</td><td>(34, 12)</td><td>(1800, 16000)</td><td>(146, 127)</td></tr><tr><td>Res4</td><td>(1500, 3000)</td><td>(71, 24)</td><td>(2200, 19840)</td><td>(178, 158)</td></tr></tbody></table>

Since this technique is based on the least square method, intrinsically it cannot provide uncertainties, but an estimation was performed in Appendix D thanks to the grid resolution. Figure C.2 shows how contour plots were performed for the Rossby number, how these contours were fitted with the best ellipse, and the definitions of radial width, Δr, and aspect ratio, χ.

However, when vortices start to split in many substructures, it is not possible to discriminate the main vortex using Rossby number contours. Therefore, we applied the same filter used in Sects. C.2 and C.3 in order to isolate the main vortex. We decided to use Rossby iso-contours, for tracking geometrical parameters, since the extensive study achieved in Section 4 showed that the only common quantity shared by all 45 vortices was the Ro=-0.08 contours, while density and the pressure maximum amplitude vary for all different structures, which makes comparison a difficult task.

## Appendix D: Uncertainties

In numerical experiments the accuracy arises from the used numerical method. In particular, for a finite volume method, it depends on the time step and on the approach used for computing fluxes between cells. Uncertainties on the computed quantities arising from finite volume method are difficult to quantify, which encouraged us to consider that uncertainties of computed quantities from raw data arise only from the grid resolution.

* **Distances:** Uncertainties related to distances (distance to the star and vortex radial width) are given by the maximum value between the resolution in the radial or azimuthal direction: $\delta_* = \text{Max}(\delta_r, r\delta_\theta)$.

* **Aspect ratios:** The aspect ratio of the fitted ellipse (see Appendix C.4) is defined as $\chi = \frac{a}{b}$, where a and b are, respectively, the semi-major and semi-minor axes. Thus, the uncertainty is: $\Delta\chi = \frac{\chi}{b} (\Delta a + \Delta b) = \frac{\delta_*}{b} (1 + \chi)$.

## Appendix E: Numerical convergency

Our simulations were tested against numerical convergency. We gathered in Table E.1 the resolutions we used for the tests related to Sections 3, 4 and 5.

### Appendix E.1. Comments for Sections 3 and 4

As an example, we show in Figure E.1 the evolution of the mass and the Rossby number for a vortex without SG for four different resolutions (case $q_0 = \infty$ of Section 3). We observe that, for the two lowest resolutions, there is a similar mass loss of 15% in only 250 orbits and a mean Rossby number increase of ~6%. If only looking at Res1 and Res2, it seems that simulations are in the numerical convergent regime. However, for Res3 there is less mass decay (only 7%) and the mean Rossby number evolution is the same as Res4 up to the 115th orbit. After
<figure><img src="image_16.png" /><figcaption>Fig. C.3: Vortex masses with respect to different density thresholds and filters.</figcaption></figure>

Top and Bottom left panels: We applied only a clip with respect to the density threshold.

Bottom right panel: We applied an additional filter that extracts only the region located at a maximal angular position of 45° from the pressure maximum (see Appendix C.2).

Top left: $\mu + 1 = 1.05$, $m_{core} = 1.69M_J$. Top right: $\mu + 1 = 1.15$, $m_{core} = 0.78M_J$

Bottom left: $\mu + 1 = 1.25$, $m_{core} = 0.28M_J$. Bottom right: $\mu + 1 = 1.25$, $m_{core} = 0.18M_J$

that the Rossby number increases. The situation is completely different for the highest resolution, Res4, where there is almost no mass loss and the Rossby number remains stationary. Even if not shown here, we highlight that when SG is included, vortices decay is stronger compared to simulations without SG, at LR. The above mass losses, with and without SG, are due to a density contrast decrease, while the vortex parameters don't change.

We conclude that: (1) working only with the lowest resolutions can be misleading, and we could get the impression that the numerical scheme is convergent when it is not; and (2) our highest resolution for Section 3 is in the convergent regime and leads to a quasi-stationary vortex. Based on these conclusions we suggest that numerical simulations on vortices should resolve at

least ~ 70 and 25 times the pressure scale height in the radial and azimuthal directions, respectively. If this condition is not satisfied, spurious decay of the vortices is likely.

## Appendix E.2. Comments for Section 5

We show in Figure E.2 the evolution of the aspect ratio, the mass and the Rossby number for the self-gravitating vortex, hosted in a non-isothermal disc and studied in Section 5. Compared to the previous subsection, we show an additional plot exhibiting the Rossby-number aspect ratio at the $R_o = -0.01$ level. This allows us to check how the vortex-core morphology is affected when
<figure><img src="image_17.png" /><figcaption>Fig. E.1: **Numerical convergency for Sections 3 and 4**</figcaption></figure>

*Top*: Convergence of the vortex mass.

*Bottom*: Convergence of the Rossby number.

All quantities were normalised with their respective value at $t = 50 t_0$.

resolution exceeds the threshold $R_H = H/3Q$. Convergence tests were performed, but were strongly time-consuming; they were carried out only during eight orbital periods. Starting from the $130^{th}$ orbit as an initial state, we performed a 2D cubic interpolation to get the desired resolution. We highlight that for *Res1*, the distance $H/3Q$ is well resolved only in the radial direction, while in the azimuthal direction it is resolved only 1.3 times. For the three other, resolutions above quantity is resolved at least five times in both directions. We observed that, for all resolutions, there is no notable difference in the vortex mass and in the evolution of the mean Rossby number. *Res2*, *Res3*, and *Res4* share the same aspect ratio evolution, with a small divergence from the $136^{th}$ orbit. The situation is different for *Res1*; even if the tendency is the same as for HR simulations, we notice from the beginning that there is a 2% difference. This suggests that vortex-core morphology is affected by SG when small scales are resolved. In conclusion: (1) vortex core morphology is affected when $H/3Q$ is well resolved in both directions; and (2) our simulations fit the convergent regime.

## Appendix F: Potential and forces due to SG

Figure F.1 shows, at $t = 140t_0$, the spatial distribution of the potential and the force components of SG for the non-isothermal run of Section 5.

<figure><img src="image_18.png" /><figcaption>Fig. E.2: **Numerical convergency for Section 5**</figcaption></figure>

*Top*: Vortex mass evolution. *Middle*: Rossby number evolution.

*Bottom*: Aspect ratio evolution.

All quantities were normalised with their respective value at $t = 130 t_0$.
<figure><img src="image_19.png" /><figcaption>Fig. F.1: **SG: potential and forces at t = 140 t₀.** Left: SG potential ($\Phi_{SG}$) and iso-contours. Middle: SG radial force ($f_{r,SG}$). Right: SG azimuthal force ($f_{\theta,SG}$).</figcaption></figure>