© ESO 2021

1 Introduction

Planets are born and grow in accretion disks around young stars. This is supported by observations of protoplanets embedded in a disk of gas and dust captured during their growth phase (e.g., Keppler et al. 2018). A protoplanet interacts with the disk around it in every stage of its growth (Kley & Nelson 2012), for example via exchange of angular momentum. This results in the launching of spiral arms (Ogilvie & Lubow 2002); if the planet is massive enough, the opening of a gap; and, in some cases, the formation of multiple rings around the planet’s orbit (Rafikov 2002). The number of spirals, gaps, and rings as well as their contrast scales with the planet mass, such that Jupiter-sized planets can have a strong impact on their environment in the right conditions, possibly resulting in multiple ring-like and nonaxisymmetric observable features (Zhang & Zhu 2020; Miranda & Rafikov 2020a). This makes the planet–disk interaction scenario a popular interpretation of such features in the numerous high-fidelity ALMA observations.

One promising scenario to explain observational asymmetries is the existence of vortices because they naturally accumulate dust at the pressure maxima in their center (see e.g., Marel et al. 2013; Bae et al. 2016; Pérez et al. 2018; Hammer et al. 2019; Barge & Sommeria 1996). Among the various ways to form vortices, the Rossby-wave instability (RWI, Lovelace et al. 1999) is particularly relevant in the vicinity of gaps. The RWI readily happens in 2D disks at the outer and inner edge of planet-opened gaps (Li et al. 2005; De Val-Borro et al. 2007). Additional mechanisms that could be relevant in this context are the subcritical baroclinic instability (SBI, Klahr & Bodenheimer 2003; Lesur & Papaloizou 2010) and the zombie-vortex instability (ZVI, Marcus et al. 2015, 2016). Vortices are then susceptible to viscous spreading as well as secondary instabilities such as the elliptical instability (Lesur & Papaloizou 2009), which cause vortex decay. The lifetime of vortices is therefore determined by a competition between vortex-forming and -decaying mechanisms.

Aside from possibly causing observable features in the disks, vortices can also affect planet migration in a stochastic fashion (Regály et al. 2013; Ataiee et al. 2014; McNally et al. 2019) and even cause temporary outward migration (Lega et al. 2021) for otherwise inwardly migrating planets. Understanding their formation pathways and lifetimes is therefore critical to the modeling of planet migration using global, low-viscosity simulations.

In previous numerical studies, vortex properties have been found to depend on various physical processes such as turbulent viscosity and disk self-gravity. Lower viscosity allows vortices to live longer (Godon & Livio 1999; De Val-Borro et al. 2007; Ataiee et al. 2013; Fu et al. 2014; Regály et al. 2017) whereas the inclusion of self-gravity tends to weaken vortices, shortening their lifespan (Lin & Papaloizou 2011; Zhu & Baruteau 2016; Regály & Vorobyov 2017; Pierens & Lin 2018).

In recent numerical studies, radiative effects have been discovered to have a significant impact on the gap-opening capabilities of planets and therefore the structure of the gaps themselves (Ziampras et al. 2020b; Miranda & Rafikov 2020b), affecting the development of the RWI and by extension vortices around their edge (Tarczay-Nehéz et al. 2020). The aim of the present study is to investigate the role of radiative effects for properties of vortices created by planets. More precisely, we explore how the thermal relaxation timescale of the gas affects the lifetime of vortices created during the growth of Jupiter-sized planets.

The impact of thermal relaxation on vortex formation and lifetime was studied for nearly inviscid disks by Les & Lin (2015), and recently Fung & Ono (2021) ran 2D shearing box simulations of RWI-induced vortices. As their simulations did not include a planet, the RWI was triggered by an artificial density bump. These authors described a baroclinic effect that spins down vortices where the decay is fastest for thermal relaxation times of the order of a tenth of the vortex turnover time.

We ran a suite of global two-dimensional hydrodynamics simulations with an embedded Jupiter-sized planet – which naturally creates vortices in the disk – for different choices of turbulent viscosity and thermal relaxation timescales, among other physical parameters. The results of these simulations are then post-processed with our newly developed pipeline for the detection and characterization of vortices.

In Sect. 2 we describe our physical model and numerical setup. We present a typical life track of a vortex in our models in Sect. 3, report the dependence of vortex properties on physical parameters in Sect. 4, and present the case of long-lived vortices in Sect. 5. We discuss and comment on our findings in Sect. 6. Finally, Sect. 7 contains a summary of our main results and our conclusions.

2 Physics and numerics

In this section, we describe the physical and numerical framework that we used in our simulations. We justify the approximations in our model, explain in detail the initialization process, and list technical parameters such as our grid setup and parameter space.

2.1 Hydrodynamics

We consider a thin disk of neutral, ideal gas with adiabatic index γ = 7∕5 and mean molecular weight μ = 2.353 that is orbiting around a star with one solar mass M_⋆ = M_⊙. The two-dimensional, vertically integrated Navier-Stokes equations in a polar coordinate system {r, ϕ} read (1a) (1b) (1c)

where u = (u_r, u_ϕ) and ε are the velocity and specific internal energy of the gas evaluated at the midplane, and Σ is the surface density. The vertically integrated pressure p is defined through the ideal gas law p = (γ − 1)Σε = R_gΣT∕μ, with R_g being the gas constant and T the gas temperature. The isothermal sound speed of the gas is then given by and relates to the adiabatic sound speed c_s as . For a disk in Keplerian motion and vertical hydrostatic equilibrium, we can also write c_s,iso = HΩ_K, where is the Keplerian orbital frequency at radius r and H is the pressure scale height of the gas.

The viscous stress tensor σ (following Tassoul 1978) appears in both the momentum Eq. (1b) and the dissipation function: (2)

where ν = αc_sH is the kinematic viscosity parametrized according to the α-viscosity model of Shakura & Sunyaev (1973). Here, α is a parameter that captures both radial angular momentum transport – that leads to accretion onto the star – and heating of the disk due to viscous friction. Numerical simulations of (magneto)hydrodynamical instabilities such as the vertical shear instability (VSI, Nelson et al. 2013) or the magneto-rotational instability (MRI, Balbus & Hawley 1991) have provided numerical estimates of α, while observations of young stellar objects surrounded by disks have constrained these estimates (Dullemond et al. 2018). To probe a wide range of diffusion regimes from practically inviscid to moderately viscous, we choose α ∈ {10⁻⁶, 10⁻⁵, 10⁻⁴, 10⁻³} for our models.

Viscous dissipation leads to the heating of the disk. An embedded planet can also deposit significant amounts of thermal energy via the dissipation of spiral shocks (Rafikov 2016; Ziampras et al. 2020a). As a cooling solution, we allow the disk to relax to a prescribed temperature profile T₀ (see, Eq. (5)) over a relaxation timescale τ_relax = β∕Ω_K (Gammie 2001).

The thermal relaxation term appears as an additional source term to the energy equation (3)

where is the heat capacity of the gas at constant volume. The parameter β controls the relaxation timescale, as well as the overall planet–disk interaction process (Miranda & Rafikov 2020b); we choose the values β ∈{0.01, 1, 100} which correspond to very fast, moderate, and very slow relaxation.

The gravity of the star and planet are included as a source term in g. We work in a star-centered coordinate system and embed a planet with mass M_p at a position r_p. Thus, the source term reads (4)

The terms g_⋆, g_p, and g_ind denote the acceleration due to the star, the planet, and the indirect term which is a correction needed because the star-centered frame is not an inertial frame. As we are considering fixed, nonmigrating planets, disk feedback on the star and planet is neglected. The planet’s gravitational pull (second term in the RHS of Eq. (4)) is smoothed using a Plummer potential with a smoothing length ϵ = 0.6H(r) that captures theeffect of the vertical structure of a more realistic 3D disk (Müller et al. 2012) and prevents singularities near the location of the planet.

For simplicity, we do not allow the planet to migrate. We chose to limit the degrees of freedom in our model to focus on the dynamics of the vortex and avoid complex and potentially chaotic interplay between the vortex and the planet (Lega et al. 2021). For the same reason of simplicity, we neglect planetary accretion in our models.

2.2 Numerics

We use two different codes for our numerical models: PLUTO 4.2 (Mignone et al. 2007), a finite-volume, energy-conserving, shock-capturing code that treats transport by solving the Riemann problem across the interfaces of adjacent cells in both directions (r, ϕ) in an unsplit fashion; and our custom FARGO (Masset 2000) version, FargoCPT (Rometsch et al. 2020), which uses a finite-difference, dimensionally split, second-order upwind method for gas advection. Both codes utilize the FARGO method (implemented into PLUTO by Mignone et al. 2012), in which orbital advection is essentially performed via the Keplerian rotation on top of which the code solves for the residual velocity deviations, significantly relaxing time-step limitations and reducing numerical dissipation in the process (Masset 2000).

