Vortex formation in protoplanetary discs induced by the vertical shear instability

6 MB·

by  Samuel Richard

Mon.Not.R.Astron.Soc.000,1–??() Printed11January2016 (MNLATEXstylefilev2.2)
Vortex formation in protoplanetary discs induced by the
vertical shear instability
Samuel Richard1(cid:63), Richard P. Nelson1 & Orkan M. Umurhan2,3
1 Astronomy Unit, Queen Mary University of London, Mile End Rd, London, E1 4NS, U.K.
6
2 Space Sciences Division, NASA Ames Research Center, Moffett Field, CA 94035, USA
1
3 SETI Institute, 189 Bernardo Way, Mountain View, CA 94043, USA
0
2

n
a
J
ABSTRACT
8
We present the results of 2D and 3D hydrodynamic simulations of idealized proto-

planetary discs that examine the formation and evolution of vortices by the vertical
]
P shearinstability(VSI).Inagreementwithrecentwork,wefindthatdiscswithradially
E decreasing temperature profiles and short thermal relaxation time-scales, are subject
. to the axisymmetric VSI. In three dimensions, the resulting velocity perturbations
h
give rise to quasi-axisymmetric potential vorticity perturbations that break-up into
p
discretevortices,inamannerthatisreminiscentoftheRossbywaveinstability.Discs
-
o withveryshortthermalevolutiontime-scales(i.e.τ (cid:54)0.1localorbitperiods)develop
r strong vorticity perturbations that roll up into vortices that have small aspect ratios
t
s (χ (cid:54) 2) and short lifetimes (∼ a few orbits). Longer thermal time-scales give rise to
a vortices with larger aspect ratios (6 (cid:54) χ (cid:54) 10), and lifetimes that depend on the
[
entropy gradient. A steeply decreasing entropy profile leads to vortex lifetimes that

1 exceed the simulation run times of hundreds of orbital periods. Vortex lifetimes in
v discs with positive or weakly decreasing entropy profiles are much shorter, being 10s
1 oforbitsatmost,suggestingthatthesubcriticalbaroclinicinstabilityplaysanimpor-
2 tant role in sustaining vortices against destruction through the elliptical instability.
9
Applied to the outer regions of protoplanetary discs, where the VSI is most likely to
1
occur, our results suggest that vortices formed by the VSI are likely to be short lived
0
structures.
.
1
Keywords: accretion,accretiondiscs-hydrodynamics-instabilities-protoplanetary
0
6 discs - turbulence
1
:
v
i
X
1 INTRODUCTION Anumberofhydrodynamicinstabilitieshavebeensug-
r gested as vortex formation mechanisms. The Rossby wave
a The presence of long-lived vortices in protoplanetary discs
instability (RWI; Lovelace et al. 1999; Li et al. 2000, 2001)
has long been considered as a means of enhancing the
has been shown to produce large scale vortices from an ini-
planet building process because of the ability of these anti-
tialaxisymmetricstatewhen anonself-gravitatingdischas
cyclonic structures to capture and concentrate dust grains
asufficientlystronglocalminimuminthepotentialvorticity
(von Weizsa¨cker 1944; Barge & Sommeria 1995; Johansen
(Umurhan2010;Lovelace&Hohlfeld2013;Yellin-Bergovoy
etal.2004).Thepresenceofvorticesindiscscanalsoleadto
etal.2015).Theextremumcanbeduetoalocalmaximum
significant transport of angular momentum through the ex-
insurfacedensity,suchascanappearattheedgeofthedead
citation of spiral density waves (Johnson & Gammie 2005).
zone where a sharp change in effective disc viscosity arises
Inspiteofthepotentialimportanceofvorticesfortheevolu-
(e.g. Varni`ere & Tagger 2006), or because a planet opens
tionofprotoplanetarydiscsandplanetformation,however,
a gap in the disc (de Val-Borro et al. 2007). Another way
the following questions do not yet have definitive answers:
to produce vortices is through baroclinic instability. Klahr
If vortices exist in protoplanetary discs, what are their for-
&Bodenheimer(2003)introducedtheconceptoftheglobal
mationmechanisms?;Whichregionsofprotoplanetarydiscs
baroclinicinstability(GBI),andsuggestedthatadiscwitha
can support the existence of vortices?; What are the life-
global negative entropy gradient could form vortices. More
times of vortices in protoplanetary discs?
recently, Petersen et al. (2007) and Petersen et al. (2007)
have shown that vortices can form when the disc has an
unstable radial stratification and undergoes thermal relax-
(cid:63) E-mail: [email protected]; [email protected]; ation (or ‘cooling’). Lesur & Papaloizou (2010) showed in
[email protected]
(cid:13)c RAS
2 S.Richard, R.P.Nelson & O.M.Umurhan
their study that this instability is actually a nonlinear in- 2 THEORETICAL BACKGROUND AND
stability (requiring finite amplitude perturbations to be ac- EXPECTATIONS
tivated) and called it the subcritical baroclinic instability
Beforepresentingoursimulationresultswediscussanumber
(SBI). The nonlinearity of the instability means that the
of theoretical results that are of relevance to this numerical
finite amplitude perturbations need to be generated by an-
study of the VSI. We make use of both spherical (r, θ, φ)
otherunspecifiedprocess.Alinearinstabilitythathasbeen
and cylindrical (R, φ, Z) coordinates in this paper. The
suggestedasapossiblesourcefortheseperturbationsisthe
cylindrical coordinates are used in the formulae giving the
‘convective overstability’, in which the growth of epicyclic
discstructure(density,temperatureandvelocityprofile)for
oscillations is powered by the same unstable stratification
convenience,whilethesphericalcoordinatesareusedforthe
and thermal relaxation required for the subcritical instabil-
simulations because they fit better with the shapes of the
ity, leading to the formation of long-lived vortices (Klahr
disc models.
& Hubbard 2014; Lyra 2014). Finally, there have been re-
cent suggestions that a zombie vortex instability may arise
in stably stratified flows by the formation of a critical layer
2.1 RWI
which rolls up into vortices, which then excite new critical
layers and vortices, leading eventually to space filling tur- The first vortex-forming instability that has been studied
bulence that is dominated by large vortices (Marcus et al. in the context of protoplanetary discs with near-Keplerian
2013, 2015). rotation profiles is the RWI. Using a combination of linear
analysisandnonlinearnumericalsimulations,Lovelaceetal.
(1999), Li et al. (2000) and Li et al. (2001) showed that a
In this paper, we examine the possibility that vortices non-axisymmetric instability may develop, leading to the
may be formed by the vertical shear instability (VSI). This formationofanumberofvortices,whenthedisccontainsa
is a linear instability that was first studied by Goldreich local extremum in the function:
& Schubert (1967) and Fricke (1968) in the context of dif- Σ(cid:16) P (cid:17)2/γ
ferentially rotating stars. The presence of vertical shear in L= , (1)
ω Σγ
z
protoplanetary discs must arise when there is a radial tem-
whereω istheverticalcomponentofthevorticity,Σisthe
perature gradient, and this can lead to the destabilization z
surfacedensityandP isthepressure(ω /Σisthepotential
ofinertial-gravitywaves(oscillationsforwhichrotationand z
vorticity and P/Σγ is related to the entropy). It has been
buoyancyproviderestoringforces)whenthermaltime-scales
observed that these vortices often merge to form a single
areshorterthanviscoustime-scales,andareoftheorderof,
vortex during the advanced stages of nonlinear evolution.
or shorter than, dynamical time-scales (Urpin & Branden-
Considered within the context of protoplanetary discs, the
burg 1998; Urpin 2003; Nelson et al. 2013; Barker & Latter
RWI is normally observed to develop when a local pressure
2015; Umurhan et al. 2015). The recent study by Nelson
or density maximum is present in the disc, such as may oc-
et al. (2013) adopted simple equations of state and cool-
curattheedgeofaplanet-inducedgap(deVal-Borroetal.
ing prescriptions, and showed that this instability can lead
2006, 2007) or at the edge of a dead zone where there is
tosustainedhydrodynamicturbulencewithaShakuraSun-
yaev angular momentum transport parameter α∼10−3 for a sharp transition in the disc viscosity (Lyra & Mac Low
2012). Vortices have also been observed to develop sponta-
veryshortcoolingtimes.Stoll&Kley(2014)performednon-
neously in global magnetized disc models that sustain the
linear hydrodynamic simulations with radiation transport
magnetorotationalinstability(Fromang&Nelson2005).Al-
and showed that the instability operates in the presence of
thoughthislatterphenomenonhasnotbeenexploredinde-
a more complete description of the gas thermal evolution,
tail,apossibleexplanationisthatthevorticesarisebecause
albeitwithareducedefficiencyofangularmomentumtrans-
theRWIfeedsofftheso-calledzonalflowsthatareobserved
port.Therequirementforveryshortcoolingtimessuggests
to arise in discs with MHD turbulence (Steinacker & Pa-
thatthisinstabilityismostlikelytooperateintheouterre-
paloizou 2002; Papaloizou & Nelson 2003; Johansen et al.
gionsofprotoplanetarydiscsbeyond∼10au(Nelsonetal.
2009; Bai & Stone 2014).
2013; Umurhan et al. 2013). An analysis of linear growth
rates in discs with energy transport in both the optically
thick and thin regimes presented by Lin & Youdin (2015)
2.2 SBI
suggests that the VSI should operate at radii in the range
10 50 au. The existence of a vortex-forming baroclinic instability op-
eratinginprotoplanetarydiscswasfirstsuggestedbyKlahr
& Bodenheimer (2003), based on a series of nonlinear sim-
Thispaperisorganizedasfollows.InSect.2,wereview ulations conducted using disc models with negative radial
the different processes and instabilities that can lead to the entropygradients.ThelinearpropertiesofthisGBIwerein-
formationanddestructionofvorticesindiscs,andinSect.3 vestigatedbyKlahr(2004)andJohnson&Gammie(2005),
we describe the disc models that are the basis of our study whofoundevidenceforonlytransientgrowth.Thenonlinear
and the numerical scheme used in the simulations. The re- evolution was investigated using shearing box simulations
sults of two-dimensional, axisymmetric simulations are pre- by Johnson & Gammie (2006), who found no instability.
sented in Sect. 4, and the results of three-dimensional runs Petersen et al. (2007) examined disc stability using global
that examine the formation and evolution of vortices are anelasticsimulationsofbaroclinicdiscswithprescribedcool-
presentedinSect.5.Finallywediscussourresultsanddraw ing and observed the growth of vortices that were found to
conclusions in Sect. 6 survive for hundreds of orbits (Petersen et al. 2007). In a
(cid:13)c RAS,MNRAS000,1–??
Vertical shear instability and vortices 3
subsequentpaper,Lesur&Papaloizou(2010)usedbothin- For a nearly inviscid disc in which perturbations are no
compressible and compressible shearing box simulations to longeradiabatic,thethermalevolutionofperturbedfluidel-
examine the growth and survival of vortices in discs with ementscanremovethestabilizinginfluencesofentropygra-
radial entropy gradients and imposed cooling, finding that dients when the cooling time is short enough. This leads to
thesediscsareunstabletoafinite-amplitudeinstabilitythat the well-known Goldreich Schubert Fricke instability (Gol-
leads to the formation of long-lived vortices – the SBI. dreich & Schubert 1967; Fricke 1968), which in the context
TheSBIisanon-linearconvectiveinstabilitythatleads of accretion discs is known as the VSI (Urpin 2003). The
to the formation and amplification of vortices when the ra- VSIdevelopswhentheflowisverticallyshearedandalmost
dialstratificationsatisfiestheSchwarschildinstabilitycrite- locally isothermal. The instability criterion is:
rion:
∂j2 k ∂j2
− R <0. (7)
N2 <0 (2) ∂R k ∂Z
R Z
and when the flow undergoes thermal relaxation. N is the Vertical shear is always present in a protoplanetary
R
radial Brunt Vaisala frequency and is defined by disc, unless the flow is isothermal or homentropic, and
|∂j2/∂R|(cid:29)|∂j2/∂Z|, so the instable modes have k (cid:29)k
R Z
1 ∂P ∂S
N2 =− (3) (radial wavelengths are much shorter than vertical wave-
R C ρ∂R∂R
p lengths). In the locally isothermal limit, the maximum
whereρisthedensity,P isthepressure,S =C ln(P1/γρ−1) growth rate of the VSI depends on the temperature profile
p
and the scaleheight :
is the entropy per unit mass, and C is the specific heat
p
capacity at constant pressure. The thermal relaxation time (cid:16)H(cid:17)
Γ ∼|q| Ω, (8)
scale must be of the order of one orbital period. Too short max R
a time scale prevents vortex formation, while too long a
where q is the temperature profile power-law index. In a
time scale does not allow vortices to be amplified. When
recent study, Nelson et al. (2013) showed that the VSI
thestratificationisstable(N2 >0)vorticesareobservedto
R can cause a disc to become highly turbulent in the locally
form,butinsteadofbeingamplifiedtheydecayanddissolve
isothermal regime, with velocity perturbations having very
into the background flow. The non-linear character of this
short radial wavelengths such that there are strong local
instability implies that sufficiently strong perturbations are
gradientsintheflow.Theonlysimulationsconductedin3D
required to trigger it.
in that work utilized a locally isothermal equation of state,
and the resulting turbulence led to the excitation of spiral
density waves in the flow but no obvious signs of long lived
2.3 Convective overstability
vortices. It therefore remains an open question whether or
The convective overstability is a linear instability that de- nottheVSIcanleadtotheformationofperturbationsthat
pends on having NR2 < 0 and thermal relaxation on ∼ dy- generate vortices when thermal relaxation is not treated as
namical time-scales (Klahr & Hubbard 2014; Lyra 2014). being instantaneous, perhaps through the RWI or the SBI
The instability involves the growth of horizontal epicyclic acting on the primary perturbations generated by the VSI.
oscillations, hence the name by which it is known, and ac-
cordingtothestudybyLyra(2014)itleadstotheformation
of long lived vortices. It has been suggested as a possible 2.5 Elliptical instability
sourceofthefiniteamplitudeperturbationsrequiredforthe
Lesur & Papaloizou (2009), using a local approach, showed
SBI to operate. Vertical buoyancy has not been included
that 3D elliptical vortices may be unstable. The vertical
in the analysis of this instability so far, so its behaviour in modes ((cid:126)k =k e(cid:126)) are the dominant modes in vortices with
z z
vertically stratified discs remains unexplored.
an aspect ratio 1.5 < χ < 4, and have a growth rate given
by:
2.4 VSI (cid:115) (cid:18) (cid:19)(cid:18) (cid:19)
2Ω χ 2Ω 1
Γ=S − − − (9)
TheRayleighcriterionforhydrodynamicstabilityissatisfied S χ−1 S χ(χ−1)
by discs with strictly Keplerian rotation profiles because
whereSistheshearandχtheaspectratioofthevortex.
dj2
>0, (4) InthecaseofaKeplerianprotoplanetarydiscS =1.5Ω,and
dR
the growth rate is:
wherej =R2Ω(R)isthespecificangularmomentum.More (cid:115)
(cid:18) (cid:19)(cid:18) (cid:19)
generally, a disc with angular velocity varying with both 3 4 χ 4 1
Γ= Ω − − − (10)
radius and height, Ω(R,Z), that is subject to adiabatic, 2 3 χ−1 3 χ(χ−1)
axisymmetric perturbations is stable according the Solberg
Thefactthatthemodeispurelyverticalandthegrowth
Hoilandcriteria(e.g.Tassoul1978).Foraccretiondiscswith
rate is independent of the wavelength makes the instability
negativeradialandverticalpressuregradientsthesestability
quite easy to capture and very high resolution simulations
criteria can be written
arenotrequiredforittoberesolved.Itshouldbenotedthat
1 ∂j2 1 (cid:16)(cid:12)∂P(cid:12) ∂S (cid:12)∂P(cid:12) ∂S(cid:17)
+ (cid:12) (cid:12) +(cid:12) (cid:12) >0 (5) thisresult,however,isvalidinthelocalapproximationand
R3 ∂R ρCp (cid:12)∂R(cid:12)∂R (cid:12)∂Z(cid:12)∂Z is not true for long wavelengths.
∂j2 ∂S ∂j2 ∂S The case of larger aspect ratio vortices is more com-
− >0. (6) plex because the instability is fully three dimensional. The
∂R ∂Z ∂Z ∂R
(cid:13)c RAS,MNRAS000,1–??
4 S.Richard, R.P.Nelson & O.M.Umurhan
instability is due to the resonance between inertial waves We adopt a radial power law for the disc temperature
inanunstratifiedfloworhighfrequencybuoyancywavesin while assuming that the disc is initially isothermal in the
a stratified flow and the turnover frequency of the vortex. vertical direction, and we also adopt a radial law for the
Asnoinertialmodescanmatchtheturnoverfrequencyina midplane density:
vortex with 4<χ<5.9, these vortices are stable when the
(cid:16) R (cid:17)q
flow is not stratified, while in a stratified flow vortices are T(R) = T
0 R
always unstable. 0
(cid:16) R (cid:17)p
However, despite the unstable character of vortices in ρ(R,0) = ρ , (14)
0 R
that case, the elliptical instability is difficult to observe be- 0
cause of the high resolution needed to resolve it. The un- where R is a reference radius. Other quantities of interest
0
stable mode will have kmax ∼ k and kmax ∼ χk , so the are determined using the equations of force balance in the
φ z r z
radialresolutionneededtoresolveitdependsontheaspect radial and vertical directions:
ratio. Moreover, the growth rate for χ>6 is about 50 time
∂P GMR
smallerthanforvorticeswithχ<4.Thesetwopointsmake =− +RΩ2 (15)
∂R r3
theellipticalinstabilitydifficulttoobserveinnumericalsim-
∂P GMZ
ulations that contain large aspect ratio vortices. =− (16)
∂Z r3
√
where r= R2+Z2 is the spherical radius.
2.6 Streaming instability The equilibrium solutions for density and angular ve-
locity give
Although the streaming instability applies to a disc com-
posed of interpenetrating gas and solids, rather than to the (cid:16) R (cid:17)p (cid:16)GM (cid:104)1 1(cid:105)(cid:17)
ρ(R,Z)=ρ exp − (17)
single fluid system considered here, we discuss it briefly for 0 R c2 r R
0 s
completeness. It arises as a linear instability when aerody-
(cid:114)
namic drag causes inwards radial drift of solid particles, (cid:16)H(cid:17)2 (cid:16) R(cid:17)
Ω(R,Z)=Ω 1+(p+q) +q 1− (18)
and the backreaction on the gas is included in the dynam- k R r
ics(Youdin&Goodman2005).Thelinearlygrowingmodes
(cid:112)
consist of particle density enhancements with growth times whereΩk = GM/r istheKeplerianvelocityandH isthe
thatliebetweenfastdynamicaltime-scalesandslowerradial scaleheight defined through:
drift time-scales, and maximal growth rates arise for parti- c
H = s . (19)
clestoppingtimescomparabletodynamicaltime-scalesand Ω
k
when the local solids-to-gas ratio is of the order of unity.
c is the isothermal sound speed defined through:
Nonlinear simulations indicate that it can lead to severe s
clumpingofsolids,suchthatparticleconcentrationsareable P
c2 = , (20)
to collapse directly to form planetesimals (Johansen et al. s ρ
2007).Assuch,thisprovidesanalternativetotheconcentra-
and we define the local disc aspect ratio h=H/R.
tionofdustinvorticesasameansofformingplanetesimals
In the disc model considered here, the Brunt-Vaisala
in protoplanetary discs.
frequency in the midplane is:
T
N2 = (p+q)(p(γ−1)−q) (21)
R γR2
3 NUMERICAL METHOD
thentheradialstratificationisstablewhenq<p(γ−1)and
Theequationsofmotionareforacompressibleandinviscid unstable otherwise.
fluid subject to a central gravitational potential: These equations are solved in spherical coordinates
∂ρ +∇(cid:126).(cid:0)ρV(cid:126)(cid:1)=0, (11) (mr,eθth,φod)u(sRinicghaarfidniettevaol.lu2m01e3c)oWdeeucshinogostehereMflUecStCinLgHbaonucnodk-
∂t
aryconditionsattheradialandverticalboundariesbecause
∂ρV(cid:126) +∇(cid:126).(cid:0)ρV(cid:126)V(cid:126)(cid:1)=−∇(cid:126)P −ρ∇(cid:126)Φ, (12) of their ease of implementation, and because Nelson et al.
∂t (2013) have shown that the boundary conditions appear to
∂ρE ρ T −T have no effect on the development of the VSI.
+∇(cid:126).V(cid:126)(ρE+P)=ρV(cid:126).∇(cid:126)Φ− 0. (13) The computational domains of our simulations are as
∂t γ−1 τ
follows.Inradiustheinnerboundaryislocatedatr=1and
Here ρ is the density, V(cid:126) the velocity, P the pressure and the outer boundary at r = 1.5, and the azimuthal domain
E the total energy per unit mass, and Φ = −GM/r the runs between φ = 0 and φ = π/4. The meridional domain
gravitationalpotentialduetothecentralstar.Thesystemis extends ±5×h above and below the disc midplane. In all
closedusingtheperfectgasequationofstate:ρE =P/(γ− simulationsexceptthosewithh=0.05,weuse500gridcells
1)+1/2ρV(cid:126)2, where we adopt γ = 7/5. The last term in intheradialdirection.Formodelswithh=0.05,wedouble
the energy equation (13) is a thermal relaxation term: the the radial resolution and use 1000 grid cells. All runs use
temperature relaxes to the initial temperature T with a 300gridcellsintheazimuthaldirectionand200cellsinthe
0
relaxation time τ. τ is assumed to be a function of R, the meridional direction. The equilibrium state is perturbed by
cylindrical radius, being a fixed multiple or fraction of the adding 10−6c amplitude white noise to each component of
s
local Keplerian orbital period. the velocity in all simulations presented in this paper.
(cid:13)c RAS,MNRAS000,1–??
Vertical shear instability and vortices 5
Figure1.Kineticenergyforadiscwithah/r=0.05anddiffer- Figure2.Kineticenergyforadiscwithah/r=0.1anddifferent
entvaluesofthecoolingtime. valuesofthecoolingtime.
4 AXISYMMETRIC SIMULATIONS
An investigation of the potential for the VSI to generate
vortices in discs clearly requires 3D simulations to be per-
formed. The fact that the unstable VSI modes have radial
wavelengthsmuchshorterthanverticalscalesindicatesthat
high resolution simulations are required, leading inevitably
tolongsimulationruntimes.Inlightofthis,mostofthesim-
ulations presented in this paper adopted density and tem-
perature power-law profiles p = −1.5 and q = −2, where
the adoption of the steep temperature profile reduces the
growth time scale for the VSI, hence allowing a suite of 3D
simulationstobeundertaken.Wenotethatthesediscmod-
elsarestableaccordingtotheSolberg-Hoilandcriteriagiven
byequations(5)and(6)andhaveimaginaryvaluesofthera-
dialBrunt-VaisalafrequencyN .Dependingontheadopted
R
thermalrelaxationtimescale,thesemodelsmaybeunstable
to the SBI and the convective overstability. It is therefore
Figure3.Kineticenergyforadiscwithah/r=0.2anddifferent
possible that any vorticity perturbations generated by the
valuesofthecoolingtime.
VSI may be amplified and sustained by the SBI, leading to
long lived vortices.
Given that the parameters used for the simulations in stable to the VSI. This critical value depends strongly on
this paper were not considered in the study presented by the scaleheight. For h = 0.05, the disc become stable for a
Nelson et al. (2013), we have undertaken a suite of 2D ax- cooling time 0.05 < τ < 0.1, for h = 0.1 it stabilizes for
isymmetric simulations to examine the growth times of the 0.3<τ <0.5,andforh=0.2thediscremainsunstablefor
VSI in these models, and to also examine the critical cool- a cooling time as large as one orbital period.
ingtimesthatallowtheVSItooperate.Weconsidermodels We notice that the level of saturation also depends on
withh=0.2,0.1and0.05withdifferentthermalrelaxation the cooling time: the total kinetic energy during the satu-
times. ratedstateishigherforashortercoolingtime,anddecreases
Following Nelson et al. (2013), the total kinetic energy for longer cooling times. In other words, higher amplitude
is defined by the volume integral of the sum of the radial velocityandvorticityperturbationsareexpectedforshorter
and meridional kinetic energies, normalized by the volume cooling times when we consider 3D simulations below.
integral of the initial azimuthal kinetic energy:
(cid:82)
ρ(v2+v2)dV
Ec= V(cid:82) θ r (22) 5 3D SIMULATIONS
ρv2 dV
V φ0
5.1 Fiducial model
Theresultsofthe2DsimulationareshowninFigs.1,2
and3,whichdisplaytheevolutionofthetotalkineticenergy Theprimarygoalofthispaperistoexaminewhetherornot
in the disc models with h = 0.05, 0.1 and 0.2 respectively. thenonlineardevelopmentoftheVSIleadstotheformation
Foreachvalueofhthegrowthratedecreaseswhenthecool- of vortices. Secondary goals include determining the range
ing time increases, in agreement with Nelson et al. (2013), of conditions under which vortices form, understanding the
andforacriticalvalueofthecoolingtimethediscbecomes natureofthesevorticesasafunctionofsystemparameters,
(cid:13)c RAS,MNRAS000,1–??
6 S.Richard, R.P.Nelson & O.M.Umurhan
Figure 4.Vorticityprofileinthemidplaneandmeridionalplaneafter34,54and64orbitalperiods,withH/R=0.2andτ =0.1.
and the possible roles of the RWI and the SBI in creating The bottom panels of Fig. 4 also show that the break-
and maintaining vortices. Given our interest in the poten- upoftheinitiallyaxisymmetricvorticityperturbationsinto
tial role of the SBI, we have chosen a disc model which in discrete vortices produces structures with relatively small
principle allows the development of the VSI and SBI. Both length scales in the vertical direction. In other words, the
instabilities are more efficient in thick discs, so we choose vorticesformedinthissimulationarenotlargescalecolum-
h = 0.2. Our 2D runs described in Section 4 indicate that narstructures,butinsteadappeartobecoherentoververti-
this model remains unstable to the VSI even for relatively cal length scales that are significantly shorter than the ver-
longcoolingtimes,whicharenecessaryfortheSBItooper- tical scaleheight (i.e ∼ 0.1H). Although we only plot the
ate (Lesur & Papaloizou 2010). vorticityinthediscmidplane,wefindthatvorticesformat
all heights in the disc.
Our fiducial model has a thermal relaxation time τ =
0.1 local orbits. The 2D simulations described in Section 4 Looking in particular at the first and last of the top
indicate that the growth rate of the VSI should be quite panels in Fig. 4, we see that the vortices have relatively
large in this case, and that this cooling time is significantly small aspect ratios when projected in the R-φ plane. De-
shorter than the critical value for which the VSI no longer tailed measurements indicate that the aspect ratio χ ∼ 2
operates. During the early phases the VSI in 3D remains for these vortices. As expected from our discussion of the
axisymmetric and so develops very similarly to 2D simula- elliptical instability presented in Section 2.5, these vortices
tions (Nelson et al. 2013). The axisymmetric velocity per- do not survive for very long, and we measure a typical life
turbations correspond naturally to axisymmetric vorticity time of between 2 and 3 orbits. Inspection of animations of
perturbations which grow with time. When these axisym- the midplane vorticity indicates that this disc model devel-
metricvorticitybandsreachacriticalamplitude,theytend opsavigorousturbulentflowinwhichvorticescontinuously
todestabilizeandvorticesareformed.Thisevolutionisillus- appear and disappear on time-scales of a few orbits.
tratedbyFig.4whichshowscontoursoftheperturbedver- Examination of Fig. 4 indicates that the formation
ticalcomponentofthevorticityinthemidplane(toppanels) mechanism of the vortices is the creation of narrow, ax-
and in a slice along the meridional plane (bottom panels). isymmetric vorticity perturbations that then break up into
AsdiscussedinNelsonetal.(2013),theVSIisfirstobserved discrete vortices when the perturbation amplitude becomes
athighlatitudesinthediscanddescendsdowntowardsthe largeenough.FromourdiscussionoftheRWIinSection2.1,
midplane, as seen in Fig. 4. we can see that perturbations to the vertical component
(cid:13)c RAS,MNRAS000,1–??
Vertical shear instability and vortices 7
Figure 5. Profiles of the quantity L demonstrating the devel-
opmentofextremainthisquantityandhencetheconditionsre-
quiredfortheRWItooperateandformvortices.
Figure 6. Kinetic energy for a disc with a h/r = 0.2, τ = 0.1
of the vorticity that are narrowly confined in radius corre- anddifferentviscosity.
spond to the creation of extrema in the quantity L defined
by equation (1). As such, it appears that vortices form in
(ν = 2.5×10−8) the vorticity profile is very similar to the
this simulation because the VSI generates narrow vorticity
inviscid case. This suggests that the numerical diffusion is
perturbationsandextremainLthatdestabilizethroughthe
of the order of 10−8. This value correspond to a Reynolds
RWI.TheSBIdoesnotappeartoplayanimportantrolein
number Re=Hc /ν ≈106 and a Shakura-Sunyaev viscous
this particular simulation. s
stressparameterα=2.5×10−7 (Shakura&Sunyaev1973).
Regarding the mechanism of roll-up observed in our
Thisvalueistoosmalltoberesponsibleoftheshortlifetimes
simulations, the axisymmetric VSI is accompanied by ra-
of the vortices, so we are confident that the disappearance
dialazimuthallysymmetricpressureperturbationswhichin-
ofthevorticeshasaphysicaloriginratherthananumerical
duces jets in the zonal (azimuthal) direction. Such jets are
one.
characterized by narrow abutting azimuthal strips of pos-
itive/negative vertical vorticity (recall vorticity anomalies
are here understood to be with respect to the background
5.3 Dependence on cooling time
Keplerian frame). When the amplitude of the jet gets large
enough, the RWI induces the roll-up of the negative vortic- Inthissection,weinvestigatehowvortexformationdepends
ity strip leaving the positive vorticity strip more or less in- on the thermal relaxation by considering the cooling times
tactsincepositivevorticityanomaliesinnonself-gravitating τ = 0.05 and 0.5, while keeping all other disc parameters
discs are stable to the RWI (Umurhan 2010). An examina- the same. Inspection of Fig. 3 shows that the τ =0.05 run
tion of the top row of Fig.s in Fig. 4 shows exactly this shouldleadtorapidgrowthoftheVSI,andthegenerationof
pattern, where the negative vorticity anomalies created by relativelylargevelocityofvorticityperturbations.Acooling
theVSIeventuallyrollupintolocalizedvorticitywhileleav- time of τ =0.5 should lead to a longer growth time for the
ing the positive vorticity strips alone. The profile of L at VSI and weaker velocity and vorticity perturbations.
different times during this simulation is shown in Fig. 5, il- Fig. 8 shows contours of the vertical vorticity pertur-
lustratingthedevelopmentofthesecondaryRWIascaused bations for the τ =0.05 run. Comparing the spectrum bars
by the VSI throughout many stages and radial locations of for the contours shown in Fig. 4 and 8 shows that the lat-
the simulation. terrungeneratessignificantlylargervorticityperturbations.
The consequence of this is that the vortices formed in the
nonlinearsaturatedstateofthissimulationhavesmalleras-
5.2 Dependence on viscosity
pectratios,χ,andshorterlifetimes.Closeinspectionofthe
In this section we investigate the role of viscosity in the toppanelsinFig.8andofsequencesofsnapshotssimilarto
development of the instability. We have performed several thesepanelsindicatesthatχ∼1.5forthisrun,withvortex
simulationsusingthesamediscmodeldescribedinthepre- lifetimes being approximately one orbit.
vious section, but adding viscous stresses and varying the Fig.9showtheresultsofthesimulationswithτ =0.5.
kinematic viscosity, ν. Fig. 6 shows the growth of the in- Inspectionofthespectrumbarthatindicatestheamplitude
stability for different values of ν. As expected, it shows of the vorticity perturbations shows that this run produces
that the viscosity has a stabilizing effect on the VSI. The significantly weaker vorticity perturbations than the runs
growth rate increases when the viscosity decreases until withτ =0.1and0.05.Consequently,thenatureoftheflowis
ν =2.5×10−8,forwhichthebehaviourisclosetotheinvis- verydifferentinthiscase,consistingofelongatedandlong-
cidcase.Alargeviscositytotallyinhibitsthedevelopmentof lived vortices. We estimate that the typical vortex aspect
the VSI. The vorticity contours for the different viscosities ratio in this run is χ∼6 and the average life time exceeds
areplottedinFig.7aftersimulationruntimesof26orbits. thesimulationruntimes.Closeinspectionofthevorticesin-
For ν = 2.5×10−6 the flow is still axisymmetric and the dicates that they maintain a turbulent core throughout the
vortices have not yet formed, while for the lower viscosity evolution, presumably due to the elliptical instability oper-
(cid:13)c RAS,MNRAS000,1–??
8 S.Richard, R.P.Nelson & O.M.Umurhan
Figure 7.Vorticityinthemidplaneafter26orbitsandforν=2.5 × 10−6,2.5 × 10−7,2.5 × 10−8 and0.
Figure 8.Vorticitycontoursinthemidplaneandmeridionalplaneafter33,41and74orbitalperiods,withh=0.2andτ =0.05.
ating in this case. The long-lived nature of these vortices to the formation of small aspect ratio (χ(cid:54)2) vortices that
suggests that their survival is due to the action of the SBI, form through the RWI and which extend only a small dis-
thatcontinuouslyattemptstoincreasetheamplitudeofthe tanceinheight(approximately10%ofthelocalscaleheight).
vortices in opposition to the elliptical instability which at- Thelifetimesofthesevorticesisfoundtobeveryshort,be-
temptstodestroythem.Note,however,thatthisstatement ing of the order of a few orbital periods. A longer cooling
is somewhat speculative since we may not be able to fully time of τ = 0.5 orbits gives rise to lower amplitude vortic-
resolve the elliptical instability in this case. ityperturbations,leadingtotheformationofelongated(i.e.
Examining the vertical structure of the vortices in this χ ∼ 6) vortices that extend in height by approximately 2
case, we note that they occupy a greater height in the disc scaleheights and have lifetimes that exceed the simulation
than those observed in the runs withτ =0.1 and 0.05. The run times. These vortices are observed to have turbulent
vortices shown in the upper panels of Fig. 9 extend above cores, presumably due to the elliptical instability, suggest-
and below the midplane by approximately one scaleheight. ing that their long lifetimes are due to the action of the
Insummary,wefindthattheVSIgivesrisetorelatively SBI maintaining the integrity of these structures. Table 1
large amplitude axisymmetric vorticity perturbations when summarize the growth rate of the instability and the vor-
thecoolingtimeisshort(i.e.τ (cid:54)0.1orbits),andthisleads
(cid:13)c RAS,MNRAS000,1–??
Vertical shear instability and vortices 9
Figure 9.Vorticityprofileinthemidplaneandmeridionalplaneafter201,317and401orbitalperiods,withh=0.2andτ =0.5.
criticalvalueforwhichtheVSIswitchesoff.Thus,wemight
Coolingtime Growthrate Vorticity Vortexaspectratio
reasonably expect behaviour that is similar to that shown
(orbits) (orbits)−1 ωz/Ω
by the run with h = 0.2 and τ = 0.5 in this case. We also
0.5 0.31 -0.3 6 presentasimulationwithh=0.05andτ =0.01.Theresults
0.1 0.9 -1.1 2 presentedinSection4indicatethatthismodelshouldexpe-
0.05 1.7 -3 1.5 rience rapid growth of the VSI as the cooling time is much
shorter than the critical thermal relaxation time scale. We
Table 1. Resultsfortheh=0.2,q=−2,p=−1.5simulations. mightreasonablyexpectresultssimilartothoseobservedfor
First column shows cooling time, second column shows growth
the runs with h=0.2 and τ =0.1 and 0.05 in this case. In
rate of the VSI, third column shows peak vorticity perturbation
order to resolve the smaller length scale features expected
(inunitsofthelocalorbitalangularvelocity),andfourthcolumn
in this run the simulation was conducted with double the
showsmeanaspectratiooftheemergingvortices.
resolution in the radial direction.
The results of the simulation with h=0.1 and τ =0.2
ticity and the aspect ratio of the resulting vortices for each are shown in Fig. 10. As expected, the growth time of the
simulations. instability is long in this case and the vorticity perturba-
tions that arise are of relatively low amplitude. The vor-
tices that form are observed to have aspect ratios χ ∼ 10,
5.4 Dependence on scaleheight
smooth non-turbulent cores, and lifetimes that exceed the
In the previous simulations the scaleheight was set to h = simulation run time of 695 orbits. The vertical heights of
0.2, which corresponds to a rather thick disc that may be the vortices formed near the disc midplane again appear to
unrealisticforatypicalprotoplanetarydisc.Here,weinves- extendapproximatelyplus-and-minusonescaleheightabout
tigate the evolution of two thinner disc models. As shown themidplane,asobservedforthelong-livedvorticesformed
inSection4,forthinnerdiscsweneedshortercoolingtimes in the run with h=0.2 and τ =0.5.
to trigger the VSI. We present a simulation with h = 0.1 The results for the simulation with h = 0.05 and τ =
andτ =0.2.InspectionoftheresultsdescribedinSection4 0.01areshowninFig.11.TheshortgrowthtimeoftheVSI
indicates that the growth time of the VSI in this case will leads to the formation of vortices with small aspect ratios
be quite long as the cooling time is relatively close to the χ∼1.5 with lifetimes of the order of one orbital period.
(cid:13)c RAS,MNRAS000,1–??
10 S.Richard, R.P.Nelson & O.M.Umurhan
Figure 10.Vorticityprofileinthemidplaneandmeridionalplaneafter331,555and694orbitalperiods,withh=0.1andτ =0.2.
Figure 11.Vorticityprofileinthemidplaneandmeridionalplaneafter100,120and130orbitalperiods,withh=0.05andτ =0.01.
(cid:13)c RAS,MNRAS000,1–??