ebook img

Large-scale N-body simulations of the viscous overstability in Saturn's rings PDF

5.6 MB·English
Save to my drive
Quick download
Download
Most books are stored in the elastic cloud where traffic is expensive. For this reason, we have a limit on daily download.

Preview Large-scale N-body simulations of the viscous overstability in Saturn's rings

Mon.Not.R.Astron.Soc.000,000–000(0000) Printed15October2012 (MNLATEXstylefilev2.2) Large-scale N-body simulations of the viscous overstability in Saturn’s rings Hanno Rein1 and Henrik N. Latter2 1InstituteforAdvancedStudy,1EinsteinDrive,Princeton,NJ08540,e-mail:[email protected] 2 2DAMTP,UniversityofCambridge,CMS,WilberforceRoad,CambridgeCB30WA,UK 1 0 2 Submitted:18August2012. t c O ABSTRACT Wepresentresultsfromlarge-scaleparticlesimulationsoftheviscousoverstabilityinSaturn’s 1 rings. The overstability generates a variety of structure on scales covering a few hundred 1 metrestoseveralkilometres,includingaxisymmetricwavetrainsandlarger-scalemodulations. SuchpatternshavebeenobservedinSaturn’sringsbytheCassinispacecraft.Oursimulations ] P model the collisional evolution of particles in a co-rotating patch of the disk. These are the E largest N-bodysimulationsoftheviscousoverstabilityyetperformed.Theradialboxsizeis . five orders of magnitude larger than a typical particle radius, and so describes a 20-50 km h radialportionoftherings.Itsevolutionistrackedformorethan10000orbits.Inagreement p withhydrodynamics,our N-bodysimulationsrevealthattheviscousoverstabilityexhibitsa - o rich set of dynamics characterised by nonlinear travelling waves with wavelengths of a few r hundred meters. In addition, wave defects, such as sources and shocks, punctuate this bed t s of waves and break them up into large-scale divisions of radial width ∼ 5 km. We find that a the wavelength of the travelling waves is positively correlated with the mean optical depth. [ In order to assess the role of the numerical boundary conditions and also background ring 1 structure,weincludesimulationsofbroadspreadingringsandsimulationswithagradientin v the background surface density. Overall, our numerical results and approach provide a tool 8 with which to interpret Cassini occultation observations of microstructure in Saturn’s rings. 5 Wepresentanexampleofsuchasyntheticoccultationobservationanddiscusswhatfeatures 3 toexpect.Wealsomaketheentiresourcecodefreelyavailable. 3 0. Keywords: instabilities–waves–methods:numerical–planetsandsatellites:rings 1 2 1 : v 1 INTRODUCTION Fine-scalestructureoriginatesinalinearinstabilityconnected i totheviscousstressofthering,termedthe‘viscousoverstability’. X Saturn’sringsexhibitaphenomenalwealthofaxisymmetricstruc- Its fastest growing modes take the form of axisymmetric density r turethatcoversavastrangeoflengthscales,from100mto100km. wavesonwavelengths∼50m(Schmit&Tscharnuter1995).Such a On the shortest scales, ultraviolet and radio occultation observa- wavesgiverisetoperturbationsintheviscousstressthatextracten- tionsrevealthatthevariationsarequasi-periodicwithwavelengths ergyfromthebackgroundorbitalshear.Thissurplusenergyisfed between150mand220m.Thesewavetrainsintermittentlypopu- intothewave,whereitoutcompetesstandardviscousdamping,and lateregionsofopticaldepth∼1-2,suchasintheA-ringandless thusleadstorunawaygrowth.Weconcentrateinthispaperonthe opaque areas in the B-ring (Colwell et al. 2007; Thomson et al. small-scaleaxisymmetricmanifestationoftheinstability,butnote 2007;Colwelletal.2009).Inaddition,theCassinicamerashave thatthesamemechanismcanleadtogrowthofaglobaleccentric uncovereddisjunctbandsofvaryingbrightnesswithwidthsof1-10 mode in both narrow and broad rings (Borderies et al. 1985; Pa- km.Theseintermediate-scalevariationsoftencoincidewiththepe- paloizou&Lin1988;Longaretti&Rappaport1995;Ogilvie2001). riodicmicrostructure(Porcoetal.2005).Ontheotherhand,broad undulationson100kmscalespervadetheinnerB-ring,whilethe The viscous overstability has been studied in the context of centralB-ringsplitsintoalternatingzonesoflowandhighoptical continuummodels(bothhydrodynamicalandkinetic)andN-body depth,alsoonscalesof100km(Colwelletal.2009).Weseekto simulations, with the early stages of its evolution receiving the understand how and why these different features develop and the greatestattention(Schmit&Tscharnuter1995;Salo2001;Schmidt nature of any relationship between them, while concentrating on etal.2001;Saloetal.2001;Latter&Ogilvie2006a,2008).The thestructureonshortandintermediatescales. hydrodynamicalstudiesrevealthatinstabilityoccursoncethering (cid:13)c 0000RAS 2 HannoReinandHenrikN.Latter viscosityincreasessufficientlysteeplywithsurfacedensity;i.e.if (M.Hedman,privatecommunication).Adirectcomparisontonu- β ≡ dlnν/dlnτ (cid:38) 1, where ν is viscosity and τ is normal opti- mericalsimulationscouldallowforthemeasurementsofparticles caldepth.InkineticmodelsandN-bodysimulationsthisinstability properties such as size, density and coefficient of restitution. We conditioncorrespondstoτ>τ ,whereτ isacriticalopticaldepth. present the first step towards this goal and present synthetic oc- c c Thelattertheoreticalformalisms,however,yieldnoupperboundon cultationobservationsfromourN-bodysimulationsthatshowthe τ,abovewhichoverstabilitydisappears,asisindicatedbyCassini feasibilityofthisapproach. observations of the B-ring. The nonlinear saturation of the insta- Thepaperisorganisedasfollows.InSection2wedescribethe bility, on the other hand, has mainly been explored with hydro- physical model and governing equations, and then the numerical dynamical models (Schmit & Tscharnuter 1999; Schmidt & Salo methodweusetosolvethem.Section3providesabriefcomparison 2003;Latter&Ogilvie2009,2010).Themostrecentstudiesshow between our code and the kinetic equilibria and linear theory of thatthelong-termevolutionoftheringisdominatedbynonlinear Latter & Ogilvie (2008), as well as a number of numerical tests. travellingwavetrainswithwavelengths∼ 200m—shorterwaves InSection4wepresentourmainresults,beginningwithadetailed arevulnerabletosecondaryinstabilities.Thesimulatedwavetrains analysis of a large-scale fiducial simulation and then moving on undergosmallchaoticfluctuationsandcanbepunctuatedbywave to the role of the main parameters and global disk structure. We defects (sources and shocks) that partition the radial domain into summariseourresultsanddiscusstheimplicationsinSection5. structures on larger scales (cid:38) 1 km (Latter & Ogilvie 2010). The similarity between the shorter wave features and the larger scale divisions, on the one hand, and the diverse Cassini observations, 2 PHYSICALANDNUMERICALMODEL ontheother,isencouraging,thoughadditionalworkisnecessary tomakethecorrespondencetighter.Inparticular,theroleofself- 2.1 Equationsofmotion gravityhasyettobeclarified. WesolvetheequationsofmotionintheHillapproximation(Hill Inthispaperwesimulatethenonlinearsaturationofthevis- 1878)whichisalocalcoordinatesystemthatisco-rotatingwitha cous overstability with large-scale particle simulations in a local particleonacircularorbit.Thegravityfromthecentralobjectis co-rotating patch of disk (the shearing box approximation). Our linearisedinlocalco-ordinatesandtheorbitalfrequencyisacon- simulationsgeneralisetheprevioushydrodynamicalresultstothe stant.Thisallows,butdoesnotrestrict,ustoshear-periodicbound- morerealisticregimeofadensegranularflow,inwhichthering’s aryconditions.InthatcasetheHillapproximationisalsoreferred pressuretensordeviatessignificantlyfromtheNewtonianprescrip- to as the shearing sheet. The x, y, and z coordinates point in the tion(Latter&Ogilvie2006a,2008).Theysimultaneouslyextend radial,azimuthal,andverticaldirectionrespectively.Theequations previous N-body simulations to the greater length and time-scale ofmotionforatestparticlecanthenbewrittenas necessarytocapturethenonlineardynamicsproperly.Thesescales arenowcomputationallyfeasiblethankstothefreelyavailablecol- x¨=2Ωy˙+3Ω2x lisionalN-bodyREBOUND(Rein&Liu2012)whichincludesasym- y¨=−2Ωx˙ (1) plecticintegratorfortheshearingsheetandanefficientcollisionde- z¨=−Ω2z, tectionalgorithm.Weomitself-gravity;itsinfluencewillbetested z in future work. A beneficial consequence of this omission is that whereΩandΩ aretheangularvelocityandverticalepicyclicfre- z wecanshortentheazimuthalextentofthebox.Thispermitsusto quency,respectively.Thesolutiontotheseequationscanbewritten simulate boxes with radial sizes of 55 km (55000 particle radii) as epicycles (e.g. Rein & Tremaine 2011). In the case where the for some 10000 orbits, orders of magnitudes more than previous centralobjectisapointmassandself-gravityisneglected,wehave simulations. Ω = Ω. We denote the radial length of the box by L and the z x Ourresultsareingeneralqualitativeagreementwiththehy- azimuthallengthby L .Typically L (cid:28) L inthesimulationswe y y x drodynamical simulations of Latter & Ogilvie (2010). Both ap- perform. proaches yield simulations that are dominated by travelling non- The only further ingredients needed besides Eqs. (1) are the linear wavetrains of comparable morphology and wavelength (∼ finite particle size r and a collision model. We treat particles as p 100m),thoughtheparticlesimulationsexhibitaprotractedearly hardspheres(theyarenotpermittedtooverlap)andtheoutcome stageintheevolutioninwhichthecounterpropagatingwavesinter- ofacollisionisdescribedusinganormalcoefficientofrestitution actinacomplicatedway,givingrisetostandingwavefiguresand (cid:15),whichcanbeeitherconstantorafunctionoftheimpactvelocity beatingpatterns(‘waveturbulence’).Inaddition,theN-bodysim- v .Theparticleshavenospin. imp ulationsyield sourcesandshocks structuresasin hydrodynamics Inthispaperwedonottreatself-gravitydirectly.Thisisdone thoughtheirsizesandmorphologydiffer.Toinvestigatethepoten- only for simplicity and clarity. However, we increase the vertical tial role of global ring structure (and to minimise that of the pe- epicyclicfrequencyinsomesimulationsΩ =3.6Ωinordertocon- z riodic boundary conditions), we also performed simulations with centrate particles in the mid-plane, thus mimicking self-gravity’s viscouslyspreadingrings(mimickingopenboundaries)andrings vertical compression. This convenient method has been used in with a background density gradient. We stress that the structures most previous ring studies (e.g. Wisdom & Tremaine 1988; Salo that emerge in our simulations never exhibit lengthscales longer etal.2001).Itallowsustoremoveonescalefromtheproblem(i.e. thanroughly5-10km.Weconcludethattheobserved100kmfea- the particle density or Hill radius) and more easily focus on the turesintheB-ringareunlikelytobegeneratedbyviscousoversta- fundamentaldetailsofthecollectiveaxisymmetricdynamics. bility. We plan to study simulations including self-gravity in a fu- NewCassiniobservationsofstaroccultationshavearesolu- turepaper.Self-gravityisimportantbecauseitmayinstigatenon- tionthatisbelowafewhundredmeterssuchthatindividualpeaks axisymmetricwakestructuresthatcouldcompetewiththeoversta- oftheoverstability-supportedwavetrainscanbedirectlyobserved bilityandalteritssaturation(Saloetal.2001).Self-gravitywillfur- (cid:13)c 0000RAS,MNRAS000,000–000 Large-scale N-bodysimulationsoftheviscousoverstabilityinSaturn’srings 3 thermorechangethepropertiesoftheoverstabilitywavesdirectly 2.3 N-bodycodeandalgorithmsused throughtheirnonlineardispersionrelation. We use the freely available N-body code REBOUND (Rein & Liu 2012)toperformallofthesimulationspresentedinthepaper.The 2.2 Diagnostics codeisideallysuitedforthiskindofstudybecauseofbothspeed andaccuracy.Weusethehighlyefficientandaccuratemixedvari- Inordertoprobethecollectivebehaviourofthegranularflow,we ablesymplecticintegratorSEI(Rein&Tremaine2011).Collisions requireanumberofaveragedquantities.Wedefinethemeannor- are detected using a so called plane-sweep algorithm. The plane- mal geometrical optical depth τ¯ as the total projected area of the sweepalgorithmscaleslinearlywiththenumberofparticles,O(N), particles on the (x,y) plane divided by the total area of the (x,y) forhighlyelongatedboxessuchasthoseconsideredhere.Bothal- plane.Inotherwords gorithmsarealreadyimplementedinREBOUND.Theseadvantages τ¯ =Nπr2/(L L ), (2) allowustouseboxeswhichare105 timeslargerthanthetypical p x y particleradius.Wethereforedonotneedtorelyonperiodicbound- whereNisthenumberofparticles.Thus,τ¯isstipulatedatthebe- aryconditionsandcanstudyeffectsnearringboundariesandden- ginningofeachrunandwillnotchange. sitygradients. We also define the radially and temporally varying optical Onecrucialaspectinsimulatingtheviscousoverstabilitycor- depth,denotedsimplybyτ,bysubdividingtheradialdomaininto rectlyistomakesurethatthereisnopreferreddirectionintroduced thinstripsofradiallengthL : S bythenumericalscheme.Infact,thesymplecticintegratorintro- τ(x)=N πr2/(L L ), (3) duced by Rein & Tremaine (2011) ensures this symmetry is not i i p S y broken, even on the scale where machine precision becomes im- where x is the radial location and N is the number of particles i i portant.However,collisionsformallybreakthesymplecticnature inthei’thstrip.Aswetypicallybeginarunfromahomogeneous inourimplementation.Therefore,werandomisetheorderinwhich state,att = 0wehaveτ(x) = τ¯ foralli.Theensuingoverstable i collisions are resolved after each timestep. This has proven to be fluctuationswillsubsequentlybecapturedbythevariationsinτ(xi). the cleanest and most efficient way to remove all spurious corre- L ischosentobesmallerthananyscaleweareinterestedinbut S lationswhichmightotherwisebeintroducedintheorderinwhich notsosmallastobeinfluencedbyPoissonnoise. collisions are resolved. Unfortunately, this complicates the paral- Wefurtherintroducethephotometricopticaldepthτˆwhichwe lelisationondistributedmemorysystems. defineasoneminusthetransmission.Thetransmission(thefrac- Forthereader’sconvenience,welistallthesimulationspre- tionoflightpassingthroughthering)isdefinedtobe1forτˆ = 0 sentedinthispaperinTable1.Thefirstcolumndescribesthesim- and0forτˆ =1.ThisiscalculatedwithaMonteCarloRaytracing ulation.Thesecondcolumnliststheparticleradius.Thethirdand algorithm.Thealgorithmshootsraysthrougharandompointinthe fourth columns list the radial and azimuthal size of the box. The simulationboxbutwithapredefineddirection(theviewingangle). fifthandsixthcolumnslisttheinitialmeanopticaldepthτ¯andthe Then,wecheckforintersectionswithringparticlesalongtheray. numberofparticles.Theseventhcolumnliststheverticalepicyclic Thisisrepeateduntiladesiredaccuracyhasbeenreached.Unless frequencyinunitsofthestandardepicyclicfrequency.Theeighth otherwisenotedweuseaviewinganglenormaltotheringplane. columnliststhecoefficientofrestitution. Thefillingfactorisdefinedastheproportionofvolumetaken up by the particles. For spherical particles it can be defined as FF = (4π/3)nr3, where n is volumetric number density. Particu- p larlyusefulisthefillingfactoratthemidplaneFF ,whichrequires 3 COMPARISONSANDNUMERICALTESTS 0 thecalculationofthenumberdensityatz=0. 3.1 Backgroundequilibriumstate Let us further introduce the mean velocity dispersion tensor whichiscomputedfrom We begin by calculating the properties of homogenous statistical W =(cid:104)x˙x˙ (cid:105), (4) equilibriawithournumericalcodeandcomparingthemtotheki- ij i j netictheoryofAraki&Tremaine(1986)asreformulatedbyLatter where(x1,x2,x3) = (x,y,z)andtheanglebracketsindicateasuit- & Ogilvie (2008). We use a simulation box with a width of only ableaverageovertheparticlesandpossiblyovertime.Thevelocity 25 particle radii. As a consequence, growing overstability modes dispersionc2isthenWii/3.Thetranslational(local)componentof cannot fit into the domain (see also Sect. 3.3). Depending on the thekinematicviscosityissimply opticaldepth,thenumberofparticlesrangesfromN = 20to600. 2 Aconstantcoefficientofrestitution(cid:15) =0.5isused.Otherthanthe νtrans= 3Wxy/Ω. (5) particleradiusrpandthefrequenciesΩandΩz,whichwesetequal to1inourunits,thereisnootherscale. Thecollisional(nonlocal)componentoftheviscosityis In Figure 1 we plot the filling factor in the midplane FF , 0 ν = 2 (cid:88)(x −x )∆p , (6) the velocity dispersion c2, the xx, yy and xy components of the coll 3ΩMT > < y velocity dispersion tensor W , and the total kinematic viscosity ij where the sum is taken over all binary collisions that occur in a ν =ν +ν .Thesequantitiesareplottedasfunctionsofthe total coll trans timeintervalT.Here Misthetotalmassofallringparticles,∆p opticaldepthτ¯.Overall,thenumericalequilibriaareingoodagree- y denotes the transfer of y momentum from the inner to the outer mentwiththekinetictheory,witherrorsofasimilarordertopre- particleinsuchacollision,andx andx theradiallocationsofthe vious comparisons (Wisdom & Tremaine 1988; Latter & Ogilvie > < twoimpactingparticles(Wisdom&Tremaine1988;Daisakaetal. 2008).Atlargerτ¯,systematicdiscrepanciesemergewhicharethe 2001).Asweneglectself-gravitythereisnogravitationalorwake resultofthelargefillingfactorsFF ≈0.4.InthisregimetheEn- 0 contributiontotheoverallmomentumtransport. skogapproximationbeginstobreakdown,andthekineticresults (cid:13)c 0000RAS,MNRAS000,000–000 4 HannoReinandHenrikN.Latter Name/Description rp[m] Lx[m] Ly[m] τ¯ N Ωz[Ω] (cid:15) equilibrium 1 25 25 0.1-3 20-600 1 0.5 convergence 1 32-8192(1024) 2-40(5) 1.64 32-21382 3.6 0.5 fiducial 1 18641 15.53 1.64 151179 3.6 0.5 varyingtau(probingwavelengthdependence) 1 3976 15.53 0.1-2.7 1967-53097 3.6 0.5 verylowopticaldepth 1 18641 15.53 0.50 46092 3.6 0.5 lowopticaldepth 1 18641 15.53 1.00 92183 3.6 0.5 highopticaldepth 1 18641 15.53 2.00 184365 3.6 0.5 densityramp 1 18641 15.53 0.04-3.24(left-right) 143070 3.6 0.5 openboundaries 1 16641 15.53 1.64 134960 3.6 0.5 fiducialsuperwide 1 55922 5.17 1.64 151162 3.6 0.5 sizedistribution 0.5-2 18641 15.53 1.64 204178 3.6 0.5 inelasticcollisions 1 18641 15.53 1.64 151179 3.6 0.0 Bridgesetal.(1984)collisionlaw 1 18641 15.53 1.64 151179 3.6 var. Table1. Parametersofsimulationspresentedinthispaper.Seetextfordetails. 0.5 kinetic theory 3.5 kinetic theory 4.5 kinetic theory 0.45 simulation simulation simulation 3 4 0.4 filling factor FF0 000 .0.0.123..52535 nviscosity total 12.. 5512 2velocity dispersion c 23 . .5523 0.1 0.5 1.5 0.05 0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3 optical depth -t optical depth -t optical depth -t 8 kinetic theory 2 kinetic theory 0.8 kinetic theory simulation simulation 0.7 simulation 7 1.8 0.6 6 0.5 Wxx 5 Wyy 1.6 Wxy 0.4 4 1.4 0.3 3 0.2 1.2 2 0.1 1 0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3 0 0.5 1 1.5 2 2.5 3 optical depth -t optical depth -t optical depth -t Figure1.Comparisonofnumericalsimulationstokinetictheory.Fromtoplefttobottomright:fillingfactor,totalviscosity,velocitydispersion,xx,yyandxy componentsofthevelocitydispersiontensor.Seetextfordetails. are less reliable. There is also a 10% error in the xy component cousoverstabilityemergesoncetheopticaldepthincreasesabove oftheviscousstressatlowerτ¯,whichhasbeenpreviouslynoted acriticalvalue.Thisτ¯ dependenceisactuallyaproxyforcentral (Latter&Ogilvie2008).Theseresultsshowthatourparticlecode fillingfactorFF ,whichmoredirectlyinfluencesthebehaviourof 0 correctlyreproducesestablishedkineticringequilibria. thering’sviscousstress.AsFF increases,sodoesthesensitivity 0 ofthestresstodensityvariations,asmediatedbytheEnskogfac- torinthekinetictheory(Araki&Tremaine1986;Latter&Ogilvie 3.2 Linearstability 2008).Ifthestressincreasessufficientlyfastwithdensity(i.e.ifβ issufficientlylarge)theninstabilityensues(Schmit&Tscharnuter Next we examine how well the code treats the linear stability of 1995;Schmidtetal.2001). theequilibriumstates.Insteadofexaminingthegrowthratesofthe viscous overstability modes, as in Schmidt et al. (2001), we de- Wecomputethenumericalstabilitycriterioninthefollowing scribethecurveofmarginalstabilityinthefullparameterspaceof way. A grid of τ¯ and (cid:15) values is set up and a sequence of simu- (cid:15)andτ¯.ThisresultextendspreviousN-bodyworkthattestedonly lationsrun,eachcorrespondingtoaunique(τ¯,(cid:15)).Thesimulations ahandfulofcases(Saloetal.2001),andcandirectlyconnecttothe use L = 1000mandareevolvedforupto104 orbits,whichen- x predictionsofthekineticmodelling(Latter&Ogilvie2008). suresthatthereissufficienttimefortheunstablemodesinthebox For given collisional properties, it has been shown that vis- toachieveappreciableamplitudes.Thelevelofactivityinthebox (cid:13)c 0000RAS,MNRAS000,000–000 Large-scale N-bodysimulationsoftheviscousoverstabilityinSaturn’srings 5 0.01 gives the overstability enough time to develop. The power spec- 0.9 tra of the radial variations in the density, the average height of a 0.8 particleabovethemid-planeandtheaveragevelocityofaparticle en 0.7 arethenmeasured100timesperorbitandaveraged.Oneofthese o uti 0.6 0.001 quantities, the power spectrum of the density, is plotted in Fig. 3 nt of restit 00..45 -2tt>-<() fsoimrdiliIaffnr.eFreign.t3d(ta,)Lwx aenfidxLLy.T=he1o0t2h4ermspaencdtrLaar=eq5umaliatantdivvealyryvethrye e x y effici 0.3 0.0001 timestep.Onecanseethatatimestepdt=2π·10−3Ω−1orsmaller co 0.2 issufficienttocapturetheoverstabilityaccurately.FromFig.3(b), 0.1 whereallparametersexcepttheazimuthalboxwidtharefixed,it isclearthatL < 4msuppressestheoverstability.ForL (cid:62) 10m 0 y y 1e-05 thespectraareindependentofL .Thustheboxsizeissufficiently 0.5 1 1.5 2 2.5 y optical depth -t largetocapturetheoverstabilityaccurately.Finally,byinspecting Fig.3(c)onecanseethataradialboxwidthofL =1024missuf- x Figure2. Deviationfromtheinitialopticaldepthasafunctionoftheinitial ficienttoaccuratelycapturethedominantwavelengthsinthenon- geometricopticaldepthandthecoefficientofrestitution.Thisisaproxyfor linearsaturatedstate. theactivityoftheviscousoverstability. Tosummarise,thedevelopmentoftheviscousoverstabilityin thesesmallerrunsisrelativelyinsensitivetothenumericalparam- eters.Ofcoursethisdoesnottelluswhetherthelarge-scalenonlin- ismeasuredby(cid:104)(τ(t)−τ¯)2(cid:105),wheretheanglebracketsdenoteabox earbehaviouriscapturedadequately,butitnonethelessprovidesa average.InordertoobtainsufficientlylargeFF0 forinstabilityat foundationofconfidencethatthisisindeedthecase. areasonableτ¯weamplifyΩzto3.6Ω.Doingsomimicsthecom- Weshouldalsopointoutthatinsimulationswithhigheropti- pressionofthediskviaself-gravitation(seeabove).Theresultof caldepth,τ¯ (cid:62)2,weneedalargerL toavoidcrystallisationofpar- y thesesimulationsareplottedinFig.2,inwhichweshowthemag- ticles.Quasi-crystallinestructuresarealowerenergystateandact nitude of the activity at the end of the run. Superimposed on the asanattractorifthevelocitydispersionistoosmalltorandomise figure in red is the stability curve calculated from the dense gas theparticlemotion.Higherdensityringsaremorepronetothisef- kinetictheory(Latter&Ogilvie2008). fectbecauseahighercollisionrateleadstoacooleranddenserring Asanyactivitycorrespondstotheonsetofviscousoverstabil- withstrongexcludedvolumeeffectsandparticlecorrelations.This ity,thepurelyblackregionsinFig.2signifystability.Allnon-black ispurelyanartefactofusingequal-sizedparticlesandwillnotoc- regions indicate instability, except for a handful of outliers influ- curinarealring1.Wecouldsimplyallowforasizedistributionto encedbynoise.Clearly,forgiven(cid:15)lessthanroughly0.9thereisa avoidthisproblem(wepresentonesuchsimulationinSect.4.2). criticalτ¯abovewhichwehaveinstability,inagreementwithprevi- However,thiswouldcomplicatetheanalysisandcomparisontoan- ousresults.Larger(cid:15) generallyrequiresalargerτ¯.Thisisbecause alyticalestimates. thevelocitydispersionanddiskthicknesstendtoincreasewith(cid:15), therebeinglessdissipation.Asaconsequence,greaterfillingfac- torsaremoredifficulttoachieveatloweropticaldepths.Therealso appearstobeahardcut-offaroundτ¯ =0.5belowwhichinstability 4 NUMERICALRESULTS canneveroccur. There is general qualitative agreement with the kinetic Inthissectionwepresentourmainresultsonthenonlinearsatu- marginalcurve.Thequantitativediscrepanciesarisebecauseinsta- rationoftheviscousoverstabilityinlargeradialdomainsandover bility occurs for larger FF0 (above roughly 0.35), the regime in longtimes.Weranasignificantnumberofsimulationsvaryingour whichthekinetictheoryislessreliable.Forlarger(cid:15) thisisespe- twomaincontrolparametersτ¯and(cid:15),butwewilldescribeonlyone ciallythecase,ashereFF0approaches0.5soastoforceinstability setofparametersindetail.This‘fiducialsimulation’conveniently intheselessdissipativeconditions. exhibitsallthemainfeatureswitnessedintheotherruns.Inaddi- In summary, both the equilibrium and linear stability results tion,wetestedhowdifferentglobalstructuresimpactedonthelong showthatourparticlecodeisbehavinginthecorrectmanner,with termbehaviour.Specifically,simulationsareshownofaspreading regardstopreviouscalculations.Thisgivesusconfidenceinitsca- ringandaringwithadensitygradient. pabilitytodescribingthenonlinearevolutionoftheviscousover- stabilitywhichisthemainfocusofthispaper. 4.1 Fiducialsimulation 3.3 Numericalconvergence Ourillustrativerunemploysanopticaldepthofτ¯ = 1.64.Weuse shearing sheet boundary conditions with an azimuthal and radial We verify how the results of the simulations depend on the size oftheboxandthetimestep.Thesearetheonlytwonumericalpa- rametersweneedtocheck.Inthesesimulationsthegeometricop- ticaldepthisτ¯ = 1.64andweuseaconstantcoefficientofresti- 1 Analogous crystallisation phenomena have been hypothesised as a re- tution(cid:15) = 0.5.AswehaveshowninSect.3.2,thissystemwillbe sult of the (poorly constrained) surface adhesion between ring particles viscouslyoverstable. byTremaine(2003).SeealsoPerrineetal.(2011);Perrine&Richardson We integrate the system for approximately 320 orbits which (2012). (cid:13)c 0000RAS,MNRAS000,000–000 6 HannoReinandHenrikN.Latter 10 10 10 1 1 1 density (power) 0 .00.11 density (power) 0 .00.11 density (power) 0.1 0.001 0.01 0.001 0.0001 100 1000 0.0001 0.001 wavelength 100 1000 100 1000 dt = 5(cid:215)10-1 dt = 2(cid:215)10-2 dt = 1(cid:215)10-3 wavelength wavelength dddttt === 215(cid:215)(cid:215)(cid:215)111000---112 dddttt === 152(cid:215)(cid:215)(cid:215)111000---233 dddttt === 521(cid:215)(cid:215)(cid:215)111000---444 LLLyyy === 234 mmm LLyyL y== = 11 055 mmm LLyy == 2400 mm LxLL xx= == 1 236824 mmm LxLL xx= == 1 025251462 mmm LLLxxx === 248001499862 mmm (a)Densitypowerspectrafordifferenttimesteps (b)Densitypowerspectrafordifferentazimuthal (c)Densitypowerspectrafordifferentradialbox dt.Theunitofthetimestepis2πΩ−1. boxwidthsLy. widthsLx. Figure3.Convergencestudy.Powerspectraofthedensityfordifferenttimestepsandboxsizes. thevelocitydispersioncaftert = 8000orbits.Wefindthatthefi- 800 nalstablenonlinearwavetrainstatepropagatesinasingledirection. y g Thefirstplotconfirmsthatthenonlinearevolutionremainsaxisym- er 400 en metric(withinthelimitedazimuthalextentofthebox),whilethe τ profile exhibits the characteristic cusp-like shape shared by all 0 0 1000 2000 3000 4000 5000 6000 7000 8000 densitywavesinthindisks,whereparticlemotionsarecontrolled time [orb] by epicycles. The system has yet to completely relax to a wave- trainofuniformamplitude,anddisplaysalarge-scalemodulation Figure4.Kineticenergyevolution(onlyxandycomponents)inthefiducial in wavelength. Ultimately this modulation decays. Note that the simulation. wavelengthvariationcorrespondstoanamplitudevariationinthe velocityfield,asshownbythehydrodynamicaltheory.Theaver- agewavelengthexhibitedissome450m.Thevelocitydispersion widthofL =15.53mandL =12.43km,respectively.Allparti- y x profile reveals spikes in the random motions at wavetrain peaks. clesareinitiallyrandomlydistributed.Thesimulationwasthenin- Itisinthesehighdensitypeaksthattheorbitalenergyisprimar- tegratedfor8200orbitswhichcorrespondstoover10Earth-years, ily thermalised, as the collision frequency (and associated colli- assumingadistancefromSaturnofa=130000km. sionalstresses)aremaximisedthere.The(x,z)cross-sectionshows Viscousoverstabilityattainsnonlinearamplitudesafterafew theseregionsofhigh-densityareaccompaniedbycoherentvertical hundred orbits, and then settles into a disordered state for some expansions (splashing) via excluded volume effects. This vertical 5000orbits.Afterroughly5300orbitsthisstateresolvesintoasin- structureisconnectedtothevertical‘breathingmode’exhibitedby gletravellingwavetrainofrelativelyuniformamplitude.Thefirst viscousoverstabilityingaseousdisks(Latter&Ogilvie2006b). disordered stage of the evolution can be arranged into two con- secutivephases,(a)a‘waveturbulence’,characterisedbynonlin- Thesource/shockstateisdescribedinFig.5(b),whichissim- earlyinteractingtravellingwaves,standingwaves,andbeatingpat- ilar to the previous figure but omits the v and c profiles. The τ y terns,and(b)acellularsource/shockstate,inwhichrelativelywell- profileshowstheexistenceoftwopatchesofdensitywavesofvary- definedpatchesofcounterpropagatingwavesareseparatedbywave ingwavelengthwiththeirdirectionofpropagationindicatedbyar- defects(shocksandsources).InFig.4wepresentthetotalfluctu- rows.Thetwowavetrainsareseparatedbyasourceregionlocated atingkineticenergyoftheboxasafunctionoftime.Initthethree atroughlyx=0m,whichgenerateswaves,andashockregionlo- phases of the evolution can be distinguished. From roughly 500 catedaroundx=−9000m,whichiswherethecounterpropagating orbits to 2000 orbits, there is a linear increase of energy associ- waves collide and destroy each other. The shock region straddles atedwithstage(a),the‘waveturbulence’.Thenfromorbits2000 the(periodic)radialboundary.Thisconfigurationofdynamicalel- to5300theenergysaturates,whileundergoinglong-termoscilla- ementsreflectsbehaviourinthecomplexGinzburg-Landauequa- tions. This corresponds to phase (b), the source/shock state. Af- tion,whichservesasagoodreducedmodelfortheoverstablesys- ter5300orbitsthesystemrelaxestoasingletravellingwavetrain tem(Aranson&Kramer2002;Latter&Ogilvie2009).Italsore- thatfillsuptheentireboxandyieldsasingleuniformenergy.This producesqualitativelythehydrodynamicalsimulationsofLatter& is the end-state of the evolution and is shared by all simulations Ogilvie(2010),yetthereexistinterestingdiscrepancies.Thewave- thatusedshear-periodicradialboundaries.Wenowdiscusseachof trainsthemselvesarefarmoredisorderedandundergolargeram- thesethreephasesinmoredetail. plitude fluctuations in their amplitude and phase. In addition, the We begin in reverse with the final state, which we plot in sourceandshockregionsarebroaderandlessdefined.Inparticular, Fig.5(a).Thisfigureshowsasnapshotoftheparticlepositionsin the shock region can extend over 1000-1500 km, i.e. many char- the(x,y)planeanda(x,z)cross-section,theopticaldepthτ,they- acteristic wavelengths of the wavetrain. This is because colliding andz-averagedazimuthalvelocities(relativetotheshear)v ,and wavetrains penetrate each other more successfully in the particle y (cid:13)c 0000RAS,MNRAS000,000–000 Large-scale N-bodysimulationsoftheviscousoverstabilityinSaturn’srings 7 5 ]y [rp 0 -5 4 2 ]z [rp - 20 -4 (cid:27) (cid:27) (cid:27) 5 h pt 4 e optical d 123 0 20 ] W [rvpy-2 00 4 ] W [rcps 2 0 -8000 -6000 -4000 -2000 0 2000 4000 6000 8000 radius x [rp] (a) Fiducialsimulationaftert=8000orbits. 5 ]y [rp 0 -5 4 2 ]z [rp - 20 -4 h 5 shock (cid:27) (cid:27) source (cid:45) (cid:45) pt 4 e optical d 123 0 -8000 -6000 -4000 -2000 0 2000 4000 6000 8000 radius x [rp] (b) Fiducialsimulationaftert=4000orbits. 5 ]y [rp 0 -5 4 2 ]z [rp - 20 -4 h 4 ept 3 d al 2 ptic 1 o 0 10 ] W [rvpy-1 00 ] 4 W [rcps 2 0 -8000 -6000 -4000 -2000 0 2000 4000 6000 8000 radius x [rp] (c) Fiducialsimulationaftert=300orbits. Figure5.Snapshotofthefiducialsimulationatdifferenttimes. (cid:13)c 0000RAS,MNRAS000,000–000 8 HannoReinandHenrikN.Latter simulationsandmustundergomultiple‘collisions’beforetheydis- sipate.Bycontrast,thehydrosimulationssupportshocksthatare 8000 verynarrowandwell-defined:wavesarealmosttotallydestroyed in a single collision. This difference in the nonlinear interaction betweentravellingwavetrainsisperhapsthekeydisagreementbe- tweenthehydrodynamicalandN-bodymodelling. InFigure5(c)weplotasnapshotoftheevolutionduringthe 7000 5 initialdisorderedwavestate.Thedynamicsinthisphaseiscompli- cated,comprisingstretchesofcounterpropagatingtravellingwaves interpenetratingeachotherandgivingrisetovariousnonlinearin- teractionssuchasbeatingandstandingwavepatterns.Itispossi- 6000 blethatsourceandshockregionsexistinthewavemaelstrombut theyaresoradiallyextendedthattheyoverlapsignificantlyanddo 4 not function in the same way as in the shock/source stage. As is clearbycomparisonwithFig.5(b),characteristicwavelengthsare much shorter in this evolutionary stage, closer to the wavelength 5000 of the most unstable mode (∼ 100m), and they steadily increase withtime(Schmit&Tscharnuter1999).This,infact,drivesthein- creaseinkineticenergywitnessedinFig.4,asthewaveamplitude h is correlated with wavelength.. This stage of the evolution bears orb] 3dept similaritieswiththeone-dimensional‘waveturbulence’thatchar- e [ 4000 al m c acterisescertainplasmas,nonlinearoptics,andotherphysicalphe- ti pti o nomena (Nazarenko 2011, for e.g.). Notably, the overstable ring shares with those systems an inverse cascade of energy to longer lengthsandageneralmovementtowardsgreaternonlinearityand coherentstructure. 3000 Thestagesofthedynamicalevolutioncanbesummarisedin 2 a stroboscopic space-time diagram, Fig. 6, in which the optical depth is plotted as a function of x and t, but only at every 10th orbit.Thisscreensoutthefasttimeassociatedwiththewavephase 2000 speedyetkeepsthesenseoftheirandthegroupvelocity’sdirec- tion, while also indicating important features such as the sources and shocks. The first clear signal of a source region emerges at 1 aroundt = 1800orbitsnear x = −5000m.Itcanbedistinguished bythelightdiagonalrayspropagatingawayoneitherside.These 1000 raysrepresenttheslowergroupvelocityofthewaves(notthevery rapidphasevelocity∼ 100mperorbit).Theshockregionisless easytopickout.Att=2000orbitsitextendsfromx=−100mto x=1000m.Asillustratedbyhydrodynamicalsimulations,source andshockregionsmigratearoundthedomainonaslowtime-scale. 0 0 Eventuallytheycollideandannihilateeachother.Thisoccursnear -8000 0 8000 t = 5300orbits.Afterthispointwehaveasinglewavetrainprop- x [r] p agatinginonedirection(alltherayspointinthesamedirection). Anotherinterestingfeatureisthatthewavelengthofthewavesgen- eratedbythesourcevarywithtimeandcanbedifferentdepending Figure6.Stroboscopicspace-timeplotofthefiducialsimulation. onwhichsideofthesourcetheyoriginate.Thetime-variablebe- haviourwasnotedinthehydrodynamicalsimulationsbutwasless prominent; the source asymmetry appears a novel feature of the hibitthedynamicsexhibitedbysimulationsonintermediatetimes, N-bodysimulations. before the influence of the boundary conditions dominate — in Inalltherunsundertakenwithshearingperiodicboundaries other words, before, the box realises that it is periodic. Viscous thelongtermsaturationofoverstabilityisasingletravellingwave- overstabilityintherealringswillprobablymanifestaseither(a)the train of uniform amplitude. However, this final state is to some disorderedwavestateor(b)thesource/shockstate.Perhapsdiffer- some extent an artefact of the shearing box’s translational sym- entregionsintheringssupportoneortheotherstate,givingriseto metry in radius. As far as the collective axisymmetric dynamics differingobservationsonintermediatescales1−10km.Asrevealed are concerned, there is no preferred radius. A real ring, on the bytheCassinicameras,thereexistsarichandirregularsetofradial other hand, exhibits radial structure in particle size, composition, variationsonthesescales(Porcoetal.2005).InSection4.3wetrial andnumberdensityonintermediatetolongscales.Moreover,res- otherglobalset-upsandboundaryconditionsinordertobreakthe onancesandwakescreatespatiallyfixedperturbations.Thesevari- radial translational symmetry and to assess how the overstability ationsshouldprohibitthetransitiontoasmoothandcoherentstate saturatesinthepresenceofsimpleglobalstructures. withonewave-train.Itisthereforelikelythattherealringswillex- Weranoneadditionalsimulationwiththesameparametersas (cid:13)c 0000RAS,MNRAS000,000–000 Large-scale N-bodysimulationsoftheviscousoverstabilityinSaturn’srings 9 2.5 0.1 To develop this idea we subsequently ran a sequence of smaller scale simulations for different τ¯, in which L = 4000m. x 0.01 Each simulation was run for 5000 orbits. Because of the smaller 2 boxsizes,afterthislengthoftimethesystemhasresolvedintoa -th 0.001 singlestablewavetrain.Fig.7summarisestheseresultsbyplotting ept 1.5 the power in each wavelength for each τ¯. As is clear, the wave- al d 0.0001 lengthofmostpowerincreaseslinearlywithopticaldepth,begin- c pti 1 ningat210matτ¯ = 1andrisingtoroughly400matτ¯ = 2.Note o 1e-05 thattheviscousoverstabilityisnotactiveforτ¯ < 0.75,consistent withtheresultsinSect.3.2.Acorrelationofwavelengthwithopti- 0.5 1e-06 caldepthmighthavebeenobservedinrecentCassiniobservations (M.Hedman,privatecommunication). 1e-07 Allsimulationspresentedsofaruseasinglesizeofringpar- 100 200 300 400 500 600 ticles. We also ran an additional simulation using a size distribu- wave length [m] tion.ThedistributionhasaslopeofdN/dr ∝ r−3 andr ranges p p p Figure7.Power(arbitraryunits)asafunctionofwavelength(horizontal from 0.5 m to 2 m, motivated by current estimates of the typical axis)andtheopticaldepthτ¯ (verticalaxis).Eachrowcorrespondstoone sizedistributioninSaturn’srings(Cuzzietal.2009).Thesystemis simulation. viscouslyoverstableandthenon-linearbehaviourisindistinguish- ablefromthefiducialsimulation.Asanillustration,wepresenta 3Drenderingofasmallpartofthe(x,z)planeinFig.8.Theplot thefiducialsimulationbutwithanextremelylargeradialdomain, showsapproximatelytwowavelengths.Clearlyvisibleistheverti- spanning over 55000 particle radii (55 km). More details of the calsplashingofparticles.Thesesimulationsmakeusconfidentthat largersimulation,aswellasastroboscopicspace-timediagram,are our results are general and are not an artefact of the single sized presentedinAppendixA.Asisdiscussedthere,preciselythesame ringparticles.Theremightbequantitativedifferencesthatdepend phenomenaappearsonsimilarscalesinthelargebox.Thisshows ontheprecisenatureofthesizedistribution.However,ourresults theboxindependenceofourresults.Italsosuggeststhatthelongest indicate that these effects will be small and we therefore do not scalesgeneratedbyviscousoverstabilitydonotexceed5-10km. investigatethemanyfurtherinthispaper. 4.2 Explorationofdifferentparameters 4.3 Explorationofdifferentglobalstructures In order to test the robustness of the results presented in the pre- Inthissubsectionweextendtheshearingboxmodeltoaccountfor vious subsection, we conducted simulations with different (cid:15) and different global structures: a large spreading ring and a ring with τ¯. In particular, we used both (cid:15) = 0 (no bouncing) and the ve- aradialdensitygradient.Wethencheckhowtheseinfluencedthe locitydependantBridgesetal.(1984)collisionlawwhilekeeping long term saturation of the overstability. The large extent of our τ=1.64.Neitherresultedinnonlinearevolutionssignificantlydif- simulations(coveringscalesseparatedbyalmost5ordersofmag- ferenttothatofthefiducialrun.Thismakesusconfidentthat,for nitude—1mto50km)allowsaconvincingexplorationofglobal sufficientlyinelasticcollisions,theprecisetreatmentofthecoeffi- structureindirectN-bodysimulations. cientofrestitutionisunimportantandthatavalueof(cid:15) =0.5yields representativeresults.Itislikelythatlarger(cid:15)mayleadtodifferent 4.3.1 Openboundaryconditions nonlinearbehaviour,astheringwillgenerallybe‘hotter’andless dense. However, experimental and theoretical evidence points to Byopenboundaryconditionswemeanaringthathasafinitera- verydissipativecollisions(Porcoetal.2008;Schmidtetal.2009), dialwidthlessthantheradialboxsizeL .Thustheringcanspread x sothishotregimeisprobablyirrelevanttoconditionsinSaturn’s viscously.Theboundariesare,infact,notreallyopen:ifparticles rings. exittheboxononeside,theyentertheboxontheotherside.How- Weundertookadditionallarge-scalesimulationswithtwodif- ever,weinitiallyplacetheboxboundariesfarfromthediskedge ferent τ¯, equal to 1 and 2. In each case (cid:15) = 0.5. In neither case, (2000m)andstopthesimulationwhentheringhasspreadtothe wasthegeneralstructureofthenonlinearevolutiondifferenttothe pointthatparticlescrosstheboundary. fiducialcase.Itshouldbenoted,however,thatthelowerτ¯ runex- We present a space-time plot of such a simulation in Fig- hibited less violent chaotic fluctuations and more ordered wave- ure9(a).Onecanseethattheringedgesareindeedspreadingvis- √ trainsduringthesource/shockphase.Whatthesesimulationsalso couslyasexpected(i.e.∼ t).Thedevelopmentoftheoverstability revealedwasthatthecharacteristicwavelengthofthedynamics(in inthemiddleoftheringisinitiallyunchangedbytheseboundary whichthepowerwasconcentrated)positivelycorrelatedwithopti- conditions:theinstabilityremainslocalonearlyandintermediate caldepth.Lowerτ¯ runsyieldshorterwavelengthbehaviour,while times.Thisisaidedbythefactthatoverstability-generatedwave- higherτ¯ givesalongercharacteristicwavelength.Thiscorrelation trainsarenotreflectedbytheringedges.Waves,oncetheyreachthe wasalsonotedindirectlyinLatter&Ogilvie(2009)intheirTable edge,arelosttothesystemandsothespreadingedgebehaveslike 3,whichshowsthatthewavelengthofthefirststablewavetrainin- anopenboundaryorthebufferedboundariesemployedinLatter& creasedwithτ¯ (asmediatedbythedimensionlessparameterβ).It Ogilvie(2010).Neart = 500orbitstwosourceshavemanifested; islikelythatinthedisorderednonlinearstate,thisfirststablewave- oneat x = −4000mandtheotherat x = 4000m.After1200or- lengthworksasanattractorandfixesthecharacteristicwavelength bits,afteracomplicatedinteraction,theinterveningshockregion onsufficientlylongtimes. annihilatesitselfandtheleftmostsourceleavingonlytherightmost (cid:13)c 0000RAS,MNRAS000,000–000 10 HannoReinandHenrikN.Latter Figure8.3Drenderingthatshowsasmallpartofthe(x,z)planeinthesimulationusingaparticlesizedistribution. 5000 7 5 5000 6 4000 4 4000 5 3000 4 3 3000 h h orb] dept orb] dept e [ al e [ al m c m c ti pti ti pti o o 3 2000 2 2000 2 1000 1 1000 1 0 0 0 0 -8000 0 8000 -8000 0 8000 x [r] x [r] p p (a) Stroboscopic space-time diagram of a spreading ring with open (b) Stroboscopicspace-timediagramofaringwithaninitialsurface boundaryconditions. densitygradient. Figure9.Stroboscopicspace-timediagramofdifferentglobalstructures. source.Thisremainingsourcestaysatthislocationsfortherestof absenceofanyotherwavedefect,whichwouldeventuallydestroy the simulation, as in the hydrodynamical simulations of buffered thesource.Howeverstatic,thesourcestillexhibitsinterestingdy- boundaries. namics. The wavelength of the waves it generates can vary, as at Thispersistenceofasinglesourceisthekeydifferencetosim- t = 1400orbits when the wavelength drops abruptly to a much shortervalue.Thecauseofthisbehaviourisunclear,asisitscon- ulations using the period boundary conditions. This is due to the (cid:13)c 0000RAS,MNRAS000,000–000

See more

The list of books you might like

Most books are stored in the elastic cloud where traffic is expensive. For this reason, we have a limit on daily download.