The inherent differences between the two numerical schemes make it worthwhile to carry out our simulations using both codes, to verify the robustness of our results, and to test for numerical convergence. Namely, the strictly energy-conserving nature of PLUTO and the necessity for artificial viscosity to stabilize FargoCPT are discussed in more detail in Sect. 6.7, among others.

2.2.1 Grid setup

Our computational domain spans the full azimuthal extent and a radial range of r ∈ [0.2, 5.0] r_p = [1.04, 26.0] au, with square cells logarithmically spaced so that the cell aspect ratio is preserved. After carrying out a thorough investigation of the effects of our numerical resolution of the recovery of both radial and azimuthal features caused by the planet, we decided to execute our simulations using a resolution of 8 and 16 cells per scale height (hereafter “cps”) in both directions (r, ϕ). Because we use a constant aspect ratio together with a logarithmically spaced radial grid, the resolution in cps is constant throughout the domain.

At this resolution, the two codes reach good convergence in terms of the presence and contrast of features shaped by the planet and results agree between the codes. This translates to a fiducial resolution of (N_r, N_ϕ) = (528, 1024) cells for 8 cps, or (1056, 2048) cells for 16 cps. In addition, using the same resolution in both directions in terms of cps ensures that the effects of numerical viscosity are isotropic (see Appendix A).

2.2.2 Initial and boundary conditions

Our disk is initially axisymmetric and in equilibrium in the radial direction, such that the initial radial velocity profile results in a constant accretion rate through the disk. The azimuthal velocity is close to the Keplerian profile, with the correction due to the radial pressure gradient. The initial surface density and temperature profiles are simple power laws such that (5)

with r_p = 5.2 au. This temperature profile translates to a disk with a constant aspect ratio h(r) = H∕r = 0.05. While the general consensus is that protoplanetary disks are flared (i.e., the aspect ratio increases with distance, see e.g., Dullemond 2000), we choose to use a constant aspect ratio because the behavior and lifetime of vortices depends on this quantity (Hammer et al. 2021). Thus, we can isolate the dependence of vortices on the physical and numerical parameters in our suite of simulations.

The radial and azimuthal velocity components at t = 0 are then (6)

Near the boundaries, within the radial extent r ∈ [0.2, 0.25] ∪ [4.2, 5.0] r_p, the surface density and velocity are both damped to their initial profiles (see Eqs. (5) and (6)) using the method of De Val-Borro et al. (2006) over a damping timescale of 0.3 periods at the respective boundary. While the radial boundary edges are closed, this minimizes the reflection of spiral waves back into the computational domain. The boundaries are periodic in the azimuthal direction.

We then embed a Jupiter-sized planet (M_p = 1 M_J = 10⁻³M_⋆) in most models, with some simulations instead containing a less massive planet of M_p = 0.5 M_J. To smoothly introduce the planet into the disk, we typically allow the planet to grow over 100 orbits at r_p using the formula by De Val-Borro et al. (2006). The importance of the growth timescale and planet mass are discussed in Sect. 4.

2.3 Vortex detection

We use the gas vortensity (7)

– where ẑ is the unit vector in the vertical direction – as a proxy to detect and track the evolution of vortices over hundreds of snapshots for every model. As these vortices consist of anticyclonic motion, the center of a vortex corresponds to a local minimum in vorticity, ω = (∇ × u) ⋅ẑ. Because vortices tend to accumulate mass towards their center and Σ is enhanced inside the vortex, the transition from the background flow to the vortex region is stronger and the vortex is more easily identified in a map of ϖ than in the case of ω alone.

More precisely, we use the gas vortensity normalized by the background vortensity from the initial conditions, ϖ₀ = (∇ × u_K) ⋅ẑ ∕ Σ₀. This eliminates the radial dependence of the Keplerian velocity and the disk’s surface density and ensures that our vortex proxy quantity, ϖ∕ϖ₀, is of order unity everywhere in the disk except for the gap region due to its very low surface density. The quantity ϖ∕ϖ₀ usually varies between −1, for strongly counter-rotating vortices, and 1 for the background flow.

We use our new Python module, Vortector, which extracts iso-vortensity contours using the computer vision library OpenCV (Bradski 2000) to detect vortex candidates and then fits a 2D Gaussian to the vortensity and surface density data. The FWHM (or 2.355σ) of this Gaussian is used to define the radial and azimuthal extent of a vortex. Using this method, we also extract information about the shape of the vortex, including its radial and azimuthal extent and the mass it encloses. A more detailed description can be found in Appendix B.

One drawback is that this automated process sometimes produces detection artifacts, as can be seen for example in Fig. 6 below (top panel, dashed orange line), such that the vortex size (and therefore its mass) is overestimated near the end of its lifetime as it blends into the disk background. While this effect is partly counteracted by using a median filter in time, we do not manually edit the output of the Vortector on a model-by-model basis.

In the following three sections, we present the results of our simulations. First, we present a typical example of vortex formation and evolution (Sect. 3). We then go on to describe the dependence of vortices on physical parameters for the group of vortices with short and intermediate lifetimes (Sect. 4). Finally, long-lived and migrating vortices are presented (Sect. 5).

3 Typical life track of a vortex

The Jupiter-sized embedded planet opens a deep gap in all of our simulations. Figure 1 shows maps of Σ (left) and ϖ (right) normalized by their initial values at five time-stamps during the vortex lifetime for a model with α = 10⁻⁵, β = 1, and a resolution of 8 cps performed with the FargoCPT code. Horizontal dotted lines at r = 1.45 r_p are superimposed as a reference marking the final location of the vortex center. Here, four small-scale vortices (top row) first merge into two slightly larger vortices (second row) and then finally into one massive vortex (middle) that will last for a little over 1100 orbits. The vortex slowly decays over time, maintaining a large size (fourth row). In the later stages, the vortex is no longer present (bottom row). The nonaxisymmetric structure still visible exists due to the planet’s spiral arm and is corotating with the planet.

Early in the gap-opening process, the outer gap edge grows Rossby-wave unstable (Lovelace et al. 1999) and several small-scale vortices form around it (top two rows). Figure 2 shows radial profiles of the Lovelace parameter, , (top) and Σ (bottom) at different time-stamps during the vortex formation up until t = 180 orbits. The Lovelace parameter is defined as (8)

with the entropy S = P∕Σ^γ. The development of a maximum in is visible, which is one condition for the onset of the RWI. Vertical lines at the center of the maxima (as determined by eye) are added to both panels to guide the eye for a comparison of the location of the maxima in and Σ at each time-stamp. The maxima are located on the slope of the gap edge slightly inward from the Σ maxima and coincide with the location of the small vortex centers. The maximum in moves outward following the maximum in Σ as the gap opens. This illustrates that the vortices form due to the RWI at the slope of the outer gap edge.

In the absence of self-gravity, these small vortices then quickly merge together (within ~ 100 planet orbits) into a single large vortex that slowly moves outwards following the gap edge as the gap deepens and widens (third and fourth rows in Fig. 1). The surviving vortex then typically decays over ~200–2000 orbits. The evolution of three vortex properties is illustrated in Fig. 3, where we show, from top to bottom, the mass M_vort enclosed within the FWHM ellipse of the 2D Gaussian fit to Σ, the location of the center of the vortex r_vort and the radial FWHM width Δr as the shaded area, and the vortensity at the vortex center normalized by the azimuthal median. The vertical dotted lines indicate the time when the planet reached its full mass (typically 100 orbits). A short phase of vortex formation is followed by a slow and steady decay process, as can be seen in the decrease of mass and radial size. Because the vortensity contribution of the anticyclonic vortex is negative, an increase in vortensity indicates a decay as well. The line in the bottom panel of Fig. 3 is continued (in orange) for another 100 orbits after the vortex has decayed – according to our criterion presented below in Sect. 4.1 – in order to illustrate the return of the curve to 1, which corresponds to an azimuthally symmetric state.

During its lifetime, the vortex can become as large as Δr = 0.4 r_p (2 au for r_p = 5.2 au) with a typical vortex aspect ratio (rΔϕ∕Δr) of 6-10. Its mass, M_vort, is typically some tenths of M_J but can be as large as one M_J, with a surface density enhanced by a factor of up to seven compared to the initial value.

The vortices form around the location where the radial Σ profile reaches 10%of its initial value (see bottom panel of Fig. 2), which we define as the gap edge in a similar way to Crida et al. (2006). During their lifetime, most vortices tend to stick to this gap edge in the sense that their inner boundary, r_vort − Δr∕2, roughly coincides with the location of the gap edge. For some models, we observe that the vortex detaches from the outer gap edge after several hundred orbits and starts migrating outward. These models are discussed in Sect. 5.

