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

5.6 MB·English

by  Hanno Rein

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