Issue 
A&A
Volume 666, October 2022



Article Number  A92  
Number of page(s)  23  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202243518  
Published online  13 October 2022 
Morphology and dynamical stability of selfgravitating vortices
Numerical simulations^{★}
^{1}
AixMarseille Univ., CNRS, CNES, LAM,
Marseille, France
email: steven.rendonrestrepo@lam.fr; pierre.barge@lam.fr
^{2}
AixMarseille Univ., CNRS, Centrale Marseille, IRPHE,
Marseille, France
Received:
10
March
2022
Accepted:
20
June
2022
Context. Theoretical and numerical studies have shown that largescale 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 submillimetre 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 selfgravity (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 selfgravitating disc must satisfy in order to permit vortex survival at long timescales. Finally, we characterise as well as possible the persistent selfgravitating vortices we have found in isothermal and nonisothermal 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 nonhomentropic and nonadiabatic flow for an ideal gas. A series of 45 runs were carried out starting from a Gaussian vortexmodel; 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 smallscale vortices seem to be more robust than largescale 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 selfgravitating vortices that we found persist during hundreds of rotation periods and look like the quasisteady vortices obtained in the nonselfgravitating case. A number of these selfgravitating 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 selfgravitating vortices in protoplanetary discs. Vortices reach a nearly steadystate more easily in nonisothermal discs than in isothermal discs.
Key words: protoplanetary disks / hydrodynamics / methods: numerical / planets and satellites: formation
Movies are only available at https://www.aanda.org
© S. R. Restrepo and P. Barge 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1 Introduction
The formation of planetesimals from the dustgrains embedded in the gas of protoplanetary discs (hereafter PPDs) remains one of the key problems of planet formation. Standard scenarios start from the sublayer that the dust forms in the midplane 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 selfgravitating dustclumps that are thought to compact into kilometresized planetesimals (Safronov 1972; Goldreich & Ward 1973). However, these models are known to stumble over a number of ‘barriers’ (metresize, 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 dustclumps to deliver selfgravitating dustclouds that are dense enough to form big planetesimals (up to Ceressized 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 runaway gas accretion before the solid material is accreted by the star and the gas is dissipated by winds and photoevaporation.
Another scenario proposes that anticyclonic 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 dusttogas ratio in the core of the vortices (this massratio ~1 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 hydrodynamical 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^{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 ValBorro 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 highresolution 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 MidInfraRed Instrument (MIRI) instrument in the thermal midIR region) or the future Extremely Large Telescope (ELT; using HIghREsolution 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 KelvinHelmoltz instability or specific resonantdrag 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 dustladen 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 backreaction of the dust particle becomes problematic at a high dusttogas ratio. 3D simulations are scarce and seem to show that vortices remain stable and behave as Taylor columns, even if the midplane region is perturbed by a thin dustlayer (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 3Dvortices (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 selfgravity (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 twodimensional studies: Lovelace & Hohlfeld (2013) have shown that SG can affect RWI only if the Toomre parameter is less than the inverse of the nondimensional 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 (TarczayNehéz et al. 2020). On the other hand, threedimensional 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, highresolution (HR) simulations were carried out to characterise as well as possible the internal structure of a selfgravitating vortex hosted in an isothermal and a nonisothermal 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 quasistationary 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 nonisothermal 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.
2.1 Disc Model and Evolution Equations
We investigate the evolution of an ideal gas (molecular dihydrogen, µ ~ 2.34) in a thin (2D), nonhomentropic disc orbiting a solartype 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: (1) (2) (3)
where P, σ, and v are the vertically integrated pressure, density, and the gas velocity, respectively, is the internal energy of an ideal gas, and γ = 1.4 is the adiabatic index for a diatomic gas. The gravitational potential has three different contributions: Φ_{⊙}, Φ_{SG}, and Φ_{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 nonactive region of the PPDs, although Stoll & Kley (2014) noted that, even in this region, verticalshear 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 α/(10^{−4}) ≪ 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 nonnegligible 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 σ_{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 β_{σ} and the temperature index ß_{T} to account for the observational constraints. With these assumptions, the unperturbed axisymmetric disc’s Toomre parameter, Q_{0}(r) = Ω_{k}c_{s0}/(π_{G}σ_{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: (4)
where r is the normalised radius at 1 AU, r_{0} = 50, σ_{0}(r) and σ_{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 power law with index ß_{p} = ß_{T} + β_{σ}. In Sects. 3 and 4, and for the nonisothermal disc studied in Sect. 5, we have set ß_{T} = −0.5. On the other hand, we have chosen ß_{T} = 0 for the isothermal disc considered in Sect. 5.
SG parameter and surface density at r_{0} = 50 AU.
2.2 Gravity Terms
In polar coordinates the SG potential reads: (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 anticyclonic 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: (6)
where and .
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} ≤ 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 σ_{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 Bifluids) 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 highperformance parallelism; it is presented in detail elsewhere (Rendon Restrepo et al. 2022) with many tests 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 & Papaloizou 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 longlived 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 = 50t_{0}, the time at which SG is plugged in. It is interesting to note that in presence of viscosity, α = 10^{−4}, and simulations with or without SG have shown that vortices decay at a rate proportional to the disc mass (Regály & Vorobyov 2017b).
Vortex parameters at t = 50 t_{0} with respect to the initial Gaussian model parameters.
2.3.1 A Smooth Activation of SG
In the presence of SG, numerical tests (at lowresolution, hereafter LR) have shown that persistent vortices can also survive, but after a transition period that is longer and more complex than in the nonselfgravitating case. This transient evolution is due to the difference between the injected and the exact vortexsolutions; 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: (7)
The time constant was chosen to be equal to the vortex relaxation time in the nonSG case; that is T_{act} ≃ 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 selfgravitating 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_{θ}) = (360, 720), during 150 orbits and for a large computational window (r ∈ [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 ∈ [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 quasisteady 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 bulkmass, 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 χ and its radial extent (scaled to the pressure scale height), δ = Δ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 isocontours of the Rossby number at the Ro = −0.08 level, where Ro = ∇ × v′/2Ω_{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 bulkvortex evolution once a nearly steadystate is reached, such as a longterm 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 longlived and robust against instabilities, and for this reason we computed the mean Rossby number in the vortex core or spin, , 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/χ), the general trend is the same with a saturation for elongated vortices: Ro ≃ −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 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 overdensity (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.
Parameter definitions.
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 from a Gaussian vortex with aspect ratio χ = 14 and radial extent δ = 1.5 which was injected at r_{0} = 50 AU and θ_{0} = π. 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 ∞ (nonSG 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 endofrun 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} = ∞ 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.
Fig. 1 Gaussian vortex model injected at t = 0. Top: Rossby number. Bottom: normalised density (σ(r, θ)/σ_{0}(r)). 
3.1 Numerical Procedure for the Run Series
Computations started at t = 0 with a LR defined by (N_{r}, N_{θ}) = (500, 720). During the first 50 orbits, the vortex morphological and global parameters did not stay constant but 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 quasistationary regime (see Table 2). Then, at t = 50 t_{0}, the resolution of the simulations was increased using a 2D cubic interpolation^{2} into an intermediate grid resolution of (N_{r}, N_{θ}) = (1500, 3000), which is equivalent to ~71 cells/H and ~24 cells/H in the radial^{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_{θ}) = (500, 720),
0 < t < 50 t_{0} : Vortex relaxes to a quasisteady state,
t = 50 t_{0} : Resolution is increased to (1500, 3000),
t = 60 t_{0} : Plugin 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 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 ∈ [0, 300 t_{0}], we found that the vortex massloss is ~15% for every 250 orbits, if the resolution is (N_{r}, N_{θ}) = (500, 720), while massloss is negligible if resolution is increased to (N_{r}, N_{θ}) = (1500, 3000).
Fig. 2 Comparison of the endofrun surface density for increasing SG (t = 300 t_{0}). Top row and from left to right: q_{0} = ∞, 4, 2; bottom row and from left to right: q_{0} = 1, 0.5, 0.25. The initial vortex parameters are: δ = 1.5 and χ = 14. The simulation window was centred in r ∈ [42.5, 57.5] AU in order to better appreciate the inner vortex structure. 
Fig. 3 Comparison of the endofrun Rossby number for increasing SG (t = 300 t_{0}). Top row and from left to right: q_{0} = ∞, 4, 2; bottom row and from left to right: q_{0} = 1, 0.5, 0.25. The initial vortex parameters are: δ = 1.5 and χ = 14. The simulation window was centred in r ∈ [42.5, 57.5] AU in order to better appreciate the inner vortex structure. 
Fig. 4 Morphological evolution for different values of the SG parameter. Top: semiminor axis with respect to the pressure scale height at the vortex core for q_{0} = 0.25, 0.5, 1, 2, 4, ∞. Bottom: aspect ratio at the Ro = −0.08 level. Initial vortex parameters are: δ = 1.5 and χ = 14. For readability, uncertainties are plotted at every fourth time steps (see Appendix D). 
Fig. 5 Bulkvortex evolution for different values of the SGparameter. 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 . The initial vortex parameters are: δ = 1.5 and χ = 14. For readability, uncertainties are plotted at every fourth time steps (see Appendix D). 
3.2 Weak SG: q_{0} = 2, 4
There is no significant difference compared to the nonSG case where the vortex is slightly stretched by the shear. Only a 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 vortex 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 5. In particular, the aspect ratio rapidly reaches a maximum at t ~ 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 accelerated to 0.15 AU in 300 orbits but, similarly to previous cases, it remains very weak. Vortex mass declines by 5%. It is important 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 χ, 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 ~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, 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 fluidforce 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 coorbital 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 noencounter state. The secondary, in gravitational interaction with the main vortex, has horseshoe oscillations of period T_{hs} ~ 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 threebody interactions with the primaryvortex, 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 maximum in the annulus reaches half the peak density of the initial vortex. At the end of the run, the mainvortex 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 ≲ q_{0}): there is no significant departure from the nonSG case. The outer regions of the vortex are slightly stretched by shear and SG. The vortex reaches a quasisteady regime with nearly constant mass and migration. The second regime is the domain of (B) selfgravitating vortices (0.5 ≲ q_{0} ≲ 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 longterm evolution, similar to the nonSG case. This suggests the existence of quasisteady 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 releasing 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 coorbital 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 secondaries 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 ~ 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 following 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.
Fig. 6 Stability domain of the selfgravitating 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 lefthand side and the righthand side of the box. χ and δ are, respectively, the vortex aspect ratio and the radial extent at t = 50 t_{0}. 
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_{σ}(t_{f}) − C_{σ}(t_{i})∣/C_{σ}(t_{i}) ≤ 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 endofrun state of the vortex in the (χ, δ, 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 χ ~ 5, vortex stability is independent of δ; the vortex is stable if q_{0} > 0.375 and unstable otherwise and (ii) if χ ~ 14, vortex stability depends on δ; the critical value of q_{0} increases with δ (from 0.375 to 0.65), and therefore large vortices are more easily destabilised than small vortices (even at lower surfacedensities). When the radialextent is constant: (i) if δ ~ 1, vortex stability is independent of the aspect ratio; the vortex is stable if q_{0} > 0.375 and unstable otherwise and (ii) if δ ~ 2.5, vortex stability depends on the aspect ratio; the critical value of q_{0} increases with χ (from 0.375 to 0.65); in this case, elongated vortices are more easily destabilised than more compact vortices (even at lower surfacedensities).
In summary, SG is able to significantly change the stability of gaseous vortices. Assuming standard disc configurations with h(r_{0}) = 0.05 and β_{σ} = −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 ~8.
4.2 Theoretical Estimate
Gravitational stability of circumstellar discs is constrained by the famous Toomre criterion c_{s}κ/(πGσ) > 1, where κ 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 ; assuming that , we get: (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 σ = σ_{0} (1 + C_{σ})^{4}. Accounting for a scaling factor of ~1.65 to better fit the numerical data, we propose an empirical condition for SGstable vortices: (9)
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 ∣Ro∣_{max}, C_{σ}, 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 (χ, δ, q_{0}) space. Applying the above criteria to the vortex of Sect. 3, we find that such a vortex is SG stable if Q_{0} ≳ 10.9, or Q ≳ 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.
Triplets (∣Ro∣_{max}, C_{σ}, and C_{P}) for simulated vortices at t = 50 t_{0}.
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 SGparameter, q_{0}, in the (χ, δ) space. Vortex stability is thus presented in the form of a simple map of the threshold value q_{0,c} as a function of δ and χ, selected in the range 0.8 ≤ δ ≤ 3 and 5.7 ≤ χ ≤ 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 selfgravitating vortices in greater detail at HR.
5 HR Simulations
We performed HR simulations to better characterise the internal structure of the selfgravitating 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 nonSG and nonisothermal 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 nonisothermal disc in order to show the differences between quasisteady 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 isocontours of the Rossby number have been plotted for a nonisothermal 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 nonisothermal 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 constantdensity core, the enclosed mass is m ~ π (Δr)^{2} σ. 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_{⊙})^{1/3}. Then, if the radius of the circular patch reaches ~R_{H}, the enclosed mass is and we get: (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_{θ}) = (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 HRstage 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 ∈ [34, 63] AU, h_{0} = 0.05, r_{0} = 50 AU, and (χ, δ) = (14, 0.92) at t = 1201). Furthermore, in order to choose the Toomre parameter most suited to obtaining selfgravitating vortices, several values of q0 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 nonisothermal cases, respectively. Finally, we found that the best critical values of the parameters are (q_{0,c} = 0.7, σ_{0} = 3.35 g cm^{−2}) and (q_{0,c} = 0.425, σ_{0} = 5.52 g cm^{−2}) in the isothermal and the nonisothermal cases, respectively. For both simulations we get Q ~ 6, which means that, in our simulations, the Hill radius of the vortex is resolved at least eight times.
We will now use our HR results to discriminate between vortices evolving in an isothermal or a nonisothermal disc. We will also use the results to better describe both the internal structure of a SG vortex and its coorbital environment.
Fig. 7 Rossby number at t = 140 t_{0} for an isothermal and a nonisothermal disc with SG. Radial zoom with r ∈ [46, 54] AU. Top: selfgravitating vortex in an isothermal disc. Middle: selfgravitating vortex in a nonisothermal disc. Bottom: reference vortex, without SG and in a nonisothermal disc. 
Fig. 8 Density and velocitydivergence at t = 140 t_{0} for an isothermal and a nonisothermal disc (with and without SG). Left column: relative density (σ/σ_{0}). The main vortex is masked to distinguish the annular overdensity and the spiral waves. Right column: normalised relative velocity field divergence (∇ · v′/(2Ω_{k})). From top to bottom: selfgravitating vortex in an isothermal disc, selfgravitating vortex in a nonisothermal disc, and the reference vortex (without SG and in a nonisothermal disc), respectively. 
Fig. 9 Radial and azimuthal profiles of the Rossby number and the relative density at t = 140 t_{0} for a nonisothermal disc. Top: radial profile of the Rossby number and the relative density (σ/σ_{0}), from left to right, respectively. Bottom: azimuth profile of the Rossby number and the relative density (σ/σ_{0}), from left to right, respectively. The dotted green line corresponds to the level Ro = −0.08 used to compute the vortex geometrical parameters. 
Fig. 10 Isocontours of the Rossby number in a nonisothermal disc at t = 140 t_{0}. Left: without SG. Right: with SG. 
5.3 Isothermal Versus NonIsothermal
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 nonisothermal cases is the presence of a rotating spirat in the Rossby number and the density distributions, which both evoke TaylorGreen vortices. This spiral pattern, which barely appears during the LRstage, tends to strengthen during the HRstage. 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 vortexcore seems to keep a stable elliptical shape. In contrast, in the nonisothermal case, the vortex keeps a smooth and steady elliptical shape with a nearly constant vorticity (or ≃ −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 nonisothermal discs.
5.3.2 Annular Bump and Waves
Gas flowing in the horseshoe region of the main vortex is one of the peculiarities of selfgravitating vortices. This is, indeed, obvious in the two upper panels of Fig. 7. The flow appears rather laminar in the nonisothermal 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 t0 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 nonisothermal case, the flow structure at r_{0} = 50 and near θ = [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 nonisothermal discs and stresses the difference between selfgravitating and gravityfree 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 isocontours 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. Selfgravitating vortices appear flatter than the gravityfree 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 isocontours for the external vortex regions suggest that gas compression is more important for selfgravitating vortices, as will be discussed in the next subsection.
Finally, it must be noted that a nonstationary 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 gravityfree 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 Uturn 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 threebody 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 nonisothermal cases, but also between gravityfree and selfgravitating vortices in the nonisothermal case. Important conclusions are reached concerning the shape, the evolution, and the coorbital environment of the vortices. When SG is taken into account: (i) vortices in an isothermaldisc 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 nonisothermal disc keep a stable elliptical shape until the end of run, and outside the vortex the flow of gas in the coorbital region is horseshoe shaped. In the case of nonisothermal discs: (i) selfgravitating vortices are flatter than gravityfree 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 largescale vortices reach a quasisteady state in nonisothermal 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 smallscale 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 selfgravitating 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), TarczayNehéz et al. (2022), and Lin & Pierens (2018), particularly where they concern vortex stability. Finally, we discuss possible observational implications.
6.1 SingleVortex 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 longlived and quasisteady vortices. This strategy was successful in the case of nonSG 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 gravityfree 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 , where δ_{*} 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 SGactivation procedure.
6.1.2 Required Time Steps for SG Activation
Moreover, following our strategy, HR should not be increased before a quasisteady 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 quasisteady state is reached and (iii) increase resolution.
6.2 InstabilityGenerated Vortices
Another way to study the evolution of selfgravitating vortices in PPDs is to assume that they are first formed following one of the numerous instabilities known to produce largescale 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 modeldependant, which is not the case in the singlevortex approach developed above.
6.3 Comparisons with Other Works
We now compare our results with those obtained in three different papers. These are TarczayNehéz et al. (2022), Pierens & Lin (2018) and Zhu & Baruteau (2016).
6.3.1 Comparison with TarczayNehé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 αdisc model with α = (10^{−4}, 10^{−5}), and their simulations had LR with (N_{r}, ν_{θ}) = (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 ≃ 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 ≲ Ro ≲ 0.2 and we get 0.6 ≲ 1.2 Q ≲ 1.2, which is consistent with our criterion in Eq. (9).
Simulation parameters for isothermal runs in Lin & Pierens (2018) at the 500th orbit.
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 βcooling at the frontier of the deadzone, where the alphaviscosity 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 ≤ Ro ≤ −0.10 and 5 ≤ χ ≤ 7, or very different with Ro = −0.05 and χ ≥ 25. These very weak and elongated structures are likely due to the artificial αviscosity and the low numerical resolution. Care should be paid to the fact that the authors named theirs studies without and with βcooling and thermal radiation as isothermal and nonisothermal, 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, χ ≥ 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 ≲ 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 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 (χ, δ) = (6, 0.85) and with our Fig. 6, we estimated that q_{0,c} ≤ 0.375^{5}. The results gathered in our Table 6 are found to be consistent with our criterion (q_{0} ≥ 0.375 ≥ q_{0,c}), except in the last line where the inequality is indefinite.
Our SG parameter, q_{0}, corresponding to the density factor, Σ_{0}, of Zhu & Baruteau (2016).
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, τ_{c}) ≥ 1, where Re is the Reynolds number and τ_{c} is the betacooling 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 TarczayNehéz et al. (2022), checked if a simple stability criterion involving power laws, , can succeed.
6.4 Theoretical Models for Steady Vortices in SGFree 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^{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^{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 (nonisothermal and gravityfree 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 nonbarotropic 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 gravityfree runs (i.e. X_{Ro} = 9, χ_{σ} = 7.7, and χ_{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_{⊙}, 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 substellar 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 largescale 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 largescale 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} ≥ 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 ~ 1400 K with a similar molecular weight as H_{2} and a vortex located at r_{0} = 0.3 AU from the star with χ ~ 7 and δ ~ 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:
Finally, if a vortex is really present, the surface density of the disc should not be larger than 62 × 10^{3} 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 VortexInduced 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, 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 background 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 stability 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 vortex remains very weak (see Figs. 8 and A.1), and the only way to obtain stronger waves is to have large amount of solid material 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 trappinginvortex 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 possible 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 nonnegligible amount of solid particles, can be considered as a tank of intermediatesized 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 vortices 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 viscosity. 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 nonSG case.
Selfgravitating vortices. vortices stretch out, contract radially, 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 smallscale vortex remnants; this can occur even if Toomre’s parameter Q_{0} ~ 10.
From these simulations it seems that for selfgravitating vortices, vortex strength is not necessarily proportional to the aspect ratio inverse. This should be confirmed by other independent 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 destabilising effects if 1.2 ∣Ro∣_{max} Q ≳ 1, or again if q_{0} ≳ q_{0,c}. This criterion is supported by comparison with the results of three independent studies (TarczayNehé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_{θ}) = (1800, 16 000), have shown that numerical resolution is key to correctly and finely describe the shape and evolution of SG vortices. Nearly uniform resolution, in r and θ, is necessary and simple estimates show that simulations should resolve distances of the order of ~H/(3Q).
Disc temperature. Comparisons between the HR simulations of Sect. 5 clearly show that SG vortices in nonisothermal 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 nonisothermal discs. The key point here is that a temperature gradient of the disc could be necessary for SG vortices to reach a nearly steadystate.
Zhu & Baruteau (2016), Regály & Vorobyov (2017b) and TarczayNehé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 initial state in our simulations), they cannot be sustained. Thus, if observed crescentshaped asymmetries are indeed vortex presence tracks, these vortices should be hosted in nonmassive discs. However, in light of selfgravitating vortices stability (as defined in this work), the definition of nonmassive disc should be clarified with respect to the other physical ingredients 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 dusttrapping simulations in permanent selfgravitating vortices. For example, these simulations should be carried out in nonisothermal discs in which q_{0} ≥ 1.2 q_{0,c} minimises the destabilising effects of SG, in anticipation of their possible enhancement by the dust backreaction; but also a q_{0} value relatively close to q_{0,c} would imply nonnegligible 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’AixMarseille 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 2020A0090411982 made by GENCI. This work was granted access to the HPC/AI resources of CINES under the allocation 2018A0050407407 made by GENCI.
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.
Fig. A.1 Compressibility waves, ∇ · v′/ (2Ω_{k}), emitted by a quasisteady vortex in its full simulation box, r ∈ [25, 72] AU, as used in Section 3. 
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 (χ, δ, q_{0}) the vortex stability with respect to SG. In order to improve results, three additional simulations were performed for the following triplets: (χ,δ, 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 (χ, δ), 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 midvalue 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.
Fig. B.1 Vortex stability in the (χ, δ, q_{0}) space. Each point is an intermediate resolution simulation : (N_{r}, N_{θ}) = (1500, 3000). χ and δ are, respectively, the vortex aspect ratio and the radial extent at t = 50 t_{0}. 
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 . We denote the vortex angular position by max.
Fig. C.1 Mass of the vortex core with respect to the µ threshold at t = 158 t_{0}. 
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 µ. In order to choose appropriately this threshold we plotted in Figure C.1 contours for different values of µ 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 µ. However, a lower threshold, such as µ = [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 µ = 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° 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 ∣θ − θ_{max}∣ ≤ 45°. In summary, the mass of a vortex core is computed numerically in our code as:
where τ = {(r, θ) ∈ ℝ^{+} × [0, 2π[ : C_{σ} ≥ 0.25 and ∣θ − θ_{max}∣ ≤ 45°}. µ + 1 = [1.25, 1.40,1.50, 1.60] (from the exterior to the vortex core). The respective vortexcore masses are 0.18, 0.1, 0.06, and 0.02 (in M_{J}). The initial vortex parameters are: (χ, δ) = (14, 1.5) (t = 0) and q_{0} = 0.5.
In the (r, θ) coordinate system, Rossby isocontours are ellipses, while in the Cartesian coordinate system, ellipses are curved.
Fig. C.2 Definition of aspect ratio, χ and radial extent, Δr, fitting the Ro = −0.8 isocontours. 
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 ≤ −0.8. The mean Rossby number is computed numerically in the extracted region:
where dr dθ is the total surface of the extracted region and Γ = {(r, θ) ∈ ℝ^{+} × [0, 2π[ : Ro ≤ −0.8 and ∣θ − θ_{max}∣ ≤ 45°}.
Fig. C.3 Vortex masses with respect to different density thresholds and filters. 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: µ + 1 = 1.05, m_{core} = 1.69 M_{J}. Top right: µ + 1 = 1.15, m_{core} = 0.78 M_{J} Bottom left: µ + 1 = 1.25, m_{core} = 0.28 M_{J}. Bottom right: µ + 1 = 1.25, m_{core} = 0.18 M_{J} 
Appendix C.4 Aspect ratios and radial width
Aspect ratios were computed thanks to the python routine Least Squares fitting of ellipses (Hammel & SullivanMolina 2020), which provides the ellipse parameters (semimajor and semiminor axes, inclination) which better fits a given set of points. 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, Ar, 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 isocontours, 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 particularly, 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: δ_{*} = Max(δ_{r}, rδ_{θ}).
Aspect ratios: The aspect ratio of the fitted ellipse (see Appendix C.4) is defined as , where a and b are, respectively, the semimajor and semiminor axes. Thus, the uncertainty is : .
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} = ∞ 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 115^{th} orbit. After 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 quasistationary 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.
Resolution for numerical convergency tests
Fig. E.1 Numerical convergency for Sections 3 and 4 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}. 
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 selfgravitating vortex, hosted in a nonisothermal disc and studied in Section 5. Compared to the previous subsection, we show an additional plot exhibiting the Rossbynumber aspect ratio at the Ro = −0.01 level. This allows us to check how the vortexcore morphology is affected when resolution exceeds the threshold R_{H} = H/3Q. Convergence tests were performed, but were strongly timeconsuming; 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 vortexcore 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.
Fig. E.2 Numerical convergency for Section 5 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}. 
Appendix F Potential and Forces due to SG
Figure F.1 shows, at t = 140 t_{0}, the spatial distribution of the potential and the force components of SG for the nonisothermal run of Section 5.
Fig. F.1 SG: potential and forces at t = 140 t_{0}. Left: SG potential (Φ_{SG}) and isocontours. Middle: SG radial force (f_{r,SG}). Right: SG azimuthal force(f_{θ,SG}). 
References
 Andrews, S. M., Huang, J., Pérez, L. M., et al. 2018, ApJ, 869, L41 [Google Scholar]
 Barge, P., & Sommeria, J. 1994, in Circumstellar Dust Disks and Planet Formation, eds. R. Ferlet, & A. VidalMadjar (Berlin: Springer), 10, 295 [NASA ADS] [Google Scholar]
 Barge, P., & Sommeria, J. 1995, A&A, 295, L1 [NASA ADS] [Google Scholar]
 Barge, P., Richard, S., & Le Dizès, S. 2016, A&A, 592, A136 [Google Scholar]
 Barranco, J. A., & Marcus, P. S. 2005, ApJ, 623, 1157 [Google Scholar]
 Baruteau, C., Barraza, M., Pérez, S., et al. 2019, MNRAS, 486, 304 [Google Scholar]
 Blum, J., & Wurm, G. 2008, ARA&A, 46, 21 [Google Scholar]
 Bracco, A., Chavanis, P. H., Provenzale, A., & Spiegel, E. A. 1999, Phys. Fluids, 11, 2280 [NASA ADS] [CrossRef] [Google Scholar]
 Cazzoletti, P., van Dishoeck, E. F., Pinilla, P., et al. 2018, A&A, 619, A161 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chavanis, P. H. 2000, A&A, 356, 1089 [NASA ADS] [Google Scholar]
 de ValBorro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529 [Google Scholar]
 de ValBorro, M., Artymowicz, P., D’Angelo, G., & Peplinski, A. 2007, A&A, 471, 1043 [Google Scholar]
 Dong, R., Liu, S.Y., Eisner, J., et al. 2018, ApJ, 860, 124 [Google Scholar]
 Dullemond, C. P., & Dominik, C. 2005, A&A, 434, 971 [Google Scholar]
 Frigo, M., & Johnson, S. G. 2005, Special issue on “Program Generation, Optimization, and Platform Adaptation” Proceedings of the IEEE, 93, 216 [Google Scholar]
 Fu, W., Li, H., Lubow, S., & Li, S. 2014a, ApJ, 788, L41 [Google Scholar]
 Fu, W., Li, H., Lubow, S., Li, S., & Liang, E. 2014b, ApJ, 795, L39 [NASA ADS] [CrossRef] [Google Scholar]
 Godon, P., & Livio, M. 1999, ApJ, 523, 350 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Ward, W. R. 1973, ApJ, 183, 1051 [Google Scholar]
 Goodman, J., Narayan, R., & Goldreich, P. 1987, MNRAS, 225, 695 [Google Scholar]
 GRAVITY Collaboration (Abuter, R., et al.) 2017, A&A, 602, A94 [Google Scholar]
 Hammel, B., & SullivanMolina, N. 2020, https://doi.org/10.5281/zenodo.3723294 [Google Scholar]
 Hammer, M., Kratter, K. M., & Lin, M.K. 2017, MNRAS, 466, 3533 [Google Scholar]
 Hayashi, C. 1981, Prog. Theor. Phys. Suppl., 70, 35 [Google Scholar]
 Huang, P., Dong, R., Li, H., Li, S., & Ji, J. 2019, ApJ, 883, L39 [NASA ADS] [CrossRef] [Google Scholar]
 Huré, J. M., Hersant, F., Carreau, C., & Busset, J. P. 2008, A&A, 490, 477 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Inaba, S., & Barge, P. 2006, ApJ, 649, 415 [NASA ADS] [CrossRef] [Google Scholar]
 Inaba, S., Barge, P., Daniel, E., & Guillard, H. 2005, A&A, 431, 365 [Google Scholar]
 Jacquet, E., Balbus, S., & Latter, H. 2011, MNRAS, 415, 3591 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., & Youdin, A. 2007, ApJ, 662, 627 [Google Scholar]
 Johansen, A., Henning, T., & Klahr, H. 2006, ApJ, 643, 1219 [Google Scholar]
 Johansen, A., Oishi, J. S., Mac Low, M.M., et al. 2007, Nature, 448, 1022 [Google Scholar]
 Kida, S. 1981a, J. Fluid Mech., 112, 397 [NASA ADS] [CrossRef] [Google Scholar]
 Kida, S. 1981b, J. Phys. Soc. Japan, 50, 3517 [Google Scholar]
 Klahr, H. 2004, ApJ, 606, 1070 [Google Scholar]
 Klahr, H. H., & Bodenheimer, P. 2003, ApJ, 582, 869 [Google Scholar]
 Klahr, H., Pfeil, T., & Schreiber, A. 2018, Handbook of Exoplanets, Instabilities and Flow Structures in Protoplanetary Disks: Setting the Stage for Planetesimal Formation, eds. D., Hans, J., & B. Juan Antonio (Berlin: Springer), 138 [Google Scholar]
 Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32 [Google Scholar]
 Lesur, G., & Papaloizou, J. C. B. 2009, A&A, 498, 1 [Google Scholar]
 Lesur, G., & Papaloizou, J. C. B. 2010, A&A, 513, A60 [Google Scholar]
 Li, H., Finn, J. M., Lovelace, R. V. E., & Colgate, S. A. 2000, ApJ, 533, 1023 [Google Scholar]
 Li, H., Colgate, S. A., Wendroff, B., & Liska, R. 2001, ApJ, 551, 874 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, M.K. 2012, ApJ, 754, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, M.K., & Papaloizou, J. C. B. 2011, MNRAS, 415, 1445 [Google Scholar]
 Lin, M.K., & Pierens, A. 2018, MNRAS, 478, 575 [NASA ADS] [CrossRef] [Google Scholar]
 Lopez, B., Lagarde, S., Jaffe, W., et al. 2014, The Messenger, 157, 5 [NASA ADS] [Google Scholar]
 Lovelace, R. V. E., & Hohlfeld, R. G. 2013, MNRAS, 429, 529 [NASA ADS] [CrossRef] [Google Scholar]
 Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805 [Google Scholar]
 Lyra, W., Johansen, A., Klahr, H., & Piskunov, N. 2009, A&A, 493, 1125 [Google Scholar]
 Manger, N., Klahr, H., Kley, W., & Flock, M. 2020, MNRAS, 499, 1841 [NASA ADS] [CrossRef] [Google Scholar]
 Manger, N., Pfeil, T., & Klahr, H. 2021, MNRAS, 508, 5402 [NASA ADS] [CrossRef] [Google Scholar]
 Marcus, P. S., Pei, S., Jiang, C.H., et al. 2015, ApJ, 808, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Meheut, H., Casse, F., Varniere, P., & Tagger, M. 2010, A&A, 516, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Meheut, H., Meliani, Z., Varniere, P., & Benz, W. 2012, A&A, 545, A134 [Google Scholar]
 Meheut, H., Lovelace, R., & Lai, D. 2013, MNRAS, 430, 3 [Google Scholar]
 Miranda, R., Li, H., Li, S., & Jin, S. 2017, ApJ, 835, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Mittal, T., & Chiang, E. 2015, ApJ, 798, L25 [Google Scholar]
 Müller, T. W. A., Kley, W., & Meru, F. 2012, A&A, 541, A123 [Google Scholar]
 Paardekooper, S.J., Lesur, G., & Papaloizou, J. C. B. 2013, Eur. Phys. J. Web Conf., 46, 05001 [CrossRef] [EDP Sciences] [Google Scholar]
 Pérez, L. M., Benisty, M., Andrews, S. M., et al. 2018, ApJ, 869, L50 [CrossRef] [Google Scholar]
 Petersen, M. R., Julien, K., & Stewart, G. R. 2007a, ApJ, 658, 1236 [NASA ADS] [CrossRef] [Google Scholar]
 Petersen, M. R., Stewart, G. R., & Julien, K. 2007b, ApJ, 658, 1252 [NASA ADS] [CrossRef] [Google Scholar]
 Pierens, A., & Lin, M.K. 2018, MNRAS, 479, 4878 [NASA ADS] [Google Scholar]
 Raettig, N., Klahr, H., & Lyra, W. 2015, ApJ, 804, 35 [Google Scholar]
 Raettig, N., Lyra, W., & Klahr, H. 2021, ApJ, 913, 92 [NASA ADS] [CrossRef] [Google Scholar]
 Regály, Z., & Vorobyov, E. 2017a, A&A, 601, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Regály, Z., & Vorobyov, E. 2017b, MNRAS, 471, 2204 [CrossRef] [Google Scholar]
 Regály, Z., Juhász, A., Sándor, Z., & Dullemond, C. P. 2012, MNRAS, 419, 1701 [Google Scholar]
 Regály, Z., Sándor, Z., Csomós, P., & Ataiee, S. 2013, MNRAS, 433, 2626 [Google Scholar]
 Rendon Restrepo, S., Barge, P., & Vavrik, R. 2022, A&A, in press, https://doi.org/10.1051/00046361/202243518 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rich, E. A., Wisniewski, J. P., Sitko, M. L., et al. 2020, ApJ, 902, 4 [Google Scholar]
 Richard, S., Barge, P., & Le Dizès, S. 2013, A&A, 559, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Richard, S., Nelson, R. P., & Umurhan, O. M. 2016, MNRAS, 456, 3571 [NASA ADS] [CrossRef] [Google Scholar]
 Robert, C. M. T., Méheut, H., & Ménard, F. 2020, A&A, 641, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rometsch, T., Ziampras, A., Kley, W., & Béthune, W. 2021, A&A, 656, A130 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Safronov, V. 1972, Evolution of the protoplanetary cloud and formation of the earth and the planets, NASA technical translation (Israel Program for Scientific Translations) [Google Scholar]
 Soon, K.L., Momose, M., Muto, T., et al. 2019, PASJ, 71, 124 [Google Scholar]
 Squire, J., & Hopkins, P. F. 2018, MNRAS, 477, 5011 [Google Scholar]
 Stoll, M. H. R., & Kley, W. 2014, A&A, 572, A77 [Google Scholar]
 Surville, C. 2013, Ph.D. Thesis, AixMarseille University, ed. P. Barge, Marseille, France [Google Scholar]
 Surville, C., & Barge, P. 2015, A&A, 579, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Surville, C., Mayer, L., & Lin, D. N. C. 2016, ApJ, 831, 82 [NASA ADS] [CrossRef] [Google Scholar]
 Tanga, P., Babiano, A., Dubrulle, B., & Provenzale, A. 1996, Icarus, 121, 158 [NASA ADS] [CrossRef] [Google Scholar]
 TarczayNehéz, D., Regály, Z., & Vorobyov, E. 2020, MNRAS, 493, 3014 [CrossRef] [Google Scholar]
 TarczayNehéz, D., Rozgonyi, K., & Regály, Z. 2022, MNRAS, 511, 6055 [CrossRef] [Google Scholar]
 Teed, R. J., & Latter, H. N. 2021, MNRAS, 507, 5523 [CrossRef] [Google Scholar]
 Terquem, C. E. J. M. L. J. 2008, ApJ, 689, 532 [NASA ADS] [CrossRef] [Google Scholar]
 Toomre, A. 1964, ApJ, 139, 1217 [Google Scholar]
 Tsukagoshi, T., Muto, T., Nomura, H., et al. 2019, ApJ, 878, L8 [Google Scholar]
 van der Marel, N., van Dishoeck, E. F., Bruderer, S., et al. 2013, Science, 340, 1199 [Google Scholar]
 Varga, J., Hogerheijde, M., van Boekel, R., et al. 2021, A&A, 647, A56 [Google Scholar]
 Varnière, P., & Tagger, M. 2006, A&A, 446, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weidenschilling, S. J. 1977, MNRAS, 180, 57 [Google Scholar]
 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 [Google Scholar]
 Whipple, F. L. 1972, in From Plasma to Planet, ed. A. Elvius (Hoboken: John Wiley & Sons Inc), 211 [Google Scholar]
 Youdin, A. N., & Goodman, J. 2005, ApJ, 620, 459 [Google Scholar]
 Zhu, Z., & Baruteau, C. 2016, MNRAS, 458, 3918 [Google Scholar]
 Zsom, A., Ormel, C. W., Güttler, C., Blum, J., & Dullemond, C. P. 2010, A&A, 513, A57 [Google Scholar]
Using definitions exposed in Table 3 we obtain the sound speed in the vortex core, , where c_{s,0} is the unperturbed sound speed.
The map of Fig. 6 is valid only in the range 0.8 ≲ δ ≲ 3.
All Tables
Vortex parameters at t = 50 t_{0} with respect to the initial Gaussian model parameters.
Simulation parameters for isothermal runs in Lin & Pierens (2018) at the 500th orbit.
Our SG parameter, q_{0}, corresponding to the density factor, Σ_{0}, of Zhu & Baruteau (2016).
All Figures
Fig. 1 Gaussian vortex model injected at t = 0. Top: Rossby number. Bottom: normalised density (σ(r, θ)/σ_{0}(r)). 

In the text 
Fig. 2 Comparison of the endofrun surface density for increasing SG (t = 300 t_{0}). Top row and from left to right: q_{0} = ∞, 4, 2; bottom row and from left to right: q_{0} = 1, 0.5, 0.25. The initial vortex parameters are: δ = 1.5 and χ = 14. The simulation window was centred in r ∈ [42.5, 57.5] AU in order to better appreciate the inner vortex structure. 

In the text 
Fig. 3 Comparison of the endofrun Rossby number for increasing SG (t = 300 t_{0}). Top row and from left to right: q_{0} = ∞, 4, 2; bottom row and from left to right: q_{0} = 1, 0.5, 0.25. The initial vortex parameters are: δ = 1.5 and χ = 14. The simulation window was centred in r ∈ [42.5, 57.5] AU in order to better appreciate the inner vortex structure. 

In the text 
Fig. 4 Morphological evolution for different values of the SG parameter. Top: semiminor axis with respect to the pressure scale height at the vortex core for q_{0} = 0.25, 0.5, 1, 2, 4, ∞. Bottom: aspect ratio at the Ro = −0.08 level. Initial vortex parameters are: δ = 1.5 and χ = 14. For readability, uncertainties are plotted at every fourth time steps (see Appendix D). 

In the text 
Fig. 5 Bulkvortex evolution for different values of the SGparameter. 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 . The initial vortex parameters are: δ = 1.5 and χ = 14. For readability, uncertainties are plotted at every fourth time steps (see Appendix D). 

In the text 
Fig. 6 Stability domain of the selfgravitating 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 lefthand side and the righthand side of the box. χ and δ are, respectively, the vortex aspect ratio and the radial extent at t = 50 t_{0}. 

In the text 
Fig. 7 Rossby number at t = 140 t_{0} for an isothermal and a nonisothermal disc with SG. Radial zoom with r ∈ [46, 54] AU. Top: selfgravitating vortex in an isothermal disc. Middle: selfgravitating vortex in a nonisothermal disc. Bottom: reference vortex, without SG and in a nonisothermal disc. 

In the text 
Fig. 8 Density and velocitydivergence at t = 140 t_{0} for an isothermal and a nonisothermal disc (with and without SG). Left column: relative density (σ/σ_{0}). The main vortex is masked to distinguish the annular overdensity and the spiral waves. Right column: normalised relative velocity field divergence (∇ · v′/(2Ω_{k})). From top to bottom: selfgravitating vortex in an isothermal disc, selfgravitating vortex in a nonisothermal disc, and the reference vortex (without SG and in a nonisothermal disc), respectively. 

In the text 
Fig. 9 Radial and azimuthal profiles of the Rossby number and the relative density at t = 140 t_{0} for a nonisothermal disc. Top: radial profile of the Rossby number and the relative density (σ/σ_{0}), from left to right, respectively. Bottom: azimuth profile of the Rossby number and the relative density (σ/σ_{0}), from left to right, respectively. The dotted green line corresponds to the level Ro = −0.08 used to compute the vortex geometrical parameters. 

In the text 
Fig. 10 Isocontours of the Rossby number in a nonisothermal disc at t = 140 t_{0}. Left: without SG. Right: with SG. 

In the text 
Fig. A.1 Compressibility waves, ∇ · v′/ (2Ω_{k}), emitted by a quasisteady vortex in its full simulation box, r ∈ [25, 72] AU, as used in Section 3. 

In the text 
Fig. B.1 Vortex stability in the (χ, δ, q_{0}) space. Each point is an intermediate resolution simulation : (N_{r}, N_{θ}) = (1500, 3000). χ and δ are, respectively, the vortex aspect ratio and the radial extent at t = 50 t_{0}. 

In the text 
Fig. C.1 Mass of the vortex core with respect to the µ threshold at t = 158 t_{0}. 

In the text 
Fig. C.2 Definition of aspect ratio, χ and radial extent, Δr, fitting the Ro = −0.8 isocontours. 

In the text 
Fig. C.3 Vortex masses with respect to different density thresholds and filters. 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: µ + 1 = 1.05, m_{core} = 1.69 M_{J}. Top right: µ + 1 = 1.15, m_{core} = 0.78 M_{J} Bottom left: µ + 1 = 1.25, m_{core} = 0.28 M_{J}. Bottom right: µ + 1 = 1.25, m_{core} = 0.18 M_{J} 

In the text 
Fig. E.1 Numerical convergency for Sections 3 and 4 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}. 

In the text 
Fig. E.2 Numerical convergency for Section 5 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}. 

In the text 
Fig. F.1 SG: potential and forces at t = 140 t_{0}. Left: SG potential (Φ_{SG}) and isocontours. Middle: SG radial force (f_{r,SG}). Right: SG azimuthal force(f_{θ,SG}). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.