Fig. 1 Multiple snapshots of the α = 10−5, β = 1, 8 cps model showcasing the vortex merging process during the early stage of gap opening, the fully grown size of the resulting vortex, and its subsequent decay. The surface density and vortensity contrast compared to their initial profiles is shown in the left and right panels, respectively. Time is quoted in units of planetary orbits. The horizontal line at r = 1.45 rp serves to highlight the outward radial movement of those structures as the gap around the planet grows wider. The planet is located at r = 1 rp and ϕ = 0.

Fig. 2 Evolution of radial Lovelace parameter (see Eq. (8)) and Σ profiles during vortex formation over the first 200 orbits of the sample case from Sect. 3. The vertical lines indicate the center of the plateau in (estimated by eye) to guide the eye to the corresponding location of the Σ profile. is calculated as the azimuthal average at each radius. The dotted horizontal line in the bottom panel marks 10% of Σ0, which we define as the location of the gap edge.

Fig. 3 Evolution of vortex properties for the showcase simulations presented in Sect. 3. The panels show, from top to bottom, the mass enclosed in the FWHM ellipse of the vortex fit Mvort in Jupiter masses, the radial location of the vortex rvort, and its FWHM Δr, indicated by the shaded area, and the ratio between minimum vortensity inside the vortex and the azimuthal median of vortensity at the radial location of the vortensity minimum. A dotted vertical line indicates the time when the planet has reached its final mass. The curves are smoothed with a median filter that spans over the preceding and following five datapoints (± 50 orbits at rp). The orange parts of the line in the bottom panel show the evolution of the vortensity prior to the “birth” and after the “death” of the vortex.

4 Dependence of vortex properties on physical parameters

Having described a typical lifetrack of a vortex in our simulations, we now present the effects of different physics and numerics on vortex lifetime, location, and impact on the overall disk structure. The model parameters are listed with the main results in Table C.1.

4.1 Vortex lifetime

We define the vortex lifetime as the time difference between its “birth” and “death” by analyzing the ratio of ϖ to the azimuthal median value as a function of time. The normalization with instead of ϖ₀ is done to eliminate the ϖ evolution of the background disk due to changes in Σ and radial pressure gradients, which affect (∇ × u) ⋅ẑ by changing the azimuthal velocity.

The “birth” is identified as the time when drops from its initial value of 1 (for an axisymmetric disk) down to lower values (see bottom panel of Fig. 3). Because drops already for small vortices, the lifetime also includes the stage where there are multiple small vortices (see Sect. 3).

The “death” of the vortex on the other hand is less obvious to identify. At this stage, usually slowly rises back to the value of the background disk. Usually, there is a “knee” visible in at or slightly after the point in time where the vortex dies and where ϖ approaches the background flow (see the orange part of the line in the bottom panel of Fig. 3 where is continued for another 100 orbits after the vortex disappeared at t = 1250 orbits). For some models, this knee is not visible, and we manually inspect the 2D contour plots of and identify the point after which no further closed iso-value lines (with spacing in of 0.05) are present. As an additional measure for less obvious cases, we analyze the gas streamlines at different time-stamps.

In our models, the drop in happens in a matter of tens of orbits. A conservative estimate for the uncertainty of this birth time measurement is 50 planetary orbits. From applying the manual method to models where the knee exists in the curve (implying the death of the vortex), we estimate a conservative uncertainty to be 100 planetary orbits. This leaves a total uncertainty of 150 planetary orbits for the lifetime of our vortices.

The lifetime of vortices in our grid of simulations is shown as an overview in Fig. 4. The left and right panels show vortex lifetimes as a function of β for 8 and 16 cps, respectively. The viscous α is encoded in color, and the symbol indicates the simulation code and the inclusion of self-gravity. For each value of α and β, we calculated the average (“avg(f,p)”) between the two codes (not including the self-gravity models) when the results are close together. The solid-colored lines connect the averages to help visualize the trends. For parameters for which the two codes showed different vortex lifetimes, we added separate lines connecting the average to the FargoCPT and PLUTO results to highlight the differences.

Lifetimes range from approximately 200–2000 orbits for the shorter-lived vortex group up to at least 15 000 orbits for the long-lived vortices discussed later in Sect. 5. The most prominent features of the distribution are the trend of decreasing lifetime with increasing α and the minimum vortex lifetime at β = 1 for low α and high resolution. In the following sections, we address the influence of our model parameters on vortex lifetime.

Fig. 4 Lifetime of vortices as a function of β for 8 cps (left) and 16 cps resolution (right). Colors encode α, and the different symbols denote the code and the inclusion of self-gravity. The solid lines help guide the eye and connect the lifetime averages between the two codes (without self-gravity) for each value of α, where thetwo codes agree sufficiently. For parameters where there is a difference between the codes, dashed and dotted lines connect to the datapoints of the FargoCPT and PLUTO runs, respectively. A “↦ ” next to a symbol marks models that were terminated due to runtime constraints but still contain an active vortex. The horizontal gray line in the right panel indicates the top of the y-axis of the left panel. A list of all vortex lifetimes shown here is provided in Table C.1.

4.2 Influence of the thermal relaxation timescale

The dimensionless thermal relaxation timescale β has a strong effect on vortex lifetime. For α = 10⁻⁴, lifetimes are of the order of several hundred to 1000 orbits with a downward trend as β increases. Vortex lifetimes are shortest for β = 1 (around 1250 orbits) and increase towards both sides to around 2000 orbits for β = 100 and to values of the order of 10 000 orbits for β = 0.01, high-resolution runs. This decrease in lifetime for nonisothermal disks is consistent with the results of Tarczay-Nehéz et al. (2020). Exceptions to this trend are the 8-cps PLUTO models for α = 10⁻⁵–10⁻⁶ and β = 100. We were not able to identify the reason why the two codes did not agree for these parameters, but we note that the two codes match well once again for 16 cps in the same configurations. Models with very long vortex lifetimes are analyzed in Sect. 5.

Fung & Ono (2021) reported a similar trend in vortex lifetime in 2D shearing-box simulations without planets, in which the vortex was introduced by initializing the simulation with a radial density bump. These authors found that vortex decay is fastest for intermediate β in the range 1–10, but their disk model assumes a constant background disk without gradients in T and Σ, which change baroclinic effects. Our results indicate that a similar mechanism might be at play in the presence of an embedded planet with strong spiral-arm shocks. However, the strong enhancement of vortex lifetime for β = 0.01 hints at the presence of an additional mechanism which keeps the vortices alive. We discuss these hypotheses further in Sect. 6.1.

For a comparison of vortex evolution at different β, see Fig. 5, where the evolution of vortex properties (analogous to Fig. 3) of three FargoCPT simulations at 8 cps resolution with β = 0.01, 1, and 100 is shown. The absolute radial location of vortices varies with β as well. This is due to the tendency of the vortices to form and subsequently stick to the outer planet gap edge and the gap-opening process being strongly influenced by β. Miranda & Rafikov (2020b) showed that “extreme” values of β (i.e., β → 0 or β→∞) result in narrower planet-opened gaps but additional gaps in the inner disk, whereas intermediate values of β ~0.1–10 lead to a single, wide gap around the orbit of the planet. In our simulations, models with β = 0.01 show the widest gaps, narrower gaps are present for β = 1, and β = 100 models showed an even slightly narrower gap. This is reflected in the vortex locations which are further in for higher β (see the center panel of Fig. 5). The difference between our results and those of Miranda & Rafikov (2020b) might be due to the presence of the vortex.

4.3 Planet growth timescale

Hammer et al. (2017) observed that the lifetime of planet-induced vortices can depend on the timescale over which the planet mass is increased in order to introduce the planet into the simulation. These authors found that vortex lifetime decreased with a longer planet growth time. In our models, increasing the planet growth timescale from 100 to 1000 orbits caused vortices to live longer by 470 orbits for β = 1, and by up to 1900 orbits for β = 0.01. Figure 6 shows the evolution of vortex quantities comparing the FargoCPT runs with a τ_ramp = 100 orbits, already presented in Fig. 5, with their respective counterparts with τ_ramp = 1000 orbits. The curves of runs with τ_ramp = 100 orbits are shifted to the right by the difference in lifetime Δt compared to the respective τ_ramp = 1000 orbits model. This shift clearly illustrates that the decay of these vortices is almost the same for both values of τ_ramp in terms of their mass, location, and vortensity curves. The only difference caused by the planet injection timescale appears in the time it takes for the vortex to reach the turnover point, after which it starts to decay.

Fig. 5 Evolution of vortex properties for varying values of the thermal relaxation β parameter. The panels are as in Fig. 3. Shown are models run with the FargoCPT code with α = 10−5 and at 8 cps resolution (orange “f8” dots in Fig. 4).

4.4 Planet mass

