Vortex Formation and Survival in Protoplanetary Disks subject to Vertical Shear Instability

Several protoplanetary disks observed by ALMA show dust concentrations consistent with particle trapping in giant vortices. The formation and survival of vortices is of major importance for planet formation, because vortices act as particle traps and are therefore preferred locations of planetesimal formation. Recent studies showed that the vertical shear instability (VSI) is capable of generating turbulence and small vortices in protoplanetary disks that have the proper radial and vertical stratification and thermally relax on sufficiently short time scales. But the effect of the azimuthal extend of the disk is often neglected as the disks azimuth is limited to $\Delta \phi \leq \pi/2$. We aim to investigate the influence of the azimuthal extent of the disk on the long-term evolution of a protoplanetary disk and the possibility of large vortices forming. To this end, we perform 3-dimensional simulations for up to 1000 local orbits using different values of $\Delta \phi = \pi/2 $ to $2\pi$ for VSI in disks with a prescribed radial density and temperature gradient cooling on short timescales. We find the VSI capable of forming large vortices which can exist at least several hundred orbits in simulations covering a disk with $\Delta \phi \geq \pi$. This suggests the VSI to be capable to form vortices or at least to trigger vortex formation via a secondary instability, e.g. Rossby Wave Instability or Kelvin Helmholtz Instability.


