Gap formation and stability in non-isothermal protoplanetary discs

6.2 MB·English

by  Robert Les

Mon.Not.R.Astron.Soc.000,000–000 (0000) Printed1April2015 (MNLATEXstylefilev2.2)
Gap formation and stability in non-isothermal
protoplanetary discs
5 Robert Les 1⋆ and Min-Kai Lin 1,2 †
1
1Canadian Institute for Theoretical Astrophysics, 60 St. George Street, Toronto, ON, M5S 3H8, Canada
0
2Department of Astronomy and Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA
2

r
a
M 1April2015

0
3 ABSTRACT
Several observations of transition discs show lopsided dust-distributions. A potential

] explanationistheformationofalarge-scalevortexactingasadust-trapattheedgeof
P
agapopenedbyagiantplanet.Numericalmodelsofgap-edgevorticeshavesofarem-
E
ployed locally isothermal discs in which the temperature profile is held fixed, but the
. theoryofthis vortex-formingor‘Rossbywave’instabilitywasoriginallydevelopedfor
h
p adiabatic discs. We generalize the study of planetary gap stability to non-isothermal
- discs using customized numerical simulations of disc-planet systems where the planet
o
opens an unstable gap. We include in the energy equation a simple cooling function
r −1
t with cooling timescale tc =βΩk , where Ωk is the Keplerian frequency, and examine
s the effect of β on the stability of gap edges and vortex lifetimes. We find increas-
a
[ ing β lowers the growth rate of non-axisymmetric perturbations, and the dominant
azimuthal wavenumber m decreases. We find a quasi-steady state consisting of one

2 large-scale,over-dense vortex circulating the outer gap edge, typically lasting O(103)
v orbits. We find vortex lifetimes generally increase with the cooling timescale tc up
9
to an optimal value of tc ∼ 10 orbits, beyond which vortex lifetimes decrease. This
7
non-monotonicdependenceisqualitativelyconsistentwithrecentstudiesusingstrictly
9
isothermaldiscs thatvarythe discaspectratio.Thelifetime andobservabilityofgap-
1
0 edge vortices in protoplanetary discs is therefore dependent on disc thermodynamics.
.
1 Key words: accretion, accretion discs, protoplanetary discs, hydrodynamics, insta-
0 bilities, planet-disc interactions, methods: numerical
5
1
:
v
i
X 1 INTRODUCTION Osorio et al. 2014), with possible evidence of companions
withinthem(e.g.Quanzet al.2013a;Reggiani et al.2014).
r The interaction between planets and protoplanetary discs
a
plays an important role in the theory of planet formation A recent theoretical development in the study of plan-
and disc evolution. Disc-planet interaction may lead to the etary gaps is their stability. When the disc viscosity is low
orbital migration of protoplanets and modify the structure and/or the planet mass is large, the presence of potential
of protoplanetary discs (see Baruteau & Masset 2013, for a vorticity (PV, the ratio of vorticity to surface density) ex-
recent review). trema can render planetary gaps dynamically unstable due
A sufficiently massive planet can open a gap to what is now referred to as the ‘Rossby wave instability’
in a gaseous protoplanetary disc (Papaloizou & Lin (RWI, Lovelace et al. 1999; Li et al. 2000). This eventually
1984; Bryden et al. 1999; Crida, Morbidelli & Masset 2006; leads to vortex formation (Li et al. 2001; Koller, Li & Lin
Fung,Shi & Chiang2014),whilelowmassplanetsmayalso 2003; Li et al. 2005; de Val-Borro et al. 2007), which can
opengapsifthediscviscosityissmallenough(Li et al.2009; significantlyaffectorbitalmigrationoftheplanet(Ou et al.
Dong, Rafikov& Stone 2011; Duffell & MacFadyen 2013). 2007;Li et al.2009;Yu et al.2010;Lin & Papaloizou2010).
Support for such disc-planet interaction have begun to
emerge in observations of circumstellar discs that reveal Vortex formation at gap edges may also have observ-
annular gaps (e.g. Quanzet al. 2013b; Debes et al. 2013; able consequences. Because disc vortices represent pres-
sure maxima, they are able to collect dust particles
(Barge & Sommeria1995;Inaba& Barge 2006;Lyra & Lin
2013). Dust-trapping at gap-edge vortices have thus been
⋆ [email protected] suggested to explain asymmetric dust distributions ob-
† [email protected] served in several transition discs (e.g. Casassus et al. 2013;
(cid:13)c 0000RAS
2 Les and Lin
van derMarel et al.2013;Isella et al.2013;Fukagawa et al. whereΣisthesurfacedensity,v=(v ,v )thefluidvelocity,
r φ
2013;P´erez et al. 2014; Pinilla et al. 2015). p is the vertically-integrated pressure, e = p/(γ−1) is the
However, studies of Rossby vortices at planetary gap- energy per unit area and the adiabatic index γ = 1.4 is
edges have adopted locally isothermal discs, where the assumed constant.
disc temperature is a fixed function of position only (e.g. The total potential Φ includes the stellar potential,
Lyraet al. 2009; Lin & Papaloizou 2011; Zhu et al. 2014; planetpotential(describedbelow)andindirectpotentialsto
Fu et al. 2014). On the other hand, the theory of the RWI accountforthenon-inertialreference frame.Inthemomen-
was in fact developed for adiabatic discs (Li et al. 2000), tum equations, f represent viscous forces, which includes
ν
whichpermitsheating.Inadiabaticdiscs,therelevantquan- artificial bulk viscosity to handle shocks, and a Navier-
tity for stability becomes a generalization of the PV that Stokesviscositywhosemagnitudeischaracterizedbyacon-
accounts for entropy variations (Lovelace et al. 1999). stant kinematic viscosity parameter ν. However, we will be
Gap-opening is associated with planet-induced spiral considering effectively inviscid discs by adopting small val-
shocks. In an isothermal disc, PV-generation across these uesof ν.
isothermal shocks leads to theRWI (Koller, Li& Lin 2003;
Li et al. 2005; de Val-Borro et al. 2007; Lin & Papaloizou
2010).However,ifcoolingisinefficientandtheshockisnon- 2.1 Heating and cooling
isothermal,thenshock-heatingmayaffectgapstability,since
In theenergy equation,the heating term H is definedas
the relevant quantity is an entropy-modified PV (described
below), and there is entropy-generation across the shocks. H≡Q+−Q+ Σ, (4)
For example, previous linear simulations of the RWI i Σi
found that increasing the sound-speed favours instability where Q+ represents viscous heating (from both physical
(Li et al. 2000;Lin 2013).In thecontext of planetary gaps,
and artificial viscosity) and subscript i denotes evaluation
however,theincreasedtemperaturemayalsoacttostabilize
at t=0. The cooling term C is definedas
the disc by making gap-opening more difficult. It is there-
foreoftheoreticalinteresttoclarifytheeffectofheatingand 1 Σ
C ≡ e−e , (5)
cooling on thestability of planetary gaps. tc (cid:18) iΣi(cid:19)
Inthiswork,weextendthestudyofplanetarygapsta-
where t = βΩ−1 is the cooling time, Ω = GM/r3 is
bility against vortex formation to non-isothermal discs. We c k k
the Keplerian frequency and β is an input parameter. This
includein thefluid energy equation an one-parametercool- p
cooling prescription allows one to explore the full range of
ingprescriptionthatallows ustoprobediscthermodynam-
thermodynamic response of the disc in a systematic way:
ics ranging from nearly isothermal to nearly adiabatic.
β≪1isalocallyisothermaldiscwhileβ ≫1isanadiabatic
Thispaperisorganizedasfollows.In§2wedescribethe
disc.
equations governing the disc-planet system and initial con-
Note that the energy source terms have been chosen
ditions. Ournumerical approach, includingdiagnostic mea-
to be absent at t = 0, allowing the disc to be initialized
sures, are given in §3. We present results from two sets of
closetosteadystate.TheC functionattemptstorestorethe
numerical experiments. In §4 we use disc-planet interaction
initialenergydensity(andthereforetemperature)profile.In
to set up discs with gaps, but study their stability without
practice,thisisacoolingtermatthegapedgebecausedisc-
further influence from the planet. We then perform long-
planet interaction leads to heating.
termdisc-planetsimulationstoexaminethelifetimeofgap-
edge vortices in §5, as a function of the imposed cooling
rate.Weconcludeandsummarizein§6withadiscussion of
2.2 Disc model and initial condition
important caveats.
The disc occupies r ∈ [rin,rout] and φ ∈ [0,2π]. The initial
disc is axisymmetric with surface density profile
2 DISC-PLANET MODELS
r −s rin
Tgahsedsiysscteomrbiistinagtwao-cdeinmteranlsiostnaarl (o2fDm)ansosn-Mse∗lf.-gWraevitaadtoinpgt Σ(r)=Σref(cid:18)rin(cid:19) (cid:20)1−rr+Hi(rin)(cid:21) (6)
cylindrical co-ordinates (r,φ,z) centred on the star. The where the power-law index s = 2, H(r) = cisoΩk defines
frame is non-rotating. Computational units are such that the disc scale-height where ciso = p/Σ is the isothermal
G = M∗ = R = µ = 1 where G is the gravitational con- sound-speed. The disc aspect ratiopis defined as h ≡ H/r
stant, R is the gas constant and µ is the mean molecular and initially h = 0.05. For a non-self-gravitating disc, the
weight. surface density scale Σref is arbitrary.
The disc evolution is governed by the standard fluid The initial azimuthal velocity vφi is set by centrifugal
equations balance with pressure forces and stellar gravity. For a thin
disc, v ≃ rΩ . The initial radial velocity is v = 3ν/r,
φ k r
∂Σ +∇·(Σv)=0, (1) where ν = νˆri2nΩk(rin), and we adopt νˆ = 10−9, so that
∂t |v /v | ≪ 1 and the initial flow is effectively only in the
r φ
∂v +v·∇v=−1∇p−∇Φ+f , (2) azimuthal direction. With this value of physical viscosity,
∂t ρ ν the only source of heating is through compression, shock-
∂e heating (via artificial viscosity) and the C function when
+∇·(ev)=−p∇·v+H−C, (3)
∂t e/Σ<e /Σ .
i i
(cid:13)c 0000RAS,MNRAS000,000–000
Gaps in non-isothermal discs 3
2.3 Planet potential where κ2 = r−3∂ (r4Ω2) is the square of the epicyclic fre-
r
quency,Ω=v /r istheangularspeed,andS ≡p/Σγ isthe
The planet potential is given by φ
entropy.Thefirstfactoristheusualpotentialvorticity(PV,
Φ =− GMp , (7) or vortensity).
p |r−r |2+ǫ2 The generalised PV appears in the description of
p p
the linear stability of radially-structured adiabatic discs
where Mpp is the planet mass and we fix q ≡ Mp/M∗ = (Lovelace et al. 1999; Li et al. 2000), where the authors
10−3 throughout this work. This corresponds to a Jupiter-
showanextremuminη˜mayleadtoadynamicalinstability,
mass planet if M∗ =M⊙. The planet’s position in the disc the RWI. In a barotropic disc where p=p(Σ), the entropy
r = (r ,φ ) and ǫ = 0.5r is a softening length with
p p p p h factor is absent and theimportant quantityis thePV.
r =(q/3)1/3r beingtheHillradius. Theplanetis heldon
h p
afixedcircularorbitwithrp =10rinandφp =Ωk(rp)t.This
also defines the time unit P0 ≡2π/Ωk(rp) used to describe
3.2.2 Fourier modes
results.
The RWI is characterized by exponentially growing pertur-
bations. Though in this paper we do not consider a formal
3 NUMERICAL EXPERIMENTS linear instability calculation, modal analysis will be useful
to analyse the growth of perturbations with different az-
Thedisc-planetsystemisevolvedusingtheFARGO-ADSGcode
imuthalwavenumbers,which isassociated with thenumber
(Baruteau & Masset 2008b,a).This is amodified version of
of vortices initially formed by theRWI.
theoriginalFARGOcode(Masset2000)toincludetheenergy
The Fourier transform of the time-dependent surface
equation. The code employs a finite-difference scheme sim-
density is
ilar to the ZEUS code (Stone& Norman 1992), but with a
modified azimuthal transport algorithm to circumvent the Σ (r,t)= 2πΣ(r,φ,t)e−imφdφ (12)
time-step restriction set by the fast rotation speed at the m
Z0
innerdiscboundary.Thediscisdividedinto(N ,N )zones
r φ
where m is the azimuthal wave number. We define the
intheradialandazimuthaldirections,respectively.Thegrid
growth rate σ of the mth component of the surface density
spacing is logarithmic in radius and uniform in azimuth.
through
dh|Σ |i
3.1 Cooling prescription m r =σh|Σ |i (13)
dt m r
Inthisworkweonlyvaryonecontrolparameter:thecooling
where h|·|i denotes the average of the absolute value over
r
time. The cooling parameter β is chosen indirectly through
aradialregionofinterest.ByusingEq.13thegrowthrates
theparameter β˜such that
of the unstable modes can be found from successive spatial
Fourier transforms overan appropriate period of time.
tc(rp+xs)=βΩ−k1(rp+xs)=β˜tlib(rp+xs), (8)
wherex isthedistancefromtheplanettoitsgapedge,and
s
3.2.3 Rossby number
tlib is the time interval between successive encounters of a
fluidelementatthegapedgeandtheplanet’sazimuth.That The Rossby number
is,wemeasurethecooling timeinunitsofthetimeinterval
between encounters of a fluid element at the gap edge and Ro= ˆz·∇×v−hˆz·∇×viφ, (14)
theplanet-induced shock. 2hΩiφ
Assuming Keplerian orbital frequencies and x ≪ r
s p isadimensionlessmeasureofrelativevorticity.Hereh·i de-
gives tlib≃4πrp/(3Ωkpxs), where Ωkp =Ωk(rp). Therefore notesanazimuthalaverage.ValuesofRo<0corresponφdto
4πr 3x anti-cyclonicrotation with respect tothebackground shear
β=β˜ p 1− s , (9)
3x 2r and thus can be used to identify vortices and quantify its
s (cid:18) p(cid:19)
intensity.
wherex ≪r wasusedagain.Weusex =2r inEq.9.For
s p s h
a planet mass with q = 10−3, Eq. 9 then gives β ≃ 23.9β˜.
In terms of planetary orbital periods, this is
β r 3/2 r 3/2 4 GROWTH OF NON-AXISYMMETRIC
tc(r)= 2π r P0 ≃3.8β˜ r P0. (10) MODES WITHOUT THE INFLUENCE OF
(cid:18) p(cid:19) (cid:18) p(cid:19) THE PLANET
3.2 Diagnostic measures In thissection, theplanet is introduced at t=20P0 and its
potential is switched on over 10P0. At t = 30P0 we switch
3.2.1 Generalised potential vorticity offtheplanetpotentialandazimuthallyaveragethesurface
density,energyandvelocityfields.(Atthispointtheplanet
The generalised potential vorticity is definedas
hascarvedapartialgapandtheRWIhasnotyetoccurred.)
η˜= κ2 ×S−2/γ, (11) Effectively,weinitialisethediscwithagapprofile.Wethen
2ΩΣ perturb the surface density in the outer disc (r > r ) and
p
(cid:13)c 0000RAS,MNRAS000,000–000
4 Les and Lin
continue to evolve the disc. We impose sinusoidal pertur-
bations with azimuthal wavenumbers m ∈ [1,10] and ran-
dom amplitudes within ±0.01 times the local surface den-
sity.Thisprocedureallows ustoanalysethegrowth ofnon-
axisymmetric modes associated with the gap, but without
complications from non-axisymmetry arising directly from
disc-planet interaction (i.e. planet-inducedwakes).
Note that these ‘planet-off’ simulations are not linear
stabilitycalculationsbecausethecoolingterminourenergy
equationrestorestheinitialtemperatureprofilecorrespond-
ingtoconstantH/r=0.05,ratherthantheheatedgapedge.
However, we will examine a nearly adiabatic simulation in
§4.3.1, which is closer to a properlinear problem.
Simulations here employ a resolution of (N ,N ) =
r φ
(1024,2048) with open boundaries at r = rin and rout =
25rin. We compare cases with β˜ = 0.1,1,10 corresponding
to fast, moderately, and slowly cooled discs.
4.1 Gap structure
We first compare the gap structures formed by planet-
disc interaction as a function of the cooling time. The
azimuthally-averaged gap profiles are shown in Fig. 1 for
differentvaluesofβ˜.Gapsformedwithlowerβ˜(fastercool-
ing) are deeper with steeper gradients at the gap edges.
Fastercoolingratesalsoincrease(decrease)thesurfaceden-
sity maxima(minima). However,a clean gap does not form
in thisshort time period.
Increasingβ˜leadstohigherdiscaspectratiosh=H/r,
i.e. higher temperatures. Heating mostly occurs at the gap
edges due to planet-induced spiral shocks. Increasing the
cooling timescale implies that this heat is retained in the
disc. In the inviscid limit the gap opening condition is
r > H or q > 3h3 (Crida, Morbidelli & Masset 2006),
h ∼ ∼
which indicates that for hotter discs (higher h), it becomes
more difficult for a planet of fixed q to open a gap. This
explains theshallower gaps in surface density when β˜is in-
creased.
Theimportantconsequenceofaheatedgapedgeisthat
thegeneralisedvortensityprofiles,η˜,becomessmootherwith
increasing cooling times, with the extrema becoming less
pronounced.Previouslocally isothermal disc-planetsimula-
tions show the RWI associated with PV minima (Li et al.
2005; Lin & Papaloizou 2010). Wecan therefore expect the
RWItobeassociatedwithminimainthegeneralisedvorten- Figure1.Azimuthallyaveragedgapprofilesatt=30P0 forthe
sity (corresponding to local surface density maxima) in the initialpartialgapopenedbeforeinstabilityemergesforfast(solid,
non-isothermalcase.Becausetheextremaarelesssharp,the β˜ = 0.1), moderate (dashed, β˜= 1), and slow cooling (dashed-
dot, β˜ = 10). The relative surface density perturbation (top),
RWIisexpectedtobeweakerandthegaptobemorestable
discaspectratio(middle)andgeneralisedvortensityperturbation
with longer cooling times.
(bottom) areshown.
However, we remark that the change in the gap struc-
turebecomeslesssignificantatlongcooling times,asFig.1
shows that the profiles with β˜= 1 and β˜= 10 are similar. to axisymmetric instabilities. The generalised local axisym-
ThisimpliesthattheeffectofcoolingtimescaleontheRWI metric stability condition is theSolberg-Hoiland criterion,
throughthesetupofthegapprofile,becomeslessimportant
for large β˜. κ2+N2 >0 (15)
where
1 ∂P 1 ∂Σ 1 ∂P
N2 = − (16)
4.2 Axisymmetric stability Σ ∂r Σ ∂r γP ∂r
(cid:18) (cid:19)
The initial planet-disc interaction form bumps and grooves is the square of the Brunt-V¨ais¨al¨a frequency. At the outer
in the gap profiles which can potentially be unstable due gap edge r = r +2.5r , where the RWI is excited (see
p h
(cid:13)c 0000RAS,MNRAS000,000–000
Gaps in non-isothermal discs 5
Table 1.Dominantmodeandgrowthratesforβ˜=0.1,1.0,10.0
(fast,moderate,andslowcooling)valuesduring‘planet-off’sim-
ulations
β˜=0.1 β˜=1.0 β˜=10.0
m 102σ/Ω(rp) m 102σ/Ω(rp) m 102σ/Ω(rp)
6 7.3 3 2.0 1 1.1
7 7.8 4 2.2 2 1.6
8 7.9 5 2.3 3 1.7
9 7.9 6 1.6 4 1.2
10 6.8 7 1.1 5 0.1
Figure3.Generalisedvortensityperturbation(relativetot=0)
forcasesofβ˜=0.1,1,10(left,middle,right)duringthegrowthof
non-axisymmetricmodes.Theplanetpotentialhasbeenswitched
off. The number of vortices decrease as β˜ increases. Note that
snapshotsaretakenlaterforincreasingβ˜becauseittakeslonger
forthevorticestogrowandbecomevisiblewithincreasingcooling
time.
Figure 2. Evolution of azimuthal Fourier modes of discsurface
density, non-dimensioqnalized by the initial axisymmetric com- 0.079 → 0.017. However, despite two orders of magnitude
ponent Σ0(t = 0), for the ‘planet-off’ simulation with β˜ = 10. increaseinthecooling time,theinstability remainsdynam-
Colourscorrespondtodifferentmvalues.Them=3component ical with characteristic growth time <∼ 10P0. Snapshots of
is the fastest growing mode during linear growth with a growth theinstability in for thedifferent β˜are shown in Fig 3.
rateofγ=0.017Ω(rp).
These ‘planet-off’ simulations show that gap edges be-
comemorestablewithlongercoolingtimes.Thisisexpected
below),wefindκ2+N2 reacheslocalminimumwithavalue because larger β˜ results in hotter gap profiles at t = 30P0
∼0.45Ω2(r ) for all β˜.The Brunt-V¨ais¨al¨a frequencyat the with less pronounced generalised vortensity minima. Stabi-
p
outergapedgeisN ∼0.1Ω(r ),decreasingmarginallywith lization with increased cooling time is therefore due to a
p
longercoolingrate.TheSolberg-Hoilandcriteriaissimilarly smoother basic state for the instability, as it is more diffi-
satisfied for the entire 2D disc throughout the simulations. cult for the planet to open a gap if the disc is allowed to
Thus for all values of β˜ the planet-induced gaps are stable heat up.
to axisymmertic perturbations. Forcompletenesswealsosimulatedalocallyisothermal
disc with h = 0.05 where the sound-speed is kept strictly
equal to its initial value. This simulation yield a most un-
4.3 Non-axisymmetric instability stablegrowthrateσ/Ω(rp)≃0.05atm=5,comparedwith
a value of 0.085 at m = 6 for a corresponding simulation
We now examine the evolution of the gap for t > 30P0, that includes the energy equation but with rapid cooling
with the planet potential switched off, but with an added β˜=0.01. However, thevortex evolution is similar.
surface density perturbation. For all three cooling times
β˜ = 0.1, 1.0, 10, we observe exponential growth of non-
axisymmetric structures. Anexample is shown in Fig. 2 for
4.3.1 Nearly adiabatic discs
β˜ = 10. We characterize these modes with an azimuthal
wavenumbermandgrowthrateσ(m)asdefinedbyEq.12— The above ‘planet-off’ simulations are not formally linear
13. Mode amplitudes were averaged over r−r ∈ [2,5]r . stability calculations, because the cooling time is compara-
p h
Table1liststhegrowthratesmeasuredduringlineargrowth ble or shorter than the instability growth time, tc <∼ γ−1.
for5valuesofmcentredaroundthatwithmaximumgrowth Thus the disc cools back to its initial temperature corre-
rate. spondingtoh=0.05beforeorduringtheinstabilitygrowth,
Table 1 show that as β˜ is increased from 0.1→10 the so we do not have a steady basic state to formulate a stan-
dominant azimuthal Fourier modedecreases from m=9→ dard linear stability problem.
3 and the respective growth rate decreases from γ/Ω(r )= In order to capture the effect of a heated gap edge, we
p
(cid:13)c 0000RAS,MNRAS000,000–000
6 Les and Lin
ran a simulation with β˜=100, corresponding to an almost
adiabatic disc. In this simulation the cooling rate is slow
enough that the gap temperature profile (e.g. middle panel
ofFig.1)changesonlymarginallyovertheinstabilitygrowth
timescale.
Wefindverysimilargapprofilesandmodegrowthrates
for β˜ = 100 as with β˜ = 10. At t = 30P0, the disc only
heats up to values h ≃ 0.06 in the nearly adiabatic case.
This is close to the original temperature of h = 0.05, so
linear growth rates are not expected tochange significantly
(Li et al. 2000).
AccordingtoLi et al.(2000),increasinghincreaseslin-
ear growth rates of the RWI because it is pressure-driven.
However,inthecaseof disc-planetinteraction,increasing h
hasastabilizingeffectthroughthesettingupthegapprofile
because it results in smoother gap edges. The fact that we
observesmallergrowthratesashisincreasedindicatesthat Figure 4. Long term simulations without the planet potential
for planetary gaps, the importance of h on the linear RWI after the gap is set up. The m = 1 surface density component,
is through setting upthe gap profile, i.e. basic state for the non-dimensionalizedbytheinitialm=0component,attheouter
instability (as opposed to thelinear response). gapedgeisshownasafunctionofthecoolingtimescale.
4.4 Long term evolution
Wealsoextendedthese‘planet-off’simulationsintothenon-
linear regime. Afterthelineargrowth phaseof thevortices,
vortexmergingtakesholdontimescalesofupto150P0,until
there is one vortex left. We find the vortex merging time is
dependentonthegrowth ratesofthemodesandsaturation
timescales,withtheslowestgrowingmodesinβ˜=10taking
thelongest to merge.
Fig.4showsevolutionofthem=1surfacedensitycom-
ponent, which represents the amplitude of the post-merger
single vortex. For completeness we also ran intermediate
cases with β˜ = 0.5 and 5.0. The amplitude of the initially
formed vortexwas foundtodecreasewith increased cooling
rate. The vortices simply decay on a timescale of O(103)
orbitswith faster decayfor strongervortices(which areob-
tained with faster cooling rates).
Thisdecayisprobablyduetonumericalviscosity.Dur- Figure 5. Evolution of the non-dimesnionalized m = 1 surface
ing the slow decay the vortex elongates (weakens) while its densitycomponentforlongtermsimulationswiththeplanetpo-
tential kept on. Notice there is vortex growth after the initial
radialwidthremainsO(H),soitssurfacedensitydecreases.
Inaddition,forβ˜=0.1, 0.5and1.0wealsoobservetheap- lineargrowth, whichcontrasts to Fig.4,wherevortices decay in
theabsenceofcontinuous disc-planetinteraction.
pearance of spiral waves associated with the vortex, which
maycontributetoitsdissipation (seebelow).Wewillseein
the next section that this decay after linear growth is very
5.1 Generic evolution
different to when the planet potential is kept on.
The linear growth of the RWI and vortex-formation is fol-
lowed by vortex merging. We now find merging timescales
independent of β˜, and by 60P only one vortex remains.
5 NON-LINEAR EVOLUTION OF GAP-EDGE o
The evolution of the amplitude of the m = 1 surface den-
VORTICES WITH FINITE COOLING TIME
sitycomponent,averagedoverr−r ∈[2,10]r ,isshownin
p h
Wenowexaminelong-termsimulationsofgap-edgevortices Fig.5fordifferentβ˜.Theinitial,post-mergervortexampli-
for β˜ = 0.1,0.5,1,5,10. (Additional cases are presented in tude is found to be weaker for longer cooling rates (which
§5.4whenexaminingvortexlifetimesasafunctionofβ˜.)The havesmaller linear growth rates).
planet potential is kept on throughout. We employ a grid In all cases the system remains in a quasi-steady state
with (Nr,Nφ) = (512,1024) in order for these simulations for >∼ 800P0 with a single vortex circulating the outer gap
to be computationally feasible. We also use a larger disc edgeatthelocalKeplerianfrequency.Fig.6showsatypical
with rout = 45rin to minimise boundary effects on vortex plotoftherelativesurfacedensityperturbationinthisstate.
evolution, and apply open boundaries at r=rin, rout. Duringthisstage,thevortexintensifies.Thisisbettershown
We comment that lower-resolution simulations with inFig.7astheevolutionofRossbynumbersmeasuredatthe
(N ,N ) = (256,512) show similar behaviour and trends vortex centres. The Rossby number increases in magnitude
r φ
as the high-resolution runsreported below. duringquasi-steady state, but themaximum |Ro| is similar
(cid:13)c 0000RAS,MNRAS000,000–000
Gaps in non-isothermal discs 7
Figure 8. Average value of the relative surface density pertur-
bation at vortex centres for various cooling times. Initial vortex
over-densitiesdecreasewithcoolingratewhilevortexgrowthrates
increasewithcoolingrate.
eration of vorticity by planet-disc interaction. This is sup-
Figure 6.Relative surfacedensity perturbation for the β˜=0.1 ported by the observation that in the previous simulations
case during quasi-steady state with a single vortex at the outer
withouttheplanet,theamplitudeofthepost-mergervortex
gap edge. The plot for other values of the cooling time β˜ are
does not grow (Fig. 4).
similar.Thedecreaseinthesurfacedensitynearr∼40rin arises
Fig.5showsthatthedurationofthequasi-steadystate
from mass loss due to the open boundary condition imposed at
rout. varies with the cooling rate: for β˜=0.1 and 10, the vortex
amplitude begins to decay around t∼800P0, while for β˜=
1, 5 the decay begins at t ∼ 1200P0. This non-monotonic
dependence suggests that there exists an optimal cooling
rate to maximise the vortex lifetime. We will discuss this
issuefurtherin§5.4.Noticethedecaytimescalecanbelong
withrapidcooling: forβ˜=0.1ittakes∼400P0 whereasfor
β˜= 10 it takes ∼ 100P0 for the m=1 amplitude to decay
significantly after reaching maximum.
5.2 Additional analysis on vortex decay
Inthis subsection we examinethevortex decay observed in
oursimulationsinmoredetail.Fig.9showsnapshotsofthe
vortexforthecaseβ˜=1.Theplotsshowthesurfacedensity
perturbationand thesurface densitygradient duringquasi-
steadystate(t=700P0),whenthem=1amplitudebegins
todecrease (t=1300P0)andjust aftertherapidamplitude
decay (t=1510P0).
Figure 7. Evolution of Rossby numbers at the centres of the Inquasi-steady(t=700P0)thevortexiselongatedwith
vortices formed in discs with different cooling rates. A negative
a vortex aspect ratio ≈4 , but becomes more compact ap-
Rossby number implies anti-cyclonic motion. Boxcar averaging
proachingaratio of2duringitsdecay(t=1310P0).Notice
wasusedtoremovecontributions fromtheplanet-inducedspiral
inFig.9theappearanceofwakesextendingfromeitherside
shock.
of the vortex at t = 1300P0. These wakes correlate with
large gradients in surface density (bottom panel), and are
for all β˜: the vortex reaches a characteristic value of Ro ≈ firstseeninthelaterhalfofthequasi-steadystate.Wefind
−0.35 for β˜=0.1 and Ro≈−0.45 for β˜=10. thetimeatwhichthevortexbeginstodecaycoincideswith
We find the vortices become significantly over-dense. theemergence of these wakes.
Fig. 8 plots the surface density perturbation measured at During the quasi-steady state the vortex orbits at r ∼
the vortex centres, showing ∆Σ/Σ0 >∼ 7 for all cases of β˜ rp+6rh. Wedonot see significant vortex migration at this
in quasi-steady state, and max(∆Σ/Σ0) ∼ 11 for β˜ = 5. stage, since the vortex is located at a surface density max-
The maximum over-density typically increases with longer imum (Paardekooper, Lesur& Papaloizou 2010). However,
coolingtimes,despitethevorticesareinitiallyweakeratfor- simultaneous with theappearance of the wakes, we observe
mation with increasing β˜. The large increase in the surface thevortex begins to migrate inwards to r∼r +5r .
p h
density is due to vortex growth as there is continuous gen- During quasi-steady state the average value of the sur-
(cid:13)c 0000RAS,MNRAS000,000–000
8 Les and Lin
Figure9.Thevortexinthecasewithβ˜=1duringquasi-steady Figure10.Machnumberrelativetothevortex,averagedovera
state(left),startofdecay(middle),andjustafterthedecayinthe regionwithin2H ofthevortexcentrewithrespecttotime.This
m=1amplitude(right).Thesurfacedensityperturbation(top) plot can be compared to the evolution of the vortex amplitude
andthe associated surfacedensitygradient (bottom) areshown. showninFig.5.
Wake-like features corresponding to large density gradients are
foundtooriginatefromthevorticesduringthelatephaseoftheir
quasi-steadystates andintodissipationtimes.Thisplotistobe
consideredinconjunctionwithFig.5.
face density gradient along the wakes is |w ∇Σ/Σ| ∼ 0.4,
s
where w ≃0.1 (code units) is a typical length scale of the
s
surface density variation across the wake. Just before the
m=1 amplitude begins to decrease, we observe this quan-
tity sharply increases to ∼ 0.6, and remains around this
value until the vortex dies out, at which point the associ-
ated Rossby number begins decreasing to zero. After the
vortex reaches small amplitudes (1 <∼ ∆Σ/Σ), it migrates
out to r∼r +6.5r .
p h
We also measured large increases in the Mach number
near the vortex as the m = 1 surface density amplitude
reaches maximum and begins to decay. Fig. 10 plots the
Mach number M = |v −vvor|/cs, where vvor corresponds
to the bulk velocity of the vortex around the disc. Values
Figure 11.Azimuthallyaveragedprofilesoftherelativesurface
in Fig. 10 have been averaged over a region within 2H of density perturbations for β˜ = 1.0 before (solid) and after vor-
the vortex centre. During the quasi-steady state the Mach
texdecay(dashed). Alongsidedrasticreductionintheoutergap
number increases steadily, and for all cases M maximizes maxima,vortexdecaysmoothoutthesharpouter gapedge.
about ∼ 100P0 after the start the m = 1 surface density
amplitude starts to decay.
Putting the above observations together, we suggest case β˜ = 1. The vortex resides around the local surface
that vortex decay (in the m = 1 surface density ampli- densitymaximum at theoutergap edge(r∼rp+6rh). We
tude) is due to shock formation by the vortex. When the see that after its amplitude has decayed (t ∼ 1500P0, Fig.
vortex reaches large amplitude, it begins to induce shocks 5),thislocalsurfacedensitymaximumisalsosmoothedout.
in the surrounding fluid, as supported by the increase in We characterize the smoothness of the outer gap edge
Mach number and the appearance of wakes with large sur- with a dimensionless gap edge gradient parameter
face density gradients. The vortex may lose energy through
shockdissipation.Inaddition,astrongvortex(orshockfor-
∂hΣ(t,r)i r
mation) can smooth out the gap structure that originally δΣ(t)= φ · (17)
gave rise to the RWI, which would oppose vortex growth. * ∂r hΣ(t=0,r)iφ+
∆r
Weexamine this below.
where∆r=r∈[r ,r +6r ]istheradialrangeofaveraging,
p p h
spanningfrom centreofthegap totheradiusofthesurface
density maximum. A larger δΣ characterizes a sharper gap
5.3 Vortex decay and gap structure
edge and larger local surface density maxima.
We find vortex decay modifies the gap structure. Fig. 11 A plot of thegradient parameter overtime for theβ˜=
shows the gap profile before and after vortex decay for the 1.0 case is shown in Fig. 12.Duringvortex decay,theouter
(cid:13)c 0000RAS,MNRAS000,000–000
Gaps in non-isothermal discs 9
Figure12.Non-dimensionalmeasureofthesurfacedensitygra-
dient at the outer gap edge β˜ = 1.0. During the vortex quasi- Figure13.Characteristictimescalesassociatedwithvortexevo-
lution,asafunctionofthecoolingparameterβ˜:timeinwhichthe
steady state the gap edge is found to have a large gradient and
sharp peak while vortex dissipation, which occurs at t ≈ 103P0 m = 1 surface density amplitude begins to decay, tdiss (black);
time at which the over-density at the vortex centre decreases to
asseeninFig.5,workstosmoothoutthegapedge.
∆Σ/Σ∼1after reaching maximum,tlife (red); andtMach isthe
time taken for the average Mach number around the vortex to
maximise(blue).
gapedgeisdrasticallysmoothedout,changingfromavalue
ofδΣ=1.2duringquasi-steadystateto0.4afterdissipation.
Thiscanbeinterpretedasthevortexprovidingaviscos- performed for comparison. In this case we did not observe
ity; and we measure a typical alpha viscosity α = O(10−2) significantvortexdecaywithinthesimulationtimescale(im-
associatedwiththevortex.Thisactsagainstgap-openingby plyingtdiss>∼2000P0). Weexpect thecorresponding vortex
theplanet,andsmoothsouttheoutersurfacedensitybump, lifetime to exceed that for β˜ = 0.01. Including an energy
sotheconditionfortheRWIbecomeslessfavourable.Inor-
equation with rapid cooling (or setting γ close to unity),
der to re-launch the RWI, the surface density bump should
could still lead to discrepancies with a locally isothermal
reform.However,thisisdifficultasthereisnomorematerial
disc. This is due to the advection-creation of specific en-
in the planet’s vicinity to clear out (Fig. 11) to form a sur-
tropy within the planet’s horseshoe region with the former
face density bump outside the gap. This may explain why
case,therebyaffectingthegeneralized vortensityandthere-
vorticesdonot reform again (at least within thesimulation
fore the instability (§4.3) and subsequent vortex evolution.
timescale).Afterfulldecaytheaspectratioattheoutergap
Nevertheless,alongervortexlifetimeinalocallyisothermal
edge is h∼0.05.
disc would be consistent with theabove trend of increasing
vortex lifetimes as β˜→0.
Intheprevioussection,weobservedthatvorticesbegan
5.4 Vortex lifetimes as a function of cooling rate
to decay when it starts to induce shocks. We thus suggest
Wenowexaminevortexlifetimesasfunctionoftheimposed that the time needed for the vortex to grow to sufficient
cooling times. For this study, additional simulations with amplitude to induce shocks in the surrounding fluid, which
β˜=0.01, 0.25, 0.75, 2.5, 7.5werealsoperformed.Wedefine maybeconsidered asthedurationofthequasi-steadystate
thevortexlifetime,tlife,asthetimeatwhichtheover-density or tdiss, is an important contribution to the overall vortex
ofthevortexreturnsto∆Σ/Σ0 ∼1after reachingmaximum lifetime.Wediscussbelowsomecompetingfactorsthatmay
(which is on the order of the initial over-density associated resultinanon-monotonicdependenceoftdiss onthecooling
with thegap formation). rate.
WeplottlifewithrespecttocoolingtimesinFig.13.We
alsoplottdiss:thetimeelapsedbeforethevortextobeginsto
5.4.1 Factors that lengthen vortex lifetimes
dissipate(whenthem=1surfacedensityamplitudebegins
to decay); and tMach: the time taken for the average Mach IthasbeenshownthattheamplitudeatwhichtheRWIsat-
numberaround thevortex to maximise. urates increases with the growth rate of the linear instabil-
For fast cooling rates (β˜ <∼ 1), the vortex lifetime is ity(Meheut, Lovelace & Lai 2013).Our‘planet-off’simula-
maximizedforβ˜→0:wefindtlife≈2100P0forβ˜=0.01and tionsyieldslowergrowthrateswithincreasingcoolingtimes,
decreases to tlife ≈ 950P0 for β˜ = 0.25. Note that for very which suggest weaker vortices are formed initially with in-
smallβ˜,thereissignificantcontributiontotheoverallvortex creasing β˜. This is consistent with the present simulations:
lifetime due to a long decay timescale. For longer cooling at the beginning of the quasi-steady state (t ∼ 100P0) we
times (β˜ >∼ 1) the vortex lifetimes maximizes at β˜ = 2.5 find the over-density at the vortex centre is ∆Σ/Σ0 = 2.5
with tlife≈1650P0. for β˜=0.1 and ∆Σ/Σ0 =1.48 for β˜=10.
Wecomment herethatalocally isothermal simulation, Thegrowthofthepost-mergersinglevortexismediated
wherethediscsound-speediskeptconstantintime,wasalso by disc-planet interaction. However, gap-opening becomes
(cid:13)c 0000RAS,MNRAS000,000–000
10 Les and Lin
more difficult in a hotter disc, and we find the generalised Inthesecondsetofcalculations,weincludedtheplanet
vortensityprofiles aresmoother with increasing β˜.This op- potential throughout the simulations and examined the
posestheRWI.Furthermore,thevortexshouldreachlarger long-term evolution of the gap-edge vortex that develops
amplitudes to induce shocks on account of the increased fromtheRWI.Thevortexreachesaquasi-steadystatelast-
sound-speed. ing O(103) orbits. Unlike the ‘planet-off’ simulations, in
These considerations suggest, with increased cooling which vortices decay after linear growth and merging, we
times,ittakeslongerforthepost-mergervortexgrowtosuf- find that with the planet potential kept on, the vortex am-
ficientamplitudetoinduceshocksanddissipate.Thisfactor plitude grows during this quasi-steady state, during which
contributestoalongerquasi-steadystatewithincreasingβ˜. no vortex migration is observed, until it begins to induce
shocks, after which thevortex amplitude begins to decay.
For our main simulations with β˜ > 0.1, the duration
5.4.2 Factors that shorten vortex lifetimes of the quasi-steady state increases with increasing cooling
timescales until a critical value, beyond which this quasi-
Notice in Fig. 5 and Fig. 8, the vortex growth during the
quasi-steady state is actually faster for β˜ = 10 than for steady state shortens again. We find the timescale for the
βh˜a=sa5l.aFrgoerreaxmamplpitlue,deatthta∼nf5o0r0βP˜0=t5h.eTvhoirsteixs awlsiothreβ˜fle=ct1ed0 vloonrgtefxortosmdeacllayβ˜,awftherichreaccohnitnrgibumtaexsitmouamloanmgpolviteurdaellcvaonrtbexe
lifetimewithrapidcooling.Wedoobservevortexmigration
in Fig. 10, where the Mach number reaches its maximum
valuesooner for β˜=10 than for β˜=5. duringits decay,which may influencethisdecay timescale.
We suggest a non-monotonic dependence of the quasi-
This observation is consistent with the RWI being
steadystateonthecoolingtimescaleβ˜canbeattributedto
favoured by higher temperatures (Liet al. 2000; Lin 2012)
the time required for the vortex to grow to sufficient am-
through theperturbations (as opposed to its effect through
plitude to induce shocks in the surrounding fluid, thereby
the set up of the gap profile discussed previously), which
losing energy and also smooth out the gap edge.
corresponds to longer cooling times in our case. While our
Forshortcoolingtimescales,theplanetisabletoopena
‘planet-off’ simulations indicate this is unimportant for the
deepergap which favourstheRWI,leading tostronger vor-
linearinstability,itmayhavecontributedsignificantlytothe
tices. For long cooling timescales, we find the vortex grows
vortexgrowthduringquasi-steadystateatverylongcooling
times (e.g. β˜=10). This effect shortens the vortex lifetime fasterduringthequasi-steadystate.Inaccordancewithpre-
vious stability calculations (Li et al. 2000), we suggest the
by allowing it to grow faster and induceshocks sooner.
latterisduetotheRWIbeingfavouredwithincreasingdisc
temperature, and that this effect overcomes weaker gap-
openingforsufficientlylongcoolingtimes.Thesecompeting
6 SUMMARY AND DISCUSSION
factorsimplyforbothshortandlongcoolingtimescales,the
In this paper, we have carried out numerical simulations vortex reaches its maximum amplitude, shock, and begins
of non-isothermal disc-planet interaction. Our simulations to decay,sooner than intermediate cooling timescales.
were customized to examine the effect of a finite cooling (However, for very rapid cooling, e.g. β˜ = 0.01, the
time on the stability of gaps opened by giant planets to quasi-steady state is also quitelong. This suggests that the
the so-called vortex or Rossby wave instability. To do so, aboveeffectsthemselvesdonothaveasimpledependenceon
we included an energy equation with a cooling term that the cooling timescale when considering β˜ → 0 and/or that
restores thedisctemperaturetoitsinitial profileon achar- otherfactorsbecomeimportantinthislimit.Thisshouldbe
acteristic timescale t . We studied the evolution of the gap investigated in future works.)
c
stability as a function of t . This is a natural extension to We remark that a non-monotonic dependence of the
c
previous studies of on gap stability, which employ locally vortex lifetime was also reported by Fu et al. (2014), who
or strictly isothermal equations of state. Weconsidered the performed locally isothermal disc-planet simulations with
inviscidlimitwhichfavorstheRWI(Li et al.2009;Fu et al. differentvaluesofthediscaspect ratio.Intheirsimulations
2014) and avoids complications from viscous heating other theoptimum aspect ratio is h=0.06. In our simulations, h
thatshockheating.However,thismeansthatthevortexlife- is a dynamical variable, but by analyzing the region where
times observed in our simulations are likely longer than in thevortexislocated(r−rp∈[2,10]rh),wefindforadimen-
realistic discs with non-zerophysical viscosity. sionless cooing timescale of β˜= 2.5, which has the longest
Weconsidered twotypesofnumericalexperiments.We vortex lifetime in the presence of moderate cooling, that
first used disc-planet interaction to self-consistently set up h≈0.058onaverage.OurresultisconsistentwithFu et al.
gapprofiles,whichwerethenperturbedandevolvedwithout (2014).
furthertheinfluenceoftheplanetpotential.Thisprocedure
isolates the effect of cooling on gap stability through the
6.1 Caveats and outlooks
set up of the initial gap profile. We find that as the cool-
ing time t is increased, the gaps became more stable, with There are several outstanding issues that needs to be ad-
c
lowergrowthratesofnon-axisymmetricmodesandthedom- dressed in future work:
inant azimuthal wavenumber also decreases. This is consis- Convergence. Although lower resolution simulations
tent with thenotion that increasing t leads to higher tem- performed in the early stages of this project gave similar
c
peratures or equivalently the disc aspect ratio h, which op- results (most importantly, the non-monotonic dependence
poses gap-opening by theplanet. This means that the gaps ofvortexlifetimesonthecoolingtimescale),wedidfindthe
opened by the planet in a disc with longer t are smoother lower resolution typically yield longer vortex lifetimes than
c
and therefore more stable to theRWI. that reported in this paper. This could be due to weaker
(cid:13)c 0000RAS,MNRAS000,000–000