From our M_p = 0.5 M_J models we cannot draw any conclusions regarding the dependence of vortex lifetime on planet mass, because for the set of parameters β = 0.01 and α = 10⁻⁵ the vortices are long-lived outliers like the ones discussed in Sect. 5. However, the location of the vortex is also influenced by the planet mass. Lower-mass planets open narrower gaps and cause the location of the vortex to be further in compared to more massive planets because this location is linked to the location of the gap edge. In our models, the vortices in the M_p = 0.5 M_J were located ~0.15 r_p closer to the star.

4.5 Viscosity

The observed vortex lifetime typically increases with lower values of α. Simulations with α = 10⁻³ show only small vortices forming. These disappear within 100 orbits, and are therefore already gone by the time the planet has grown to its full mass. For models with α = 10⁻⁴, we observe vortex lifetimes of up to around 1000 orbits.

Simulations with a lower viscosity (α = 10⁻⁶–10⁻⁵) show even longer lifetimes, usually in the range between 1000 and 2000 orbits, excluding the outliers that we discuss later in Sect. 5. For this range of α, vortices usually have similar lifetimes for simulations sharing the same β value. An example of this is shown in See Fig. 7, which shows, from top to bottom, the evolution of the mass enclosed in the region of the vortex (FWHM), the location and radial extent (in FWHM) of the vortex as determined by the surface density fit, and the ratio of normalized vortensity to the azimuthal median of the latter at the location of the vortex.

The vortex location is not influenced by viscosity. Although the gap-opening time is according to the estimate in Kanagawa et al. (2017), the bulk of the gas in the vicinity of the planet is cleared within the first few hundred orbits. During this time, Σ is lowered by two orders of magnitude within the gap region, and the radial gradient of Σ becomes steep enough to facilitate vortex formation.

Fig. 6 Influence of the planet introduction time on the evolution of vortex properties. The panels are as in Fig. 3. Solid and dashed lines show models with τramp = 100 orbits and 1000 orbits, respectively. The τramp = 100 orbits curves are shifted to the right (see the horizontal lines) to illustrate that the curves have the same shape in the decay phase, independent of τramp. We note that the final evolution of the vortex, after it has reached its minimum in vortensity, is the same independent of planet introduction time.

4.6 Self-gravity

Several studies showed that vortices in weakly or strongly self-gravitating disks might not grow as large because small vortices do not merge into one large vortex (Lin & Papaloizou 2011) and dissipate more rapidly because of stretching in the azimuthal direction (Lovelace & Hohlfeld 2013; Regály & Vorobyov 2017; Zhu & Baruteau 2016). This can be the case even for low-mass disks as long as the Toomre stability parameter Q is lower than 50 or . For the choice of parameters in our models, the Toomre parameter is (), dropping under 5 at roughlyr∕r_p = 2.8. To check the effect that disk self-gravity has in our models, we ran additional simulations with FargoCPT with self-gravity activated for all three values of β (0.01, 1, and 100) and α = 10⁻⁵.

The lifetimes of vortices in these simulations are shown in Fig. 4 as the rightmost datapoint in each column (models “f8sg”). An example evolution of their properties is shown in Fig. 7.

Self-gravity inhibits the merging of the small initially formed vortices into one large vortex. Instead, two smaller vortices usually remain until they decay. This leads to a significantly shorter lifetime compared to the analogous simulations without self-gravity, which is consistent with the findings of the studies mentioned above. However, this does not apply to the long-lived migrating vortices discussed in Sect. 5.

Figure 7 shows that the center of the vortex in a model with self-gravity and β = 1 is further in compared to its nonself-gravitating counterpart. This is due to the smaller radial extent of the vortex in the run with self-gravity and the tendency of the inner edge of each vortex to coincide with the gap edge. Because self-gravity does not noticeably change the radial disk profile for the mass regime of our models, the inner edge of the vortices is at the same location, independent of whether self-gravity is included or not. The same effect is also observed for β = 0.01 and β = 100.

Fig. 7 Evolution of vortex properties for varying values of α. Panels are shown as in Fig. 3. Shown are models run with the FargoCPT code with β = 1 and at 8 cps resolution. The α = 10−3 run is excluded because no vortex forms. In addition, a run with disk self-gravity enabled is added for the α = 10−5 case. The similarity between simulations with α = 10−5 and 10−6 is apparent.

5 Long-lived and migrating vortices

In some of the cases, a much longer lived vortex is observed. In these models, vortices stay close to their peak mass for several thousand orbits, and in some cases migrate outwards after having stayed at the planet gap edge. This happens only for very low viscosities (α ≤ 10⁻⁵) and β = 0.01 or locally isothermal simulations (β → 0). For our standard M_p = 1 M_J planets, the long-lived outliers appear only at the highest resolution of 16 cps but not at 8 cps. For the corresponding M_p = 0.5 M_J model, the long-lived vortex also appeared at 8 cps. Spiral arms launched by the vortex are clearly visible for these long-lived large vortices (see Fig. B.2). They are more pronounced for lower values of α.

Figure 8 shows the evolution of vortex properties for a selection of models to highlight the observed behavior. The most prominent example is the model with α = 10⁻⁵, β = 0.01, and a 16 cps resolution. The vortex in those runs lived for 15 100 orbits before we terminated the two simulations because of their long runtime. Both codes, PLUTO and FargoCPT, agree well for the long-lived cases. Specifically, they are in exceptionally close agreement for α = 10⁻⁵ and only differ at later stages for α = 10⁻⁶ (see orange and blue lines in Fig. 5). We do not currently fully understand the mechanism that allows these long-lived vortices to sustain themselves for such long timescales. We attempt to provide a speculative explanation in Sect. 6.1.

For β ≠ 1, a secondary radial density and pressure bump is observed in the outer disk. This is the result of the vortex generating spiral arms which transport angular momentum. Radially outwards, this results in the accumulation of mass in a second bump (see panels for β ≠ 1 in Fig. 9). This does not happen for β = 1 because of the less efficient angular momentum transport by spiral arms for this intermediate value of β (Miranda & Rafikov 2020b). For β = 0.01, some models show vortices migrating radially outwards (e.g., the α = 10⁻⁶ models in Fig. 8). This is likely related to the formation of the secondary bump outside of the vortex location (see Fig. 9) and the fact that vortices typically migrate towards pressure bumps (Paardekooper et al. 2010). For vortices that migrate far enough outside, which only happens for β = 0.01, a weaker secondary vortex appears between them and the edge of the planet-generated gap (see Fig. B.2). These secondary vortices then decay over a few hundred orbits, already having decayed by the time the primary vortex disappears. While they are treated as independent entities, these secondary vortices are not included in Fig. 4 or the discussion above; their occurrence is likely the result of a multistage process which begins with the secondary bump forming and the primary vortex migrating radially outwards towards it and meanwhile supplying mass towards the edge of the planet-generated gap. This then feeds the emerging “secondary” vortex.

Fig. 8 Selection of models with long-living and migrating vortices at 16 cps resolution for the two different codes. Both codes agree remarkably well for the blue and orange cases. The panels are as in Fig. 3. In models shown here, β = 0.01. The values of α, the resolution, and the code used are indicated in the legend.

Fig. 9 Azimuthally averaged surface density profiles as a function of different physical (α, β) and numerical (cps) parameters at two different time-stamps. The peak around r∕rp = 1.5–1.6 corresponds to the pressure bump formed by the planet as the latter pushes material away, forming a gap around its orbit. The smaller, secondary peak at around r∕rp = 2.1 is caused by the vortex that forms near the primary, planet-generated bump. Top: radial profiles at t = 1000 orbits. At this stage, all models pictured feature a vortex near the primary bump. We note the absence of a secondary bump for the models with β = 1. Bottom: same profiles at t = 5000 orbits. Here, the primary bump has moved radially outwards as the planet’s gap gets deeper and wider. We highlight the depletion of gas near the primary pressure bump for the second panel from the left (β = 0.01). This is caused by the combination of a vortex migrating outwards to the secondary bump, and the inability of the planet to resupply that zone with material from its now-depleted gap region. We also note the difference between resolutions of 8 and 16 cps (dashed and solid lines), especially for the β = 100 models and the β = 0.01, α = 10−5 model (the evolution of this model is shown as the orange line in Fig. 8).

6 Discussion

In this section, we address some ways in which our results could be interpreted and their relevance in explaining observations. We also underline some caveats of our models.

6.1 The conditions needed to form and sustain a vortex

To form a vortex, one needs to create a local vortensity extremum. In the absence of nonconservative forces, the evolution equation for the vortensity in a 2D flow reads (9)

where is the baroclinic term and describes viscous diffusion of vortensity which can lead to vortex decay.

As outlined in the introduction, several instabilities have been discovered that provide a mechanism to form or destroy large-scale vortices, but they all fundamentally rely on Eq. (9) to change the vortensity of the flow. The mechanism responsible for the formation of the vortex in our simulations is most likely the RWI which is triggered during the gap-opening process, as we demonstrated in Sect. 3 and Fig. 2.

