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