High Resolution Parameter Study of the Vertical Shear Instability
MMNRAS , 1–13 (2019) Preprint August 21, 2020 Compiled using MNRAS L A TEX style file v3.0
High Resolution Parameter Study of the Vertical ShearInstability
Natascha Manger, , (cid:63) Hubert Klahr , Wilhelm Kley , and Mario Flock Max Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Ave, New York, NY 10010, USA Institut für Astronomie und Astrophysik, Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Theoretical models of protoplanetary disks have shown the Vertical Shear Instability(VSI) to be a prime candidate to explain turbulence in the dead zone of the disk.However, simulations of the VSI have yet to show consistent levels of key disk turbu-lence parameters like the stress-to-pressure ratio α . We aim to reconcile these differentvalues by performing a parameter study on the VSI with focus on the disk densitygradient p and aspect ratio h := H/R . We use full 2 π
3D simulations of the disk forchosen set of both parameters. All simulations are evolved for 1000 reference orbits, ata resolution of 18 cells per h. We find that the saturated stress-to-pressure ratio in oursimulations is dependent on the disk aspect ratio with a strong scaling of α ∝ h . , incontrast to the traditional α model, where viscosity scales as ν ∝ αh with a constant α . We also observe consistent formation of large scale vortices across all investigatedparameters. The vortices show uniformly aspect ratios of χ ≈ and radial widths ofapproximately 1.5 H . With our findings we can reconcile the different values reportedfor the stress-to-pressure ratio from both isothermal and full radiation hydrodynam-ics models, and show long-term evolution effects of the VSI that could aide in theformation of planetesimals. Key words: planets and satellites: formation, protoplanetary discs, hydrodynamics,turbulence.
The mechanisms governing angular momentum transport inprotoplanetary disks are still not well understood. The gasmolecular viscosity is orders of magnitude to weak to facil-itate the angular momentum transport required to explainthe timescales found in observations. To allow neverthelessfor a simple description, Shakura & Sunyaev (1973) intro-duced a parametric turbulent viscosity prescription wherethe strength of the viscosity depends only on a dimensionlessparameter α relating the viscous stresses linearly to thermaland magnetic pressure, see eq. (5) below.Many physical explanations yielding sufficiently largevalues for the α -parameter have since been proposed. Thefirst strong candidate was the Magneto-Rotational Insta-bility (MRI, Balbus & Hawley, 1991), an instability act-ing in rotating disks where the equations of ideal Magneto-Hydrodynamics (MHD) are applicable. Since then, furtherresearch has shown that the ideal MHD equations are onlysatisfied in close proximity to the protostar ( < AU) or (cid:63)
E-mail: nmanger@flatironinstitute.org high up in the disk atmosphere where gas densities are lowand ionisation is created by stellar irradiation. In the gasrich midplane of the disk, non-ideal MHD processes sup-press the MRI (e.g. Turner et al., 2014; Dzyurkevich et al.,2010). Additionally, Lenz et al. (2019) showed that the tur-bulent α ∼ − generated by the MRI is too high to ex-plain the solid mass distribution in the solar system. How-ever, recent models suggest that magneto-centrifugal windslaunched from the disk could facilitate angular momentumtransport on the level of α ∼ − (Béthune et al., 2017).In recent years, pure hydro-dynamical models haveemerged to explain turbulent angular momentum transportin the absence of strong coupling to magnetic fields. Amongthese are the Convective Overstability (COS, Klahr & Hub-bard, 2014; Lyra, 2014) and its non-linear cousin SubcriticalBaroclinic Instability (SBI, Klahr & Bodenheimer, 2003; Pe-tersen et al., 2007a,b; Lesur & Papaloizou, 2010), the ZombieVortex Instability (ZVI, Marcus et al., 2015, 2016; Barrancoet al., 2018) and the Vertical Shear Instability (VSI Gol-dreich & Schubert, 1967; Fricke, 1968; Nelson et al., 2013).All these instabilities rely on the fact that baroclinic disksthat cool on a certain timescale can become unstable de- c (cid:13) a r X i v : . [ a s t r o - ph . E P ] A ug N.Manger et al. spite them being stable according to the Solberg-Hoilandcriteria (Rüdiger et al., 2002; Urpin, 2003; Arlt & Urpin,2004). Each of these instabilities works in a specific rangesof cooling times: The ZVI requires the disk to cool on verylong time-scales ( τ cool (cid:29) / Ω ), whereas the COS and SBIneed the cooling time to be on the order of τ cool ∼ / Ω ,and the VSI requires the cooling of the disk to be very fast( τ cool (cid:28) / Ω ), or even instantaneous (Lin & Youdin, 2015).Both the VSI and the COS have been shown to facilitate an-gular momentum transport at the level of α ∼ − − − (Raettig et al., 2013; Lyra, 2014; Stoll & Kley, 2014, 2016;Flock et al., 2017; Manger & Klahr, 2018).Richard et al. (2016) showed the VSI to be capableto spawn small, short lived vortices in localized disk mod-els, identified by Latter & Papaloizou (2018) as generatedby a parasitic Kelvin-Kelmholtz instability. In Manger &Klahr (2018), we showed the VSI to be also capable to gen-erate large scale, long lived vortices via the Rossby-Wave-Instability (RWI, Lovelace et al., 1999; Li et al., 2000, 2001).These vortices are high pressure maxima and are thereforeable to trap dust particles efficiently (Weidenschilling, 1977;Barge & Sommeria, 1995), making them prime locationsfor planetesimal formation via gravo-turbulent mechanisms.Ricci et al. (2019) showed that these vortices should be de-tectable with ALMA and ngVLA.In this work, we present the results of the first high res-olution 3D parameter study conducted for protoplanetarydisks with active Vertical Shear Instability. Because full 3Dsimulations are resource consuming, we focus specifically onthe parameters p and h controlling the radial density gradi-ent and the disk aspect ratio, respectively. Because we areinvestigating the behaviour of the VSI within a controlledparameter set, we choose to employ a simplified cooling lawin favour of a more physically correct radiative model (e.g. asin Stoll & Kley, 2016; Flock et al., 2017), because it allowsus the freedom to set uniform cooling parameters throughoutthe disk. To make the simulations comparable, we choose thenumerical setup to achieve a resolution of approximately 18grid cells per scale height in all directions in all simulations.Our goals in this study are three-fold:(i) Investigation on whether the results obtained in pre-vious 2D axisymmetric studies of the VSI (Nelson et al.,2013; Richard et al., 2016) are applicable to full 3D disks.We will show that the α -value generated by the VSI turbu-lence scales with the disk’s aspect ratio as α ∝ h . , whilebeing otherwise consistent with previous works.(ii) Reconciliation of differences in turbulence parametersfound when comparing recent work on full 3D simulationsutilising different disk parameters (e.g Stoll & Kley, 2014;Flock et al., 2017; Manger & Klahr, 2018).(iii) Further investigation into vortex formation withinthe framework of the VSI and it’s relation to the disk pa-rameters. We show that vortex formation is ubiquitous inVSI turbulent disks and that the vortex size is related tothe disk aspect ratio.This paper is structured as follows: In section 2 webriefly describe the methods and initial conditions used inthis study. In section 3 we present our simulation results,and in section 4 we take a closer look at the vortices formedin the disk. Section 5 presents theoretical explanations toour findings and connections to previous studies. Finally, in section 6, we summarise our results and present conclusionsto our work. As in our previous work Manger & Klahr (2018, hereafterMK18) we study a three-dimensional section of a protoplan-etary disk numerically solving the equations of ideal hydro-dynamics with a specific cooling prescription as given below.We use a setup similar to the one used in MK18. The ini-tial conditions are defined in force equilibrium, where thedensity is given by ρ = ρ (cid:18) RR (cid:19) p exp (cid:18) − Z H (cid:19) , (1)and the initial angular velocity of the disk by v φ = Ω K R (cid:34) q − q R √ R + Z + ( p + q ) (cid:18) HR (cid:19) (cid:35) . (2)While in the simulations we use spherical polar coordinates, R and Z in the above equations denote cylindrical coor-dinates, and H ( r ) is the half thickness (vertical pressurescale height) of the disk. We use an ideal equation of state ρe = Pγ − with the pressure defined as P = c ρ . Here, theisothermal sound speed, c = c (cid:16) RR (cid:17) q , is given as a func-tion of radius, where we choose q = − in all our simulations.All other quantities used are defined in table 1.Because we are interested in the dynamics of the VSIwith changing cooling law, we use a simple prescription with dPdt = − P − ρ c , init τ relax (3)where, to force the disk to be approximately isothermal, wedefine the relaxation time τ relax equal to the simulation timestep dt ∼ − π Ω K . The influence of longer cooling times onthe VSI will be explored in a future publication.For our computations we use the multi-purpose MHDGodunov code PLUTO (Mignone et al., 2007) with the hllcsolver (Toro, 2009). We use the piecewise parabolic method(PPM) by Mignone (2014) for the spatial reconstruction anda 3rd-order Runge-Kutta time integration. All simulationsperformed are listed in table 2 with their parameters. Wechoose the grid sizes listed in column 3 of table 2 to ensureall simulations run share a common resolution of circa 18cells per H in all directions.We use periodic boundary conditions in azimuthal di-rection and modified outflow conditions in vertical direction,where we ensure zero inflow into the domain and addition-ally extrapolate the Gaussian density profile into the ghostzones. In radial direction, we use reflective boundaries com-bined with buffer zones, where we relax all variables to theirinitial values. These buffer zones are excluded in our analysisof the simulations. MNRAS , 1–13 (2019)
SI parameter study Symbol Definition Description
R, φ, Z cylindrical coordinates r, θ, ϕ spherical coordinates R , Z cylindrical reference coordinates ρ density ρ reference density H = c s Ω K disk pressure scale height h = HR = c s v K disk geometric scale height P = c ρ pressure c s = c (cid:16) RR (cid:17) q/ isothermal sound speed c = kT µm H reference sound speed p = d log ρ d log R radial density slope q = d log T d log R = − radial temperature slope e specific internal energy density γ = 1 . adiabatic index v φ azimuthal velocity v K = (cid:113) GM (cid:63) R Kepler azimuthal velocity Ω K = v K R − Kepler angular frequency τ relax temperature relaxation time T r,φ eqn. 4 r φ component of viscous stress tensor α = T r,φ P turbulence parameter v rms eqn. 6 root mean squared velocity E eqn. 8 spectral specific kinetic energy ˜ v i eqn. 9 velocity in i direction in fourier space ω = ∇ × (cid:126)v vorticity L eqn. 16 RWI criterion Table 1.
List of all symbols used within this work.
In this section, we present the results of our parameter studyon the influence of radial density gradient and disk aspectratio on the vertical shear instability, focusing particularlyon the angular momentum generating stresses and gas rms-velocities.
To measure the strength of the VSI turbulence in the disk,we calculate the T r,φ component of the Reynolds tensor T r,φ ≡ (cid:104) ρδv r δv φ (cid:105) = (cid:104) ρv r v φ (cid:105) − (cid:104) ρv r (cid:105)(cid:104) v φ (cid:105) , (4)where (cid:104) (cid:105) denotes an average in azimuth. T r,φ measures thestrength of angular momentum transport generated by thedisk turbulence. To present the values in a non-dimenionalfashion, we normalise T r,φ by the azimuthally averaged pres-sure to obtain α as defined by Shakura & Sunyaev (1973):
100 200 300 400 500 𝑡 / (2𝜋𝛺 −10 ) −8 −7 −6 −5 −4 −3 −2 𝛼 h = 0.1h = 0.05 p = -0.66p = -1.5
50 100 150 200 250 300 𝑡 / (2𝜋𝛺 −10 ) −11 −10 −9 −8 −7 −6 −5 −4 −3 𝛼 h = 0.1h = 0.07h = 0.05h = 0.03 Figure 1.
Cumulative space and time average of the stress-to-pressure ration for simulations with different density slopes p (top) and aspect ratios h at constant p = − . (bottom). α = T r,φ (cid:104) P (cid:105) . (5)In figure 1 we present the cumulative space and timeaveraged values of α . The top panel of figure 1 comparessimulations with two different density slopes p = − . and p = − . at two different values of the disk aspect ratio h = 0 . and h = 0 . , whereas the bottom panel comparessimulations with four different values of h and a constantvalue of p = − . .The time evolution of the simulations with h = 0 . shows a rapid growth of the α value in the first few tens oforbits of the simulation, after which a slower growth phaseleads to growth to the final saturated phase of the turbu-lence after around 100 orbits, where values of α = 9 · − arereached. These values are comparable the ones we reportedfor our simulations in MK18 for simulations with lower az-imuthal but similar radial and meridional resolution. A sim-ilar behaviour is observed for the simulations with h = 0 . ,which also show a first strong growth phase up to around50 orbits, after which slower growth phase is observed untilaround 300 orbits, where a steady state value of α = 1 · − is reached.Comparing all four simulation runs, we find that thedensity slope p does not significantly influence the aver-age value of the turbulent stresses. The comparison however MNRAS , 1–13 (2019)
N.Manger et al. model r in , out /R grid size ( N r × N θ × N ϕ ) p h (cid:104) α (cid:105) / − Γ [2 π/ Ω] p0.6h0.1 . − . × × − .
66 0 . . ± . . p1.5h0.1 . − . × × − . . . ± . . p1.5h0.07 . − . × × − . .
07 2 . ± . . p0.6h0.05 . − . × × − .
66 0 .
05 1 . ± . . p1.5h0.05 . − . × × − . .
05 1 . ± . . p1.5h0.03 . − . × × − . .
03 0 . ± . . Table 2.
List of simulation and model parameters. From left to right: model name, radial domain size, numerical grid size, densityslope parameter, disk aspect ratio and space and time averaged stress to pressure value. The vertical and azimuthal domain sizes for allsimulations are θ = ± . H and φ = 0 − π , respectively. −2 −1 0 1 2 Z/H 𝛼 h = 0.1h = 0.05 p = -0.66p = -1.5 −2 −1 0 1 2 Z/H −4 −3 𝛼 h = 0.1h = 0.07h = 0.05h = 0.03 Figure 2.
Stress-to-pressure ratio as a function of height. Thetop panel shows the dependence on the density slope p for twovalues of h . The bottom panel shows the dependence on diskaspect ratio h for a slope p = − . . The values are averaged inradial and azimuthal direction and from 600-1000 orbits. shows a clear correlation of the disk aspect ratio with theturbulent α values of the disk. This result has been expected,as the disk aspect ratio is proportional to the disk tempera-ture and therefore a larger value of h leads to a larger overalltemperature in the disk and to stronger turbulent velocitiesas visible in figure 3. We therefore ran additional simulationswith aspect ratios h = 0 . and h = 0 . , presented in thelower panel of figure 1 along the simulations discussed above.In comparing these four simulation, we see a clear trend with disk aspect ratio emerging: Simulations with larger h show overall stronger turbulent angular momentum trans-port, with α = 9 · − for h = 0 . going down to · − for h = 0 . . A full list of averaged saturated α -values is listedin 2, where the errors listed are calculated for fluctuations intime only. The additional simulations also confirm the trendfor later onset of turbulent growth and longer times untilsaturation.Figure 2 shows the dependence of α on height above themidplane, where the top panel again compares two differentdensity slopes at h = 0 . and h = 0 . , and the bottompanel for different h at p = − . . The values are averagedin radial and azimuthal direction and between 600 and 1000reference orbits. We again find no evidence that the initialdensity gradient p influences the generated stresses. For both h = 0 . and h = 0 . the deviations between the curvesrepresenting p = − . and p = − . are minor and can beexplained with statistical effects. In both the top and bottompanel the trend of overall increasing α with increasing h is observed. We also find that higher h leads to a largerdifference between the α values in the midplane and theupper disk layers, though the general dependence of α withheight described in MK18 is found in all simulations. This isalso in accordance with simulations presented by e.g. Stoll& Kley (2014), Stoll et al. (2017) and Flock et al. (2017).In Appendix A, we present a resolution study of theazimuthal direction. We find that the choice of resolutionin this direction is not influencing the time evolution of the α -value, even at one-eighth the fiducial resolution. This islikely due to the axisymmetric nature of the initial linearVSI growth. However, at the lowest resolution the verticalprofile of α becomes steeper, as the resolution is only 2 cellsper H . We now take a look at the initial growth phase of the VSIby calculating the rms-velocity defined as v rms = (cid:113) ( v r − (cid:104) v r (cid:105) ) + ( v θ − (cid:104) v θ (cid:105) ) + ( v ϕ − (cid:104) v ϕ (cid:105) ) , (6)with brackets representing spatial averages in φ direction.The averages (cid:104) v r (cid:105) and (cid:104) v θ (cid:105) are assumed equal to zero.In figure 3 we plot the spatial average of the rms-velocityas a function of time for the first 200 orbits of the simula-tions. We find that the simulations for h = 0 . with different MNRAS , 1–13 (2019)
SI parameter study 𝑡 / (2𝜋𝛺 −1 ) −3 −2 −1 𝑣 r m s / 𝑐 s , h = 0.1h = 0.05 p = -0.66p = -1.5 𝛤 = 0.36𝛤 = 0.2 𝑡 / (2𝜋𝛺 −10 ) −5 −4 −3 −2 −1 𝑣 r m s / 𝑐 s , h = 0.1h = 0.07h = 0.05h = 0.03 𝛤 = 0.38𝛤 = 0.3𝛤 = 0.2𝛤 = 0.08
Figure 3.
Rms-velocities normalised by reference sound speedmeasured for the initial growth phase of the instability. The toppanel compares the disks with p = − . and p = − . at h =0 . and h = 0 . , while the bottom panel compares different h from 0.03 to 0.1 at a common p value. We additionally plot curveswith exponential growth rates Γ for comparison. values of p start their growth at a similar time and with iden-tical growth rates of Γ = 0 . per orbit. We also observe thesecondary growth phase seen already in the time evolutionof the α value. A similar qualitative behaviour is found forthe simulation with h = 0 . . The simulations both settleto a steady state value of v rms ≈ . c s . The onset of theVSI growth is later than for the h = 0 . cases in agreementwith the results from the analysis of the α stresses. The sec-ond growth phase is however not observed in the simulationswith h = 0 . , instead we find a plateau in the v rms valuesafter the initial growth phase and further growth is only ob-served after additional 50 orbits. For both values of p wefind the initial growth rate to be Γ = 0 . per orbit.In the bottom panel of figure 3 we present the rms-velocities for the simulations with p = − . and differentvalues of h . For the simulation with h = 0 . we find agrowth rate of Γ = 0 . which is in good agreement with thelinear scaling predicted by analytical calculations (Nelsonet al., 2013). For h = 0 . we however get Γ = 0 . , whichis lower than we would expect based on the growth rates ofthe other simulations in combination with a linear scaling,from which a growth rate of Γ = 0 . would be expected.A possible explanation is that the simulation with h = 0 . was run using the FARGO scheme to increase computational −2 −1 0 1 2 Z/H 𝑣 r m s / 𝑐 s , h = 0.1h = 0.05 p = -0.66p = -1.5 −2 −1 0 1 2 Z/H 𝑣 r m s / 𝑐 s , h = 0.1h = 0.07h = 0.05h = 0.03 Figure 4.
Rms-velocities as a function of distance from the diskmidplane. The averages are calculated for the remaining dimen-sions and over the simulation time. Simulations shown in eachpanel ar the same as in figure 3. efficiency, which resulted in a time step time step twice thanin the other simulation. As the cooling time in our setup islinked to the simulation time step, the simulation thus hasa slightly longer cooling time, possibly influencing the VSIgrowth.In figure 4 we display the vertical dependence of therms-velocities. In the top panel, we again find the curves cor-responding to the same h to agree well with each other, sup-porting our claim that the initial density gradient does notinfluence the simulation result. For h = 0 . the rms-velocityat the midplane is v rms = 0 . c s and rises to v rms = 0 . c s inthe upper layers. These results are consistent with the re-sults obtained in MK18 taking into account that we includedthe φ component of velocity in the calculation of v rms pre-sented in this work. Comparing the results for different h value shown in the bottom panel of figure 4, we find thatthe overall shape of the v rms stays similar for all cases, al-though the difference between the velocity at z = 0 and z = 2 H increases with increasing h . A similar correlationexists between h and the v rms value at z = 0 , which is largerfor higher h . This is expected, as a warmer disk has moreshear energy available to be converted into turbulence.Because of the vertical variation of the velocity fluctu-ations and the fact that the growth of an instability is gov-erned by the largest and not the root-mean-square velocity, MNRAS , 1–13 (2019)
N.Manger et al. 𝑡 / (2𝜋𝛺 −10 ) −3 −2 −1 𝑣 m a x / 𝑐 s , h = 0.1h = 0.07h = 0.05h = 0.03 𝛤 = 0.4𝛤 = 0.28𝛤 = 0.2𝛤 = 0.075
Figure 5.
Maximum perturbed velocity as a function of time forthe initial growth phase for different values of h. We again showexponential growth rates as a comparison. we calculate the maximum velocity as v max = (cid:113) v r + v θ . (7)Here, we consider only the poloidal motion and neglect the φ component of velocity in this case as the VSI has beenshown to grow axisymmetrically. The results are shown infigure 5 for the simulations with p = − . , where we againplot growth rates for comparison. We find the growth ratesobtained using v max to be in good agreement with the onesobtained using v rms . To gain additional insight into the turbulent behaviour ofthe VSI we look at the specific kinetic energy spectrum de-fined as E ( m ) = (cid:104)| ˜ v R | + | ˜ v θ | + | ˜ v φ | (cid:105) R (8)where m is the azimuthal wave number and ˜ v k is the Fouriertransform of the k component of velocity in the azimuthaldirection defined by ˜ v k ( m ) = 12 π ∞ (cid:90) −∞ v k ( φ ) e − imφ d φ . (9)We plot the results as a function of azimuthal wavenumberin figure 6. We find that the azimuthal kinetic energy spec-trum is not consistent with a 3D Kolmogorov turbulent en-ergy cascade, which would be ∝ m − / . This was expected,as the VSI produces weak turbulence that is always and atall scales dominated by the disk rotation. This means, if theinjection scale has a Rossby-number Ro (defined as the ra-tio of inertial to coriolis forces) of smaller than unity, thenthe turbulence will always be in the Ro < regime at allscales and will not develop isotropic 3D turbulence at anyscale (Sharma et al., 2018). Instead, the results show thatthe spectrum is proportional to m − with a turn-over intoa shallower power law at low m , which is consistent with anupward cascade transporting energy to larger scales com-bined with an downward cascade transporting enstrophy, aquantity related to the square of vorticity, to smaller scales 𝑚 −16 −14 −12 −10 −8 −6 −4 −2 ⟨ ∑ 𝑖 ( ̃ 𝑣 𝑖 ( 𝑚 )) ⟩ 𝑟 ∼ 𝑚 −5 ∼ 𝑚 −5/3 p = -0.66 h = 0.1p = -1.50 h = 0.1p = -0.66 h = 0.05p = -1.50 h = 0.05p = -1.50 h = 0.07p = -1.50 h = 0.03 Figure 6.
Specific kinetic energy spectrum as a function of az-imuthal wavenumber. For comparison, we show both the m − / behaviour expected for 3D isotropic Kolmogorov-like turbulenceand the closest matching slope m − . −1 L/H −16 −14 −12 −10 −8 −6 −4 ⟨ ∑ 𝑖 ( ̃ 𝑣 𝑖 ( 𝑚 )) ⟩ 𝑟 p = -0.66 h = 0.1p = -1.50 h = 0.1p = -0.66 h = 0.05p = -1.50 h = 0.05p = -1.50 h = 0.07p = -1.50 h = 0.03 Figure 7.
Specific kinetic energy spectrum as a function of tur-bulent length-scale in disk scale heights. The dashed line marksthe scale L = 3 H , in close proximity of which the spectrum of allsimulations shows a break in the power-law. (Lyra & Umurhan, 2019). This behaviour is characteristic of2D turbulence (Kraichnan, 1971) and is consistent with find-ings from our previous work presented in MK18, althoughfor lower h this behaviour is less clearly visible. We also findthe maxima in the spectrum for m = 2 − reported pre-viously, although these are also less clearly visible for lowerh. We also observe h = 0 . differ from this behaviour, herethe energy is deposited at m ≥ .In figure 7 we plot the same spectrum as a functionof the turbulent length scale L = 1 /k associated with thewavenumber m in units of the pressure scale height H . Wefind that the spectrum flattens to a shallower power law forall cases at or around L = 3 H , indicating a universal energyinjection scale. We also observe that the energy is depositedat scales ranging from L ≈ − H in all simulations,explaining the shift in preferred wavenumber of the RWI atdifferent h/r . MNRAS , 1–13 (2019)
SI parameter study 𝜙 𝑅/𝑅 𝜙 𝑅/𝑅 𝜔 𝑧 / 𝛺 K Figure 8.
Midplane vorticity for different values of p and h.The top row shows simulations with p = − . and the bottomrow with p = − . . For the left column, the aspect ratio of thesimulations is h = 0 . , resulting in less extended structures thanseen in the right column where h = 0 . . To identify vortices in our simulations, we use the verticalcomponent of the vorticity, defined as the rotation of thevelocity vector (cid:126)v : ω z = ( ∇ × (cid:126)v ) z . (10)In this section, we use the vorticity scaled by the orbitalfrequency Ω to identify and characterise the vortices formingin the disks at different values of p and h. Figures 8 and 9 show the midplane value of the vertical vor-ticity for the simulations with different p and h after 900reference orbits. We find large vortices forming in all oursimulations irrespective of the values assumed for p and h .Each simulation has formed between 2 and around 10 vor-tices at this point in the simulation, but there is no correla-tion of the number of vortices with the chosen parameters.The strength of the vortices however shows a correlationwith the disk aspect ratio, best seen in figure 9. The vor-tices in the bottom left panel corresponding to the simula-tion p1.5h0.1 appear to have lower absolute vorticity thanthe simulations in the top row corresponding to simulationsp1.5h0.03 and p1.5h0.05. Because the disk has overall vor-ticity ω z = Ω K , the lower absolute vorticity in the centre ofthe vortices in the lower right panel corresponds to a largerrelative vorticity and therefore vortex strength.In figure 9 it can also be seen that in both simulations inthe bottom row two or more vortices are in close radial prox-imity to each other. Therefore the vortices influence eachother and each will eventually merge into one vortex. 𝜙 𝑅/𝑅 𝜙 𝑅/𝑅 𝜔 𝑧 / 𝛺 K Figure 9.
Midplane vorticity for different values of the disk as-pect ratio h . The top left panel shows the simulation with thesmallest h=0.03, increasing to h = 0 . in the top right, h = 0 . in the bottom left and h = 0 . in the bottom right panel. Allsimulations have initially p = − . . Note the increase in vorticityand decrease of overall structure with increasing h . We also find that the run p1.5h0.03 (top left panel)shows a ordered band structure in vorticity, which does notappear in the other simulations at this stage, although thesimulations with h = 0 . show axisymmetric bands at someradii. The band structure however appears for all simula-tions during the growth phase of the VSI, but it breaks downsoon after. We discuss the implications of this in section 5.In Appendix A we show the influence of the azimuthalresolution on the vorticity. For a resolution of n φ ≥ we find large vortices forming, corresponding to a min-imal resolution of 4 cells per H needed to resolve thenon-axisymmetries of the non-linear stage of the VSI andthe generation of the vortices through the Rossby-Wave-Instability mechanism as discussed in section 5. Because we are also interested in the size and shape of thevortices formed, we choose one vortex from each simulationfor a close up inspection. Figures 10 and 11 show the vorticestaken from figures 8 and 9, respectively, in a local coordinatesystem centred on the vortex. We express the radial andazimuthal coordinate as a function of the local pressure scaleheight for easier comparison. The centre coordinates of eachpanel are listed in the corresponding figure caption.We find that the vortices share a common size of 1-1.5 local scale heights in radial diameter and an azimuthalsize between 10 and 20 H , leading to aspect ratios χ = r ∆ φ/ ∆ r in the range of 8.5 - 20. Figure 11 again shows thedependence of the relative vorticity on the disk scale height,but no correlation of vortex size or aspect ratio with diskaspect ratio h can be found. Figure 10 seems to indicatethat the vortices are larger for simulations with p = − . , MNRAS , 1–13 (2019)
N.Manger et al. −1.0−0.50.00.51.0 𝑅 − 𝑅 c 𝐻 𝑐 −15 −10 −5 0 5 10 15 (𝜙 − 𝜙 c )𝑅 c 𝐻 c −1.0−0.50.00.51.0 𝑅 − 𝑅 c 𝐻 𝑐 −15 −10 −5 0 5 10 15 (𝜙 − 𝜙 c )𝑅 c 𝐻 c 𝜔 𝑧 / 𝛺 K Figure 10.
Vortex extent for different values of p and h . The coordinate systems are positioned at the centres of the vortices, R c , φ c . Thesimulations and the respective centre coordinates used are: p0.66h0.05 with R c , φ c = 1 . , . (top left); p0.66h0.1 with R c , φ c = 1 . , . (top right); p1.5h0.05 with R c , φ c = 1 . , . (lower left) and p1.5h0.1 with R c , φ c = 1 . , . (lower right). −1.0−0.50.00.51.0 𝑅 − 𝑅 c 𝐻 𝑐 −15 −10 −5 0 5 10 15 (𝜙 − 𝜙 c )𝑅 c 𝐻 c −1.0−0.50.00.51.0 𝑅 − 𝑅 c 𝐻 𝑐 −15 −10 −5 0 5 10 15 (𝜙 − 𝜙 c )𝑅 c 𝐻 c 𝜔 𝑧 / 𝛺 K Figure 11.
Vortex extent for different values of h and p = − . . The simulations shown, and the respective centre coordinates used,are: p1.5h0.03 with R c , φ c = 1 . , . (top left); p1.5h0.05 with R c , φ c = 1 . , . (top right); p1.5h0.07 with R c , φ c = 1 . , . (lowerleft) and p1.5h0.1 with R c , φ c = 1 . , . (lower right). but as the vortices depicted from the simulation p1.5h0.1 arecurrently merging, this could be artificial. This is supportedby the fact that the top left panel of figure 8 shows runp0.66h0.05 to also form vortices with a sizes comparable tothe ones found in p1.5h0.05. To track the radial position of the vortices over time, we usethe same technique as employed in MK18. Therein, we cal-culate the radial positions of the vortices in each timestepby first applying a box filter to the vertical vorticity to elim-inate all structures smaller than 1 H in radius and 6 H inazimuth. Then the azimuthal average of the vorticity is sub-tracted from this to exclude possible zonal flows and thenthe minimum value in azimuthal direction is calculated.The results are shown in figures 12 and 13. In all cases,we find the formation of long lived vortices at the radial po-sitions of the vortices depicted in figures 8 and 9, thoughthe time after which the stable vortices appear seems notdirectly correlated to neither density gradient nor aspect ra-tio. For example, the simulations p1.5h0.05 and p0.66h0.1show stable vortex formation early on, while p1.5.h0.1 firstshows intermittent large vortices before forming long timestable vortices. The same can be found upon inspecting ad- ditional frames of the run p1.5h0.07 which in figure 13 onlyshows strong stable vortices emerging well after 750 orbits,but vortices can be found much earlier in the simulation,though their constant interaction makes it hard to detectthem with this method. An exception in both the time ofonset of vortex formation and the number of vortices formedis again run p1.5h0.03 which shows only one stable vortexforming after around 600 orbits, though this could be againattributed to the longer evolution timescale of the VSI itself.Upon inspection of the time series of the midplane vor-ticity, we also find that even if a vortex has established itself,new vortices can be formed at the same radial distance tothe star, of which one example can be seen in the bottom leftpanel of figure 9. The newly formed vortex catches up withthe dominant vortex after a few tens to hundred orbits andis eventually absorbed into the stronger vortex. We also findevidence of this in other runs, e.g. p0.66h1.5, but as bothexamples occur for vortices located close to the reference ra-dius we cannot exclude that this happens for other vorticesalso. The effect is most easily observed close to R becausewe take one snapshot after each completed orbit at R , lead-ing the azimuthal position of the vortices with R > R and R < R to drift due to the radial dependence of Ω .Different to our results from MK18, we observe vortexdestruction in our simulations with h = 0 . . In both sim- MNRAS , 1–13 (2019)
SI parameter study 𝑅 / 𝑅 𝑡/(2𝜋/𝛺 K ) 𝑅 / 𝑅 𝑡/(2𝜋/𝛺 K ) −0.10−0.08−0.06−0.04−0.020.00 m i n 𝜙 ( 𝜔 𝑧 − ⟨ 𝜔 𝑧 ⟩ 𝜙 ) Figure 12.
Evolution of the vortex radial position as a functionof time. The panels are sorted as in figure 8 with columns corre-sponding to aspect ratio and rows to density gradient. The colourshows the azimuthal minimum of the vertical vorticity subtractedby the azimuthally averaged vorticity. To extract only larger scaleminima, we apply an image filter prior to the calculation. 𝑅 / 𝑅 𝑡/(2𝜋/𝛺 K ) 𝑅 / 𝑅 𝑡/(2𝜋/𝛺 K ) −0.10−0.08−0.06−0.04−0.020.00 m i n 𝜙 ( 𝜔 𝑧 − ⟨ 𝜔 𝑧 ⟩ 𝜙 ) Figure 13.
Same as figure 12, but now for the different diskaspect ratios as in figure 9. ulations, vortices formed early and close in in the disk aredestroyed after ca. 500 and 300 orbits for p = 0 . and p = − . , respectively. For the case p=-0.66 this is due to avortex forming radially close to the vortex and interactingwith it until it is destroyed. For the case p = − . some-thing similar seems to happen, but we cannot exclude theinfluence of the boundary in this case.We also observe the migration of the vortices over timein all cases except p0.66h0.05 and p1.5h0.03. For the formerthis is likely due to the large amount of vortices formed in the disk which prevents radial migration due to the interactionof the vortices with each other. For the latter, it is simplythe late time of the formation of the vortices that prohibitsus from detecting migration, though it is possible that itoccurs later on. For the migrating cases, the migration isstronger for simulations with larger aspect ratio. This can beexplained by the fact that the vortices formed in disks withlarger h are stronger, but we cannot exclude the influenceof the disk surface density gradient, which changes most forthe cases with large h due to overall mass loss. In this section, we discuss our results presented 3 and com-pare them to other models. We start with the turbulentvelocities and then present an interpretation of the relation-ship between the stress-to-pressure ratio α and h . In section 3, we showed the growth rate of the VSI in oursimulations with h = 0 . to be Γ ≈ . and Γ = 0 . forthe simulations with h = 0 . . These values of Γ are ingood agreement with the values obtained in MK18 and Stoll& Kley (2014), respectively. The observed linear scaling ofmost of the reported values with aspect ratio is expectedby linear theory. Nelson et al. (2013) showed that for nearlyisothermal disks the growth rate can be to first order ex-pressed as Γ ∼ | q | h Ω . (11)To compare our results with the theoretical ones pre-sented in Lin & Youdin (2015) we calculate the h indepen-dent growth rate σ = Γ h Ω . (12)For the simulations with h = 0 . , . and 0.05 we obtain σ = 0 . , which is in good agreement with the results pre-sented by Lin & Youdin (2015) for a disk simulation with z max = 3 H and radial wavenumber k=30. Because we ob-tained a smaller than expected growth rate for h = 0 . ,this case also has a lower σ and does not align well withtheir largest expected growth rate.We additionally observed that the onset of VSI growthis increasingly delayed with decreasing aspect ratio. The oc-currence of this trend is consistent with the results of Nelsonet al. (2013), who observed this behaviour for the kinetic en-ergy in their 2D axisymmetric simulations. Our simulations with h = 0 . show a saturated α of · − .This value compares well with the values we obtained inMK18. It is however about one order of magnitude largerthan the values reported in both Flock et al. (2017) andFlock et al. (2020), who reach a similar scale height in theirradiation hydrodynamics simulation. This is likely due tothe cooling time varying between τ cool = 10 − − − or-bital periods in their simulations in contrast to our fixedcooling time of approximately − orbits. An investigation MNRAS , 1–13 (2019) N.Manger et al. of the influence of the cooling time on the VSI strength ishowever beyond the scope of this work. For our simulationswith h = 0 . we find a saturated α of − , which is com-parable to the values of a few times − reported for thehigh resolution cases in Stoll & Kley (2014).When we plot the α values in our 3D simulations as afunction of the pressure scale height h (see Fig.14) we finda scaling of 2.6, which seeks for an explanation. For thatpurpose we measure the strength of the VSI, i.e. the r.m.s.velocity of the gas as a function of scale height h (see Fig.15)and find it to scale linearly, which is easy to explain.The VSI is driven by vertical shear S = | ∂ z v φ | = R | ∂ z Ω | (see Eq. 2) over a scale height H : S · H = R Ω K (cid:18) HR (cid:19) q = c s q h, (13)thus if v r . m . s . ∝ SH , then we find: v r . m . s . ∝ c s q h, (14)exactly what we find in Fig. 15. We ran additional 2D simu-lations (listed in table 3) to confirm this trend over an evenwider range of h . Although α as defined in equation 5 is notstrictly applicable to 2D simulations, we used a time averageto determine the steady state v φ velocity component and theresults align well to the linear scaling law retrieved from the3D data for the simulations with h < . . The higher val-ues of v rms and α for the 2D simulations with h > . canbe explained by the non-linear state of the 2D simulationsdeveloping large scale motions not present in 3D.As now α measures the Reynolds stresses scaled withthe speed of sound, we find α ∝ v . m . s . c s ∝ h , (15)clearly leading to a steep rise in α as function of h . Thisestimated scaling is not as steep as the one we measured inFig.14, which probably needs additional effects. One sugges-tion would be the decreased pitch angle of the spiral den-sity waves that actually transport the angular momentum(Rafikov, 2002) with h . Thus h modifies the correlation be-tween radial v r and azimuthal v φ velocity fluctuations. Inother words, the fewer pressure scale heights fit into thecircumference of the disk, the less tightly the spirals arewrapped and the stronger the angular momentum transportbecomes. In MK18 we put forward the hypothesis that the Rossby-Wave-Instability (RWI) is working as a secondary instabilityin our simulations. (Lovelace et al., 1999) showed that theRWI develops in disks that form regions with an extremumin L , defined as L = Σ2 ω z (cid:18) P Σ γ (cid:19) /γ . (16)Within a region around this extremum, standing Rossbywaves develop and amplify. Eventually multiple vorticesform within this region which subsequently merge into onelarge vortex. Looking at L as a function of radius supportsthis: Figure B1 shows L after 200 and 500 orbits, where the model r in , out /R N r × N θ × N ϕ h (cid:104) α (cid:105) / − p1.5h0.01 . − . × × . . ± . p1.5h0.02 . − . × × .
025 0 . ± . p1.5h0.05 . − . × × .
05 1 . ± . p1.5h0.1 . − . × × . ± p1.5h0.2 . − . × × . ± Table 3.
List of additional 2D simulation parameters. From leftto right: model name, radial domain size, numerical grid size, diskaspect ratio and space and time averaged stress to pressure value.The vertical and azimuthal domain sizes for all simulations are θ = ± . H and φ = 0 − π , respectively, and we assume p = − . throughout. −2 −1 ℎ −5 −4 −3 −2 𝛼 𝛼 ∝ ℎ Figure 14.
Time and space averaged α -values as a function ofdisk scale height. Clearly visible is the power law relationshipbetween both, which we determine to be α ∝ h . for the 3Dsimulations. Although α is not strictly defined in 2D, we find goodagreement with the 3D data using a time averaged to determine (cid:104) v φ (cid:105) . −2 −1 ℎ −1 𝑣 r m s / 𝑐 s , fit 𝑣 𝑟𝑚𝑠 /𝑐 𝑠 ∝ ℎ data 2Ddata 3D Figure 15.
Time and space averaged v rms -values as a function ofdisk scale height. Clearly visible is the linear relationship betweenboth. The outliers at high h can be explained by the non-linearstate of the VSI in 2D. MNRAS , 1–13 (2019) SI parameter study profile of L changes significantly in amplitude for the caseswith h = 0 . and . where the RWI has formed vorticeswell before 500 orbits, whereas in the remaining runs RWI isstill developing, and therefore L retains its large amplitude.For those cases, it is more instructing to look at L after 500and 700 orbits, respectively. There we can see that for thecase with h = 0 . the amplitude of L at the position ofthe vortex is declining between the two times, and the vortexappears after around 600 orbits. This finding is in line withMeheut et al. (2010), who showed that the extremum in L diminishes after the instability saturates. For the case with h = 0 . the vortices form very late in the simulation, sothe extremum grows as time progresses, and we see a strongextremum only at time t = 700 orbits, which will eventuallydevelop into a vortex about 50-100 orbits later.Additionally, we find small vortices forming in all sim-ulations with h ≥ . , which can be seen in figures 10 and11. These smaller vortices have been previously identified inRichard et al. (2016) and are in line with the work of Latter& Papaloizou (2018), who showed that a parasitic Kelvin-Helmholtz instability (KHI) develops off the VSI modes,leading to the formation of small vortices. We however arguethat these KHI-vortices are transient and that the instabil-ity is not responsible for the larger long-lived vortices wefind in our simulations.In our simulations with h = 0 . and h = 0 . presentedin section 4, we also find azimuthal bands of low vorticityin addition to the vortices forming in the disk, generatedby the vertical structure of the body modes of the VSI. Weinterpret these azimuthal bands as zonal flow like structuresforming from the VSI body modes in the disk, representingthe axisymmetrical nature of the VSI. If these modes indeedact as zonal flows, they likely still exist also in the simula-tions with large h as an underlying construct obscured bythe high turbulence level present. This would allow for theformation of dust rings in the presence of e.g. magnetic fieldsinhibiting the long term survival of the vortices.The vortices shown in figures 10 and 11 also show inter-nal turbulent structure. As our simulations are conducted in3D, the vortices are in principle susceptible to the EllipticInstability (EI, Lesur & Papaloizou, 2009). But the EI iseasily suppressed by numerical diffusion, especially for vor-tices with χ (cid:38) . We therefore argue that the EI is notsignificantly effecting the vortices in our simulations. Theintrinsic turbulent structure of the vortices can more eas-ily be explained with the intrinsic turbulence generated bythe VSI and the small Kelvin-Helmholtz parasite vorticesgenerated by it (Latter & Papaloizou, 2018). In this work we present the first 3D high resolution param-eter study on the Vertical Shear Instability. We focus ourattention on how the initial density slope p of the disk andthe disk aspect ratio h impacts the VSI and vortex forma-tion. The dependence on the cooling time and temperaturegradient will be the the focus of future work. For now, wehave assumed τ ≈ − and q = − .We find that the VSI is capable to support angular mo-mentum transport with α values up to α = 10 − for thelargest scale height of h = 0 . . The α values we find in our simulations are comparable but slightly lower than the val-ues we reported in MK18 and the values found by Stoll &Kley (2014) for the respective h used. The values reporteddiffer however from values obtained from radiation hydro-dynamics simulations of the VSI (Flock et al., 2017, 2020),which will need to be be investigated in the future.With decreasing scale height, the α value decreasesdown to a value of a few times − for the lowest scaleheight investigated. Overall, we find α to scale with h roughly as h . . We explain this behaviour with the driv-ing shear of the VSI, which increases as h , generating a de-pendence of α ∝ h . The additional h dependence can beattributed to the formation of spiral density waves, whichalso transport angular momentum (Rafikov, 2002). The de-pendence of α on h could have a profound influence on thedisk structure which should be investigated further.The growth rates of the VSI we find are in good agree-ment with previous studies of MK18 and Stoll & Kley(2014). In the vertical direction, the quantitative behaviourwe observed in MK18 is recovered. Contrary to the aspectratio of the disk, the density gradient does not have any in-fluence on the disk turbulence. This result is expected, asthe VSI itself is not sensitive to the density gradient.We find that the VSI is able to seed multiple vorticesin all considered parameter combinations and the vorticesin all simulations live for hundreds of orbits, reinforcing thefindings of MK18. The vortices generally have a radial diam-eter between 1 and 1.5 local scale heights, which is close tothe diameter of H allowed by the disks radial shear. Mostvortices have aspect ratios of χ (cid:39) but values of up to 20 arefound. The time at which the first stable vortex appears isnot correlated with any of the investigated parameters, andit is likely that such a time is random. We also do not finda correlation between the number of vortices found simulta-neously in the disk and any of the investigated parameters.We also observe that the turbulence sustained by the VSIsteadily creates new vortices even at radii at which a largescale vortex has already been established for a longer periodof time. This could indicate that the VSI constantly worksto replenish the vorticity gradient that drives the Rossby-Wave-Instability, which once replenished sufficiently createsnew vortices. These new, weaker vortices are then eventuallyabsorbed by the larger vortex already present.In this work, we considered the disk to be comprisedpurely of gas. In reality, the disk contains about 2% solids,which have been shown to change the buoyancy of the diskand can therefore have an impact on the growth of the VSIespecially near the midplane (Lin, 2019; Schäfer et al., 2020).Future simulations should therefore investigate whether vor-tices also emerge in a dusty disk when the VSI is suppressedclose to the midplane. We also neglect the influence of mag-netic fields in this work. Cui & Bai (2020) showed thatthe VSI can in principle coexist with magnetically launchedwinds at the surface, but future work should address whetherthe non-axisymmetric effects found in this work are affectedby non-ideal MHD processes. ACKNOWLEDGEMENTS
This research has been supported by the DeutscheForschungsgemeinschaft Schwerpunktprogramm (DFG
MNRAS , 1–13 (2019) N.Manger et al.
References
Arlt R., Urpin V., 2004, A&A, 426, 755Balbus S. A., Hawley J. F., 1991, Astrophysical Journal, 376, 214Barge P., Sommeria J., 1995, A&A, 295, L1Barranco J. A., Pei S., Marcus P. S., 2018, ApJ, 869, 127Béthune W., Lesur G., Ferreira J., 2017, A&A, 600, A75Cui C., Bai X.-N., 2020, ApJ, 891, 30Dzyurkevich N., Flock M., Turner N. J., Klahr H., Henning T.,2010, Astronomy and Astrophysics, 515, A70Flock M., Nelson R. P., Turner N. J., Bertrang G. H.-M.,Carrasco-González C., Henning T., Lyra W., Teague R., 2017,ApJ, 850, 131Flock M., Turner N. J., Nelson R. P., Lyra W., Manger N., KlahrH., 2020, ApJ, 897, 155Fricke K., 1968, Zeitschrift für Astrophysik, 68, 317Goldreich P., Schubert G., 1967, Astrophysical Journal, 150, 571Klahr H. H., Bodenheimer P., 2003, ApJ, 582, 869Klahr H., Hubbard A., 2014, The Astrophysical Journal, 788, 21Kraichnan R. H., 1971, Journal of Fluid Mechanics, 47, 525Latter H. N., Papaloizou J., 2018, MNRAS, 474, 3110Lenz C. T., Klahr H., Birnstiel T., 2019, ApJ, 874, 36Lesur G., Papaloizou J. C. B., 2009, A&A, 498, 1Lesur G., Papaloizou J. C. B., 2010, A&A, 513, A60Li H., Finn J. M., Lovelace R. V. E., Colgate S. A., 2000, ApJ,533, 1023Li H., Colgate S. A., Wendroff B., Liska R., 2001, ApJ, 551, 874Lin M.-K., 2019, MNRAS, 485, 5221Lin M.-K., Youdin A. N., 2015, ApJ, 811, 17Lovelace R. V. E., Li H., Colgate S. A., Nelson A. F., 1999, ApJ,513, 805Lyra W., 2014, ApJ, 789, 77Lyra W., Umurhan O. M., 2019, PASP, 131, 072001Manger N., Klahr H., 2018, MNRAS, 480, 2125Marcus P. S., Pei S., Jiang C.-H., Barranco J. A., HassanzadehP., Lecoanet D., 2015, ApJ, 808, 87Marcus P. S., Pei S., Jiang C.-H., Barranco J. A., 2016, ApJ, 833,148 Meheut H., Casse F., Varniere P., Tagger M., 2010, A&A, 516,A31Mignone A., 2014, Journal of Computational Physics, 270, 784Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O.,Zanni C., Ferrari A., 2007, The Astrophysical Journal Sup-plement Series, 170, 228Nelson R. P., Gressel O., Umurhan O. M., 2013, Monthly Noticesof the Royal Astronomical Society, 435, 2610Petersen M. R., Julien K., Stewart G. R., 2007a, ApJ, 658, 1236Petersen M. R., Stewart G. R., Julien K., 2007b, ApJ, 658, 1252Raettig N., Lyra W., Klahr H., 2013, ApJ, 765, 115Rafikov R. R., 2002, The Astrophysical Journal, 569Ricci L., Flock M., Blanco D., Lyra W., 2019, arXiv e-prints, p.arXiv:1902.01897Richard S., Nelson R. P., Umurhan O. M., 2016, MNRAS, 456,3571Rüdiger G., Arlt R., Shalybkov D., 2002, A&A, 391, 781Schäfer U., Johansen A., Banerjee R., 2020, A&A, p. A190Shakura N. I., Sunyaev R. A., 1973, Astronomy and Astrophysics,24, 337Sharma M. K., Verma M. K., Chakraborty S., 2018, Physics ofFluids, 30, 115102Stoll M. H. R., Kley W., 2014, A&A, 572, A77Stoll M. H. R., Kley W., 2016, A&A, 594, A57Stoll M. H. R., Kley W., Picogna G., 2017, A&A, 599, L6Toro E. F., 2009, Riemann Solvers and Numerical Methods forFluid Dynamics. Springer-Verlag Berlin HeidelbergTurner N. J., Fromang S., Gammie C., Klahr H., Lesur G., WardleM., Bai X.-N., 2014, Protostars and Planets VI, pp 411–432Urpin V., 2003, A&A, 404, 397Weidenschilling S. J., 1977, Monthly Notices of the Royal Astro-nomical Society, 180, 57 MNRAS , 1–13 (2019)
SI parameter study model grid size ( N r × N θ × N ϕ ) (cid:104) α (cid:105) / − n128 256 × × . ± . n256 256 × × . ± . n512 256 × × . ± . n1024 256 × × . ± . Table A1.
List of parameters for the resolution study. The simulation n1024 is identical to p1.5h0.1 from table 2. 𝑡 / (2𝜋𝛺 −10 ) −8 −7 −6 −5 −4 −3 −2 𝛼 𝑛 𝜙 = 128𝑛 𝜙 = 256𝑛 𝜙 = 512𝑛 𝜙 = 1024 −2 −1 0 1 2 Z/H 𝛼 𝑛 𝜙 = 128𝑛 𝜙 = 256𝑛 𝜙 = 512𝑛 𝜙 = 1024 𝑡 / (2𝜋𝛺 −10 ) −3 −2 −1 𝑣 r m s / 𝑐 s , 𝑛 𝜙 = 128𝑛 𝜙 = 256𝑛 𝜙 = 512𝑛 𝜙 = 1024 −2 −1 0 1 2 Z/H 𝑣 r m s / 𝑐 s , 𝑛 𝜙 = 128𝑛 𝜙 = 256𝑛 𝜙 = 512𝑛 𝜙 = 1024 Figure A1.
Turbulence analysis of the simulations performed for the resolution study. The top row shows the α value as a function oftime on the left and the the dependence on height on the right. The bottom row shows the same for the rms-velocity. APPENDIX A: AZIMUTHAL RESOLUTION STUDY
To determine the minimum grid resolution in the azimuthal direction required to achieve a converged nonlinear state of theVSI and form vortices, we performed a small resolution study. For all simulations we used the parameters of the simulationrun p1.5h0.1 described in table 2 and varied only n φ . The simulations performed are listed in table A1.Figure A1 shows the results of the parameter study. We find that both the spatial average of α and v rms show convergedvalues which is the expected outcome as the VSI modes grow axisymmetrically. Therefore the formation of vortices does notseem to influence the overall strength of the turbulent angular momentum transport.We also find that the cases n256, n512 and n1024 give similar results when looking at the vertical profiles of α and v rms .Contrary to this, the simulation n128 shows steeper profile of both values with height than the other simulations. Looking atfigure A2, we find that n128 is the only case in which no vortex formation is observed. We find a similar behaviour in ouranalysis presented in MK18 for the α value for different spatial sizes in the azimuthal direction, where the profile is also steeperfor the case where we exclude that vortices form in the disk. Therefore, the formation of vortices could have an influence onwhere the bulk of the angular momentum transport is facilitated with respect to the midplane of the disk. APPENDIX B: SECONDARY ROSSBY WAVE INSTABILITY
We plot the Rossby Wave Instability criterion from equation 16 as a function of radius in figure B1 for the simulations with p = − . after 200 and 500 orbits. We show the existence of multiple extrema, some of which align well with the radiallocation of the vortices highlighted in figure 11 (dashed red lines). MNRAS000
We plot the Rossby Wave Instability criterion from equation 16 as a function of radius in figure B1 for the simulations with p = − . after 200 and 500 orbits. We show the existence of multiple extrema, some of which align well with the radiallocation of the vortices highlighted in figure 11 (dashed red lines). MNRAS000 , 1–13 (2019) N.Manger et al. 𝜙 𝑅/𝑅 𝜙 𝑅/𝑅 𝜔 𝑧 / 𝛺 K Figure A2.
Midplane vorticity after 500 orbits for the simulations considered in the resolution study. Clockwise from top left the phiresolutions are 128, 256, 1024 and 512. MNRAS , 1–13 (2019)
SI parameter study ℒ / ℒ ℎ = 0.03 ℎ = 0.05 ℎ = 0.07 ℎ = 0.1 𝑅/𝑅 ℒ / ℒ 𝑅/𝑅 𝑅/𝑅 𝑅/𝑅 ℒ / ℒ Figure B1.
Instability criterion for the RWI (Lovelace et al., 1999) (see eq. 16 in the text) for different h . The values are averaged invertical and azimuth, and show the simulations after 200 (top row), 500 orbits (middle row) and 700 orbits (bottom row). The red linesindicate the positions of the vortices shown in figure 11.MNRAS000