To check whether vortices can only form during the gap-opening process and not in the quasi-steady state after the bulk of the gas has been pushed out of the gap region, we removed the long-lived vortex from the α = 10⁻⁵, β = 0.01, 16 cps model by replacing the velocities and Σ with their azimuthal median values for r > r_p during the peak of its activity (t = 1880 orbits). The fact that there is no vortex forming again is an indication that the formation of vortices in our simulations depends on the gap-openingprocess to produce conditions that can trigger the RWI. This is also supported by the observation that the peak in is strongest for an intermediatetime, t = 70 orbits, during the gap opening process, after which the maximum disappears and a plateau in forms.

Vortex decay happens due to at least two mechanisms. Viscous spreading attacks the vortices for high α = 10⁻⁴–10⁻³, as illustrated by the trend of lower vortex lifetime for higher α, and vortex stretching due to self-gravity effects additionally limits the vortex lifetime, if it is considered (Lin & Papaloizou 2011; Zhu & Baruteau 2016; Regály & Vorobyov 2017). For sufficiently low α, another process that depends on β starts to be dominant. We do not fully understand the mechanism but we observed some similarities to the recent work by Fung & Ono (2021). These latter authors find that, in their simulations, vortices decay the fastest for β = 1–10 and decay slower for both smaller and larger β. Vortex lifetime in their simulations changes by up to an order of magnitude depending on β. We also find a minimum in vortex lifetime for β = 1 with lifetimes increasing as β ≠ 1.

Fung & Ono (2021) explained the decay mechanism by asymmetries in the structure of around the vortex center, which they found to be quadrupolar (see their Fig. 6) and to change with β. We also find asymmetries in the structure of , but our simulations differ from theirs in some fundamental aspects. Our simulations are global with radially varying Σ and T profiles and include a planet that continually perturbs the disk, whereas the simulations of Fung & Ono (2021) consider a local shearing sheet with a constant background Σ and T, with only an initial perturbation in the form of a density bump. As a consequence of the radially varying T in our simulations, the structure of around the vortex center is dipolar in the azimuthal direction, as can be expected for a Gaussian-like density maximum. Additionally, the planetary spiral arms strongly influence . Figure 10 shows a 2D map of for two simulations with α = 10⁻⁵ and 16 cps resolution. The left panels show a short-lived vortex with β = 1 at t = 1000 orbits and the right panels show the long-lived vortex model which exhibits the secondary vortex (see Sect. 5 for a description and Fig. B.2 for ϖ and Σ maps at the same time). The actual shapeof the perturbation of inside and around the vortex varies in time because it depends on the phase with respect to the spiral arm. It is currently not clear to us how the changes in structure of lead to the change in vortex decay and how this proposed mechanism depends on the various parameters in our system.

The long-lived group of vortices for low β (see Sect. 5) indicates that there might be another vortex formation mechanism at play. Given that the RWI already caused finite perturbations in the disk and our disks exhibit a radial entropy gradient, the SBI (Klahr & Bodenheimer 2003; Lesur & Papaloizou 2010) seems to be a natural candidate. However, we verified that the SBI is not active in our disks byanalyzing the Richardson number, the ratio of the buoyancy (also called Brunt-Väisälä) frequency to the shear rate, which needs to be negative in a radially extended region over the full azimuth of the disk for the SBI to operate. The Richardson number in our simulation is positive, except for narrow stripes following the spiral arms, which rules out that the SBI is active.

To rule out that the difference in lifetime is a result of the initial vortex formation during gap opening, we took the long-lived vortex out of the α = 10⁻⁵, β = 0.01, 16 cps model and inserted it into the α = 10⁻⁵, β = 1, 16 cps model. Although this artificial vortex has the same structure as in its original β = 0.01 model, it decays over nearly the same time as the standard β = 1 vortex. Thisis an additional indication that the difference in lifetime is caused by the dependence of the decay process on β or a possible additional vortex formation channel that sustains the vortex at low β.

This leaves us with the hypothesis that the interaction of the spiral arms with the vortices might play a major role in either slowing down vortex decay or providing an additional vortex formation channel. This hypothesis is motivated by the strong impact of the spiral arms on and the dependence of spiral-arm properties on β (Ziampras et al. 2020b; Miranda & Rafikov 2020b). Another contribution might be the vortensity jump across the spiral-arm shock, which was recently illustrated to be important for the evolution of vortensity in the case of subthermal-mass planets (Cimerman & Rafikov 2021). Providing an analysis of both mechanisms in our context is unfortunately out of the scope of the present explorative study.

Fig. 10 Baroclinic term (RHS of Eq. (10)) in the outer disk for a short-lived model (β = 1 at t = 1000 orbits) and a long-lived vortex (β = 0.01 at t = 7150 orbits, see also Fig. B.2) with α = 10−5 and 16 cps resolution on the left and right side, respectively. The top row shows maps of the baroclinic term with the detected vortices indicated with green ellipses as obtained from the Σ fit. The bottom row shows the radial Σ profile in orange, the azimuthally averaged baroclinic term in blue and the region between its minimum and maximum shaded in gray.

6.2 Effect of in-plane radiation transport

It has been shown that parametrizing radiative effects with β while omitting the effects of in-plane radiation transport can result in a potentially inaccurate radial surface density structure, mainly in the inner disk and around the gap, due to the impact of β on the capability of a planet to open multiple secondary gaps at r < r_p (Miranda & Rafikov 2020b). Here, we are not interested in the annular structures of the inner disk, and so we chose to ignore in-plane radiation transport. Nevertheless, to check for possible effects of in-plane radiation transport on the vortex dynamics, we repeated the α = 10⁻⁵, β = 1 model at 8 cps. This time we included a flux-limited diffusion (FLD) approach (Levermore & Pomraning 1981) similar to Ziampras et al. (2020a), but by parametrizing the diffusion coefficient D_rad following Eqs. (12)–(14) of Flock et al. (2017): (10)

where l_thin is the photon mean free path. We found that including FLD slightly changes the radial surface density structure in the inner disk as predicted by Miranda & Rafikov (2020b) and reduces the vortex lifetime from 1200 to 900 orbits. Studying the effect of in-plane radiation transport in more detail requires further investigation.

6.3 The assumption of a 2D disk

One of the main limitations of our models is the 2D assumption, which was chosen because of runtime constraints in our rather wide exploration of the parameter space. It is entirely possible that various 3D effects can result in quantitative differences in vortex properties. Three-dimensional vortices can be susceptible to the elliptical instability (Lesur & Papaloizou 2009) which would lower their lifetime. On the other hand, the vertical modes of the SBI could provide an additional channel to sustain the vortices,and vertical gas circulation due to the VSI might interfere with vortex growth and decay (Flock et al. 2020).

To estimate the impact of including full-3D effects, we ran one 3D simulation with FARGO3D (Benítez-Llambay & Masset 2016) using a setup analogous to our 2D setup. We chose β = 2π, α = 0 and a resolution of 8 cps in all three directions. The simulation assumed symmetry about the midplane and covered four scale heights in the vertical direction. Similar to our 2D models, a large vortex formed at the outer gap edge and lived for 7000 orbits. This illustrates that while there are differences, large vortices can survive in 3D disk simulations for a long time, even longer than in 2D for our example. We limited the 3D runs to this one test because its runtime at 8 cps resolution was close to 4 months with the simulation performed on four NVIDIA K80 GPUs.

6.4 Observability of vortices at large radii

Section 4.2 illustrates that vortex lifetime is affected by the thermal cooling timescale β. The latter is expected to vary with radius in a disk, with values of 1–10 at 5 au, 0.1 at ~10 au, and below 0.1–0.01 at ~50 au (Ziampras et al. 2020b). Thus, we expect vortices to be in the short-lived regime close to the star and in the long-lived regime far from the star. From Fig. 4 we can estimate the lifetime of vortices in disks with α ≤ 10⁻⁴ to be between 500 and 3000 orbits for β ≥ 1 and between 1000 and 15 000 orbits for β < 1 for α ≤ 10⁻⁴. Assuming a solar-mass star, this yields estimated lifetimes for a planet-induced vortex in the range of 6–30 kyr at 5 au, 175–700 kyr at 50 au, and 1–15 Myr at 100 au. On the basis of a simple lifetime-centered argument, our results therefore suggest that planet-induced vortices are more likely to be observed at larger radii.

It should be noted that planet growth timescales of 100 and 1000 planetary orbits are at the very low end of the spectrum of physically expected planet growth times. Hammer et al. (2017) provided estimates for more realistic planet growth times of several thousand to tens of thousands of orbits. It remains to be seen whether the effects observed in this study still appear for longer, more realistic planet-growth timescales. However, simulating the disks at the required resolution of at least 16 cps for longer planet-growth times along with the additional vortex evolution time is still computationally expensive.

