Table Of Content

DRAFTVERSIONJANUARY23,2013 PreprinttypesetusingLATEXstyleemulateapjv.5/2/11 APARAMETERSTUDYFORBAROCLINICVORTEXAMPLIFICATION NATALIERAETTIG1,2 WLADIMIRLYRA2,3,4,ANDHUBERTKLAHR1 DraftversionJanuary23,2013 ABSTRACT Recent studies have shown that baroclinic vortex amplification is strongly dependent on certain factors,namely,theglobalentropygradient,theefficiencyofthermaldiffusionand/orrelaxationas wellasnumericalresolution. Weconductacomprehensivestudyofabroadrangeandcombination ofvariousentropygradients,thermaldiffusionandthermalrelaxationtime-scalesvialocalshearing 3 sheet simulations covering the parameter space relevant for protoplanetary disks. We measure the 1 Reynolds stresses as a function of our control parameters and see that there is angular momentum 0 transport even for entropy gradients as low as β=−dlns/dlnr=1/2. Values we expect in proto- 2 planetarydisksarebetweenβ=0.5−2.0Theamplification-rateoftheperturbations,Γ,appearstobe n proportionalto β2 andthusproportionaltothesquareoftheBrunt-Väisäläfrequency(Γ∝β2∝N2). a The saturation level of Reynolds stresses on the other hand seems to be proportional to β1/2. This J highlightstheimportanceofbarocliniceffectsevenforthelowentropygradientsexpectedinproto- 2 planetarydisks. 2 Subjectheadings:accretion,accretiondisks,circumstellarmatter,hydrodynamics,instabilities,turbulence,methods: numerical,solarsystem: formation,planetarysystems ] P E 1. INTRODUCTION tionsofbaroclinicdisks,e.g.witharadialentropygradi- . h entandthusverticalshear,whichtheyassumedtobea Angular-momentum transport and turbulence are p kindofbaroclinicinstability(BI)modifiedbytheKeple- important issues concerning protoplanetary disks. - rianshearprofile. Observedprotoplanetarydiskshave o Magneto-hydrodynamic turbulence brought about by a non-zero radial entropy gradient β = −dlns/dlnr, r themagnetorotationalinstability(MRI,Balbus&Haw- t where s is the entropy and r the radial distance to the s ley 1991), is a reliable way to achieve a sufficient a angular-momentumtransportandwiththisalsoanac- star. Withβ=q−(γ2D−1)pΣ,whereq=−dlnT/dlnr [ cretion rate fitting observations (Andrews et al. 2009) andpΣ=−dlnΣ/dlnrarethetemperaturesurfaceden- and playing an important role in planet formation (Jo- sity gradient respectively and γ2D is the 2D adiabatic v2 hansen et al. 2007; Lyra et al. 2008; Dzyurkevich et al. index, we see that disks that fulfill pΣ < q/(γ2D −1) 4 2010; Flock et al. 2011; Uribe et al. 2011; Johansen et al. indeed have a negative entropy gradient with values 6 2011). However, for MRI to be active the gas has to be fromAndrewsetal.(2009)ofq≈0.3−0.5and pΣ=0.9. 4 sufficiently ionized. This is only the case in the outer Therefore protoplanetary disks are not barotropic but 4 regions, upper layers of the disk, and in regions close rather baroclinic which means that planes of constant . tothestar. Theotherpartsofthediskaretoocoldand pressure and constant density are misaligned, creating 2 dust-richforsufficientionizationandthusthemagnetic a thermal wind, e.g. vertical shear. In a linear stability 1 fieldscannotcoupletothegas. Becauseofthis,theMRI analysis that followed (Klahr 2004) it was shown that 2 cannot operate in this region, which is therefore called thisinstabilitycanonlybeofnon-linearnature(seealso 1 "deadzone"(Gammie1996;Turner&Drake2009). Cabot1984;Knobloch&Spruit1986). : v Since the precise ionization structure is still under ThermalrelaxationturnedouttobecrucialwhenPe- i debate (Turner & Drake 2009) as is the interplay be- tersenetal.(2007a,b)studiedbaroclinicvortexamplifi- X tween active and dead-zones (Lyra & Mac Low 2012) cation using an incompressible approximation. In fact r we want to assess the precise hydrodynamic behavior thermalrelaxationordiffusion,besidestheentropygra- a ofdeadzones,becauseaccretionhastoproceedthrough dient,arekeyingredienttoestablishbaroclinicfeedback it somehow and it is where planets form. Therefore it thatkeepstheinstabilitye.g. vorticesinbaroclinicdisks isofinteresttostudypurelyhydrodynamicturbulence growing. in circumstellar disks. Klahr & Bodenheimer (2003) While both effects e.g. the baroclinic instability and foundsuchahydrodynamicinstabilitycreatingvortices baroclinic vortex amplification are a result of the su- inthree-dimensionalradiationhydrodynamicalsimula- peradiabaticradialstratificationofadisktheyarenotto be confused. An operating linear baroclinic instability (compare Cabot 1984; Knobloch & Spruit 1986) would 1Max-Planck-InstitutfürAstronomie,Königstuhl17,69117,Hei- beabletocreatevorticesindisksfrominfinitesimalper- delberg,Germany;[email protected],[email protected] turbations, whereas the baroclinic vortex amplification 2Department of Astrophysics, American Museum of Natural dealswiththegrowthofexistingvorticalperturbations, History,79thStreetatCentralParkWest,NewYork,NY,10024,USA 3Jet Propulsion Laboratory, California Institute of Technol- forwhichLesur&Papaloizou(2010)usedtheterm"sub- ogy, 4800 Oak Grove Drive, Pasadena, CA 91109, USA; criticalbaroclinicinstability"(SBI). [email protected] TheoccurrenceofaclassicalBIinthediskinitsgeo- 4NASACarlSaganFellow physicaldefinitionisstillunderdebateandshallbedis- 2 Raettigetal. cussedelsewhere. Therearethreepossibilities: 1.) there Vorticity is conserved in quasi-incompressible isaclassicalBIworkinginprotoplanetarydiskscreating barotropic simulations, but in flows with density the initial vortices, 2.) there is an other instability op- and pressure as independent quantities vorticity is erating (see the discussion in Klahr 2004) for instance producedviathesocalledbaroclinicterm creating vortices via Kelvin-Helmholz instability from (cid:18) (cid:19) ∂ω 1 1 vorticity maxima in sheared waves of baroclinic disks =∇× − ∇p = ∇ρ×∇p∝β∂ ρ. (1) or 3.) small vortical perturbations are triggered from ∂t ρ ρ2 y othereffects,e.g. wavesfromtheMHDactiveregionof Here ρ is the gas density, p the gas pressure, and β is thedisksormaybefromthewavesemittedbyvortices the global radial entropy gradient. The ground state at other radii. In any case the vortices are then grow- of a disk is geostrophic, e.g. all centrifugal forces and ingasdescribedbytheBVAuntiltheyreachasufficient gravity are in balance with the strictly radial pressure sizetoinfluencetheevolutionofthedisk,andthisisthe gradient. Ifanentropyperturbationisintroducedwith- physicsbeingsubjectofthepresentpaper. out perturbing the pressure, then this entropy pertur- Recently, Lyra & Klahr (2011) have examined the in- bationwillefficientlycreatevorticityinthepresenceof terplay of baroclinic vortex amplification and MHD. theglobalentropyandpressuregradients. Thiseffectis Theyfoundthatassoonasmagneticfieldsarecoupled basicallyradialbuoyancybecauseofsuperadiabaticra- tothegas,theMRItakesoverandthussuperseedsvor- ticeswhichwerepreviouslyamplifiedbyvortexampli- dialstratification5. IndeedtheradialBrunt-Väisäläfre- fication. Thisisevidencethatthevortexamplificationis quency(Tassoul2000) aphenomenonrestrictedtothedead-zone. (cid:18) (cid:19) 1 ∂p ∂ p All the above mentioned (lower resolution) studies N2=− ln (2) hadtoapplyentropygradients2-4timesstrongerthan γρ ∂r ∂r ργ to be expected in protoplanetary disks (Andrews et al. is imaginary, which would lead to radial convection. 2009, Klahr 2013 submitted) to drive BVA. We show in However,shearstabilizesnon-axisymmetricmodesand thecurrentpaper,throughhighresolutionrunsthatre- forthedynamicstabilityoftheaxisymmetricsystemthe alisticentropygradientsinprotoplanetarydisksaresuf- Solberg-Høilandcriterium(Tassoul2000;Rüdigeretal. ficientforBVA. 2002) RecentlyPaardekooperetal.(2010)haveinvestigated the effect of radial vortex migration. They discov- 1 ∂j2 1 eredthatvorticesmigratequicklyradiallyinwardonce − ∇p∇S>0 (3) R3 ∂R C ρ growntotheirfullsize. Whilethiseffectwillbeofma- p jorimportancetounderstandthelife-cylceofavortex,it ∂p(cid:18)∂j2 ∂s ∂j2 ∂s(cid:19) − <0 playsaweakerroleforthesmall/stillgrowingvortices ∂z ∂R∂z ∂z ∂R inthepresentpaper. Ofcoursemigrationwillinfluence theeffectiveangularmomentumtransportgeneratedby has to be considered. If one re-writes Eq. (3) for local the vortices via the emission of waves, but this is be- approximation(seee.g.Balbus&Hawley1998)thesta- yondthescopeof2Dlocalsimulationsasinourstudy. bilizingactionofthespecificangularmomentumshows Weshallreturntovortexmigrationandhaveabetteres- upasthevalueofOort’sconstantintheCoriolisterm.If timateforangularmomentumtransportoncewereturn alsotheverticalstratificationinvelocityistakenunder toglobalsimulations. consideration,asitwilloccurinrealthree-dimensional We carry out local, compressible shearing sheet sim- accretiondisks(Fromangetal.2011),thenthecombined ulationsatvariousresolutions. Weshowthataswego action of radial buoyancy and Coriolis forces lead to a tohigherresolutionsonecanexcitethenonlinearinsta- thermalwind,e.g. averticalshearinrotationalvelocity. bility and achieve Reynolds stresses with the low en- Thisispreciselytheinitialstateasbaroclinicinstability tropy gradients deduced for observed accretion disks. inrotatingstarsandplanetaryatmospheres. Yet, insta- We conduct an extensive parameter study for entropy bilityinthesesystemsisnotobstructedbyradialshear, gradients(β),resolution,thermalcooling(τ )anddif- whereasinaKepleriandiskradialscaleswouldhaveto cool fusion times (τ ) respectively. Section 2 gives a brief beontheorderoftheverticalpressurescale-heigth(H) diff overview of the physical background of the instability. (Knobloch & Spruit 1986) to be linearly unstable with InSection3wepresentthenumericalsetupofoursim- respecttobaroclinicinstability. ulations. InSection4weexaminetheamplificationand Beforeweexplainthemotionofagasparcelinavor- decay-timesofvaluessuchasenstrophyω2=(∇×u)2 tex we want to explain the cooling and heating pro- z z and α-stresses. Here α=(cid:104)ρu u (qp )−1(cid:105) with ρ being cessesinadiskastheyprovedtobecrucialtomaintain x y 0 the baroclinic feedback (Petersen et al. 2007a,b). Dust thegasdensity, uthegasvelocity, q=1.5theshearpa- particlesabsorbphotonswhichheatsthemup. Tocool rameter,and p theinitialmeanpressure. Wealsoana- 0 theyradiatephotonsintheinfrared. Thisradiationcan lyze the saturation values, e.g. how quantities like the beabsorbedby other particles. Thishappenson atyp- entropy gradient, cooling processes in the disk or the ical length-scale. A convenient parametrization for the size of the simulated domain influence the strength of angular momentum transport. Finally we summarize 5Notethatsimilarsituationscanbefoundinsubadiabaticconfig- ourresultsandgiveaconclusioninSection5. urations. Infact,inanynonbarotropicdisk,anentropyperturbation willleadtoavorticityfluctuation.Butwithouttheglobalpressureand entropygradientpointinginthesamedirectiontheseperturbations willquicklydecay(shearaway)astheyarelackingthemechanismof 2. PHYSICALBACKGROUND vortexamplification. AParameterStudyforBaroclinicVortexAmplification 3 diffusiontimeinourvortexsystemisτ =a2/Kwhere gradient β. Note that in our approximation the gradi- diff aistheradiusofthevortexandKthediffusionconstant. entsforentropy(s)andpressure(p)arethesame.There- The diffusion constant can be approached using a flux forewedonotdistinguishbetweentheminournotation limiteddiffusionapproachasinKleyetal.(2009).There andcallboth β. However,inrealdisksbothmayeasily K = λc4a T3(ρκ)−1 where λ is the flux limiter, c the differ. R Thetotalpressurep =p¯+pconsistofalocalfluctu- speedoflight,a theradiationconstant,Tandρthegas tot R ationpandatime-independentpartthatfollowsalarge temperature and density, respectively and κ the opac- scaleradialpressuregradientβ ity. Since K is constant and the vortex grows τ will diff change over time. Thermal relaxation is the other pro- cessbywhichdustcandepositheatintothegas. When p¯= p0(r/R0)−β, (4) adustparticlehasacertaintemperatureotherthanthe equilibriumtemperatureitwillexchangeheatwiththe where r is the cylindrical radius. The full set of lin- ambient medium until it reaches the background tem- earizedequationsusedinoursimulationsis peratureagain. τ isthetimeneededtoachievethis. Thistime-scaleafcfoeocltsvorticesofallsizesequally. Dρ+(u·∇)ρ=−ρ∇·u+ f (ρ) (5) Thebaroclinicfeedbackitselfwasexplainedindetail Dt D by Petersen et al. (2007b). A nice description of the Du 1 mechanism can also be found in Lesur & Papaloizou Dt +(u·∇)u=−ρ∇p−2Ω0(zˆ ×u) (2010). Inabaroclinicflowentropyisafunctionofpres- (cid:18) (cid:19) sureanddensity, s(p,ρ). Pressureontheotherhandis +3Ω u yˆ + βp0 1 − 1 xˆ + f (u,ρ) (6) onlyafunctionofradius. Thevortexinteriortransports 2 0 x R ρ ρ ν 0 0 highentropymaterialfromsmallradiitolargeradii.Af- Ds 1 (cid:26) (T−T ) ter thermalization low entropy material is transported +(u·∇)s= ∇·(K∇T)−ρc 0 to small radii. Since the pressure variations, especially Dt ρT v τcool fromweakvortices,arenegligibleincomparisontothe βp u (cid:27) global radial pressure gradient and much smaller than + R0 (γ−x1) + fK(s). (7) the azimuthal entropy gradient, pressure can be seen 0 asapproximatelyazimuthallyconstant(Klahr&Boden- Here ρ is the gas density, u is the deviation of the gas heimer2003;Klahr2004;Petersenetal.2007a). Tokeep velocity from the Keplerian value, T the temperature, the pressure constant an azimuthal density gradient is c the specific heat at constant volume and, K the heat established, e.g. outflowing material has a lower den- v conductivity. Tthermal diffusion time-scale is denoted sityasinflowingmaterial. Thusthevortexfeelstheef- byτ . Thesymbol fectofdifferentialbuoyancywhichestablishedthepos- cool itivebaroclinicfeedback(Eq. (1)). D ∂ ∂ Ifcoolingistoofast(shorttime-scales)thenthefluid = +u(0) (8) parcel adapts the background temperature slope too Dt ∂t y ∂y quickly. The vortex becomes locally isothermal and no entropy transport is possible. Conversely, if cooling is represents the Keplerian derivative where u(0) = y too slow (long time-scales) then gas will not be ther- −3/2Ω x. 0 malizedfastenough. Thevortexgasbecomesadiabatic For a more thorough derivation of these equations withconstantentropyacrossthevortex.Inbothextreme andthelinearizationoftheglobalpressuregradientwe cases, isothermal or adiabatic, the azimuthal entropy refertoLyra&Klahr(2011)andtheappendixtherein. gradient across the vortex vanishes. As shown in Eq. In order to keep the numerical scheme stable we (1) the vorticity source ceases to amplify the vortex, or add sixth-order hyperdiffusion f (ρ), hyperviscosity D atleaststabilizesitagainstlossesfromnumericalviscos- f (u,ρ), and hyperconductivity f (s) (Lyra et al. 2008, ity from radiating vorticity perturbations, e.g. Rossby ν K 2009;Oishi&MacLow2009). waves. Therefore it is important that thermal cooling The radiation processes in the disk are implemented anddiffusiontimesareintherightregime. throughthefirst(thermaldiffusionasanapproximation Wemodelboththermalrelaxationandthermaldiffu- for flux limited diffusion of radiation energy density) sion separately because, dependent on the vortex size, andsecond(thermalrelaxationtomimicheatexchange either one or the other dominates thermalization. Al- with the surface of the disk and thermal equilibration ways the process with the shorter time-scale sets the withtheirradiationfromthecentralobject)termsonthe heatexchangebetweenvortexandambientgas. righthandsideoftheentropyequation. Asmentioned in the last chapter we keep the diffusion coefficient K, 3. NUMERICALSETUP which is defined as in (Kley et al. 2009), constant and Our simulations were conducted with the PENCIL defineitsvaluevia τ = H2/K. Soifthevortexhasa CODE6. Weuseatwo-dimensional,localshearingsheet radiusofH,thepressduifrfescale-hightofthedisk,thedif- approach. Weconsiderasheetinthemid-planethatco- fusiontime τ hasthevaluewequoteine.g. Table1. diff rotates with the co-rotational radius R0. This is a 2D If the vortex is smaller than H relaxation will be much version of the model used in Lyra & Klahr (2011). To faster. includethebaroclinictermtheydefineaglobalentropy To clarify that it is indeed the global entropy gradient that produces the vorticity we take the curl of the 6Seehttp://www.nordita.org/software/pencil-code/ Navier-StokesEq. (6)andassumeanequilibriumstate, 4 Raettigetal. TABLE1 SIMULATIONSETUPANDRESULTS run βP τcool τdiff ωz2 α uRMS ρRMS x-res x-domain (2πΩ0−1) (2πΩ0−1) (Ω02) (cs) (gridcellsH−1) H A 2.0 10 10 0.056 1.05×10−2 0.33 0.22 144 4 A2 10 10 2.15×10−2 3.09×10−3 0.19 0.19 144 8 B 10 10 0.060 1.21×10−2 0.33 0.22 288 4 C 1.0 10 10 0.051 8.67×10−3 0.31 0.22 144 4 C2 10 10 4.63×10−3 8.2×104 0.08 0.06 144 8 D 10 10 0.059 9.63×10−3 0.31 0.21 288 4 E 30 10 0.022 4.33×10−3 0.23 0.15 144 4 F 30 30 0.022 3.72×10−3 0.23 0.15 144 4 G 100 10 0.017 2.61×10−3 0.14 0.08 144 4 H 100 30 0.013 2.22×10−3 0.15 0.08 144 4 I 100 100 0.010 1.36×10−3 0.14 0.08 144 4 J 0.5 10 10 5.25×10−3 6.38×10−4 0.35 0.04 144 4 J2 10 10 1.77×10−3 8.91×10−5 0.03 0.03 144 8 K 10 10 4.30×10−3 4.30×10−4 0.57 0.05 288 4 L 30 10 0.021 3.89×10−3 0.23 0.15 144 4 M 30 30 0.021 3.01×10−3 0.23 0.15 144 4 N 100 10 6.00×10−3 1.38×10−3 0.14 0.10 144 4 O 100 30 6.00×10−3 1.38×10−3 0.14 0.10 144 4 P 100 100 8.63×10−3 1.18×10−3 0.15 0.10 144 4 u =0,and∇P=0sothat theterminEq. (9)thatcreatesthedevelopmentofnon x laminarflowstructure. Dω βp z = 0 ∂ ρ. (9) All our simulations are done in dimensionless code- Dt ρ2R0 y units. Sothat R0=Ω0=1,γ=1.4,andcs =0.1,which means that H = 0.1. All time-quantities are given in Hereweseethatthenegativeazimuthaldensitygradi- entacrossthevortexisthesourceforvorticityproduc- 2πΩ−1 which is one local orbit at the co-rotational ra- 0 tionproportionaltotheglobalentropygradient. diusR . 0 ShearingsheetsimulationswithZeus7 likefinitevol- The individual setups are given in Table 1. The ther- ume codes without explicit viscosity, e.g. the TRAMP mal cooling times and thermal diffusion times are de- code, have shown a weak amplification of kinetic en- rivedfromstandarddiskmodelslikeinBelletal.(1997), ergy for the pure adiabatic case, i.e. infinite cooling alsoseeKlahr2013submitted. time (see Klahr 2013 ApJ submitted). This numerical We explored different resolutions in our simulations, artifact does not occur with simulations performed by namely 2882, 5762 and 11522. The unusual non power the PENCIL CODE. See Appendix A for a 1D radial of2resolutioncomesfromourcomputationalplatform test/comparisonsimulation. with6coreprocessors.Typicallyweusedupto24CPUs Initiallyweapplyafiniteperturbationinthedensity totaling144coresforourlargestgrids. Stillweneeded sothat about1200hoursperrun. Thegridcovers±2Haround ρ(x,y)=ρ0+ρ(cid:48) (10) R0 in the radial and [0H,16H] in azimuthal direction. with ρ theconstantbackgrounddensityand ρ(cid:48) theac- This leads to an effective resolution of 72 (2882), 144 0 (5762)and288(11522)grid-pointsperscalehightinra- tualperturbationoftheform dial direction and 18 (2882), 36 (5762) and 72 (11522) ρ(cid:48)=ρ0Ce−(x/2σ)2× ∑kx ∑ky sin(cid:26)2π(cid:26)i x +j y +φij(cid:27)(cid:27), gnreicde-spsoarinytstopceormHprionmaizsiembuetthwaelednirreecstoiolunt.ioInt aisndalwcoamys- L L i=−kxj=0 x y putationaltime. Lowerresolutionsimulationsarecom- putationally less expensive but might not resolve the whereC describesthestrengthoftheperturbation. We necessaryscales. perturb the density in a way that ρ = 5% for β = rms 1.0,2.0 (runs A-I) and ρ = 10% for β = 0.5 (runs J- rms P). To achieve a random perturbation we apply an ar- 4. RESULTS bitrary phase φij between 0 and 1. The initial state is 4.1. SaturationValuesandConvergence non-vortical. Again, this is the identical initial condi- We show the time-developement of α-stresses in tionasusedinLyra&Klahr(2011)aswellasthesame amplitude,C,forsimulationswithβ=2.0,aswasused Fig. 1. The green line shows the resolution of 2882, intheirsimulations. black of 5762 and red 11522 for β = 2.0 (top), β = 1.0 Notethatwiththisinitialperturbationwedonotper- (middle) and β =0.5 (lower panel). In all simulations turbthepressurebuttheentropy. Thusitisreallyonly τdiff=τcool=10localorbits. Weseethatforβ=1.0and0.5andaresolutionof2882 7http://www.astro.princeton.edu/˜jstone/zeus.html theperturbationdecaysrightaway.Higherresolutionis AParameterStudyforBaroclinicVortexAmplification 5 10−1 100 10−2 10−2 10−3 10−4 α 10−4 α 10−6 β=2.0 10−5 β=1.0 res=2882 α res=5762 10−6 τ=70 10−8 res=11522 10−7 10−10 0 500 1000 1500 100 time [2π/Ω] 0 10−2 FIG.2.—Timeevolutionoftheα-valuesandenstrophyforβ=1.0 and a resolution of 5762 (run C). The red slope marks exponential 10−4 amplification with a amplification-time τ =70Ω2π. For larger en- 0 α tropy gradients (smaller entropy gradients) we get faster (slower) 10−6 β=1.0 amplification-times. res=2882 res=5762 10−8 res=11522 4×10−3forentropygradientsaslowasβ=0.5. Infact, 10−10 inSection4.5weshowthatthereisonlyaweakdepen- 100 denceofαonβasα∝β0.5. Fig. 1showsthatthesatura- tionvaluesofαdonotdependstronglyon β,butaswe 10−2 willseeinthenextsectiontheamplificationratesdo. 10−4 4.2. Amplification-andDecayRates α We analyze the amplification timescales of the vor- 10−6 β=0.5 tices,meaninghowfastavortexgrowsduetothebaro- res=2882 clinic feedback. Thus it is independent of the precise res=5762 10−8 shape of the initial condition as long as the amplitude res=11522 islargeenoughforthegivenReynoldsnumbertohave 10−10 vortexgrowth. Infact,theinitialstrongkickneededto 0 200 400 600 800 1000 getthevortexgoingdecaysratherquicklyascanbeseen time [2π/Ω] ine.g. Fig.1. Here,theα-valuesstartoutintheorderof 0 10−5 then drop to around 10−8 as the initial perturba- tiondecays. Assoonasthebaroclinicfeedbacksetsin, FIG.1.—Timeevolutionofα-stressesforthethreedifferentresolu- the values rise again. The timespan that follows is the tionsof2882,5762and11522withanentropygradientofβ=2.0(green onewherewemeasuretheamplificationtime. line),β=1.0(blackline)andβ=0.5(redline). Forallthesemodels tτhdeiffre=foτrceooaln=gu1l0ar·m2πo/mΩe0n.tuFomrtarlalnrespsoolruttciaonnsbveosreteenxfaomrsptlriofincgateionntroanpdy weInfianndatlhyaztintghethineitaimalpalmifipcalitfiiocant-iroante-rsaotef othfethinesαta-sbtirleitsys gradients(β=2.0). Forlowerentropygradientshigherresolutionis (Γ(α)), as can be seen in Fig. 2 for run C, can be fit- neededtoseethedevelopmentofvortices.Thedashedlinesshowthe ted as exponential amplification α = α exp(t/τ) with saturationsvalues(β=2.0and β=1.0)andvalueattheendofthe 0 simulation(β=0.5)respectively. τ ≈70β−2. The proportionality to β−2 is not what one wouldnaivelyexpectfromalinearconvectiveorbuoy- ancydriventurbulence. requiredtoincreasetheReynolds-numberofthesystem For a linear buoyancy driven turbulence one would andhavelessdissipationonthesmallerscalesandthus expect an amplification rate proportional to the Brunt- excitetheinstabilityagain. Väisäläfrequency,N Wetakeastrongerinitialperturbationforβ=0.5than (cid:18) (cid:19) forthehigher β. Theperturbationinentropyresultsin 1 ∂p ∂ p N2=− ln (11) aperturbationinvorticity. Thisperturbationispropor- γρ ∂r ∂r ργ tionaltoβ. Forsmallβwehavetoapplyastrongerper- turbation to get the same effect on the vorticity. How- whichinourparameterslookslike ever, we expect that if we go to even higher resolution it is possible to keep the initial density perturbation at 1 (cid:18)H(cid:19)2 N2=−β β Ω2∝−β2. (12) ρrms=5%(Petersenetal. 2007). p sγ R Ifwecomparethesaturationvaluesofrunswithdif- ferent resolution, we see that they differ by only 10 % Here we explicitly wrote βp and βs to make clear that fromoneanother(seeTable1). theBrunt-Väisäläfrequencydependsontheproductof It is important to note that the instability is excited entropyandpressuregradientwhichcanbedifferentin and we measure α-values in the converged runs up to globalsimulations. 6 Raettigetal. mentumtransport,beforetheydriftinward. 10−1 To study the numerical dissipation effects even fur- β=1.0 therwenowassesshowthevorticesdecayifbaroclinic res=5762 10−2 drivingisswitchedoff(Fig.3).Todothiswefirstevolve runs C and D with β =1.0 and the two resolutions of α 5762 and 11522 for 800 orbits and then turn off the en- & 10−3 tropygradientsothat β=0.0. Weobservethatthevor- 2ωz ticesgetsmallerandthatallrelevantquantitieslikevor- ticity, ω2, or α-stresses decay with exponential behav- 10−4 ωτz=21000 ior. Godzon&Livio(1999)sawthesameexponentialde- ω α cay of vorticity when they analyzed longevity of anti- 10−5 τα=400 cyclonicvorticesinprotoplanetarydisks. Theirdissipa- 10−1 tion was proportional to the effective viscosity applied β=1.0 in their numerical experiment. Here we find the same res=11522 10−2 decay-rateforbothresolutions,highlightingthatthede- cay of vortices is no longer through numerical effects, α butduetotheradiationofwavesasinKorotaev(1997). & 10−3 2ωz 4.3. SaturationValues Wehaveestablishedthatevenshallowentropygradi- 10−4 ωτz=21000 entsleadtovorticesbutwestillhavetoshowthatsuffi- ω α cientangularmomentumtransportcanbereachedwith 10−5 τα=400 these shallow gradients. The saturation values of en- 500 1000 1500 2000 2500 strophy, ω2, or u are of interest as well. Note that z rms time [2π/Ω] wetalkaboutsaturationvaluesofour2Dlocalsimula- 0 tions, where certain restrictions apply, see a more de- FIG.3.—Inthisrunwithβ=1.0aresolutionof5762(runC,upper tailed discussion in the conclusions. In the next sec- panel) and 11522 (run D, lower palnel)and we turn off the entropy tionswediscussthemeasuredsaturationvaluesandan- gradientafter800localorbits(indicatedbytheblackdashedline)and alyzehowthedifferentcontrollingparametersinfluence seehowtheinstabilitydecays.Enstrophyisshownwiththeblackline amplification-phaseandfinalvalues. andα-stresseswiththeblueline.Ourfitisgiventhroughtheredand fgorreethnedeansshterodplhinyeasnrdesτpαe=cti−ve4l0y0. fWoreαfi.tadecaytimeofτωz2 =−1000 4.3.1. InfluenceofEntropyGradient In Fig. 4 we compare runs A, C and J (at a resolu- All quantities in Eq. (12) are positive. Thus the tionof5762 and τ =τ =10)whichdifferonlyre- diff cool Brunt-Väisälä frequency is imaginary and therefore garding the value of β. There is an initial exponential a linear buoyancy driven turbulence would have a amplification-phaseof α,E and ω2 thatis shorterfor amplification-rateΓ∝ iN∝β. However,wefoundthat high β, followed by a satukrianted statze. We also see that Γ∝β2 providesabetterfit. Thisonceagainreflectsthat for lower β the saturation values are lower. We want thebaroclinicvortexamplificationisanon-lineareffect. to stress that we did not reach saturation for simula- Inlinearconvectiveinstabilityadisplacedparcelofgas tions J and K (at a resolution of 5762 and 11522 and feels a buoyancy force and thus accelerates propotion- τ = τ = 10). Even after 3000 local orbits vortex diff cool allytoβ. Butinthediskbaroclinicinstabilityfirstavor- amplificationwasstillongoing. Here,τ =10ismuch diff texhastoformwithanazimuthalentropygradientpro- shorter than the amplification-rate we estimated in the portionaltoβ(andτ )andinasecondstepthisvortex previous section (τ ≈300). As we will see in the next cool feelsatorqueproportionaltoβ.Thereforetheamplifica- sectiontheamplification-phaseisshortestifthosetime- tionisproportionalto β2. The β2 andτcool dependance scales are comparable, because τdiff also defines how hasalsobeenderivedbyLesur&Papaloizou(2010),see fast pressure perturbations are damped. Although we theirEq. (23)foranorderofmagnitudeestimateofthe expectthesaturationvaluesofsimulationJandKtobe growthrate. higher than what they are right now, it is possible that The amplification behavior in Fig. 1 already shows theywillstillstaybelowthesaturationvaluesobtained convergencefor576gridcellsresolution,e.g. 144/H in insimulationswithhigherβ. radialdirection. Thevorticitycanbeseenasameasureofthestrength If we compare our amplification timescales for the ofthevortex. Thehighertheabsolutevalueofthevor- lowest entropy gradients with the migration times ob- ticity the stronger the vortex. The only stable vortices tained by Paardekooper et al. (2010) we see that they indisksareanticyclonic8andthereforethevorticityhas areofthesameorderofmagnitude. Whichmeansthat negative values. So the minimum value of vorticity the vortex could have drifted into the central star be- (ω )showshowstrongavortexis.Toexplainthebe- z,min foreitreachesstrongα-values. However,Paardekooper havior of ω (3rd panel in Fig. 4), cooling processes z,min etal.(2010)alsostatethattheirtimescalesrefertofully havetobetakenintoaccount. Duringtheearlyphases grown vortices of size H. Smaller vortices drift signif- icantly slower. This gives them enough time to reach 8 Cyclonicvorticesarealsopossible,butarequicklydestroyedby a size, with which they provide sufficient angular mo- shear(Godon&Livio1999). AParameterStudyforBaroclinicVortexAmplification 7 100 10−2 10−2 10−4 10−4 Ekin 10−6 β=2.0 α 10−6 ττcool==1300,,ττdiff==1100 β=1.0 τcool=30,τdiff=30 10−8 β=0.5 10−8 ττccooooll==110000,,dττiffdiff==1300 τcool=100,τdiff=100 10−10 10−10 cool diff 100 100 10−2 10−1 10−4 c]s α 10−6 β=2.0 u [rms10−2 τττccooooll===133000,,,τττddiiffff===113000 β=1.0 τcool=100,dτiff =10 10−8 β=0.5 ττccooooll==110000,,ττddiiffff==31000 cool diff 10−10 10−3 1.0 0 500 1000 1500 time [2π/Ω] β=2.0 0 0.5 β=1.0 2]π β=0.5 nFuImGb.e5r.s—)fCorosmampaeriβso=n1o.0f(dRifufenrseCnt-Iτ)d.Tiffh(eritgohptpnaunmelbsehros)wasntdheτcαo-ovla(lleufet /Ω0 0.0 andthebottomone urms. Onecanseethattheearlyamplification- [ phaseisdeterminedbythediffusiontimesincetheheatingacrossthe min vortexismoreimportantthenverticalheattransport. Wegetfaster ωz, amplificationforhigherτ .Oncethevortexgrowslargerheattrans- −0.5 diff portgetsmoredifficultandthermalrelaxationdominates. Therefore thesaturationvaluesaredeterminedthroughτ .Saturationvalues cool arehigherforshorterτ −1.0 cool. 0 1000 2000 3000 diffusion and relaxation influence the saturation val- time [2π/Ω] ues and the amplification-phases. As long as τ = 0 diff(l) l2/K < τ , τ will dominate the heat exchange cool diff(l) from the inside of the vortex to the ambient disk. As FIG.4.—TimeevolutionofkineticenergyEkin (top),α-value(mid- thevortexgrowsτ willincreaseandwiththatonly dle)andminimumvorticityωz,min (bottom)foraresolutionof5762 diff(l) andτ =τ =10butdifferententropygradientes:β=2.0(green), contribute to the heat exchange at the outskirts of the diff cool β=1.0 (black) and β=0.5 (red) (runs A, C, J). Saturation is first vortex. τ willthen dominatetheinteriorof thevor- cool reachedforhighβalreadyafter300orbits,thenforβ=1.0Forβ=0.5 tex. nosaturationisreachedevenafter3000orbits. Theincreaseinωz,min Forthesimulationswherewesetτ =τ ,τ will afterthepointintimewhensaturationisreachedcanbeexplained diff cool cool throughtheheattransportacrossthevortex. Sinceithasreachedits take over when the vortex has reached a size of H. In finalandlargestsizeheattransporttakeslongerduetothelargersize radialextendthishappensoncethevortexhasgrownto ofthevortex. itsfinalsize. ThisisconsistentwithwhatweseeinFig.5. During thermalization is dominated by thermal diffusion (Pe- the early amplification-phase simulations with equal tersenetal.2007b). Asmentionedbeforethistime-scale τ behaveexactlythesame.Eventuallyτ takesover diff cool isshorterforsmallervortices. Thereforeheatexchange sothatthesaturationvaluesaredeterminedbyτ .For cool betweenthevortexgasandtheambientgasismoreef- longerτ saturationvaluesarelowerthanforshorter cool ficient than in later stages. Once the vortex has grown τ . cool to its final size, thermal relaxation takes over. How- everheatexchangeinthecenterofthevortexislesseffi- 4.3.3. InfluenceofPhysicalDomain cientthanintheearlierstages. Thebaroclinicfeedback, A problem with local shearing sheet simulations is e.g. theazimuthalentropygradientacrossthevortex,is that eventually vortices grow to box-size. We cannot lessefficient, thevortexgrowsweaker, and ω rises z,min saywhethertheyhavereachedtheirfinalsizeorjustdo again,creatingaflatyetextendedvortex. nothaveanymoreroomtogrow. Anotherproblemthat ariseswiththeperiodicboundaryconditionsisthatthe 4.3.2. InfluenceofThermalDiffusionandCoolingTimes vortices potentially interact with themselves and thus We take a closer look at simulations with β =1 and forcing (shaking) them to shed more waves and there- differentcombinationsofKandτ toseehowthermal foreincreasethe α-values. Todealwiththat, were-did cool 8 Raettigetal. ω ω z z −1.5 −0.5 0.5 −1.5 −0.5 0.5 10−2 β=0.5 3.0 3.0 10−3 10−4 2.5 2.5 α 10−5 2.0 2.0 10−6 small 10−7 large y y 1.5 1.5 1.5 1.5 10−8 10−2 1.0 1.0 1.0 1.0 y y 0.5 0.5 0.5 0.5 ω2z 10−3 0.0 0.0 0.0 0.0 −0.4−0.2 0.0 0.2 0.4 −0.2 0.0 0.2 −0.4−0.2 0.0 0.2 0.4 −0.2 0.0 0.2 x t=100 x x t=500 x ω ω z z small −1.5 −0.5 0.5 −1.5 −0.5 0.5 large 10−4 3.0 3.0 0 500 1000 1500 2000 2500 time [2π/Ω] 0 2.5 2.5 FIG.7.—Timedevelopmentofαandωz2withβ=0.5forsmall(black) andlarge(red)physicaldomain(runsJandJ2).Saturationvaluesare lowerinthelargeboxthaninthesmallerbox. 2.0 2.0 direction,andtheythuspasseachotherlessfrequently y y 1.5 1.5 1.5 1.5 duetotheextendedazimuthaldomain.Eventuallythey can merge as Godon & Livio (1999) saw, but the larger theboxthelongerittakes.Wedonotwanttodiscussthe 1.0 1.0 1.0 1.0 mechanism of how the process of vortex merging hap- y y pensexactly. Thishasbeenexplainedextensivelyinthe fieldoffluiddynamics(seee.g.Cerretelli&Williamson 0.5 0.5 0.5 0.5 2003). Themergingprocessitselfisnotthefocusofour study, because a) the vortex merging is strongly influ- 0.0 0.0 0.0 0.0 encedbytheboxdimensionsinashearingsheetsimu- −0.4−0.2 0.0 0.2 0.4 −0.2 0.0 0.2 −0.4−0.2 0.0 0.2 0.4 −0.2 0.0 0.2 x t=1000 x x t=1500 x lationandb)2Dflatvorticesmergedifferentlythanfull scale3Dvortices.Theimportantthingisthatvorticesdo FIG.6.—Snapshotsofthez-componentofthevorticity,ωzafter100, indeedmergeifthearesufficientlyclosetooneanother, 500,1000,1500localorbitsforthetwodifferentphysicaldomainswith β=0.5. Initiallybothrunshavevorticesofequalsize. Sincethereis butconserveωintheprocess. lessspacebetweenvortices,theycanmergesoonerinrunswiththe Anotherunphysicalprocessthatcanoccurinlocalpe- smallphysicaldomain.Thevorticesinthelargephysicaldomaintake riodic simulations is that when the vortex approaches longertogrow. Thedashedwhiteboxinthelastplotindicatesthe areaofthesmallphysicaldomain. theintegralscaleitinteractswithitself,theouteredges of the one side of the vortex almost touches the other sideofthesamevortex. Wedonotseethisfortheruns simulationsA,CandJwithadoubledphysicaldomain with the larger physical domain. Since the vortices in (simulations A2, C2, J2 in Table 1). The resolution is the larger domain do not interact with themselves, the thesame. Insteadofx=[−0.2,0.2]andy=[0.0,1.6]we saturation values are lower. However, they are still in switch to x = [−0.4,0.4] and y = [0.0,3.2]. We did not thesameorderofmagnitude(seeTable1). adjusttheinitialperturbationinanyway. Thereforethe InFig.6weshowsnapshotsofthevorticityforβ=0.5 initialstateisperturbedatsmallerwavenumbersthan (simulationsJandJ2).Initiallythereareseveralvortices. inthesmallerdomain. Ifwegotoevenlargerboxesthe Thelargeronessweepupthesmallervorticesandthus initial condition has to be adjusted so the the effective growfurther. At1500localorbitsthereisonlyonevor- perturbationinthedensityisofthesamestrengthasin tex left for the small physical domain, whereas in the thesmallerphysicaldomain. largerphysicaldomaintherearestillthreevortices. If we compare the time development of runs with a Ifwelookattheα-valueandenstrophyforthesetwo different physical domain (see Fig. 6), we see that vor- simulations (see Fig. 7) we see that the value seems to ticesinfactdonotmergeasfastinthelargedomainbe- decay in the larger box at the end of the run. How- cause there now is more space between them in radial ever this does not mean that the vortices die out. It AParameterStudyforBaroclinicVortexAmplification 9 ω uu [c2] z x y s 0.0015 −1.5 −0.5 0.5 −0.002 0.0 0.002 0.0010 3.0 3.0 0.0005 c]2s u [y ux 0.0000 2.5 2.5 −0.0005 −0.0010 2.0 2.0 −0.4 −0.2 0.0 0.2 0.4 x y y −1.5 1.5 1.5 −1.0 1.0 1.0 ωz −0.5 0.5 0.5 0.0 2.2 2.0 1.8 −0.4 −0.2 0x.0 0.2 0.4 −0.2 0x.0 0.2 y 1.6 1.4 1.20.10 0.12 0.14 0.1x6 0.18 0.20 0.22 t=1000 FIG.8.—Vorticityprofile(left)andα-stress(middle)forβ=1.0andthelargephysicaldomain(runC2).Yellowandredareasdenotepositiveα- vaueswhereasblueareasshownegativeα-stresses.Ingreenareasα=0.Onecanseethewavesexcitedbythevortex.Thosewavesareresponsible fortheangularmomentumtransport. Itisalocalizedprocess. Sincethevortexandthevorticity-wavesfilloutasmallerareaoftheboxinthe largebox(largegreenareaswherethereisnoangularmomentumtransport)andourcalculationofthesaturationvaluesaveragesovertheentire areaofthebox,thesaturationvaluesseemtobelower. Theplotintherightupperpanelshowsanazimuthalaverageovertheuxuy. Insidean idealvortexα-stresseswouldsumuptozero. However,asindicatedinthelowerrightplot,thevortexhasacomplexstructurewhichleadsto deviationsfromtheidealizedshape. more so reflects fluctuations in the vortex interaction, netα-valueacrossthevortex. modulating α, as also can be seen in the small domain Toproperlycomparethevaluesofαforbothphysical case at high frequency. We calculate the values as a domains,theboxaveragehastobetaken.Iftheaverage meanovertheentireboxbutespeciallytheangularmo- overanequalphysicalsizecenteredaroundavortex,as mentum transport is a very localized process as can be indicated by the white dashed lines in Fig. 6, is taken, seen in Fig. 8 (this time for β = 1.0 after 1000 orbits). then the α-values agree again. The α-values are gener- Here we show the product u u at each location in the atedonlyinthevicinityofvortices. x y box. Most areas of the box have an u u -value close x y to zero. However, one can clearly see bands excited 4.4. Correlations by the vortex with positive u u -values. These bands x y It is a feature of baroclinic instability that the satu- areinertia-acousticwaveswhichareresponsibleforthe angular momentum transport (Klahr & Bodenheimer ration values of urms, ωz2, ρrms seem to correlate with 2003; Mamatsashvili & Chagelishvili 2007; Heinemann each other. In Fig. 9 we plot the different quantities as &Papaloizou2009; Tevzadzeetal.2010). Ifwe hadan a function of α. Figure 9 shows the dependencies on α idealvortexwithasmoothsurfacewewouldexpectthat for all our simulations. The colors represent the differ- u u sums up to zero within the vortex. However the ententropygradients: β=2.0(black), β=1.0(red)and x y β=0.5(green). Thedifferentcombinationsofdiffusion vortex has a more complex structure as can be seen in andcoolingtimesarerepresentedthroughthedifferent thelowerrightplotofFig.8. Thisleadstoannegative symbols. We find that the following relations are good 10 Raettigetal. 10−1 0.07 0.06 0.05 10−2 ~β0.5 0.04 α τ =10 2ωz τcool=10, τdiff=10 τcool=30 0.03 ττccooooll==1300,,τ τdidfiff=f=1100 10−3 (()) ττccooooll==11000 0.02 τcool=30, τdiff=30 τdiff=30 0.01 ττcool==110000,, ττdiff==1300 10−4 τddiiffff=100 cool diff 0.00 τcool=100, τdiff=100 0.1 1.0 10.0 0.25 β 0.20 FIG.10.—Saturationvaluesofαforallourrunswiththesmallerbox dependingonβ. Runswithparenthesesaroundthemwerenotsatu- 0.15 τcool=10, τdiff=10 ratedattheendofthesimulations.Thereforewedonottaketheminto τ =10,τ =10 accountwhenwefittheα−β-relation.Thesymbolsshowthedifferent ρrms τccooooll=30, τdidfifff=10 combinationsofτcool(symbols)andτdiff(colors). 0.10 τ =30, τ =30 cool diff τ =100, τ =10 τcool=100, τdiff=30 dependencies (varying initial conditions) before we do 0.05 τcool=100, τdiff=100 three-dimensionalsimulations. cool diff 0.00 4.5. Dependenceonβ 0.4 InSection4.2weshowedthatamplificationofvortices for low entropy gradients is computationally demand- 0.3 ingintermsofevolutiontime. Thusitisdifficulttoex- tractsaturationvaluesforentropygradientsevenshal- [c]s 0.2 ττcool==1100,,τ τdif=f=1100 lowerthan β=0.5withthecomputationalresourcesat u rms τccooooll=30, τdidfifff=10 hand. τ =30, τ =30 In Fig. 10 we plot the α-stresses as a function of the cool diff 0.1 τcool=100, τdiff=10 entropygradient. Notethatwechooseadifferentcolor- τ =100, τ =30 cool diff codingthaninFig.9. Heresymbolsrepresentthether- τ =100, τ =100 0.0 cool diff malcoolingtimeswhereascolorsrepresentthermaldif- fusion times. The dashed black line illustrates a slope 0.000 0.005 0.010 0.015 ∝ β0.5 which is a reasonable fit for the set of points α with τ = 30,τ = 10 (black triangles) and τ = cool diff cool 100,τ =30(orangex). Wecannotpredictα-valuesfor diff specific entropy gradients and thermal cooling and re- FIG.9.—Saturationvaluesofωz2,ρrmsandurmsasafunctionofsat- laxationtimes. urated (value at the end of the simulation for β=0.5) α-value and allourrunswiththesmallphysicaldomain(runsA-P).Thesymbols The key issue is less a strong correlation between α showthedifferentcombinationsofτcool(leftnumbers)andτdiff(right and βbutratherthelackthereof. Thestrengthoftheα- numbers). Whereredarerunswithblackβ=2.0,β=1.0andgreen stressesreflectsthesizeandtheamplitudeofthelargest β=0.5.Theblackdashedlineshowsthedependencythatwefit. vortex. Its size is defined by H only and not by any of the other τ and β parameters. As long as τ and fitstooursimulationresults β are sufficient to replenish vorticity at the loss-rate, √ u =3 αc (13) the α-stresses should be independent of τ and β. The rms √ s loss time-scale via generation of waves and Reynolds ρ =2 αρ (14) rms 0 stresses is rather long, see Section 4.2 and Fig. 3. Thus ω2=5αΩ2. (15) aslongastheamplification-ratesarefasterthandecay- z 0 rates one should always obtain roughly the same α- We can derive the typical length-scale of angular mo- values. mentumtransportL,ofthesystemifEq. (13)isinserted intothegeneralαformalisms(Shakura&Sunyaev1973) 5. SUMMARYANDCONCLUSION ν=αc H=u Lsothat s rms Inthispaperwehaveconductedanextensiveparam- √ αH eter analysis for the baroclinic vortex amplification. In L= , (16) particular we analyzed the influence of the global en- 3 tropygradient,thermalrelaxationandcoolingaswellas indicating smaller structures than the vortices in our numericalparameterssuchasresolution, boxsize, and simulations and also smal√ler than the vorticity in stan- amplification-ratesforvorticesandsaturationvaluesof dardα-modelswhereω∝ αwithadifferentcoefficient α. (Cuzzietal.1994). Themostimportantresultofourstudyisthatwefind We do not perform a more exact analysis of these vortexgrowthevenforentropygradientsaslowasβ=