© S. R. Restrepo and P. Barge 2022

1 Introduction

The formation of planetesimals from the dust-grains embedded in the gas of protoplanetary discs (hereafter PPDs) remains one of the key problems of planet formation. Standard scenarios start from the sub-layer that the dust forms in the mid-plane of the disc after settling under the vertical component of the star’s gravity. They are based on the collisional growth of the solid particles through coagulation and/or sticking mechanisms, generally followed by the formation of self-gravitating dust-clumps that are thought to compact into kilometre-sized planetesimals (Safronov 1972; Goldreich & Ward 1973). However, these models are known to stumble over a number of ‘barriers’ (metre-size, fragmentation, bouncing; Whipple 1972; Weidenschilling 1977; Dullemond & Dominik 2005; Blum & Wurm 2008; Zsom et al. 2010) that prevent them from giving a complete picture of the problem. More recent models are based on resonant drag instabilities such as the streaming instability (Youdin & Goodman 2005; Jacquet et al. 2011), in which pressure bumps can marry with dust-clumps to deliver self-gravitating dust-clouds that are dense enough to form big planetesimals (up to Ceres-sized planetesimals; Johansen et al. 2007). Once planet cores have grown due to the accumulation of planetesimals or the sweeping of residual pebbles (Lambrechts & Johansen 2012), another important requirement for models is to reach the critical mass for run-away gas accretion before the solid material is accreted by the star and the gas is dissipated by winds and photoevaporation.

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

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

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

From a theoretical point of view, strong dust concentrations in PPDs are known to be the location of instabilities due to the dust feedback. Similar situations also occur in the strong dust concentrations produced by a vortex, as observed in numerical 2D simulations (Fu et al. 2014b; Surville et al. 2016). In this case various dust and gas instabilities are at work, such as the Kelvin-Helmoltz instability or specific resonant-drag instabilities (Squire & Hopkins 2018). This evolution is complex and seems to end in a dissipation of the vortex. Of course, a more realistic description of the problem should include the vortex generation mechanism; for example, Miranda et al. (2017) claim that dust-laden vortices could last for thousands of orbits, which would be sufficient for dust growth or gravitational collapse. In any case, in 2D models, the dissipation timescales are difficult to estimate since the back-reaction of the dust particle becomes problematic at a high dust-to-gas ratio. 3D simulations are scarce and seem to show that vortices remain stable and behave as Taylor columns, even if the mid-plane region is perturbed by a thin dust-layer (Raettig et al. 2021). These authors also find that dust densities in these vortices can be higher than Roche’s density, possibly leading to a gravitational collapse. Another issue of the vortex scenario is the possible growth of the elliptical instability, which is known to destroy 3D-vortices (Lesur & Papaloizou 2009). However, vortices with a large enough aspect ratio could survive due to a weaker growth rate of this instability (Richard et al. 2013), and furthermore, the exact role of compressibility on the growth rate also appears unclear.

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

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

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

2 Theoretical Context

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

2.1 Disc Model and Evolution Equations

We investigate the evolution of an ideal gas (molecular dihydrogen, µ ~ 2.34) in a thin (2D), non-homentropic disc orbiting a solar-type star. Compressibility must be taken into account since the vortex size often exceeds the scale height of the disc. On the other hand, viscosity is neglected since we are working in the magnetically dead zone of the protoplanetary discs. The equations governing gas dynamics in a 2D protoplanetary disc are thus the continuity equation, Euler’s equation, and the conservation of energy equation, which are: $${\partial }_t\sigma +\nabla \cdot \left(\sigma \nu \right)=0,$$(1) $${\partial }_t\nu +\left(\nu \cdot \nabla \right)\nu =-\frac{\nabla P}{\sigma }-\nabla \left({\text{Φ}}_{\odot }+{\text{Φ}}_{\mathrm{SG}}+{\text{Φ}}_{\text{ind}}\right),$$(2) $$\sigma \left[{\partial }_{t}e+\left(\nu \cdot \nabla \right)e\right]=-P\nabla \cdot \nu ,$$(3)

where P, σ, and v are the vertically integrated pressure, density, and the gas velocity, respectively, $e=\frac{1}{\gamma -1}\frac{P}{\sigma }$ is the internal energy of an ideal gas, and γ = 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 non-active region of the PPDs, although Stoll & Kley (2014) noted that, even in this region, vertical-shear instability can generate turbulence associated with an effective viscosity. According to Manger et al. (2020, 2021), Formula (14), our simulations correspond to an infinite cooling time and to α/(10⁻⁴) ≪ 1, so that the viscosity is very small and compatible with our inviscid disc assumption. In the case of a finite cooling time, we should account for a non-negligible viscosity in the computations.