6.5 Using the lifetime of vortices in simulations to explain observations

In the suite of simulations we carried out, the lifetime of vortices in models with identical physical parameters varies significantly with resolution. This was the case for low values of the viscous α parameter (α = 10⁻⁵, 10⁻⁶). We argue that the numerical viscosity of our simulation codes is comparable to α_num ≲ 10⁻⁵. This suggests that simulations with a prescribed viscosity of the order of the numerical viscosity cannot be used as a controlled numerical experiment, at least as far as the occurrence and persistence of vortices is concerned. For prescribed viscosities well above the estimated numerical viscosity (α = 10⁻⁴, 10⁻³ in our case), the consistency of vortex lifetimes between the two codes and numerical choices supports the idea that the numerical experiment is indeed a controlled one.

Recent observations of molecular line broadening (e.g., Flaherty et al. 2018) and numerical studies of VSI turbulence (e.g., Flock et al. 2017) and planet–disk interaction (e.g., Zhang et al. 2018) point to low α values. The requirement of a numerical viscosity lower than the physical viscosity necessitates high resolution, which poses a challenge for simulations of vortices in protoplanetary disks.

6.6 Resolution and numerical viscosity

Vortex evolution in “inviscid” disks is often studied using very high-resolution grids to minimize the effects of numerical viscosity (Li et al. 2005; Paardekooper et al. 2010; Lin & Papaloizou 2011; Zhu & Baruteau 2016; Hammer et al. 2017, 2021; McNally et al. 2019; Fung & Ono 2021). While the resolution of 8 and 16 cells per scale height is likely enough to resolve planet-generated features such as the gap shape and spiral arms (see Appendix A), the numerical viscosity also needs to be low enough not to interfere with vortex decay.

An estimation of the numerical viscosity, valid for first-order schemes, is , with a representative cell size Δx and the time-step Δt. For our choices of parameters and assuming this corresponds to α_num ~ 10⁻²–10⁻¹. Clearly, we see substantial changes in dynamics down to much lower values of the prescribed α. Because we employ a higher-order scheme, this simple estimate is not applicable. To our knowledge, there exists no formula to estimate the numerical viscosity for the higher-order schemes employed in this study, and so we attempt to estimate it by comparing the results of our simulations at different values of α.

In general, we observe similar behavior between models with α = 10⁻⁵ and 10⁻⁶, in terms of both the behavior of vortices during their lifetime (size, mass, migration patterns) and the overall lifetime itself (see Fig. 4). This is also true across both codes that we used in this study, with the exception of the 8 cps models for β = 100. We attribute the similarity to the numerical diffusion inherent in the different advection schemes of the two codes and expect that this translates to an effective α_num between 10⁻⁶ and 10⁻⁵ for our given choices of grid resolution. This implies that our experiments with α = 10⁻⁶ are most likely not controlled ones, and for this reason, we typically group models with α ≤ 10⁻⁵ together.

Nevertheless, we still observe a different behavior for some models with α = 10⁻⁶ when comparing them to those with α = 10⁻⁵, such as the migration of the long-lived models presented in Sect. 5 (see the different tracks of r_vort(t) in Fig. 8), most of which have a 16 cps resolution. This hints at a lower numerical diffusion for 16 cps of α_num ~ 10⁻⁶. Because the numerical viscosity in the 8 cps models might interfere with the prescribed α ≤ 10⁻⁵, our 8 cps simulations might not be as trustworthy as our higher resolution 16 cps, α ≤ 10⁻⁵ runs.

6.7 Effects of the different numerics of the two codes

We used two codes (PLUTO and FargoCPT) with fundamentally different numerical properties. The fact that the two codes agree in terms of results (see the orange lines in Fig. 8 for one striking example) is reassuring, but it is worth discussing their differences nonetheless.

FargoCPT requires an artificial viscosity prescription to stabilize the upwind method near regions of strong compressionsuch as shocks. This provides additional dissipation which could affect the evolution of vortices whenever they interact with the spiral shocks induced by the planet. With the exception of the 8 cps models for β = 100, we found no significant differences in vortex lifetimes between the two codes. The one case for which the codes disagreed might be a result of insufficient resolution because the differences disappear for 16 cps.

On the other hand, PLUTO’s strictly energy-conserving nature means that the evolved quantity in the energy equation is the sum of kinetic and thermal energy. As kinetic energy dominates over thermal in typical Keplerian flows (for our setup, ), numerical errors in the calculation of total energy could affect the thermal energy budget of the disk due to subtractive cancellation error. In order to check this effect, we reran our fiducial model using the ENTROPY_SWITCH option of PLUTO, which ensures entropy conservation outside of the vicinity of shocks (which by definition do not conserve entropy, but are captured accurately by the Riemann solver). We find that this does not affect the life track of the generated vortex.

Finally, we also reran the fiducial model with PLUTO using a third-order solver and parabolic reconstruction instead of the standard second-order solver and linear reconstruction setup. We find no differences in vortex evolution or lifetime. On the basis of our tests and the agreement of the codes for high resolution, we conclude that the vortex dynamics and effects we observed in our simulations are not numerical artifacts but are indeed physical.

7 Summary

We studied vortices created by planets in protoplanetary disks using 2D viscous hydrodynamics simulations. The equation of state was assumed to follow an ideal gas, turbulence was included following the α parametrization, and thermal processes were considered by prescribing a thermal relaxation timescale using the β formalism. Our focus was on vortices exterior to the gap opened by the planet. In order to verify our results, we carried out the simulations with both the FARGO and PLUTO codes which use different numerical schemes. The planet was treated as a nonaccreting point mass with a smoothed gravitational potential and kept on a fixed circular orbit. Properties of vortices were automatically extracted using our newly developed Vortector Python tool, which identifies and characterizes vortices. Vortex identification was performed by looking for elliptical shapes in iso-vortensity lines in the r-ϕ plane, and characterization was performed by fitting a 2D Gaussian to the vortensity and surface density.

Vortices form during the gap-opening process as the embedded Jupiter-mass planet is introduced into the simulation. At the outer gap edge, multiple small vortices form that usually merge into a single large vortex that lives, depending on parameters, between 200 and several thousand orbits. These vortices have a FWHM (as determined by the fitted 2D Gaussian) of up to 0.4 r_p (several au for a planet at r_p = 5.2 au). The mass enclosed in this vortex area is up to one planetary mass (one Jupiter-mass in our models) for our choice of disk mass.

Vortex lifetime depends on the thermal relaxation timescale such that vortices live shortest for intermediate cooling times (β = 1), a result also found by Fung & Ono (2021). We find two regimes for the lifetimes of vortices. A short-lived regime with vortex lifetimes of up to 3000 orbits is observed for slowly cooling disks (β ≥ 1), in which the vortices decay faster than expected from viscous dissipation alone. In the long-lived regime, which is observed for fast cooling (β ≪ 1) with the isothermal assumption as an extreme, vortices live for a much longer time and do not decay rapidly. Vortex lifetimes are considerably longer in this regime, with a lower bound on the maximum lifetime being 15 000 orbits (the model was terminated while the vortex was still alive due to runtime constraints). From our analysis, we suspect that the long lifetime for small β is connected to the interaction of the vortex with the spiral arms, which are a source of vorticity. Details are left to future studies.

Additionally, including the disk’s self-gravity in our models with a Toomre parameter Q ≈ 25 usually shortens the lifetime of vortices and stops the small initial vortices from merging into one large vortex. Typically, in our models, two smaller vortices remain after the initial gap opening process, which then decay faster compared to those in models where disk self-gravity is not accounted for. This finding that self-gravity is detrimental to vortex survival is in line with previous studies (Lovelace & Hohlfeld 2013; Zhu & Baruteau 2016; Regály & Vorobyov 2017; Pierens & Lin 2018).

Outward migration of the vortex is observed in some of the models with β ≪ 1 and β ≫ 1. In those cases, a second density (and thus pressure) bump forms outside of the vortex location, towards which the vortex then migrates (Paardekooper et al. 2010). In some β = 0.01 models, a small, short-lived, secondary vortex forms between the planet gap and the primary vortex.

Concerning the dependence of vortex lifetime on viscosity, we find the expected behavior that this lifetime is shorter for higher viscosity (Godon & Livio 1999; De Val-Borro et al. 2007; Ataiee et al. 2013; Fu et al. 2014; Regály et al. 2017). For the highest viscosity of α = 10⁻³, practically no vortices are observed. For α = 10⁻⁵ and 10⁻⁶ we find nearly identical results, suggesting that the numerical viscosity in our models with a resolution of 8 and 16 cells per scale height is of the order of α_8cps ≲ 10⁻⁵ and α_16cps ≈ 10⁻⁶.