INTRODUCTION
Turbulence in disks around young stars is still one of the most interesting questions in modern astrophysics (Turner et al., 2014). Balbus & Hawley (1991) introduced the magneto-rotational instability (MRI) as a promising source of turbulence with an alpha viscosity (Shakura & Sunyaev, 1973) large enough to explain angular momentum transport on timescales set by observations. But more recent work shows the MRI to be hampered by non-ideal magnetic effects such as resistivity or ambipolar diffusion (Lesur, Kunz & Fromang, 2014). They show that the instability can be damped efficiently in parts of the disk by low ionization fractions, where then other sources of turbulence can and have to be considered (Lyra & Klahr, 2011).
Several possible mechanisms for pure hydrodynamic turbulence have since been proposed, acknowledging the fact that a Keplerian flow with radially increasing angular momentum profile is hydrodynamic stable as can be seen in a most general way from the Solberg-Hoiland criteria, which is derived for no thermal relaxation (Rüdiger, Arlt & Shaly-E-mail:manger@mpia.de bkov, 2002). The global baroclinic instability (aka Subcritical Baroclinic Instability) (Klahr & Bodenheimer, 2003; and convective overstability (COS) (Klahr & Hubbard, 2014;Lyra, 2014) rely on radial temperature and density stratifications introducing convective cells in disks with moderate cooling times on the order of τ = 1/(γΩ). The Vertical Shear Instability (VSI) (Nelson, Gressel & Umurhan, 2013;Richard, Nelson & Umurhan, 2016;Stoll & Kley, 2014) relies on short cooling times to remove the stable vertical stratification to tap into the energy of the vertical shear.
In all of those turbulence models, the formation of vortices has been observed. Raettig, Lyra & Klahr (2013) showed them in local simulations of subcritical baroclinic instability, (Flock et al., 2015) in a global context at the edge of the MRI dead zones. The Rossby Wave Instability (Papaloizou & Pringle, 1984, 1985Lovelace et al., 1999;Li et al., 2000Li et al., , 2001 has been shown to break axisymmetric rings into large vortices. Richard, Nelson & Umurhan (2016) showed the possibility of small vortex formation within disks susceptible to VSI and Latter & Papaloizou (2018) found the Kelvin Helmholtz Instability acting as a parasitic instability on the VSI modes and forming small vortices in the disk.
Because vortices are naturally identified with local pressure maxima in the context of PPDs, they act as particle traps (Barge & Sommeria, 1995). They are found to enhance the dust to gas ratio locally, aiding planetesimal formation via triggering the streaming instability (Raettig, Klahr & Lyra, 2015) and accelerate core growth for giant planets via pebble accretion (Klahr & Bodenheimer, 2006). Large vortices have also been discussed as explanation for features found in (sub-)mm observations of protoplanetary disks with ALMA (van der Marel et al., 2013) and VLA (Carrasco-González et al., 2016).
In this work we investigate vortex formation and survival in protoplanetary disk simulations undergoing Vertical Shear Instability, reexamining the work of Richard, Nelson & Umurhan (2016). We find that once we loosened the restriction to small azimuthal domains (large azimuthal wavenumbers m) and allowed for simulation domains of ϕ = 180 • and ϕ = 360 • , we find large vortices forming in the disks after a few hundred orbits. We also find these larger vortices to survive for hundreds of local orbits, making them excellent particle traps and candidates for planetesimal formation sites. With this, we stress again that non-axisymmetry plays an important role in assessing disk turbulence features (Klahr, Henning & Kley, 1999).
In section 2, we will shortly revisit the theoretical background to our simulations. Section 3 lays out the numerical set-up used in all computations. In Section 4 we present our results and in section 5, we discuss them in context of recent literature. Finally, section 6 summarizes our findings and presents an outlook on future work.

THEORETICAL BACKGROUND
In this work, three instabilities have to be considered. We implement our simulation with a steady-state model susceptible to the growth of the Vertical Shear Instability, which we identify for this paper as the main driver of turbulence in the short cooling time regime. The Rossby-Wave-Instability is important in the subsequent formation of giant vortices in the disk. The elliptic instability then has to be considered acting inside the formed vortices as a driver of turbulent vortex substructures and possible vortex destruction.

Vertical Shear Instability (VSI)
Protoplanetary disks with radial gradients in temperature and entropy do not rotate on cylinders but have an angular frequency depending on height Ω = Ω(R, z) and therefore exhibit vertical shear and a non-zero vertical epicyclic frequency, i.e. κ 2 z = ∂z R 2 Ω 2 R 4 = 0. The stability of these disks against axisymmetric, adiabatic perturbations can be in general described through the Solberg-Hoiland criteria (Rüdiger, Arlt & Shalybkov, 2002), which in general predict stability due to the vertical stratification of the disk. For nearly locally isothermal disks, i.e. disks with a short cooling time, or vertically adiabatic disks (Pfeil & Klahr, subm.) however, the stabilizing effect of vertical buoyancy against axisymmetric perturbations is diminished. This leads to the Vertical Shear Instability, first described for radiative zones of differentially rotating stars by Goldreich & Schubert (1967) and Fricke (1968). The instability has been more recently suggested to also operate in protoplanetary disks by Urpin (2003) and shown by Nelson, Gressel & Umurhan (2013). In the limit of thin, locally isothermal disks ((H/R) 1), instability is determined by the criterion: with j the specific angular momentum of the disk and kR and kZ the radial and vertical wavenumber. Instability then occurs for nearly vertical modes with kR/kZ > R/H tapping into the free energy provided by the vertical shear (Nelson, Gressel & Umurhan, 2013). Nelson, Gressel & Umurhan (2013); Stoll & Kley (2014, 2016 and Flock et al. (2017) showed turbulent stresses for VSI active disks to be in the range α ≈ 5 · 10 −5 − 5 · 10 −4 . Richard, Nelson & Umurhan (2016) reported small vortex formation to be possible for certain disk parameters.

Rossby Wave Instability (RWI)
The Rossby Wave Instability, first described in Lovelace et al. (1999); Li et al. (2000) and Li et al. (2001), is a global instability occuring in rotating flows that exhibit a local extremum in pressure or vortensity. Such an extremum can be archived e.g. at the inner edge of a dead zone Lyra & Mac Low (2012). The criterion for instability is an extremum in the function: with Σ denoting the column density, ωz = (∇ × v)Z the zcomponent of vorticity, P the vertically integrated pressure and γ the adiabatic index. In a region around the extremum, Rossby waves are trapped in a standing pattern amplifying with time, which can be observed as a vortex at the extremum. Outside, the Rossby waves emit density waves due to gradients in vorticity (Meheut et al., 2010).

Elliptic instability
The elliptic instability is a linear parametric instability. It destabilizes elliptic streamlines once the vortex turnover time matches one or more inertial frequencies of the underlying flow field by creating a positive resonance. Lesur & Papaloizou (2009) studied the instability in the context of accretion disks, showing the instability to be most effective for vortex aspect rations χ=semi-mayor axis/semiminor axis to be in the range of 1 < χ < 4 for 3 dimensional disks. A weaker second range for the instability occurs for χ > 6, the region in between is stable for disks without vertical stratification. When vertical stratification is introduced, vortices become unstable for all χ, the growth rate however depends on the strength of the stratification for vortices with aspect ratios χ > 4. Vortices with smaller aspect ratios are always highly unstable independent of the strength of the stratification.

SIMULATION SETUP
We conduct 3 dimensional simulations using the magnetohydrodynamics code PLUTO 1 . In this work, we use 2 coordinate systems. The simulations are carried out on a spherical grid (r, θ, ϕ), the model setup and analysis are presented in cylindrical coordinates (R, φ, Z), We implement the disk setup following Nelson, Gressel & Umurhan (2013). The gas density is defined by with the disk scale height H = cs Ω K . We chose a ideal equation of state ρe = P γ−1 with the specific internal energy e,the adiabatic index γ = 5/3 2 and P = c 2 s ρ using a radially changing isothermal sound speed c 2 s = c 2 0 R R 0 q . Note that the temperature is related to the isothermal sound speed via c 2 s = kT /µmH with k denoting the Boltzmann constant, µ the mean molecular weight and mH the atomic mass of hydrogen, thus defining a radial temperature gradient in the disk. For all our simulations, we choose − 2 3 and −1 for p and q respectively, satisfying the requirements set for the VSI.
The q = −1 is close to the maximum of q = −1.5 that one can expect in a viscously heated disk (see Eq. 11 in Bell et al. (1997)) in regions which are dominated by icy grains. The value for p = − 2 3 is chosen to be consistent to models in which we try to simulate the COS, by using an identical setup as described here, but using a longer cooling time. This particular value of p has no significant impact on the VSI, but helps to maximise the radial buoyancy for the given p expressed in the radial Brunt-Vaissala frequency N 2 r (Klahr & Hubbard, 2014).
The geometrical scale height is then constant throughout the disk with H/R = c0/v Kepler = 0.1, which is also nice to evaluate the simulation. H/R might actually be smaller in a real circumstellar disk (Bell et al., 1997) and we tested a value of H/R = 0.05, which was also the value choice in the work by Stoll & Kley (2014), yet (Flock et al., 2017) again uses values of H/R = 0.1 for the outer disk. H/R = 0.05 simulations are computationally more expensive, as the pressure scale height has to be resolved by at least as many cells as in the H/R = 0.1 case, leading to 8 times more cells and reducing the time step by a factor of two. Once our H/R = 0.05 simulations are finished, we will publish them in comparison to the H/R = 0.1 cases, but we claim that smaller values of H/R will show similar results for runs with comparable resolution per scale height.
The scale height can be expressed as: and the initial angular velocity of the disk is given by whereas the radial and vertical velocities are set to zero. All velocity components are initially seeded with a white noise perturbation of 10 −4 cs.
The cooling in our model is described by where Tinit(R) is the initial Temperature profile and τ relax is the relaxation time scale. It is set to τ relax = dt and is about 2 · 10 −3 /Ω0 in all simulations. This is the shortest cooling time we can realise in our explicit cooling scheme, effectively leading to an almost locally isothermal disk.
In our numerical computations we use the hllc method (Toro, 2009) with a ppm reconstruction scheme (Mignone, 2014) for spatial integration and a 3rd order Runge-Kutta method for time integration. The mesh extends from (0.5r0, π/2 − 0.35, 0) to (2r0, π/2 + 0.35, φmax) with r0 being the radius at r = 1 and φmax given in table 1. We use a logarithmic grid in radial direction and a uniform grid otherwise to preserve the aspect ratio of the individual grid cells. This gives us a resolution of 18 cells per H in vertical and radial direction and 12 cells per H in azimuthal direction. For the largest run (p360), we used about 1.4 million cpu-hours.
We employ outflow boundary conditions in radial direction, reflective boundaries in vertical direction and periodic boundaries in azimuthal direction. To minimise mass loss in radial direction and generally wave reflection at the boundaries, we add damping layers at the inner and outer boundaries in the radial and polar direction with ∆R = 1.0 H and ∆θ = 0.05. Inside the damping layers, we damp the velocities to their initial values if the velocity component normal to the boundary points inside the domain. For the damping we use:   The alpha values are averaged over the analysis domain an as running average in time. We observe the turbulence to saturate rapidly within around 100 orbits to alpha values of around 10 −3 .
The alpha values for the simulation with smaller azimuthal domain are slightly larger than for the larger domain. The inset highlights the evolution for the last 200 orbits.
all simulations enables simulation times close to 1000 orbits (which for H/R = 0.05 needs 16 times the computation cost) enables us to determine the lifetime of the formed vortices and possible influences on disk evolution an planetesimal formation. The simulation parameters are summarized in Table  1. In this section, we present the results of our simulations, focusing first on the flow properties and later on vorticity.

Transport properties
We analyse the transport properties of the VSI for all runs in a subdomain of the simulation grid. The subdomain is . This is done to avoid possible influence of the imposed boundary conditions on the results. We look at the Reynolds stresses generated in the disk following the prescription by Klahr & Bodenheimer (2003) and compute the α-parameter of the disk with P = ρc 2 s φ,t . This guarantees a mass weighted α and filters out angular momentum flux associated with mean mass transport ρvr . We plot the evolution of α over time in figure 1. The values are averaged over the whole analysis domain in each direction and up to the given point in time. We find α-values of approximately 10 −3 in agreement with Stoll & Kley (2014) and (Nelson, Gressel & Umurhan, 2013), but significantly higher than those reported by Richard, Nelson & Umurhan (2016). The equilibrium alpha values (table 1) can be seen to decrease with increasing azimuthal extent. This effect is caused by the use of periodic boundary conditions in phi direction, leading to a large pitch angle for the tightly wound spiral pattern for smaller phi extents.
To look at the vertical structure of the disk, figure 2   compares the meridional α profiles of the simulated disks averaged between 500 and 700 orbits and over the radial and azimuthal subdomain. We observe different profiles for each simulation. We find the steepest profile for the simulation p45, showing high alpha values at the disk surface and negative values at the midplane, which suggest moderate inward angular momentum transport in the midplane of the disk and outward transport in the surface layer. The profile flattens for the case p90, with positive alpha values throughout the disk and lower alpha values at the surface compared to run p45. This trend continues for the cases p180 and p360, though the latter show to be in good agreement with each other. This also suggest angular momentum is transported outward at all heights, albeit with stronger transport in the upper layers of the disk. The vertical profile of the Reynolds-stresses of the disk shows the largest angular momentum transport to occur at z ≈ ±1.6H for simulation p45 and decreasing to z ≈ ±1.3H for the full 2π disk.
To investigate this behaviour further, we calculate the radial mass flow of the disk ρ vr as a function of height above the midplane in all our simulations, shown in figure 3. We again average over the radial and azimuthal domain of the simulation and over 200 snapshots taken between 500 and 700 local orbits. We find radial mass inflow for all models in the midplane, which aligns with the findings of Stoll, Kley & Picogna (2017), who find the same flow reversal applying an anisotropic viscosity model with a heightened z-viscosity component. Therefore, although there is only small to no outward radial angular momentum transport present in the midplane from the VSI, mass can be accreted efficiently. The angular momentum however is transported vertically from the midplane to the upper layers of the disk, where it is then transported radially outward. This mechanism also gives an explanation to the observed shallower vertical profiles for the models p180 and p360. We argue that due to the large azimuthal extent of the disk, the anisotropy manifests in different magnitudes and therefore leads to an overall shallower profile. This supports our view that treating the disk as too high m leads to incomplete results.
To prove that the vertically averaged angular momentum transport (i.e. our measured mean α-value) is nevertheless sufficient to prescribe the mean radial mass-accretion, we calculate the average radial velocity given from steadystate viscous accretion theory and compare this to the simulation values. Integrating the vertically-averaged steadystate angular momentum equation one obtains with Σ the disk column density, ν = αH 2 Ω the viscosity and dΩ/dR = −3/2Ω. For the case p360 with α = 10 −3 we get ΣvR/(Σv k ) = −1.5 · 10 −5 . From the corresponding simulation run we get from integrating over figure 3 an average value of ΣvR/(Σv k ) = −2.4 · 10 −5 , which agrees well to the predicted value. This shows that angular momentum transport driven by VSI is well described in the picture of α-viscosity, even if transport is most likely realised by travelling spiral waves rather than local Kolmogorov-like 3D turbulence. This also supports our claim that angular momentum is transported mainly vertically from the midplane to higher layers and then transported outward, as the vertically averaged angular momentum matches the value needed for the occurring mass transport despite the low α values measured in the midplane of the disk. Whether the turbulence is truly local, i.e. the dissipation of kinetic occurs also proportional to the measured α-viscosity, was not possible to be determined in our simulations, but should be goal for future setups. This would need the monitoring of local heating and cooling.
We also look at the time evolution of the rms-velocities in the disk, defined as We consider only the radial and vertical component of velocity because the disk rotation profile is not keplerian in this work and spatially varying systematic deviations cannot be taken into account. The results are presented in figure 4. We find the values of all runs to agree with each other for the first few orbits as expected, as the limit in azimuthal  wavenumber does not influence the onset of the instability. The overall onset of the instability is observed earlier than in other studies. This is explained by the specific choice of parameters in our setup, which allows for earlier onset of the instability due to the very short cooling time. The saturated values for cases missing the lower azimuthal wave-numbers p45 and p90 are higher than for the other cases. We measure the growth rate of the velocity perturbations as 0.4 per orbit, about double the value reported by Stoll & Kley (2014), which is in excellent agreement with theory as growth rate is proportional to the pressure scale height σ ∼ qΩ H R (Nelson, Gressel & Umurhan, 2013), which was only half the value we adopted. Note that Stoll & Kley (2014) give their value for the growth rates for the kinetic energy. Comparing the time evolution of the rms-velocity to the time evolution of the alpha value in figure 1, we find that the rms-velocity saturates after about 20 orbital periods, while α, the quantity related to angular momentum transport, saturates only after 100 orbital period. This behaviour is linked to the transport of energy to larger scales in the disk and we will discuss this further in section 5.1. Figure 5 shows the vrms-values as a function of height averaged over the radial and azimuthal subdomain and between 500 and 700 orbits. We find a similar shape for the vertical structure of all runs with low velocities in the midplane and rising with z/H. We find a systematic positive offset for all runs compared to p360 with p45 having the largest offset. This compares to our findings for the α viscosity parameter in figure 2, suggesting a smaller φmax systematically overestimates the turbulence strength.
In spectral line observations of protoplanetary disks the total rms velocity can however not be measured. The most easily accessible quantity is the vertical component of the velocity, which can be measured by determining the line broadening of face-on disks. We separately calculate the height profile for this quantity in figure 6. We again see the large offset for the simulations p45 and p90 compared to the case p360, the case p180 shows similar values to the 360 • case.  Figure 5. Meridional profile of the vrms value. The vrms values are averaged over the radial and azimuthal subdomain and between 500 and 700 orbits. We again observe significantly lower values for the run p360 than for the other runs.  Figure 6. Vertical profile of the vertical rms velocity. We average the values over the radial and azimuthal subdomain and between 500 and 700 orbits. We find a large offset for the cases p45 and p90 compared to p180 and p360, which are in good agreement with each other.
The height profile for all simulations follows a similar shape as the total vrms values with low values in the midplane growing with height z. Cuzzi et al. (2001) related the rms-velocity to the turbulent viscosity parameter α via (their equation 2) if the largest eddies have a rotation frequency comparable to the orbital frequency. To check the applicability of this relation to the turbulence induced by the VSI, we plot the ratio of the height dependent z component of the rms-velocity to the square-root of the simulation averaged alpha value (see table 1) in figure 7. For all simulations, the ratio of vrms/cs to √ α is above unity. Therefore the angular momentum transport in our simulations of the VSI is weaker than one would expect from the measured velocities. The deviation depends on height with values closer to  1). We average the values over the radial and azimuthal subdomain and between 500 and 700 orbits. All simulations show an offset from unity, so estimating the α-values from vz,rms leads to an overestimation. We again find a large positive offset for the cases p45 and p90 compared to p360, whereas p180 has a smaller negative offset with respect to p360.
unity in the midplane. The overall deviation is largest again for the case p45, which also showed the highest values for vrms/cs and α and decreases with φmax. For the case of a full disk we find however values larger than for φmax = 180 • . For protoplanetary disks subject to the VSI, the α values calculated from measured turbulent velocities should therefore be treated with caution. Depending on the height above the midplane where the measurement is taken, the values for α could be overestimated by up to one order of magnitude. Also, when comparing the alpha values of simulations and observations, the influence of the azimuthal extent of the simulation should be taken into account.

The influence of φmax on the disk structure
To asses the influence of φmax on the outcome of the numerical simulation, we plot the z component of the vorticity in the disk midplane for three different times: after 200, 500 and 700 local orbits. The results are presented in figure 8. For the case of p45, we find small vortices with aspect ratios χ ≈ 4 and zonal flows. All structures show strong variations in time, as can be seen by the differences in the time frames shown in figure 8. These structures also appear in the case of p90, where also larger structures similar to vortices emerge in the first frame but are destroyed again in the second frame of figure 8, also pointing to high variability with time. This changes for the simulations p180 and p360. Here we observe small, unstable vortices forming quickly in the beginning and additionally two larger vortices with aspect ratios χ ≈ 8 − 10 after a few hundred orbits, seen in the left frame in figure 8. The large vortices continue to appear both after 500 and 700 orbits ( fig. 8, middle and right columns), suggesting stability over larger times. The middle and right column also show additional large vortices appearing at later times. This In contrast to this we see large vortices in the lower rows for the second and third snapshot.
change in overall structure for azimuthal extents larger than 180 • suggests that the approximation of large m restricts the development of large, long lived structures in VSI disks.
Looking at the large vortices specifically, we find them to have an inner turbulent structure, similar to the one found in (Raettig, Lyra & Klahr, 2013;Lyra, 2014). Figure 9 illustrates this, showing the structures inside the outer vortex of the simulation p360 at (Rc, φc) = (1.4 R0, 3.5) after 500 orbits ( fig. 8, middle column). We again plot the z component of vorticity, now for 4 different heights z = 0, 0.5 H, 1 H and 2 H above the midplane of the disk. We multiply the azimuth with the central radius and scale both axis with the local pressure scale height to show the true size of the vortex.
The structure inside the vortex shows changes on small scales, forming a highly turbulent substructure. The turbulent structure is visible throughout the upper and the lower left panel of figure 9. In these 3 panels we also find a consistent outer boundary for the vortex as an azimuthally elongated radially narrow region of higher vorticity relative to the inside of the vortex, forming a shell around the vortex. This structure is consistent with the predictions of Lesur & Papaloizou (2009) for the elliptic instability, which we will discuss in section 5.2. The lower right panel of figure 9 shows the vortex structure fading into the turbulent back-ground structure of the disk, suggesting the vortex extends to about 2 scale heights above and below the midplane.
To asses the influence of the forming vortices on the flow structure reported for the VSI we plot the vertical velocity times density in the R-Z plane in figure 10 after 500 orbits, corresponding to the middle column of figure 8. Figure 10 (b) and (c) show the models p45 and p90 at ϕ = π/8 and ϕ = π/4 respectively and (d) shows the model p180 at ϕ = 2.5, the center position of both vortices at that time.  Figure 10 (a) shows the 2D axisymmetric simulation with the same r,z resolution from the resolution study shown in appendix A after 80 orbits. The snapshot from the 2D run is taken at an earlier time to avoid the non-linear phase of the axisymmetric VSI, which differs from the full 3D simulations due to the occurrence of the RWI in 3D.
We find a similar vertical velocity profile as reported by Nelson, Gressel & Umurhan (2013), Stoll & Kley (2014) and Flock et al. (2017) with strong vertical motions present over the whole height of the disk. We find the run p45 agreeing well with the early stage of the 2D comparison simulation, showing a pattern of disk annuli of ∆R ≈ 0.5H with alternating positive and negative vertical velocity and symmetry about the midplane of the disk. In (c) we see the overall pattern still preserved for the case p90, but the annuli now have a larger radial width of ∆R ≈ H. This picture changes once the disk is able to form vortices. Fig. (d) shows the case of p180 at the centre position of the inner vortex. The ordered pattern seen in (a)-(c) is broken at the radii of the vortices, although it is still partly visible at other radii. The size of the annuli of alternating velocity also decreases in this simulation. Both vortices have negative vertical velocities, although the velocity is greatly reduced with respect to the surrounding, especially for the outer vortex. Figure  (f) shows a similar behaviour for the case p360. The overall structure of the disk is less ordered than for the axisymmetric case and the vortex has a lower vertical velocity as the surrounding disk, although a similarity with the VSI velocity pattern is retained. This is also true for slices at ϕ positions outside the vortex, shown in (e). There we also find the pattern of positive and negative vz annuli, although their strength and width again is diminished with respect to the axisymmetric case. Interestingly, we also find a turbulent ring structure at the raidal position of the outer vortex in (e). This could indicate the vortex is embedded in a larger flow structure or the occurrence of small vortices which are elongated in z-direction. Comparing (e) and (f) in the inner part of the disk, we find the radial positions of the annuli described above do not coincide, indicating the annuli are not axisymmetric anymore for these larger disks.

Vortex lifetime
For planetesimal formation the lifetime of a vortex is a critical factor determining the trapping efficiency. Therefore we are interested in the lifetime of the vortices generated by the instability. A distinctive signature of a large anticylonic vortex is a minimum in local vorticity stretching out over a fraction of a disk annulus. We therefore apply a box filter in radial and azimuthal direction to the midplane z-vorticity values. In azimuth, we apply a filter width of 48 and 96 grid cells for p45 and p90 respectively and 200 grid cells for both p180 and p360. In radial direction, the filter width is 30 grid cells,equalling 2H. The filter damps fluctuations on scales smaller than the filter width, enhancing the visibility of large scale structures. To outline the radial position of the large vortices, we then find the minimum in the deviation of the smoothed ωz from the azimuthal average in the respective annulus. The calculation of the deviation from the azimuth ensures we exclude the large azimuthal flow structure occurring at the outer boundary in figure 8. Figure 11 shows the evolution of this value as a function of radius and time, outlining the radial position of the vortices in the disk over the run of the simulation.
We observe multiple long living vortices emerging in run p180 and two in run p360, seen in the lower panel in figure  11. They start with distances of 1 to 5 pressure scale heights and are observed to migrate inwards with a migration rate lower than 0.0005 R0/Ω. Paardekooper,  find similar migration rates for vortices in laminar 2D vertically integrated disks. Though their results are not directly transferable to the case presented here, we think the underlying mechanism is the most plausible explanation of the phenomenon observed. The outward migration of the inner vortex in p360 can be explained by the occurrence of a surface density maximum moving outward. The vortex is trapped in the maximum and is dragged outward again, as it cannot drift across it (Paardekooper, Lesur & Papaloizou, 2010). Figure 11 in conjunction with fig.8 also shows vortices interacting and after some time merging once one vortex migrates within one pressure scale height radial distance to another.
In contrast, we don't observe the same in the simulations p45 and p90, shown in the upper row of figure 11. We do not see any sign of larger vortices in simulation p45, which is in agreement with (Richard, Nelson & Umurhan, 2016). The case p90 shows a vortex forming after around 150 and another after 550 and existing at least intermittently. This suggests the phi range of a simulation to be crucial in the development and sustaining of large, long lived vortices.

RWI as secondary Instability
In this section, we now focus on evidence in our simulations that suggests the triggering of the Rossby-Wave-Instability (RWI) to initiate vortices. We therefore calculate the values of the critical function L as given in equation 2 from an azimuthally averaged density, pressure and velocity structure of our simulation. We also averaged the field quantities over θ before calculating L . Figure 12 shows L /L0 as a function of radius after 200, 500 and 700 orbits. In the top row we show the results from p45. All plots show multiple strong global maxima, suggesting the disk to be in principle susceptible to the RWI. The amplitude of the extrema decreases with time but they stay as distinctive features throughout the simulation time.
In the bottom row we show L /L0 for p360. Here we also find extrema, but fewer in number and not as distinct as for p45. Also, the overall amplitude of the extrema is lower as in the other case.
Taking a closer look at the lower right plot in figure 12, we find the most distinctive peaks at R=1.0,1.3 and 1.45. Comparing these to the initial locations of the vortices we find in the lower left subplot of figure 11, we find them to be in good agreement. This is clear evidence that the RWI is triggering the large vortices we observe in our simulations. The decrease in amplitude of the extrema after saturation of the instability has been shown for artificial RWI vortices by (Meheut et al., 2010) and can explain the decrease we witness in the p360 case.
To investigate this, we also look at the azimuth of the kinetic energy spectrum in the midplane of the disk.
F denotes the Fourier transform of the respective velocity component. (Li et al., 2000) showed the RWI has a maximum growth rate for the azimuthal wave number m in the range of m = 3 − 6. These wavenumbers can only be accessed for simulations with an azimuthal extent equal or larger than 90 degree. Therefore the RWI is expected to only grow inefficiently in the p45 simulation. Figure 13 shows the radial average of the azimuthal component of the kinetic energy spectrum after 700 orbits.
We find a broken power law for E(m) in all our simulations. At smaller wavenumbers m, the slope falls as m − 5 3 , as expected for the upward Kolmogorov cascade in a rotationally dominated flow above the Rhines scale (Rhines, 1975), i.e. where the eddy turn over time is longer than the rotational period of the system. At large wavenumbers, the spectrum behaves ∝ m −5 , a scaling in agreement with a 2D enstrophy downward cascade. This result indicates that our simulations are predominantly 2D/rotation dominated due to the fast rotating flow in the disk.This is supported by our finding that the Reynolds stresses saturate at a later time than the rms-velocity. The rms-velocity directly traces the turbulence generated at smaller scales induced by the VSI, but αRey traces the angular momentum transported at larger scales. The difference in saturation time is then explained with the observed inverse energy cascade in the disk, which has to transport the energy generated at small scales to the larger scales on which angular momentum is transported in the disk. This is supported by the fact that we find from figure 13 that the energy injection of the VSI occurs at about m = 40 − 60, which is about 1 − 1.5 pressure scale heights.
At the smallest wavenumbers m, we find a maximum for E(m) at m = 3 and m = 4 for the cases with φmax = 180 • and 360 • , whereas for the other cases the energy piles up at m = 4 and m = 8 for φmax = 90 • and 45 • , the largest respective wavenumber accessible. For the cases with φmax ≤ 180 • , the large vortices formed by the RWI extract energy from the flow, whereas the absence of large vortices in the other simulations forces the energy to be deposited in the largest mode accessible in the simulation. This can also explain the systematically higher rms-velocities we found for φmax = 45 • and 90 • in figures 4 and 5.
To confirm we are in the rotation dominated regime in our simulations, we plot the Rossby number with u being the flow velocity at length scale l. The Rossby number gives the relative importance of Coriolis forces vs. inertial forces and is smaller unity if Coriolis forces are nonnegligible for the flow in the system. We can express this number also as a function of wavenumber m using the kinetic energy spectrum to define the velocity spectrum.
We plot Ro(m) as a function of azimuthal wavenumber m in figure 14. We find the Rossby number stays below unity at all scales accessible in our simulations, confirming we are indeed rotationally dominated. We also find Ro ≈ const. at the largest wave numbers m, indicating a change of slope for E not prominently visible in figure 13. Of course we 0.5 1.0  only investigate one component of the turbulent spectrum, because neither the radial nor the vertical wave numbers are as easy measurable as the azimuthal component. VSI is known to favour high radial and low vertical wave-numbers, yet this investigation was not feasible with the current data sets, but shall be attempted in the future. In any case the low values for Ro are an indication that all our turbulence is strongly rotational dominated and thus most likely not isotropic anyway. Because Ro < 1 in our simulation, we now revisit equation 12 . Cuzzi et al. (2001) developed the equations under the assumption that the disk turbulence is isotropic and follows a Kolmogorov dissipation law, which implies that turbulence has to be in the Ro > 1 regime, which is clearly not the case in our work. From equation 15 we can estimate the wavenumber at which our system should satisfy this condi-  Figure 14. Rossby number as a function of azimuthal wavelength at 700 orbits. We find Ro to rise at large scale and later fall towards larger scales, with a maximum at the energy injection scale. The Rossby number stays below unity at all scales resolved in our simulation.
tion. For a typical velocity vrms = 0.01v k,0 we get a value much larger than the Nyquist-wavenumber of our setup. This suggests that the rms-velocities in the disk act on much larger scales than the turbulence scales, and the ansatz presented in Cuzzi et al. (2001) is not always valid in protoplanetary disks. We however consider only the azimuthal direction in our analysis, leaving the possibility of a different scaling in the other spatial directions. Nevertheless, even if there might be more energy at large radial and vertical wavenumbers kR, kz ≡ 600/R, this still means that turbulence is also on small scales highly anisotropic, asking for revisiting the Kolomogorov picture for small scale turbulence in protoplanetary disks. Also, the results of section 4 show the ansatz works quite well despite vrms being generated on different scales than assumed for the turbulent α.

Influence of elliptic instability
In section 4 we described a turbulent substructure occurring inside the large vortices in our simulation runs p180 and p360. Lesur & Papaloizou (2009) showed a similar effect for vortices with aspect ratio χ > 8 influenced by elliptic instability. The larger vortices in our simulations have aspect ratios 8 < χ < 11 and fit well with the predictions. We do not observe the predicted decay of these vortices by the elliptic instability, similar to the reports of Lesur & Papaloizou (2009). In how far this an effect of an inactive elliptic instability (EI), maybe suppressed by low resolution, or an overlap with the VSI is open for future investigations. At χ = 10 the minimum growth time of the elliptic instability should be about 17.6 Orbits. But maybe the vortices receive additional driving from absorbing smaller vortices (as observed for Jupiter's Red Spot), which so far counteract destruction from the EI. The occurrence of the EI also naturally explains the fast destruction of the small vortices we observe in all simulation runs. These vortices have χ < 4 and are strongly influenced by the fast growing modes of the elliptic instability reported for this size regime.

SUMMARY & CONCLUSIONS
We have performed full 3D hydrodynamical simulations of protoplanetary disks undergoing Vertical Shear Instability. We used four different azimuthal extents ranging between π/2 and 2π to asses the influence of non-axisymmetry on the development of the instability.
We summarize our main findings as follows: • We find the Vertical Shear Instability to be capable of seeding vortices with large (χ > 8) aspect ratios using the disk parameters p = − 2 3 , q = −1 and H R = 0.1. This has to our knowledge not been reported for simulations of this instability. The vortices we observe are long lived (lifetimes larger than 500 local orbits) and can aid in the growth of planetesimals in protoplanetary disks via particle trapping as proposed by Barge & Sommeria (1995).
• We find the angular momentum transport in VSI to be sufficiently efficient and in agreement with latest assumptions for protoplanetary disks and constraints on planetesimal formation therein (Drążkowska & Alibert, 2017). Interrestingly, despite most angular momentum is transported outward at large z, the radial mass flux is also outward at large z but inward close to the midplane. This is due to additional strong vertical transport of angular momentum in agreement with (Stoll, Kley & Picogna, 2017). This leads to transport of angular momentum away from the midplane into upper layers, where it is then transported outwards.
• As a direct consequence, depending on the height z above the midplane, the ration of vrms(z) to √ α ranges between 2.0 and 3.5. This height depending ratio has to be taken into consideration when inferring disk α values from rms-velocity measurements.
• The choice of the size of the simulation domain, especially of φmax has a significant impact on the outcome of the simulation. We find for φmax < 90 no indication for larger, long lived vortices forming. Furthermore we find a system-atic increase of α-values and rms-velocities with decreasing φmax.
Our findings for the VSI on smaller scales are therefore consistent with the results reported by Nelson, Gressel & Umurhan (2013), Stoll & Kley (2014, 2016, Richard, Nelson & Umurhan (2016) and Flock et al. (2017) for the low azimuthal extent cases, while identifying their limitations on more global disk scales. We propose future work on instabilities in protoplanetary disks include global 360 • simulations to identify possible low m effects suppressed in current simulations of disks sections.
Furthermore, we find resolution as an important factor, both for the development of the RWI and the EI in global disk simulations. We suggest the lower radial resolution used in Stoll & Kley (2014) as a possibility why they did not find vortex formation in their work. Although they use different parameters in their work (p = −1.5, q = −1 and H R = 0.05), we do not find significant differences to their work in preliminary results of a parameter study we conduct. The results of this study will be the subject of a follow up publication. We also propose further high resolution simulations should be performed to investigate the influence of the EI on the longevity of the vortices generated.  Figure A2. Rms velocities for the 2D resolution study in vertical direction. We find convergence for the domain sizes larger or equal N θ = 128. For comparison, we again also plot the values from the 3D model presented in figure 4.

APPENDIX A: 2D RESOLUTION STUDY
We performed a resolution study in 2D to confirm we resolve the instability in the radial and meridional direction. We conducted 3 runs with radial domain sizes of Nr = 128, 256 and 512 and N θ = 128 and 3 additional runs with Nr = 256 and N θ = 64, 128 and 256. In figures A1 and A2 we show the RMS-velocities for the initial growth phase. We find convergence for radial domain sizes Nr ≥ 256 and N θ ≥ 128, proving the values we chose for the 3D simulations to be sufficient to capture the relevant physics in the early growth phases of our simulations. We also show the slope for a growth rate of Γ = 0.42 per orbit as a guidance value (black dashed curve).The growth rates found for the 2D simulations are in good agreement with the results presented above for the 3D case, supporting the validity of our choice of simulation domain parameters.