At equilibrium the background flow is assumed to be axisymmetric with radial profiles for the temperature T₀(r), the surface density σ₀(r), and the pressure P₀(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₀(r) = Ω_kc_s0/(π_Gσ₀), 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: $$\begin{array}{llll}{\sigma }_0\left(r\right)\hfill & =\hfill & {\sigma }_0{\left(r/r_0\right)}^{{\beta }_{\sigma }}\hfill & \text{with\hspace{0.17em}}{\beta }_{\sigma }=-1.1,\hfill \\ T_0\left(r\right)\hfill & =\hfill & T_0{\left(r/r_0\right)}^{{\beta }_T},\hfill & \hfill \\ H/r\hfill & =\hfill & 0.05{\left(r/r_0\right)}^{{\beta }_H-1},\hfill & \text{with\hspace{0.17em}}{\beta }_H\frac{{\beta }_T+3}{2},\hfill \\ Q_0\left(r\right)\hfill & =\hfill & 0.0279\left(\frac{1700\text{\hspace{0.17em}g\hspace{0.17em}}{\text{cm}}^{-2}}{{\sigma }_0}\right){\left(r/r_0\right)}^{{\beta }_Q},\hfill & \hfill \\ \hfill & \hfill & \hfill & \text{with\hspace{0.17em}}{\beta }_Q={\beta }_H-3-{\beta }_{\sigma },\hfill \end{array}$$(4)

where r is the normalised radius at 1 AU, r₀ = 50, σ₀(r) and σ₀ are the surface density (in g cm⁻²) at radial positions r and r₀, respectively, and H is the pressure scale height. T₀(r) and T₀ = 6.2 K refer to the temperature at radial positions r and r₀, 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 non-isothermal 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.

2.2 Gravity Terms

In polar coordinates the SG potential reads: $${\text{Φ}}_{\text{SG}}\left(r\right)=-G{\displaystyle \underset{\text{disc}}{\int }\frac{{\sigma }_0\left(r\prime \right)}{\left|r-r\prime \right|}{d}^{2}r\prime -G{\displaystyle \underset{\text{disc}}{\int }\frac{\sigma \left(r\prime ,\theta \prime \right)-{\sigma }_0\left(r\prime \right)}{\left|r-r\prime \right|}{d}^{2}r\prime ,}}$$(5)

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

where $A={\displaystyle \underset{\text{disc}}{\int }\frac{G\cos \theta \prime \text{d}m\prime }{r{\prime }^2}}$ and $B={\displaystyle \underset{\text{disc}}{\int }\frac{Gsin\theta \prime \text{d}m\prime }{r{\prime }^2}}$.

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

2.3 Numerical Method and Physical Parameters

The system of Eqs. (1)–(3) was solved numerically thanks to the latest version of the RoSSBi (Rotating Systems Simulation for Bi-fluids) code, which was specifically developed to study the evolution of PPDs (Inaba et al. 2005; Inaba & Barge 2006; Surville 2013). This new version has been extended to 3D and structured for high-performance parallelism; it is presented in detail elsewhere (Rendon Restrepo et al. 2022) with many tests 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 long-lived vortices, which can then be studied in detail (Surville & Barge 2015). Here, the same strategy is used to address the problem of vortices in a disc that has SG. Thus, instead of generating vortices through RWI or baroclinic instability, we preferred to inject a Gaussian solution as modelled in Surville & Barge (2015) since this allows a large variety of vortices to be explored by tuning their control parameters (radial width, Rossby number, and aspect ratio). This method allows the most appropriate vortex parameters to be selected, but also to save computational time since the timescale necessary to generate a steady vortex from a hydrodynamic instability generally exceeds hundreds of orbits. Since the geometrical parameters of the vortex can change during the relaxation phase, we gathered in Table 2 all of the parameters used to initialise our numerical simulations, as well as their value at t = 50t₀, the time at which SG is plugged in. It is interesting to note that in presence of viscosity, α = 10⁻⁴, and simulations with or without SG have shown that vortices decay at a rate proportional to the disc mass (Regály & Vorobyov 2017b).

2.3.1 A Smooth Activation of SG

In the presence of SG, numerical tests (at low-resolution, hereafter LR) have shown that persistent vortices can also survive, but after a transition period that is longer and more complex than in the non-self-gravitating case. This transient evolution is due to the difference between the injected and the exact vortex-solutions; of course it reduces to the relaxation already discussed in the previous section for vanishing SG. For these reasons, and to keep the same strategy as previously used, it was decided to gradually introduce SG during the computations, in order to avoid spurious vortex bursts, which occur when SG is introduced with a Heaviside time step function. To this end we have implemented in the code the time function: $$g\left(t\right)=\{\begin{array}{ccc}0& \text{if}& t\le 60t_0\\ 1-\exp \left(-{\left(\frac{t-60t_0}{T_{\text{act}}}\right)}^2\right)& \text{otherwise}& \end{array}.$$(7)

The time constant was chosen to be equal to the vortex relaxation time in the non-SG case; that is T_act ≃ 15 t₀, where t₀ = 353 years is the orbital period at r₀ = 50 AU. With this procedure, SG reaches 98% of its final value in ≃30 orbital periods.

2.3.2 Choice of Computational Box

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

Since in this work we are mainly interested by variations in vortices’ structure and not in their type I migration, the initial vortex should keep its radial position during the whole run. This is possible provided that the computational box includes the same number of inner and outer Lindblad resonances. Therefore, prior to the main study, for each Gaussian vortex structure, we ran a preliminary study at LR, (N_r, N_θ) = (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 quasi-steady vortex.

2.3.3 Suitable Physical Parameters

A number of physical parameters have been chosen to characterise as well as possible the coherent gaseous vortices and to make easier the discussion. It is, for example, convenient to distinguish between morphological and global parameters such as the radial position, the bulk-mass, and the spin. Since elliptical streamline models (e.g. Kida 1981a; Goodman et al. 1987; Chavanis 2000) approximately describe the central region of 2D vortices, the most natural morphological parameters are the ellipse parameters such as the vortex aspect ratio χ 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 iso-contours 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 bulk-vortex evolution once a nearly steady-state is reached, such as a long-term dynamical evolution. In particular, we numerically computed the distance from the vortex to the star, tracking the pressure maximum as a function of time. A debated question is whether vortices can be long-lived and robust against instabilities, and for this reason we computed the mean Rossby number in the vortex core or spin, $\overline{\text{Ro}}$, 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: $\underset{\chi \to \infty }{\lim }$ 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 over-density (see Sect. 5), and, eventually, some secondary eddies (see Sects. 3.4 and 5), which leads to a prudent calculation. We encourage the interested reader to refer to Appendix C, where we explain in detail the way we define the mass of a gaseous vortex and also its mean Rossby number. In Appendix D we estimated uncertainties. The definition of all the parameters and key quantities used in this paper are gathered in Table 3.

3 Vortex Evolution

Six numerical simulations were run to address the problem, starting from a single reference vortex that was followed over 300 orbital periods. All runs performed in this section started from a Gaussian vortex with aspect ratio χ = 14 and radial extent δ = 1.5 which was injected at r₀ = 50 AU and θ₀ = π. We exhibit in Fig. 1 the initial Gaussian vortex. Disc SG was quantified thanks to the SG parameter, q₀, whose values were selected among six different values: 0.25, 0.5, 1, 2, 4 and ∞ (non-SG case). The results of the simulations are presented in Figs. 2, 3, 4, and 5, and the temporal evolution can be found in online movies 1 and 2. The end-of-run snapshots of the surface density and the Rossby number, obtained for the six values of q₀, 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₀ = ∞ and values at t = 50 t₀ 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 quasi-stationary regime (see Table 2). Then, at t = 50 t₀, the resolution of the simulations was increased using a 2D cubic interpolation² 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³ and azimuth directions, respectively. At t = 60 t₀, 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₀                : Vortex relaxes to a quasi-steady state,

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

t = 60 t₀                        : Plug-in of SG,

60 t₀ < t < 300 t₀      : 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₀], we found that the vortex mass-loss is ~15% for every 250 orbits, if the resolution is (N_r, N_θ) = (500, 720), while mass-loss is negligible if resolution is increased to (N_r, N_θ) = (1500, 3000).

Fig. 2 Comparison of the end-of-run surface density for increasing SG (t = 300 t0). Top row and from left to right: q0 = ∞, 4, 2; bottom row and from left to right: q0 = 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 end-of-run Rossby number for increasing SG (t = 300 t0). Top row and from left to right: q0 = ∞, 4, 2; bottom row and from left to right: q0 = 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: semi-minor axis with respect to the pressure scale height at the vortex core for q0 = 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 Bulk-vortex evolution for different values of the SG-parameter. Top: radial distance to the star. Middle: bulk of the vortex mass; Mi is the vortex mass at t = 0. Bottom: absolute value of the mean Rossby number $\overline{\text{Ro}}$. 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₀ = 2, 4

There is no significant difference compared to the non-SG 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₀ = 1

If q₀ = 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₀, 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.5

At t = 60 t₀, 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, $\left|\overline{\text{Ro}}\right|$ sharply varies as a function of time. This likely indicates that the release of a secondary vortex results from the changes introduced by SG in the fluid-force equilibrium (pressure gradient, star attraction, and fictive forces); in the present case, ejection of secondaries could result from the predominance of the centrifugal forces. In particular, this secondary vortex picks up mass to the primary which, at t = 90 t₀, can lose more than 25% of its initial mass.

After the first release, the secondary vortex remains in the co-orbital region and crosses the main vortex at ~117 t₀ and ~145 t₀. During each encounter, the main vortex can pick up part of the secondary mass and, at the end of the run, reaches 60% of its initial value. After the 220th orbit the flow stabilises in a no-encounter state. The secondary, in gravitational interaction with the main vortex, has horseshoe oscillations of period T_hs ~ 65t₀. 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.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₀, the vortex is strongly deformed and loses more than 60% of its mass during the ejection of two secondaries. The evolution of the secondaries is also quite complex; it is governed by three-body interactions with the primary-vortex, but also by successions of splitting, merging, or dissipation processes. After the 190th orbit, the vortex strongly decays and finally disintegrates at the 300th orbit in an annular overdensity and two residual vortices that are located at [20°, 180°]. At this stage the density maximum in the annulus reaches half the peak density of the initial vortex. At the end of the run, the main-vortex has lost 75% of its mass.

3.5 Summary of Vortex Evolution

Three different regimes have been identified in our simulations. The first evolution is dominated by (A) shear (1 ≲ q₀): there is no significant departure from the non-SG case. The outer regions of the vortex are slightly stretched by shear and SG. The vortex reaches a quasi-steady regime with nearly constant mass and migration. The second regime is the domain of (B) self-gravitating vortices (0.5 ≲ q₀ ≲ 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₀. Our simulations show that such vortices have a stable and long-term evolution, similar to the non-SG case. This suggests the existence of quasi-steady vortex solutions of the Euler and Poisson equations, but further investigations are needed for firm conclusions on this point. We must stress that regime (B) includes the case of vortices releasing a single secondary vortex, such as observed when q₀ = 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.5) is performed using the methodology presented in Sect. 5. The last occurrence refers to (C) splitting by SG (q₀ < 0.5): SG strongly affects the vortex shape, up to destruction, but also the global flow. The initial vortex releases more than one secondary vortex, which remain trapped in the co-orbital region of the primary. In this regime the global parameters of the primary significantly change as a function of time. For smaller values of q₀, 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₀ = 5 whereas, in the vortex core, the mean value of Toomre parameter is Q ~ 10. These values of Q₀ 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 self-gravitating vortices as a function of q0. Left: when estimated from numerical simulations and right, when deduced from Eq. (9). Vortices are stable above the surface and unstable below (blue and orange, respectively). Both surfaces were projected on the left-hand side and the right-hand side of the box. χ and δ are, respectively, the vortex aspect ratio and the radial extent at t = 50 t0.

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₀), (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₀ and t_f = 300 t₀. A total of 45 different runs were carried out and provide a set of dots characterising the end-of-run state of the vortex in the (χ, δ, q₀) 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.375 and unstable otherwise and (ii) if χ ~ 14, vortex stability depends on δ; the critical value of q₀ increases with δ (from 0.375 to 0.65), and therefore large vortices are more easily destabilised than small vortices (even at lower surface-densities). When the radial-extent is constant: (i) if δ ~ 1, vortex stability is independent of the aspect ratio; the vortex is stable if q₀ > 0.375 and unstable otherwise and (ii) if δ ~ 2.5, vortex stability depends on the aspect ratio; the critical value of q₀ increases with χ (from 0.375 to 0.65); in this case, elongated vortices are more easily destabilised than more compact vortices (even at lower surface-densities).

In summary, SG is able to significantly change the stability of gaseous vortices. Assuming standard disc configurations with h(r₀) = 0.05 and β_σ = −1.1, we find that SG destabilises all vortices if q₀ < 0.375 (or again Q₀ < 7.5) and all vortices remain stable for moderate SG when q₀ > 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₀ 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 $\omega =2\left|\overline{Ro}\right|{\text{Ω}}_{\text{k}}$; assuming that $\left|\overline{\text{Ro}}\right|\text{~}{\left|\text{Ro}\right|}_{\max }$, we get: $$2\frac{c_{\text{s}}{\left|\text{Ro}\right|}_{\max }{\text{Ω}}_{\text{k}}}{\pi G\sigma }\gtrsim 1.$$(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₀ (1 + C_P) and σ = σ₀ (1 + C_σ)⁴. Accounting for a scaling factor of ~1.65 to better fit the numerical data, we propose an empirical condition for SG-stable vortices: $$\begin{array}{cc}{q}_0\gtrsim q_{0,c}=\frac{0.83h{\left(1+C_{\sigma }\right)}^{3/2}}{{\left|\text{Ro}\right|}_{\max }{\left(1+C_{\text{P}}\right)}^{1/2}}& \left(\text{hosting\hspace{0.17em}disc\hspace{0.17em}criterion}\right),\end{array}$$(9)

or again, $$\begin{array}{cccc}1.2{\left|\text{Ro}\right|}_{\max }Q\gtrsim 1& & & \left(\text{vortex\hspace{0.17em}core\hspace{0.17em}criterion}\right),\end{array}$$(11)

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₀). 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₀) space. Applying the above criteria to the vortex of Sect. 3, we find that such a vortex is SG stable if Q₀ ≳ 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.

4.3 Summary

The simple theoretical estimate we used to study the stability of the vortex core against SG allowed us to define a simple criterion based on a modified form of Toomre’s criterion. The new criterion was formulated in two different ways (see Eqs. (9) and (10)) to make the comparison easier. The widened numerical experiment we carried out also allowed us to build up a mapping of the SG-parameter, q₀, 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 self-gravitating vortices in greater detail at HR.

5 HR Simulations

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

5.1 Numerical Resolution in the Vortex Core

Numerical resolution is important for relevant descriptions of SG and gas evolution in the inner region of vortices. In a circular vortex patch around a constant-density core, the enclosed mass is m ~ π (Δr)² σ. Thus, the gravitational influence of this region of the vortex extends to a mean distance given by Hill’s radius, R_H = r₀(m/3 M_⊙)^1/3. Then, if the radius of the circular patch reaches ~R_H, the enclosed mass is $m~\pi R_{\text{H}}^2\sigma$ and we get: $$R_{\text{H}}=\frac{H}{3Q}.$$(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₀                           : SG is smoothly activated,

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

120 t₀ < t < 140 t₀       : HR simulation.

The final resolution is equivalent to 146 cells/H and 127 cells/H in the radial and azimuthal directions, respectively. Such a resolution was chosen so that both directions are almost equally resolved, which is important to correctly account for the SG. For obvious reasons of computational costs, the last HR-stage was stopped after only 20 orbits; despite this drawback, results are sufficient to distinguish substructures and to guess a possible steady state in the evolution of the main vortex. The outline for our HR simulations are the following: r ∈ [34, 63] AU, h₀ = 0.05, r₀ = 50 AU, and (χ, δ) = (14, 0.92) at t = 1201). Furthermore, in order to choose the Toomre parameter most suited to obtaining self-gravitating 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 non-isothermal cases, respectively. Finally, we found that the best critical values of the parameters are (q_0,c = 0.7, σ₀ = 3.35 g cm⁻²) and (q_0,c = 0.425, σ₀ = 5.52 g cm⁻²) in the isothermal and the non-isothermal 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 non-isothermal disc. We will also use the results to better describe both the internal structure of a SG vortex and its co-orbital environment.

Fig. 7 Rossby number at t = 140 t0 for an isothermal and a non-isothermal disc with SG. Radial zoom with r ∈ [46, 54] AU. Top: self-gravitating vortex in an isothermal disc. Middle: self-gravitating vortex in a non-isothermal disc. Bottom: reference vortex, without SG and in a non-isothermal disc.

Fig. 8 Density and velocity-divergence at t = 140 t0 for an isothermal and a non-isothermal disc (with and without SG). Left column: relative density (σ/σ0). The main vortex is masked to distinguish the annular over-density and the spiral waves. Right column: normalised relative velocity field divergence (∇ · v′/(2Ωk)). From top to bottom: self-gravitating vortex in an isothermal disc, self-gravitating vortex in a non-isothermal disc, and the reference vortex (without SG and in a non-isothermal disc), respectively.

Fig. 9 Radial and azimuthal profiles of the Rossby number and the relative density at t = 140 t0 for a non-isothermal 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 Iso-contours of the Rossby number in a non-isothermal disc at t = 140 t0. Left: without SG. Right: with SG.

5.3 Isothermal Versus Non-Isothermal

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

5.3.1 Vortex Core

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

5.3.2 Annular Bump and Waves

Gas flowing in the horseshoe region of the main vortex is one of the peculiarities of self-gravitating vortices. This is, indeed, obvious in the two upper panels of Fig. 7. The flow appears rather laminar in the non-isothermal case, but a more complex evolution is possible in the isothermal case, such as observed in the top panel of Fig. 7. This snapshot at t = 140 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 non-isothermal case, the flow structure at r₀ = 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 non-isothermal discs and stresses the difference between self-gravitating and gravity-free vortices. Comparisons are based on Figs. 7 and 8 (the two lower panels), and on Figs. 9 and 10.

5.4.1 Vortex

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

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

5.4.2 Annular Bump and Waves

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

5.5 Summary on HR Results

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

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

6 Discussion

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

6.1 Single-Vortex Strategy

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

6.1.1 Learning from HR runs

In fact, the numerical procedures used in Sects. 3 and 5 manage two different effects: a slow variation of SG as a function of time (which does not correspond to a true physical situation), and an increase in the numerical resolution, which necessarily introduces noise in the computations. Of course an increase in the numerical resolution also corresponds to a decrease in the numerical viscosity, which at the small sizes scales as $\text{~}{\delta }_{*}^2\text{Ω}$, 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 SG-activation procedure.

6.1.2 Required Time Steps for SG Activation

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

6.2 Instability-Generated Vortices

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

6.3 Comparisons with Other Works

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

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

The authors performed a similar study to our own using a wider sample (1980 different runs). In contrast to our work, their outline is different since vortices are generated by RWI at the edge of a magnetically dead zone in a viscous α-disc model with α = (10⁻⁴, 10⁻⁵), 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 $\left|\overline{\text{Ro}}\right|$ Q ≲ 1.2, which is consistent with our criterion in Eq. (9).

6.3.2 Comparison with Pierens & Lin (2018)

As in the previous paper, the authors study Rossby vortices in a disc with SG. These vortices are generated with or without β-cooling at the frontier of the dead-zone, where the alpha-viscosity jumps from 10⁻⁴ inside the dead zone to 10⁻² 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 non-isothermal, respectively. For meaningful comparisons and possible testing of Criterion 10, we focus on their isothermal case and excluded the cases f = [4, 8] because in our simulations, we never encountered vortices with such low Rossby numbers and such high elongations (Ro = −0.05, χ ≥ 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₀ 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⁵. The results gathered in our Table 6 are found to be consistent with our criterion (q₀ ≥ 0.375 ≥ q_0,c), except in the last line where the inequality is indefinite.

6.3.4 Comparison Summary

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

6.4 Theoretical Models for Steady Vortices in SG-Free Discs

Simple vortex solutions of the incompressible Euler equations at steady state have been studied by a number of authors (Goodman et al. 1987, hereafter GNG; Kida 1981b; Chavanis 2000). GNG’s model provides a rough solution of the Euler and continuity equations for a polytropic⁶ 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⁷ vortex solution in a sheared flow should satisfy the requirement that the surface density, pressure, and vorticity are functions of the stream function, which indirectly implies that these three quantities must have the same aspect ratios.

Nonetheless, the nearly steady state reached after the relaxation of the Gaussian model in the present work (non-isothermal and gravity-free run studied in Sect. 5) and in Surville & Barge (2015) suggests the existence of an analytical solution of the steady, compressible Euler equations in the case of 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 gravity-free 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 sub-stellar companion, and an inner disc misalignment with the line of sight. Following their conclusions, the most likely mechanism to explain such a dynamic asymmetry is growing dust grains being captured by a gaseous vortex. Indeed, the proximity of the star could produce axisymmetrical step jumps in the density and the pressure that induce a local extrema in the generalised potential vortensity distribution (Li et al. 2000); this is known to be favourable to RWI and the generation of a large-scale vortex. The presence of a vortex seems to be supported by the shadow variations revealed by the Hubble Space Telescope (Rich et al. 2020), since vortices are also pressure bumps whose vertical extent easily exceeds the mean scale height.

The results we get in this paper have potential implications on the disc density that observational predictions could infer. If a large-scale vortex is definitely present in the HD 163296 disc and at r₀ = 0.3 AU from the star, it should also be stable against SG forces; thus, following our criterion, we should have q₀ ≥ 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₂ and a vortex located at r₀ = 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: $${\sigma }_{0,max}\left(0.3\text{AU}\right)=\frac{c_{\text{s}}{\text{Ω}}_{\text{k}}}{0.375\cdot 4\pi G{Q}_{\text{c}}}=62\times {10}^3\text{g\hspace{0.17em}}{\text{cm}}^{-2}.$$

Finally, if a vortex is really present, the surface density of the disc should not be larger than 62 × 10³ g cm⁻² at 0.3 AU to avoid being destroyed by SG splitting. This is, of course, a very crude approach to the problem, but our goal was only to stress that new constraints can emerge when SG is taken into account.

6.5.2 Vortex-Induced Spiral Waves

In the absence of SG, vortices are known to excite spiral waves that cause them to migrate inwards towards the star. These waves, 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 trapping-in-vortex stage. Indeed, dust concentrations in this region could be observational markers of the vortices. As noted in Sect. 5, this annular region could be a simple annular pressure and density bump in a nearly Keplerian shear, for weak SG, or a more complex region with a horseshoe flow and 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 non-negligible amount of solid particles, can be considered as a tank of intermediate-sized particles, which possibly ‘feed’ the vortex during its evolution.

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

7 Conclusions and Perspectives

In this paper we reported the evolution, in an inviscid disc, of 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 non-SG case.

Self-gravitating 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 small-scale vortex remnants; this can occur even if Toomre’s parameter Q₀ ~ 10.

From these simulations it seems that for self-gravitating vortices, vortex strength is not necessarily proportional to the aspect ratio inverse. This should be confirmed by other 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₀ ≳ q_0,c. This criterion is supported by comparison with the results of three independent studies (Tarczay-Nehéz et al. 2022; Zhu & Baruteau 2016; Pierens & Lin 2018). In addition, the aforementioned criteria were also applied to a possible vortex candidate in HD 163296 disc (Varga et al. 2021) which allowed an estimation to be made of the upper limit for the gas density.

Numerical resolution. In the absence of viscosity, we found that simulations with at least (71, 24) cells/H in the radial and azimuthal directions, respectively, are needed to accurately describe vortices. If such a requirement is not fulfilled, vortices lose 15% of their mass during 250 orbits (see Appendix E.1) and this trend is accentuated in the presence of SG. On the other hand, simulations with (N_r, N_θ) = (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 non-isothermal discs: (i) evolve more slowly than vortices in isothermal discs, and (ii) keep a smooth elliptical shape throughout the run. On the other hand, vortices in isothermal discs are distorted by the SG forces and vortical substructures are generated in the horseshoe region; they also excite more compressive waves than in the case of non-isothermal discs. The key point here is that a temperature gradient of the disc could be necessary for SG vortices to reach a nearly steady-state.

Zhu & Baruteau (2016), Regály & Vorobyov (2017b) and Tarczay-Nehéz et al. (2022) found that vortices cannot form in massive discs since lower mode numbers are inhibited by SG. Although the strategy in this paper was different, our results are consistent with their findings: not only can vortices not form in massive discs, but even if they are well established (as the initial state in our simulations), they cannot be sustained. Thus, if observed crescent-shaped asymmetries are indeed vortex presence tracks, these vortices should be hosted in non-massive discs. However, in light of self-gravitating vortices stability (as defined in this work), the definition of non-massive disc should be clarified with respect to the other physical 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 dust-trapping simulations in permanent self-gravitating vortices. For example, these simulations should be carried out in non-isothermal discs in which q₀ ≥ 1.2 q_0,c minimises the destabilising effects of SG, in anticipation of their possible enhancement by the dust back-reaction; but also a q₀ value relatively close to q_0,c would imply non-negligible effects from SG in the gas and dusty phases.

Acknowledgements

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

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 quasi-steady 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₀) the vortex stability with respect to SG. In order to improve results, three additional simulations were performed for the following triplets: (χ,δ, q₀) = [(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₀ direction, from the stable to the unstable region.

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

• A blue cross in the q₀ direction indicates that the transition between stable and unstable vortices is the lowest blue cross.

Fig. B.1 Vortex stability in the (χ, δ, q0) space. Each point is an intermediate resolution simulation : (Nr, Nθ) = (1500, 3000). χ and δ are, respectively, the vortex aspect ratio and the radial extent at t = 50 t0.

Appendix C Computed Quantities

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

Appendix C.1 Main Vortex Position and Distance to the Star

We found that even in presence of secondary vortices the position of the main vortex is given by the pressure maximum. For the above reason, the distance to the star and the angular position of the vortex are computed by tracking the position of the maximum of pressure contrast $C_{\text{P}}=\frac{P}{P_0}-1$. We denote the vortex angular position by max.

Fig. C.1 Mass of the vortex core with respect to the µ threshold at t = 158 t0.

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.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: $$m_{\text{core}}={\displaystyle \underset{\tau }{\int }\sigma \left(r,\theta \right)rdrd\theta }$$

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 vortex-core 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.5.

In the (r, θ) coordinate system, Rossby iso-contours 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 iso-contours.

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: $$\overline{Ro}=\left({\displaystyle \underset{\text{Γ}}{\int }Ror\text{\hspace{0.17em}d}r\text{\hspace{0.17em}d}\theta }\right)/S,$$

where $S={\displaystyle \underset{\text{Γ}}{\int }r}$ 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, mcore = 1.69 MJ. Top right: µ + 1 = 1.15, mcore = 0.78 MJ Bottom left: µ + 1 = 1.25, mcore = 0.28 MJ. Bottom right: µ + 1 = 1.25, mcore = 0.18 MJ

Appendix C.4 Aspect ratios and radial width

Aspect ratios were computed thanks to the python routine Least Squares fitting of ellipses (Hammel & Sullivan-Molina 2020), which provides the ellipse parameters (semi-major and semi-minor axes, inclination) which better fits a given set of points. 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 iso-contours, for tracking geometrical parameters, since the extensive study achieved in Section 4 showed that the only common quantity shared by all 45 vortices was the Ro=−0.08 contours, while density and the pressure maximum amplitude vary for all different structures, which makes comparison a difficult task.

Appendix D Uncertainties

In numerical experiments the accuracy arises from the used numerical method. In 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 $\chi =\frac{a}{b}$, where a and b are, respectively, the semi-major and semi-minor axes. Thus, the uncertainty is : $\text{Δ}\chi =\frac{\chi }{b}\left(\frac{\text{Δ}a}{\chi }+\text{Δ}b\right)=\frac{{\delta }_{*}}{b}\left(1+\chi \right)$.

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₀ = ∞ 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 quasi-stationary vortex. Based on these conclusions we suggest that numerical simulations on vortices should resolve at least ~ 70 and 25 times the pressure scale height in the radial and azimuthal directions, respectively. If this condition is not satisfied, spurious decay of the vortices is likely.

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 t0.

Appendix E.2 Comments for Section 5

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

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 t0.

Appendix F Potential and Forces due to SG

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

Fig. F.1 SG: potential and forces at t = 140 t0. Left: SG potential (ΦSG) and iso-contours. Middle: SG radial force (fr,SG). Right: SG azimuthal force(fθ,SG).

References

All Tables

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 end-of-run surface density for increasing SG (t = 300 t0). Top row and from left to right: q0 = ∞, 4, 2; bottom row and from left to right: q0 = 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 end-of-run Rossby number for increasing SG (t = 300 t0). Top row and from left to right: q0 = ∞, 4, 2; bottom row and from left to right: q0 = 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: semi-minor axis with respect to the pressure scale height at the vortex core for q0 = 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 Bulk-vortex evolution for different values of the SG-parameter. Top: radial distance to the star. Middle: bulk of the vortex mass; Mi is the vortex mass at t = 0. Bottom: absolute value of the mean Rossby number $\overline{\text{Ro}}$. 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 self-gravitating vortices as a function of q0. Left: when estimated from numerical simulations and right, when deduced from Eq. (9). Vortices are stable above the surface and unstable below (blue and orange, respectively). Both surfaces were projected on the left-hand side and the right-hand side of the box. χ and δ are, respectively, the vortex aspect ratio and the radial extent at t = 50 t0.

In the text

	Fig. 7 Rossby number at t = 140 t0 for an isothermal and a non-isothermal disc with SG. Radial zoom with r ∈ [46, 54] AU. Top: self-gravitating vortex in an isothermal disc. Middle: self-gravitating vortex in a non-isothermal disc. Bottom: reference vortex, without SG and in a non-isothermal disc.
In the text

Fig. 8 Density and velocity-divergence at t = 140 t0 for an isothermal and a non-isothermal disc (with and without SG). Left column: relative density (σ/σ0). The main vortex is masked to distinguish the annular over-density and the spiral waves. Right column: normalised relative velocity field divergence (∇ · v′/(2Ωk)). From top to bottom: self-gravitating vortex in an isothermal disc, self-gravitating vortex in a non-isothermal disc, and the reference vortex (without SG and in a non-isothermal disc), respectively.

In the text

Fig. 9 Radial and azimuthal profiles of the Rossby number and the relative density at t = 140 t0 for a non-isothermal 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 Iso-contours of the Rossby number in a non-isothermal disc at t = 140 t0. Left: without SG. Right: with SG.
In the text

	Fig. A.1 Compressibility waves, ∇ · v′/ (2Ωk), emitted by a quasi-steady 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 (χ, δ, q0) space. Each point is an intermediate resolution simulation : (Nr, Nθ) = (1500, 3000). χ and δ are, respectively, the vortex aspect ratio and the radial extent at t = 50 t0.
In the text

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

	Fig. C.2 Definition of aspect ratio, χ and radial extent, Δr, fitting the Ro = −0.8 iso-contours.
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, mcore = 1.69 MJ. Top right: µ + 1 = 1.15, mcore = 0.78 MJ Bottom left: µ + 1 = 1.25, mcore = 0.28 MJ. Bottom right: µ + 1 = 1.25, mcore = 0.18 MJ

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 t0.
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 t0.
In the text

	Fig. F.1 SG: potential and forces at t = 140 t0. Left: SG potential (ΦSG) and iso-contours. Middle: SG radial force (fr,SG). Right: SG azimuthal force(fθ,SG).
In the text