Allowing the planet to grow over a longer time, 1000 instead of 100 orbits, leads to longer vortex lifetimes in all the cases we tested. This disagrees with the findings of Hammer et al. (2017), who found reduced vortex lifetimes for longer planet-growth times. In our models, vortices take longer to form in the case of the more slowly growing planet. During their decay, however, their evolution is very similar, independent of planet introduction time (see Fig. 6), which in total increases their lifetime. The fact that vortex lifetime increases for longer planet-growth timescales could be an indication that the effects presented in this study, including the long-lived vortex regime, are also applicable to longer, and arguably more realistic planet-growth timescales of around 10 000 orbits.

Estimating vortex lifetime from our results, vortices are expected to live much longer at larger distances away from their host star. The increase in expected lifetime is firstly due to the longer orbital period at large radii, but also because the expected β values – the thermal relaxation timescale compared to the orbital timescale – are much lower and vortices then likely belong to the long-lived regime (see Sect. 5). From order-of-magnitude calculations, we find that large planet-induced vortices exterior to the planet at 50–100 au might live for up to several million years for low-viscosity disks (α ≲ 10⁻⁴). Considering the sensitivity of instruments like ALMA at these distances from the star, this suggests that these vortices should be observable more easily than planet-induced vortices at smaller radii.

Acknowledgements

T.R. and W.K. acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG) research group FOR 2634 “Planet Formation Witnesses and Probes: Transition Disks” under grant DU 414/22-1 and KL 650/29-1, 650/30-1. W.B., W.K., and A.Z. acknowledge support by the DFG-ANR supported GEPARD project (ANR-18-CE92-0044 DFG: KL 650/31-1). The authors acknowledge support by the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant INST 37/935-1 FUGG. Plots in this paper were made with the Python library matplotlib (Hunter 2007).

Appendix A Grid resolution and numerical convergence

A grid resolution of 8 cells per scale height (cps) is often adopted in models of planet–disk interaction in the literature. While it is widely agreed upon as sufficient, we test this statement by performing a series of test simulations with both PLUTO and FargoCPT using the same physical parameters as our locally isothermal models (α = 10⁻⁵), but using varying grid resolutions with 1, 2, 4, 8, and 16 cps in both directions (always maintaining square cells).

We find that we achieve numerical convergence on large-scale features such as the gap width and pressure bumps for a resolution of 4 cps. Convergence on more numerically sensitive features such as gap depth and vortex formation is reached for a resolution of 8 cps, with 16 cps affecting the picture relatively weakly. This was observed across both codes, with the two showing very good agreement with each other both in terms of the resolution at which different features converge and the physical properties of said features across codes.

Fig. A.1 Results for our resolution study using FargoCPT. The overall shape of the gap is resolved with around 4 cps, while it takes 8 cps to properly resolve the gap depth and the contrast of most pressure bumps far from the gap. We are interested in the region between 0.5–2.0 rp. Extending the outer boundary to r = 10 rp in the “8-8-ext” model practically made no difference. It should be noted that the 16 × 16 cps model develops some small-scale vortices in the inner disk, which causes these differences around 0.7 rp. Interestingly, a model that resolves the radial and azimuthal directions with 4 and 1 cps, respectively, captures these radial features almost as well as one with 4 cps in both directions.

Appendix B The Vortector

A major task in this study was the identification and characterization of vortices in simulation data. For this purpose, we developed a Python package, the Vortector, which automates the process for relatively generic 2D hydrodynamics planet–disk simulations.

The Vortector package allows one to visualize the vortex detection results (an example is shown in Fig. B.2), and includes information about the location, extent, and mass of a vortex along with various statistics related to the contour.

The package is publicly available on GitHub¹. We hope to make the detection and characterization of vortices in simulation data easier for other members of the community and facilitate quantitative comparison of vortices between studies by providing a common detection pipeline.

To search for possible vortex candidates, a simple search for the location of minimum vortensity is sometimes enough to find the location of a vortex. Then, the value of the vorticity ω =(∇×u) ⋅ẑ can be used to learn how strongly the vortex rotates and the local surface density can be used as an indication for the mass enclosed in the vortex. However, this method fails for many simulations, e.g., when the vortensity in a very small region close to a spiral arm of the planet is lower than inside a vortex candidate, or when the gap region intrudes into the outer disk, which can induce strong anticyclonic motion at the outer gap edge.

To get around these issues, the Vortector uses the geometrical shape of vortices as they appear in a face-on image of the disk. Looking down on the surface of a disk, vortices appear as crescent-shaped objects. In the r-ϕ plane, which is more suitable for this task, large vortices appear as elliptical objects (see also Fig. 1 of Lesur & Papaloizou 2009). In fact, contour lines of the vortensity closely resemble ellipses in the r-ϕ plane. We can therefore identify vortices in a disk by finding closed contour lines that closely resemble ellipses. To solve this task programmatically, we can make use of the computer vision library OpenCV (Bradski 2000).

Our strategy to extract vortex candidates from a simulation snapshot can be then subdivided into three tasks:

• 1.

  Extract contour lines in the r-ϕ map of the vortensity,

• 2.

  identify nearly elliptical contours as vortex candidates, and

• 3.

  fit 2D Gaussians to ϖ and Σ for characterization.

The algorithm step by step

This section describes the vortex detection process using the model presented in Sect. 5 (FargoCPT, α =10⁻⁵, β =0.01, 16 cps) which shows the emergence of a secondary vortex. The data used for this analysis corresponds to a time t = 7150 orbits.

Before the analysis is performed on ϖ, the map is periodically extended in the ϕ direction in order to be able to identify vortices that intersect the periodic azimuthal boundary. The resulting image is shown in the left panel of Fig. B.1. There, the original domain is indicated by the thick solid rectangle, which spans the azimuthal range from 0 to N, where N is the image size in pixels in the azimuthal direction. The top and bottom areas of the domain (orange and green) are repeated at the lower and upper boundaries, respectively.

Fig. B.1 Periodically continued iso-vortensity line image (left) used for extracting contours for vortex candidates and an example contour illustrating the ellipse fit (right). The snapshot shown is at time t = 7150 orbits of the model with the secondary vortex that was discussed in Sect. 5. The left panel shows how the data array is mirrored in order to allow the detection of vortices that overlap with the periodic boundary. Areas with the same color are copies of one another. The red line indicates the outline of the grey area in the right panel. The original size is marked by the black rectangle ranging from 0 to N on the vertical axis. The area shaded in blue in the right panel illustrates the definition of the deviation from the ellipse that is used to select the vortex candidates from the closed contours. For the example shown, the ratio of difference area (blue) to the total contour area is 0.122, which is below the 0.15 threshold.

Task 1: Contour lines

Contour lines are extracted for a range of ϖ values ranging from 0 to 1 in increments of 0.05. For each value ϖ_crit, a binary image is produced by setting each cell with ϖ ≥ ϖ_crit to 1 and 0 otherwise. The binary image is then analyzed using findContours from OpenCV. Only closed contours are retained. This step usually results in up to a few thousand contours, depending on the dynamical state of the disk and the choice of increments in ϖ.

Task 2: Find closed contours resembling ellipses

Next, the fitEllipse function from OpenCV is used to fit an ellipse to each closed contour. One example of this is shown in the right panel of Fig. B.1, where the ellipse is visible as an orange line in the zoom-in.

The difference in area between the contour and fit is used as a measure of deviation. The deviation from an ellipse is then defined as the ratio of this difference and the area enclosed by the contour. We only keep contours for which the deviation is smaller than 0.15. The example contour in Fig. B.1 has a deviation of 0.122.

Finally, all the contours that are contained within the largest contour that satisfies this criterion are discarded, which leaves the example red contour in Fig. B.1 as the selected vortex candidate (see also the white contour line in Fig. B.2).

We only retain contours that enclose at least two other contours. With this restriction, we make sure that ϖ∕ϖ₀ changes by a value of 0.1 from the outside to the inside of the vortex candidate. This has proven to be useful to filter out small fluctuations in the disk that otherwise appear as small transient vortices.

At this point, it becomes clear that the extent of the vortex and derived quantities such as the mass contained within are influenced by the choice of the levels used to produce the contour lines and the choice of the maximum relative ellipse deviation. The properties of the contour give an order-of-magnitude estimate nonetheless.

Task 3: Fit a 2D Gaussian

To remove the influence of the threshold parameters in the detection of the contour, a process that does not depend on our parameter choices but on the underlying data is needed.

Upon inspection of the curves of vortensity and surface density along a cut through the vortex, either radial or azimuthal, it becomes clear that these lines resemble Gaussian functions (see curves in Fig. B.2 around the 2D maps) (B.1)

Here, σ_r and σ_ϕ provide a measure for the vortex size and can even be used to give a definition of the vortex region that does not depend on additional parameters. In combination with the center coordinates r₀ and ϕ₀, σ_r and σ_ϕ can be used to define the vortex as the disk material contained within the ellipse given by (B.2)

where and denote the half width at half maximum of the 2D Gaussian function defined in Eq. (B.1). We usually use the values obtained from the surface density fit because these are less time-sensitive compared to the vortensity fit and because the shape of Σ curves more closely resemble Gaussians (see Fig. B.2).

Fig. B.2 Overview of the results produced by the Vortector package for a model showing a secondary vortex discussed in Sect. 5 at t = 7150 orbits, with a 2D map of the vortensity on the left and surface density on the right. All detected vortex candidates are indicated in the 2D plots. The extracted contour (shown in Fig. B.1) is marked with a white line, the ellipse of the vortensity fit is shown in blue and the ellipse of the surface density fit is shown in green. Note that these ellipses are defined by σr and σϕ from the fit of Eq. (B.2) and are different from the ellipse used to fit the contour. The ellipses of the most massive vortex include a crosshair indicating the center of the fit. Each 2D plot is accompanied by 1D plots of slices through the main vortex. The plots also show the values of the respective Gaussian fit in blue for vortensity and green for surface density. In this figure, the planet is located at r = 1 and ϕ = 0.

Appendix C Data table

The lifetimes and parameters of all models mentioned in Sects. 4 and 5 are listed in Table C.1.

References

All Tables

All Figures

Fig. 1 Multiple snapshots of the α = 10−5, β = 1, 8 cps model showcasing the vortex merging process during the early stage of gap opening, the fully grown size of the resulting vortex, and its subsequent decay. The surface density and vortensity contrast compared to their initial profiles is shown in the left and right panels, respectively. Time is quoted in units of planetary orbits. The horizontal line at r = 1.45 rp serves to highlight the outward radial movement of those structures as the gap around the planet grows wider. The planet is located at r = 1 rp and ϕ = 0.

In the text

Fig. 2 Evolution of radial Lovelace parameter (see Eq. (8)) and Σ profiles during vortex formation over the first 200 orbits of the sample case from Sect. 3. The vertical lines indicate the center of the plateau in (estimated by eye) to guide the eye to the corresponding location of the Σ profile. is calculated as the azimuthal average at each radius. The dotted horizontal line in the bottom panel marks 10% of Σ0, which we define as the location of the gap edge.

In the text

Fig. 3 Evolution of vortex properties for the showcase simulations presented in Sect. 3. The panels show, from top to bottom, the mass enclosed in the FWHM ellipse of the vortex fit Mvort in Jupiter masses, the radial location of the vortex rvort, and its FWHM Δr, indicated by the shaded area, and the ratio between minimum vortensity inside the vortex and the azimuthal median of vortensity at the radial location of the vortensity minimum. A dotted vertical line indicates the time when the planet has reached its final mass. The curves are smoothed with a median filter that spans over the preceding and following five datapoints (± 50 orbits at rp). The orange parts of the line in the bottom panel show the evolution of the vortensity prior to the “birth” and after the “death” of the vortex.

In the text

Fig. 4 Lifetime of vortices as a function of β for 8 cps (left) and 16 cps resolution (right). Colors encode α, and the different symbols denote the code and the inclusion of self-gravity. The solid lines help guide the eye and connect the lifetime averages between the two codes (without self-gravity) for each value of α, where thetwo codes agree sufficiently. For parameters where there is a difference between the codes, dashed and dotted lines connect to the datapoints of the FargoCPT and PLUTO runs, respectively. A “↦ ” next to a symbol marks models that were terminated due to runtime constraints but still contain an active vortex. The horizontal gray line in the right panel indicates the top of the y-axis of the left panel. A list of all vortex lifetimes shown here is provided in Table C.1.

In the text

	Fig. 5 Evolution of vortex properties for varying values of the thermal relaxation β parameter. The panels are as in Fig. 3. Shown are models run with the FargoCPT code with α = 10−5 and at 8 cps resolution (orange “f8” dots in Fig. 4).
In the text

Fig. 6 Influence of the planet introduction time on the evolution of vortex properties. The panels are as in Fig. 3. Solid and dashed lines show models with τramp = 100 orbits and 1000 orbits, respectively. The τramp = 100 orbits curves are shifted to the right (see the horizontal lines) to illustrate that the curves have the same shape in the decay phase, independent of τramp. We note that the final evolution of the vortex, after it has reached its minimum in vortensity, is the same independent of planet introduction time.

In the text

	Fig. 7 Evolution of vortex properties for varying values of α. Panels are shown as in Fig. 3. Shown are models run with the FargoCPT code with β = 1 and at 8 cps resolution. The α = 10−3 run is excluded because no vortex forms. In addition, a run with disk self-gravity enabled is added for the α = 10−5 case. The similarity between simulations with α = 10−5 and 10−6 is apparent.
In the text

	Fig. 8 Selection of models with long-living and migrating vortices at 16 cps resolution for the two different codes. Both codes agree remarkably well for the blue and orange cases. The panels are as in Fig. 3. In models shown here, β = 0.01. The values of α, the resolution, and the code used are indicated in the legend.
In the text

Fig. 9 Azimuthally averaged surface density profiles as a function of different physical (α, β) and numerical (cps) parameters at two different time-stamps. The peak around r∕rp = 1.5–1.6 corresponds to the pressure bump formed by the planet as the latter pushes material away, forming a gap around its orbit. The smaller, secondary peak at around r∕rp = 2.1 is caused by the vortex that forms near the primary, planet-generated bump. Top: radial profiles at t = 1000 orbits. At this stage, all models pictured feature a vortex near the primary bump. We note the absence of a secondary bump for the models with β = 1. Bottom: same profiles at t = 5000 orbits. Here, the primary bump has moved radially outwards as the planet’s gap gets deeper and wider. We highlight the depletion of gas near the primary pressure bump for the second panel from the left (β = 0.01). This is caused by the combination of a vortex migrating outwards to the secondary bump, and the inability of the planet to resupply that zone with material from its now-depleted gap region. We also note the difference between resolutions of 8 and 16 cps (dashed and solid lines), especially for the β = 100 models and the β = 0.01, α = 10−5 model (the evolution of this model is shown as the orange line in Fig. 8).

In the text

Fig. 10 Baroclinic term (RHS of Eq. (10)) in the outer disk for a short-lived model (β = 1 at t = 1000 orbits) and a long-lived vortex (β = 0.01 at t = 7150 orbits, see also Fig. B.2) with α = 10−5 and 16 cps resolution on the left and right side, respectively. The top row shows maps of the baroclinic term with the detected vortices indicated with green ellipses as obtained from the Σ fit. The bottom row shows the radial Σ profile in orange, the azimuthally averaged baroclinic term in blue and the region between its minimum and maximum shaded in gray.

In the text

Fig. A.1 Results for our resolution study using FargoCPT. The overall shape of the gap is resolved with around 4 cps, while it takes 8 cps to properly resolve the gap depth and the contrast of most pressure bumps far from the gap. We are interested in the region between 0.5–2.0 rp. Extending the outer boundary to r = 10 rp in the “8-8-ext” model practically made no difference. It should be noted that the 16 × 16 cps model develops some small-scale vortices in the inner disk, which causes these differences around 0.7 rp. Interestingly, a model that resolves the radial and azimuthal directions with 4 and 1 cps, respectively, captures these radial features almost as well as one with 4 cps in both directions.

In the text

Fig. B.1 Periodically continued iso-vortensity line image (left) used for extracting contours for vortex candidates and an example contour illustrating the ellipse fit (right). The snapshot shown is at time t = 7150 orbits of the model with the secondary vortex that was discussed in Sect. 5. The left panel shows how the data array is mirrored in order to allow the detection of vortices that overlap with the periodic boundary. Areas with the same color are copies of one another. The red line indicates the outline of the grey area in the right panel. The original size is marked by the black rectangle ranging from 0 to N on the vertical axis. The area shaded in blue in the right panel illustrates the definition of the deviation from the ellipse that is used to select the vortex candidates from the closed contours. For the example shown, the ratio of difference area (blue) to the total contour area is 0.122, which is below the 0.15 threshold.

In the text

Fig. B.2 Overview of the results produced by the Vortector package for a model showing a secondary vortex discussed in Sect. 5 at t = 7150 orbits, with a 2D map of the vortensity on the left and surface density on the right. All detected vortex candidates are indicated in the 2D plots. The extracted contour (shown in Fig. B.1) is marked with a white line, the ellipse of the vortensity fit is shown in blue and the ellipse of the surface density fit is shown in green. Note that these ellipses are defined by σr and σϕ from the fit of Eq. (B.2) and are different from the ellipse used to fit the contour. The ellipses of the most massive vortex include a crosshair indicating the center of the fit. Each 2D plot is accompanied by 1D plots of slices through the main vortex. The plots also show the values of the respective Gaussian fit in blue for vortensity and green for surface density. In this figure, the planet is located at r = 1 and ϕ = 0.

In the text