Collisional Simulations of Neptune's Ring Arcs

ICARUS

126, 1–27 (1997) IS965613

ARTICLE NO.

Collisional Simulations of Neptune’s Ring Arcs J. HA¨ NNINEN AND C. PORCO Lunar and Planetary Laboratory, University of Arizona, Tucson, Arizona 85721–0092 E-mail: [email protected] Received November 10, 1995; revised August 28, 1996

merical test particle simulations, and the exact positions of the arcs and their confining resonances were revised by 0.2 km (Horanyi and Porco 1993). Visible photometry derived from Voyager imaging data indicated that the arcs consist predominately of dust-sized (i.e., radii s p 1em) particles comparable in abundance to Saturn’s F ring (Porco et al. 1995). It was suggested (Smith et al. 1989) and then numerically verified (Colwell and Esposito 1990) that collisions among putative bigger particles are the dominant source of dust in the narrow Neptune rings and arcs. [Impact velocities of p1 m/sec are required to reproduce the observed dust fractions. Velocities of this magnitude are expected from the Porco (1991) model.] However, collisional velocities of order 1 m/sec should be high enough to remove the colliding (i.e., source) particles from the arcs (Porco 1991). If the arc particles are indeed undergoing disruptive collisions, as they must be to explain the abundance of dust, then why are the arcs still there? The point of this study is to investigate the collisional environment within the arcs and address two questions: (1) Do collisions occur within the arcs that are energetic enough to produce dust? (2) If so, what mechanisms, if any, keep the arcs stable against disruption by these collisions over long time scales? To answer these questions, we have modified a numerical N-body simulation code used previously to study planetary rings (e.g., Ha¨nninen and Salo 1992, 1994) to reproduce the Neptunian ring arcs/ Galatea geometry. Because of computational constraints, we have concentrated in this paper on a dynamical study of particles ranging in size from p10 to p250 m. (For such particles, radiation pressure and other nongravitational perturbations are insignificant and their dynamics can be modeled with gravitational forces and interparticle impacts alone.) We save to a later paper the problem of arc dust creation and dynamics.

The currently accepted model for Neptune arc confinement relies on the radial and azimuthal confining perturbations due to the nearby satellite, Galatea. This model calls for arc particle orbits exhibiting a negative eccentricity gradient and crossing at quadrature, a configuration that paradoxically leads to collisions energetic enough to disrupt arc confinement. We confirm with numerical collisional N-body simulations that the confinement mechanism relying on a 42:43 corotation–inclination resonance and a 42:43 outer Lindblad resonance with Galatea is indeed capable of confining a large population of 10-m-size and bigger particles over short time scales. Moreover, we find that an 84:86 outer vertical resonance, also due to Galatea, falling within 20 m of the arcs’ radial position, effectively reduces the collision frequency and relative collisional velocities and consequently stabilizes the arcs over long time scales against the disruptive effects of collisions.  1997 Academic Press

1. INTRODUCTION

Earth-based stellar occultation observations during the 1980s led to expectations of arclike rings around Neptune, which were confirmed in August 1989 by the Voyager 2 flyby of Neptune. Considerable interest in possible dynamical confining mechanisms of these features was generated on their discovery. For example, it was suggested that ring arcs might be azimuthally confined at the L4 or L5 Lagrangian point of one satellite, while another satellite on a different orbit could confine the ring material in the radial direction (Lissauer 1985). Another model invoked corotation and Lindblad resonances of the same satellite for azimuthal and radial confinement, respectively (Goldreich et al. 1986). By analyzing Voyager 2 imaging and stellar occulation data, Porco (1991) found that the geometry and kinematic behavior of the arcs were best explained by forcing due to the Voyager-discovered nearby inner Neptunian satellite Galatea. She proposed the arcs to be confined to particular longitudinal locations defined by Galatea’s 42:43 corotation–inclination resonance; radial confinement was provided by its nearby 42:43 outer Lindblad resonance. The shepherding mechanism was later confirmed by simple nu-

2. SIMULATION METHOD

In our direct N-body simulation method Aarseth’s N-body integrator (Aarseth 1972, 1985) is used. The integrator uses fourth-order force polynomial and individual 1 0019-1035/97 $25.00 Copyright  1997 by Academic Press All rights of reproduction in any form reserved.

¨ NNINEN AND PORCO HA

2

time-step scheme in the integration of particle orbits. The present code is basically the same as used by Ha¨nninen and Salo (1992, 1994), except for the one important modification. The essential requirement for the Neptune ring arc mechanism, as proposed by Porco (1991), is displacement of resonance locations due to the oblateness of Neptune. To be able to simulate system of ring arcs, gravitational harmonics J2 and J4 have been embedded into the calculation of planet’s gravitational force. In the simulation code, ring particles, satellites, and the central mass are all treated similarly. When the masses of ring particles and satellites are denoted by mi and MSj , and the mass of the central body by Mp , while ri stands for a particle position, RSj and Rp for positions of a satellite and the central body, the equations of motion for the ring particles in the inertial frame are given by

r¨i 5 2G

O

2G 2 J4 2G 1 J4

k51 k ?i

mk(ri 2 rk ) 2G uri 2 rk u 3

S S DS S S S DS S

Mp (ri 2 Rp ) uri 2 Rp u

rp Dri

3

4

70

4

i

Sj

i

Sj

4

2

2 15

2

15 Dzi 2 Dri

2

Dzi Dri

2

1

2

2

3 2

15 8

2 J4

2

15 2

,

where G is the gravitational constant, rp is the equatorial radius of the planet, and Nparticle and Nmoon denote the number of particles and satellites. We have also used Dri 5 uri 2 Rp u and Dzi 5 uzi 2 Zp u, where zi and Zp are vertical components of a particle and central body positions, respectively. For the satellites and the central mass the equations are

¨ S 5 2G R j

2G

2 J4

O

Nparticle

mi(RSj 2 ri )

i51

uRSj 2 ri u

3

O MuR(R2 2R Ru )

Nmoon

2G

Sl

l51 l?j

Sj

Sj

uRSj 2 Rp u 3

rp DRSj

4

70

1 2 J2

DZSj

DRSj

4

rp DRSj

2 15

2

15 DZSj 2 DRSj

DZSj

DRSj

2

1

2

15 2

,

rp Dri

S

(2)

mi(Rp 2 ri ) rp 1 2 J2 3 uRp 2 ri u Dri 4

70

mi(Zp 2 zi ) uRp 2 ri u

i51

rp Dri

rp DRSj

1 J4

j51

4

70

2

4

2

2

DZSj

4

DRSj

uRp 2 RSj u 3 4

rp Dri

2G

35 DZSj 2 DRSj

2

2

2 15

3J2

2

2

1

2

2

3 2

15 8

2

15 2

15 DZSj 2 DRSj

MSj(Zp 2 ZSj )

rp DRSj

Dzi Dri

2 15

3J2

3

35 Dzi 2 Dri

4

rp DRSj

Nnoon

2G

Dzi Dri

15 Dzi 2 Dri

2

2

MSj(Rp 2 RSj )

j51

uRp 2 RSj u 3

3 2

DZSj

DRSj

rp DRSj

15 2

Nmoon

2

1

15 8

2

,

(3)

where we have used DRSj 5 uRSj 2 Rp u and DZSj 5 uZSj 2 Zp u. A softening parameter is also used in the force calculations, typically set equal to p0.6 m, i.e., much less than the particle radius. Actually, the softening parameter is not needed in the simulations of massless particles. Even though the code is thus capable of simulating self-gravitating particles, we are only using massless particles in the following simulations. Our impact model is the standard one: colliding bodies are assumed to be spherical and frictionless. Hence, the perpendicular component of the relative velocity is reversed and reduced by a factor of « in each impact, (4)

3

S S DS S D D S D S S D S D DD

Mp (RSj 2 Rp )

2

c ? Vpost 5 2«c ? Vpre ,

Sl

Sl

i51

Nparticle

2G

2 J4

(1)

35 DZSj 2 DRSj

4

2

S DS S D D S D S S D S D DD O S S D S D S S D DD O S S DS S D D S D S S D S D DD O S S D S D S S D DD O

Nparticle

¨ p 5 2G R

1 2 J2

3

S DS S D D D S D DD S D D DD

Dzi Dri

35 Dzi 2 Dri

j

j5 1

rp Dri

1 2 J2

Mp(zi 2 Zp) rp 3J2 uri 2 Rp u 3 Dri

rp Dri

O mur(r22RRu )

Nmoon

rp DRSj

rp DRSj

3J2

uRSj 2 Rp u 3

1 J4

1 J4 Nparticle

S S D S D S S D DD

Mp(ZSj 2 Zp)

2G

15 8

2

2

3 2

where Vpost and Vpre are the post- and precollisional relative velocities, and c is the unit vector in the direction joining the particle centers. For the coefficient of restitution, the Bridges et al. (1984) velocity-dependent « is used. It can be written in the form «(v) 5 (vn /vc )20.234,

(5)

3

SIMULATIONS OF NEPTUNE’S RING ARCS

where vn 5 uc ? Vpreu is the perpendicular component of the relative velocity, and the constant vc has value vc 5 0.01 cm/sec. In our simulations we have used a slightly different value, vc 5 0.063 cm/sec. We will later discuss how the different values of vc affect the stability of the arcs. The collisional change in particle i’s velocity and the conservation of momentum lead to the expression vpost 2 vpre 5

mj (1 1 «)Vpre ? cc, m i 1 mj

(6)

where vpost and vpre are post- and precollisional velocities of particle i. Some precautions have been taken to prevent colliding particles from artificially sticking together. Otherwise significant amounts of CPU time would be lost in meaningless collisions actually corresponding to sliding of particles with respect to each other. First, if the perpendicular component of the relative impact velocity is less than Vs /100, where V is particle’s angular velocity around the central mass, the coefficient of restitution is set to « 5 1. Second, there is a minimum time limit for accepting subsequent collisions for a recently collided pair, typically 0.001/2f times the orbital period. In the numerical simulations extra care has to be taken to ensure that numerical integration errors are so small that they do not affect the results. All calculations are carried out in double precision. In a test particle run of 1.355 3 106 orbital periods of Galatea (about 1590 years), the observed relative error in angular momentum is DIz / Iz Q 21.2 3 1024 and, in kinetic energy, is DE/E Q 2.4 3 1024. Because the numerical integration errors grow linearly with respect to time, they are considerably less in the following simulations than those quoted above due to the very much shorter time spans investigated here, the longest of which is 100 years. On the other hand, systematic uncertainties in the relative positions of simulated particles from analytically calculated resonance locations can be as great as several meters, owing to the uncertainties in the latter. 3. ARC–SATELLITE GEOMETRY

To be able to simulate the complicated dynamics of Neptunian ring arcs, the physical parameters, size and mass of the ring particles, as well as their initial kinematic parameters, position and velocity, have to be calculated carefully. We examine only the strongest gravitational resonances due to Galatea falling in the vicinity of the arcs: a 42:43 outer Lindblad resonance (OLR), a 42:43 corotation– inclination resonance (CIR), and an 84:86 outer vertical resonance (OVR). We also briefly dwell on another resonance falling within this system—a form of parametric resonance—which was referred to as a Lindblad inclined

resonance (LIR) by Foryta and Sicardy (1996) in an analytical and numerical study of Neptune arc dynamics. The condition for the first-order OLR can be written as mLns 5 (mL 1 1)n 2 Ã ˙,

(7)

where mL is the wavenumber associated with the OLR (in this case mL 5 42), and ns and n are the mean motions of satellite and particle, respectively. Ã ˙ is apsidal precession rate of the particle. Mean motion of the particle (and also that of the satellite) n is related to geometric semimajor axis a by n2 5

S

SD

rN GMN 3 1 1 J2 a3 2 a

2

2

S DD

rN 15 J4 8 a

4

(8)

in the case of oblate planet. G is gravitational constant, MN is Neptune’s mass, and rN is equatorial radius of Neptune. Only the coefficients of the second ( J2 ) and fourth ( J4) gravitational harmonics have been taken into account. We have used Ms /MN 5 2.1 3 1028 and rN 5 25225 km for the mass and equatorial radius of Galatea and J2 5 341.05 3 1025 and J4 5 23.47 3 1025 for the coefficients of gravitational harmonics (Owen et al. 1991). The condition for the CIR is ˙ s, mcn 5 (mc 2 1)ns 1 V

(9)

˙ s is nodal regression rate of the satellite, and mc where V is the wavenumber of CIR (mc 5 43 for the Neptunian arcs). There are 2mc 5 86 corotation sites distributed around the ring. The condition for the OVR that falls p20 m inward of the exact CIR radius may be written as ˙ 5 (mv 2 1)ns 1 V ˙ s, (mv 1 1)n 2 V

(10)

˙ is nodal precession rate of the particle, and where V mv 5 85 is the wavenumber of the OVR. All the critical arguments for the resonances above are given in Table I, in which l and ls refer to the longitudes of a ring particle and satellite, respectively. Equations (7), (9), and (10) apply only at the exact resonance radius. In practice, however, resonant forcing can occur at some finite distance from the resonance radius if satellite perturbations act to keep the orbits periodic. For example, in the case of the Lindblad resonance, the condition (mL 1 1)

dn dà ˙ Da 5 Da da da

(11)

¨ NNINEN AND PORCO HA

4

TABLE I Critical Arguments and Locations for the Resonances in the Neptune Ring Arc–Galatea System Resonance

Critical Argument C

m

Corotation–inclination (CIR) Outer vertical (OVR) Outer Lindblad (OLR) Parametric resonance

Cc 5 2[mcl 2 (mc 2 1)ls 2 Vs] Cv 5 (mv 1 1)l 2 (mv 2 1)ls 2 V 2 Vs CL 5 (mL 1 1)l 2 mLls 2 Ã Cp 5 2[mpl 2 (mp 2 1)ls 2 V]

mc 5 43 mv 5 85 mL 5 42 mp 5 43

must be satisfied at a distance Da from the exact resonance radius. Now the radial gradients include the contributions from other significant gravitational perturbers; in the Neptune arc case, these are the satellite and the planet. The 2mc 5 86 stable corotation sites around the ring are potential maxima which are located where the CIR critical argument Cc 5 f. Small-amplitude librations of particles about the potential maximum can be written (Goldreich et al. 1986) as Cc 5 f 1 Dc cos u,

acaDc sin u, 3mcnc

Ac . a 2c

Ms b 2mc11 (ns as )2 b 3/2 (b), m . 0, Mp 8

0.13 0.13 0.13 0.13

At the beginning of each simulation u and Dc are selected randomly from the uniformly distributed values of 6[08, 1808] and [0, 1.4], respectively. The forced eccentricity due to the 42:43 OLR has to be carefully set for the particle orbit to be periodic in the reference frame of the satellite. Following Horanyi and Porco (1993), we chose the forced eccentricity at a given distance from the resonance to yield periodic orbits, and then calculated the eccentricities at other distances so that they obeyed the relation

eY

1 , a 2 aL

(16)

where aL is the radius of the OLR. The initial pericenter of the particle is set so that f Dc cos u. Ã 5 Vs 1 1 2 2

(14)

In all the simulations, we have set the satellite inclination is 5 0.038. This is in fact the value of the mutual inclination between the arcs and Galatea; however, there are indications that the pole of Neptune as determined from Voyager tracking data may be in need of a revision (Porco et al. 1995) that would put the arcs in their Laplace plane and would make the mutual inclination between the arcs and the satellite virtually identical to the inclination of the satellite relative to its Laplace plane. In this case, the configuration selected for these simulations is the correct one. Finally, the coefficient Ac can be written with the help of Laplace coefficients b si (b) Ac 5

6 6 6 6

(13)

and the libration frequency a is a 2 5 12m 2c i 2s

62,932.68 62,932.66 62,931.03 62,932.64

(12)

where Dc is the libration amplitude and u is the phase of the libration. The semimajor axis of a librating particle can be written a 5 ac 6 Da, where ac is the semimajor axis of the CIR site, and where Da 5

Location (km)

(15)

where Ms /Mp is the satellite–planet mass ratio, and b 5 as /a , 1 (Goldreich et al. 1986).

(17)

Initially, the satellite is set at its vertical maximum (ls 5 Vs 1 f/2) and the particles are distributed around a single CIR site, the center of which is separated from the satellite by f in longitude: l 5 ls 1 f. As a result the particle is at its apocenter at superior and inferior conjunction with the satellite and is thus in a closed orbit in the satellite’s frame. It is not possible, however, to set a particle on a closed orbit in the vertical direction with respect to the satellite. The libration of the particle’s guiding center around the CIR causes the vertical forcing to change during the libration period. Thus, the equilibrium value of inclination is expected to oscillate during a time equal to the CIR libration period of approximately 10 years. Because of the coupling between different resonances, it was empirically found that the relation (in the case of isolated vertical resonance) from the linearized theory

iY

1 , a 2 av

(18)

SIMULATIONS OF NEPTUNE’S RING ARCS

where av is the radius of vertical resonance, does not give sufficiently stable orbits with constant inclination. The best results for the arc particles gave an empirical relation

i 5 0.028 2

Dc 0.0058. 1.5

(19)

Thus the inclination will decrease linearly when moved away from the CIR point. The ascending nodes of the particles were initially set so that the nodes were exactly the same for every particle and they are aligned with the descending node of the satellite, V 5 Vs 1 f.

(20)

Comparison with Eq. (17) shows that for the particles exactly centered in the CIR zone (i.e., Dc 5 0), the nodal line is 908 away from the pericenter. We call attention to the presence of an additional resonance that falls among those mentioned above. This fourth resonance satisfies the condition ˙ 5 (mp 2 1)ns , mpn 2 V

(21)

where mp 5 43 is the resonance wavenumber. Because this particular resonance occurs when the forcing frequency is equal to twice the natural vertical oscillation frequency of the test particle, it is a form of parametric resonance. Such a disturbance is characterized by a slowly increasing amplitude of oscillation (in this case, in the forced inclination) very close to the natural vertical frequency of the particle; the amplitude growth time is p104 years (Goldreich, personal communication). Examination of the critical argument given in Table I shows that the nodal line of the satellite (and, therefore, the magnitude of its inclination) is irrelevant. Physically, this resonance arises once the particle has acquired an inclination through some other mechanism; then, the natural vertical movement of the particle in orbit around the planet brings about a resonant forcing by the satellite, whether it is on an inclined orbit or not. We show that among the strongest resonances located in the vicinity of the arcs which involve the inclination of either the satellite or the particle, this resonance is the least dynamically important.

5

and setting the gravitational harmonic coefficients to zero, i.e., J2 5 J4 5 0. In this case the 42:43 outer Lindblad resonance falls at RL 5 62932.1839 km. A test particle was given a semimajor axis a 5 RL 1 1.6 km, an eccentricity ae 5 28.481 km, and an orbital phase that keeps the orbit of the particle closed and periodic in the reference frame of Galatea. That is to say, we have set up the orbital shape and orientation so that the forcing of Galatea on the apse precession rate, which depends on eccentricity, keeps the orbit closed and nonlibrating in the satellite’s frame. (In selecting the initial boundary conditions for the test particle, we deliberately ignored the presence of the vertical resonance; i.e., the particles were given zero inclinations. These inclinations grow as a result of vertical forcing by the satellite, an effect we ignore in the first part of this paper.) The results are illustrated in Fig. 1 in which the behavior of the three planar osculating orbital elements is shown: the semimajor axis a, eccentricity e, and longitude of pericenter à are plotted during the encounter of Galatea. Comparing Fig. 1 with Fig. 3 of Horanyi and Porco (1993), it is clear that both sets of calculations using different numerical techniques yield the same results. To test the arc shepherding process, we included the nonzero J2 and J4 for Neptune. The gravitational harmonics (due to the oblateness of Neptune) will modify the mean motion and apsidal precession of the ring particles and satellites, shifting the resonance locations. In this case the outer 42:43 Lindblad resonance will be located at RL 5 62931.03 km and the 42:43 (CIR) will be located at Rc 5 62932.68 km (Horanyi and Porco 1993). We have set three particle pairs initially so that the particles start at the opposite side of the CIR with a 5 Rc 6 0.10, 60.15, and 60.20 km, respectively. In this simulation particles have zero free eccentricities. Figure 2 shows the evolution of azimuthal distance between particles in the same pair. The distance grows infinitely when the particles are not confined. Only the particles in the pair a 5 Rc 6 0.10 km are confined. Thus, the semimajor axis dispersion in the arcs falls between 200 and 300 m, confirming the result by Horanyi and Porco (1993) (compare with their Fig. 6). It is evident that our simulation method is capable of reproducing the simple dynamical features of the arc/Galatea interaction observed with other methods. 5. CONFINEMENT AND COLLISIONS

4. VERIFICATION OF THE RESONANCE CONFIGURATION

We began our numerical studies by first verifying the accuracy of the simulation code against the results of previous numerical studies of the arc/Galatea dynamics. Following Horanyi and Porco (1993), we isolated the dynamical effects of Galatea by ignoring the oblateness of Neptune

There are three different relevant time scales related to the arcs/Galatea system: the libration period of the arc particles around the corotation site, p10 years, or 9000 orbital periods of Galatea; the synodic period, p43 orbital periods of Galatea, about 18 days; and the orbital period of the arc particles and of Galatea, roughly 10.5 hr. Figure 3 illustrates the kinematics of the arc particles over the

6

¨ NNINEN AND PORCO HA

FIG. 1. Osculating orbital elements of a test particle during the encounter of the satellite (Galatea). In this simulation J2 5 J4 5 0. The elements are (from top to bottom) semimajor axis a, eccentricity e, and longitude of pericenter Ã. A time unit is Galatea’s orbital period Torb , i.e., Q0.428745 day. The test particle is in a closed periodic orbit with respect to the satellite, and at t 5 0Torb , the satellite and the test particle were 1808 apart. The data points from the test particle simulation are plotted with diamonds only during the closest encounter.

shortest relevant time scale, one orbital period of Galatea. The positions of the particles are plotted in the coordinate system corotating with the CIR site, and the origin of the coordinate system is fixed to the longitude and radius of the center of the site, taken to be f 5 0 and r 5 0. At the beginning of this simulation, the arc and its constituent particles are deliberately located near their apocenters and close to 1808 from the satellite. (Those particles with f 5 0 are exactly 1808 from the satellite at the beginning of the simulation.) The inclinations are taken to be zero, and eventually grow with time due to the nearby outer vertical resonance. The effects of this resonance turn out to be important for arc stability. For the moment, however, we concentrate only on the planar motions of the ring particles and defer discussion of the vertical motions until the latter part of this paper.

Several effects can be seen by close examination of Fig. 3. First, because of Galatea’s Lindblad perturbation, the arc particles, which are moving along an elliptic orbit around the planet, move radially in and out and we see the phased movement of the particles with respect to the particle exactly at the CIR location, taken to be r 5 0 km in Fig. 3. Thus, the origin of coordinate system is also moving radially, because it is following an elliptic orbit. Second, because particles are confined by the 42:43 corotation–inclination resonance, they do not longitudinally stray from the corotation zone which has a maximum angular width of p48. Third, because the Lindblad resonance falls at p 21.5 km from the semimajor axis of the CIR, particles with smaller semimajor axes (i.e., those closer to the Lindblad resonance) have larger forced eccentricity. This negative eccentricity gradient forces streamlines to cross at quadrature, which is seen after

SIMULATIONS OF NEPTUNE’S RING ARCS

7

FIG. 2. Azimuthal separation of test particles in the vicinity of the CIR location (Rc 5 62932.68 km). The particles in a pair were located initially at the opposite side of the CIR, a 5 Rc 6 0.10, 0.15, and 0.20 km.

one-fourth, and again after three-fourths, of an orbital period. The distortion in the shape of the arc, from tadpoleshaped to narrow ring and back again, is a result of the m 5 42 Lindblad pattern moving through the arcs with the (faster) angular speed of the satellite and the streamline crossing brought about by the negative eccentricity gradient. Notice the different radial scale in the plots. After onehalf orbital period particles are naturally close to their pericenters, and one-fourth orbital period after that the streamlines cross again at quadrature. This streamline crossing is another critical factor in arc stability and we will return to it later. To study collisional effects on the confinement of the arcs, we need to examine the behavior of particles over the longest relevant time scale: a CIR libration period. This period will vary with distance from the CIR: for the particles examined here, it varies from p7 to p10 years. The evolution of 800 massless colliding particles with radius s 5 15.5 m is followed in Figs. 4–7 for one-half libration period. The initial optical depth of the arc was t0 Q 2 3 1024. (The actual arc t p 0.1, but this is largely due to dust. The t in meter-sized particles is unknown.) In Figs. 4a and 5a, the particle radial positions and semimajor axes,

respectively, with respect to the particle exactly at the CIR location, are plotted against their longitudinal locations. (Figures 4b and 5b show the collisionless simulations as a reference.) The first four plots in each of Figs. 4a and 5a demonstrate the collisional effects during the first few synodic periods of Galatea (Tsyn 5 43Torb ). The plot t 5 3400Torb shows the arc after about one-fourth of libration period and the plot t 5 6795 Torb is after about one-half libration period. It is clear, especially comparing the collisional and collisionless simulations, that collisions deplete the arc particle population. Only a small change in semimajor axis is sufficient for removal, and in a few collisions the observed displacement in a has been several kilometers, i.e., of the order of 10 times the width of the CIR zone. Looking more closely at the last two plots in each of Figs. 4 and 5, however, shows that a few particles do survive the collisions and remain confined by the satellite, even though most of the particles have escaped the CIR zone. Because the collision frequency decreases with decreasing particle number in the arcs, the time scale for the removal of the remaining confined particles can be expected to be longer than a half of the libration period. Furthermore, not every collision leads to a removal of a colliding particle; we have observed

8

¨ NNINEN AND PORCO HA

FIG. 3. Kinematics of arc particles, with a short, one-orbital-period simulation of 800 particles. The particle positions are plotted in the reference frame that corotates with the exact CIR location. The origin of the coordinate system follows elliptic orbit around Neptune; i.e., it also moves radially. The units of time refer to the orbital period of Galatea, not to the orbital period of the arc.

particles to have several collisions in which displacements in semimajor axis are not large enough to remove a particle from the CIR zone. The 42:43 OLR is much stronger, and wider, than the 42:43 CIR: the width of the CIR is only about 2Da Q 300 m; that of the OLR is p 9 km. Consequently, even if a particle is removed from the CIR, it can remain trapped by the Lindblad resonance if the eccentricity after the collision remains close to the forced eccentricity (Foryta and Sicardy 1996). In this circumstance, it continues to follow the Lindblad resonance distortion, as can be verified from the last two plots: all particles are observed to have the Lindblad signature, even though they have acquired a bigger or smaller amount of free eccentrity in the collisions. The last two plots also clearly demonstrate why it is necessary to examine the dynamics over a period at least as long as half of the libration period: it is over this time scale that the dispersive effects of collisions are observed.

Figure 6 demonstrates the evolution of the semimajor axis distribution for only those particles whose azimuth remains within the resonance zone. The three different snapshots have been taken from the previous simulation. Initially, the semimajor axes are confined within the resonance zone (notice the difference in the scale of y axis in the first and following snapshots). After one-fourth of the libration period, collisions have caused the spread to be larger than the resonance width. However, during the next one-fourth libration period, all but a few unconfined particles leave the vicinity of the resonance site due to differential rotation; the vast majority of those that remain are clearly confined. This confinement is better illustrated in Fig. 7, in which the azimuthal profile of the arc has been plotted. Initially the particles are located within p28, about half of the length of the resonance zone. After one-fourth libration period, we see mainly an increased spread in azimuth, but after

SIMULATIONS OF NEPTUNE’S RING ARCS

9

FIG. 4. Collisional effects on the arc dynamics: (a) collisional simulation, (b) collisionless simulation. The particle positions are plotted with respect to the reference frame corotating with CIR location. The first four plots are after consequent synodic periods (roughly, the synodic period is 43 orbital periods of Galatea), and the last plot shows the arc after half of the libration period. The 800 simulation particles are massless with a radius s 5 15.5 m.

twice that time, almost all of the remaining particles near the arc are confined. Simulations with different particle sizes yield roughly similar results. Figure 8 compares dispersion in semimajor axes after one-half libration period in three different simulations with particle radii s 5 62, 31, and 15.5 m. The initial optical depth has been taken to be the same value (t0 Q 2 3 1024) in all these simulations by adjusting the number of particles in each simulation: 50, 200, and 800 particles, respectively. Figure 9 shows the azimuthal density profile for the same simulations. To study more carefully the evolution of the particle density in the arcs, we have plotted the number of particles in the arc as a function of time in Fig. 10. We have used simulations with four different particle sizes (s 5 12.4, 15.5, 31, and 62 m), but with the same initial optical depth. As collision frequency is proportional only to optical depth, or density, for identically sized particles, we should be able to see particle size effects.

The evolution of particle density can be divided into two separate groups on the basis of Fig. 10. In the simulations with s 5 31 m or less the density of the arc is decreasing slower than in the simulation with s 5 62 m; since the optical depth is initially similar, the difference between these evolutions must be due to a difference in relative collision velocities. The particle size s 5 62 m is apparently big enough to add a significant extra component to the velocity dispersion. The basic velocity dispersion is due to streamline crossing at quadrature. Finite size effects seem to be insignificant for particles s p 30 m and smaller; the critical particle size at which finite size effects become important is thus somewhere between s p 30 m and s p 60 m. We also tested for possible dependency of the density evolution on the initial optical depth. We did a simulation with s 5 62 m but with half the initial optical depth (longdashed line in Fig. 10). The impact frequency should be half of that in the other simulations, resulting in relatively

10

¨ NNINEN AND PORCO HA

FIG. 5. Similar to Fig. 4, except that the particle semimajor axes are plotted instead of the radial positions. (a) Collisional simulation, (b) collisionless simulation.

twice as many remaining particles at the end of simulation. This is exactly what we observe in Fig. 10. 6. COLLISIONS AND PARTICLE ORBITS

As the CIR site is very narrow, with a width of only 2Da Q 6300 m, even very small relative collisional velocities may be able to remove a particle from the resonance site. A small tangential impulse as low as 0.5 m/sec at quadrature could result in a displacement of 6 km in the particle’s semimajor axis (Porco 1991). In Fig. 11 the changes occurring in particles’ semimajor axes in every collision during the simulation run are plotted as a function of time of impact. The mirror symmetry in the plot is due to the lack of particle size distribution: both the colliding particles are equally massive and therefore experience a displacement equal in magnitude but opposite in direction. Most of the collisions occur between t 5 1Torb and t 5 100Torb . The dramatic drop in the frequency of impacts is due to the strong decrease in the particle density of the arc (see the previous figure). Most of the collisional changes in

a are smaller than the width of the resonance site, implying that these particles do not necessarily leave the confinement zone; however, a considerable number of particles have uDau p 100 m or larger. Furthermore, even though a particle may have uDau p 10 m in all impacts in which it participates, a collisionally induced random walk in a could result in a particle getting outside the boundaries of the CIR zone. Arc particles can be expected to collide preferentially at the quadratures of their orbits where streamline crossing occurs (Porco 1991). This is demonstrated in Fig. 12, in which the phase distribution of colliding particles is plotted: most of the collisions happen at 908 before or after the pericenter, even though these maxima are not symmetric. Surprisingly, there are still other, smaller maxima in the distribution. These maxima are due to the fact that particles initially orbiting in the equatorial plane will feel small kicks in the z direction at every satellite encounter, acquiring small inclinations. After a while the growth in the dispersion in the ascending nodes around their initial alignment with quadrature will effectively reduce the number of colli-

SIMULATIONS OF NEPTUNE’S RING ARCS

11

FIG. 6. Evolution of the distribution of semimajor axes for the arc particles. The three snapshots correspond to the initial state and roughly to one-fourth and one-half of the libration period. The snapshots are from the simulation illustrated in the previous figures.

FIG. 8. Comparison of semimajor axis distribution at the end (t 5 6795Torb Q one-half of the libration period) of three collisional simulations. The simulations were similar, except for the particle size: s 5 62, 31, and 15.5 m. Arc optical depth was kept constant by adjusting the number of particles: N 5 50, 200, and 800, respectively. Only azimuthally confined particles are used when calculating the distribution.

FIG. 7. Evolution of the azimuthal density profile for the ring par-

FIG. 9. Azimuthal density profile compared at the end of the simulations with different particle sizes. See the previous figure.

ticles.

12

¨ NNINEN AND PORCO HA

FIG. 10. Number of particles confined in the arc plotted as a function of time for several collisional simulations. With one exception, the simulations have different particle sizes, but the optical depth of the arc was initially set to the same value by adjusting the number of particles.

sions occurring at quadrature. This same dispersion causes particle orbits to cross closer to pericenter and apocenter. Finite particle size effects then lead to collisions at these locations also. And this is what is observed. The effect of an impact depends on the magnitude and direction of the relative collisional velocity and these in turn may depend on the orbital phase and size of the colliding particles. Particle size effects aside for the moment, the radial component of the relative impact velocity at quadrature will be Dvr Q 0.6 m/sec, if the colliding particles are on orbits near the edges of the corotation zone, i.e., one orbit 100 m outward from the exact CIR location and the other 100 m inward from the CIR. The azimuthal component of the relative impact velocity is only Dvw Q 2 cm/sec in this case. (For collisions occurring at pericenter, the azimuthal component of relative velocity would be Dvw Q 2 cm/sec and the radial component would be zero.)

Histograms of both the observed radial and azimuthal components of collisional relative velocity are shown in Fig. 13. Even though both distributions have a high end tail with values up to more than 2 m/sec, the peak is concentrated on the small values (Dvr,w , 10 cm/sec). There are also clear differences: the radial distribution is wider than the azimuthal distribution. There are several reasons why the observed collisional velocities differ from our crude estimate above. First of all, streamlines with differing semimajor axes do not cross exactly at the quadrature. For example, two streamlines with semimajor axis 100 m outward and inward from the exact CIR radius and with proper forced eccentricity cross at w Q 938. Furthermore, finite particle size effects will move the collisional azimuth. For streamlines 200 m apart in semimajor axis, collisions between particles of size s 5 30 m will shift collisions about Dw Q 0.98 farther from quadrature; collisions between particles with s 5 60 m will

SIMULATIONS OF NEPTUNE’S RING ARCS

13

FIG. 11. Collisional change in the semimajor axis of colliding particles plotted as a function of time of collision. The simulation with particle size s 5 15.5 m is used.

FIG. 12. Phase distribution of colliding particles. The simulation with particle size s 5 15.5 m is used.

14

¨ NNINEN AND PORCO HA

FIG. 13. Velocity distribution of azimuthal and radial components of relative impact velocity. The collisional data are collected from the simulation of 800 particles with particle size s 5 15.5 m.

shift it twice that amount. If streamlines are only 100 m apart in semimajor axis (50 m inward and outward from the CIR radius), the shifts will be doubled. As a result the azimuthal component of the relative impact velocity will increase. The radial component will decrease because the collisions no longer occur exactly at quadrature where the Dvr is at maximum. Furthermore, as there is a change in a phase of colliding particles, there will be a commensurate effect on the other orbital elements as well. Figure 14 illustrates the effects of the azimuthal and the radial components of the relative impact velocity on the semimajor axes of the colliding particles. Only the azimuthal component has a correlation with a change in semimajor axis, as is expected from Gauss’ equations. Using Gauss’ equations we can estimate that Dvw 5 0.5 m/sec should have an effect of Da Q 6 km, which is actually very close to what we observe in Fig. 14. Collisions also transform energy from orbital motions into random motions; i.e., particles will acquire free eccentricities. This may be caused by a change in particle’s eccentricity and/ or longitude of pericenter. Figure 15 shows these changes as a function of the azimuthal and radial components of the relative impact velocity. The correlation between a change in either orbital element and either of the compo-

nents of the impact velocity is not as clear as in the case of the changes in semimajor axis. Gauss’ equations also confirm this: both the radial and azimuthal components will have an effect of the same order of magnitude, the relative importance of the effect depending on the actual phase of colliding particles in their orbits. There is still one source of free eccentricities we have not discussed. Even without any change in eccentricity or in longitude of pericenter, the change in a particle’s semimajor axis means that a particle does not have a proper forced eccentricity with respect to the forcing of 42:43 OLR. A 100-m shift inward in a would leave a particle with an eccentricity p7% too small for its radial location, resulting in a free eccentricity if the particle is still confined in the arc by CIR. 7. EFFECTS OF THE 84:86 VERTICAL RESONANCE

Until now we have ignored the nearby 84:86 OVR when calculating the initial conditions; i.e., all the particles begin with zero, instead of the proper value of forced, inclination. In this circumstance, the particles are not in periodic orbits with respect to the vertical resonance, their inclinations will grow from zero, and both their inclinations and their

15

SIMULATIONS OF NEPTUNE’S RING ARCS

FIG. 14. Collisional change in the (geometric) semimajor axis plotted as a function of azimuthal and radial components of the relative impact velocity.

nodal lines will oscillate around their equilibrium values. The period of oscillation will depend on the distance from the OVR. (We refer to this as the OVR libration period, to distinguish it from the CIR libration period.) To study the effects of the vertical resonance on the behavior of arc particles, it is necessary to run simulations that extend over time scales equal to the OVR libration period, since it is over this time scale that the vertical forcing will change substantially. The 84:86 OVR is located p20 m inward from the 42:43 CIR, and because the width of the OVR, WOVR Q 400 m (Franklin et al. 1984), is larger than the width of the CIR, WCIR Q 300 m, all the arc particles are affected by the OVR. Since only a 0.0018 difference in inclination corresponds to p1-km difference in the vertical direction, the OVR has important effects on the arc stability as we show below. In the vicinity of an isolated OVR, a particle that is locked ‘‘in resonance’’ will have an entirely forced inclination (i.e., it will have no free inclination) whose magnitude depends on the distance from resonance. The particle’s critical argument will be constant in time and equal to its equilibrium value: either 08 for a , a0 or 1808 for a . a0 , where a0 is the location of the resonance. For a particle not locked into resonance but undergoing resonant pertur-

bations, the critical argument will librate around its equilibrium value and the libration frequency, «, will change with distance from the resonance, «5

˙ dC Da, da

(22)

where C is the critical argument and Da 5 a 2 a0 , the distance between the particle semimajor axis and the resonance. Nonperiodic Orbits: The 2D Case Figure 16a illustrates the inclination behavior of three particles—one placed at the center of the CIR zone and two placed 100 m on either side of center—in a 100-year simulation in which the initial inclinations were set to zero. (As we show later, the ‘‘centered’’ particle is actually slightly displaced from, and undergoes a small-amplitude libration around, the CIR site.) All test particles exhibit a periodic variation in inclination with periods ranging from p38 to p45 years. (Though we have plotted the inclinations of the particles, the periodic behavior applies to the variation of mutual inclination of particles and satellite as well,

¨ NNINEN AND PORCO HA

16

FIG. 15. Collisional change in the (geometric) eccentricity of colliding particles plotted as a function of azimuthal and radial components of relative impact velocity. Also, the collisional change in the (geometric) pericenter is shown as a function of both the impact velocity components.

since the inclination of the satellite remains constant.) The center particle (solid line) exhibits the largest and most regular oscillation, with a period of t Q 43.5 years. Using the form of the critical argument for this resonance given in Table I, we find ˙v ˙ dC dn dV 5 (mv 1 1) 2 da da da

(23)

and

S

S DD

˙v n0 RN dC 7 3 1 1 J2 5 2 (mv 1 1) da 2 a0 4 a0

2

1

˙0 7V . 2 a0

(24)

With mv 5 85, J2 5 0.0034, RN 5 25225 km, and a0 5 62932.66 km, we find a libration period of 43.5 years for a Da 5 13.5 m. To within the uncertainties, this Da is

consistent with the analytically derived difference in the locations of the OVR and the CIR, the latter being the presumed semimajor axis of the arcs. The variation in inclination of arc particles was also observed in the Neptune arc investigation of Foryta and Sicardy (1996) but was given an entirely different interpretation. These authors claimed that the period of variation in the mutual inclination is p27 years (which in fact is not consistent with their Fig. 9b), that it is due to the natural differential nodal regression between the arcs and the satellite, and therefore that the vertical forcing by Galatea is too weak to have any effect on the vertical motion of the particles. However, these authors do state that the parametric resonance, which they refer to as a Lindblad inclined resonance, is responsible, along with the CIR and the Lindblad resonance, for the motion of the arcs on a long (500-year) time scale. One has only to examine the behavior of the nodes of

SIMULATIONS OF NEPTUNE’S RING ARCS

17

FIG. 16. Evolution of inclination (a) and azimuthal distance from the CIR (b) plotted as a function of time. This test particle simulation consists of three particles that are initially at the same longitude (1808 from the satellite) but have different semimajor axes. One particle is placed at the CIR point (a 5 62932.68 km, solid line), and the two others are 100 m inward (dotted line) and 100 m outward (dashed line). The particles have zero inclinations at the beginning of the simulation.

the satellite and particle to see that this is not the correct interpretation. In Fig. 17a, we plot V 2 Vs (i.e., the longitude of the particle ascending node with respect to the longitude of ascending node of the satellite) versus time for the same three particles displayed in Fig. 16a. It is clear that the phasing of the vertical motion of the center particle (solid line) is strongly influenced by Galatea and, in fact, that instead of circulating around 3608, the particle’s ascending node librates around the descending node of the satellite with the OVR libration period of p43.5 years. One can examine the behavior of all the critical arguments involving inclination (given in Table I) to understand the coupling among all the resonances and the effect of the OVR. In Fig. 18 we plot the critical arguments for the CIR (Cc), the OVR (Cv), and the parametric resonance (Cp) versus time for the particle centered on the CIR site. Several noteworthy features are immediately obvious.

First, there are three different time scales in all the plots. The shortest period, p7 months, is caused by a forced variation in the particle eccentricity vector due to the OLR and the fact that the particle is initially not strictly on a periodic orbit with respect to the OLR. The appearance of this period in the critical arguments involving inclination is probably due to coupling among all the resonances (Foryta and Sicardy 1996). The p7-year period is due to the libration of particle’s guiding center around the CIR. (Recall that this particle is not exactly situated on the CIR site, but is displaced probably a few meters from it: hence, the low-amplitude libration of its critical argument.) This motion also causes the distance from the vertical resonance to change periodically. The longest p43-year period is apparently the effect of the OVR, as graphically demonstrated in Fig. 17a. Second, among all the critical arguments, Cv is the most

¨ NNINEN AND PORCO HA

18

FIG. 17. Evolution of the longitudes of ascending node with respect to the longitude of ascending node of the satellite is plotted for the particles in Fig. 16 (a) and in Fig. 19 (b). Note the different scale in the plots.

constant in time; i.e., the particle is apparently most strongly influenced by the OVR. (This remains true until the inclination of the particle nears zero at which point the nodal line difference between particle and satellite quickly resets from 2708 back to 908 after the OVR libration period of 43.5 years.) To the contrary, the parametric resonance critical argument circulates through 3608, indicating that it plays no role at all; moreoever, it does so with OVR libration period, as well. The CIR critical argument exhibits an intermediate behavior. Though obviously perturbed by the OVR into a libration period of 43.5 years, Cc librates around a value of 1808, as it should, for arc particles locked into the corotation–inclination resonance. Figure 18 can be easily mathematically verified by noting the relationships among the critical arguments involving inclination: Cc 5 Cv 1 (V 2 Vs),

(25)

Cp 5 Cv 2 (V 2 Vs).

(26)

Here we have written the critical arguments of the CIR and the parametric resonance as a function of the critical argument of the OVR; i.e., the addition/subtraction of Fig. 17a to/from the middle panel of Fig. 18 explains the behavior of Cc and Cp . We now focus on the two particles initially 100 m inward and outward from the CIR site, which are respectively 80 m inward and 120 m outward from the OVR. Figure 16a shows their inclination behavior. Averaged over a CIR libration period of roughly 10 years, these particles spend more time farther from the OVR than does the CIR-centered particle. Hence, their inclination variations have a smaller amplitude, though still nearing zero approximately every 40 years. Their nodal lines (upper panel, Fig. 17a) are clearly influenced by Galatea’s perturbations, and show a tendency to librate around that of the satellite.

SIMULATIONS OF NEPTUNE’S RING ARCS

19

FIG. 18. Critical arguments of the CIR (Cc), the OVR (Cv), and the parametric resonance (Cp) have been plotted for the particle located at the CIR site in the two-dimensional case; i.e., the particle that is plotted with a solid line in Figs. 16 and 17a.

In Fig. 16b is plotted the azimuthal distance between each of these particles and the particle located at the CIR site. While the inner particle (dotted line) is confined in the arc (librating irregularly around the CIR point), the outer particle (dashed line) will clearly leave the CIR site after p40 years; the azimuthal separation of CIR sites is about 48. However, this particle is later trapped into another CIR site (not the neighboring, but the next one) for p20 years. After this it will still experience yet another transition into another CIR site during the first 100 years. We have not observed this kind of migration into neighboring CIR zones as a result of collisional changes in the orbital elements, only in this long-term integration of test particles with initially zero inclinations. Periodic Orbits: The 3D Case Figure 19 shows the evolution of inclinations in a test particle simulation that is otherwise identical to that of Fig. 16, but the initial inclinations and nodal line orientations for the test particles are set as described in Section

3. The inclination stays nearly constant, as does DV, for the particle positioned at the center of the CIR zone because its distance from the OVR changes only slightly in time; it is close to being in a periodic orbit with respect to both the CIR and OVR. The equilibrium nodal line difference, DV p 1808, has the added benefit of maximizing the vertical distance between satellite and particle during conjunction, which is the most stable configuration for the CIR. This is not the case for particles located some distance from the CIR point. As these particles librate, their distance from the OVR varies, causing the strength of the vertical forcing to change periodically. Furthermore, their longitudes of ascending node exhibit a larger amplitude of libration around the equilibrium value (DV 5 1808), as demonstrated in Fig. 17b, than the centered particle. The lower plot in Fig. 19 shows the azimuthal distance of the 100-m-displaced test particles from the particle at the CIR point. The test particles are observed to librate regularly around the CIR point. In the behavior of the critical arguments in the threedimensional case (plotted in Fig. 20) compared with the

20

¨ NNINEN AND PORCO HA

FIG. 19. Similar to Fig. 16, except that now the particles have properly set inclinations and nodal lines for the 84:86 OVR.

earlier two-dimensional case, there are two clear differences. While the 7-month and 7-year periods are still present, the long-period (43 years) oscillation has disappeared. This is expected because the particle now has the stable equilibrium inclination, and consequently, the nodal line oscillation amplitude around the value DV 5 1808 has decreased to few degrees. The other difference is that now all the critical arguments are well confined, and they are behaving very similarly. In Fig. 21 the phase of the colliding particles is plotted in the simulation of 800 particles with s 5 15.5 m. The simulation is otherwise identical to that used in Section 6, except that the particle inclinations and nodes are initially set as described in Section 3. When comparing with Fig. 12, it is evident that the effect of forced inclinations is to concentrate the collisions more effectively on the quadratures of particle orbits. The previously observed third collisional frequency maximum near pericenter has now totally disappeared, and there is only a residual number of collisions occurring outside quadrature. Furthermore, the

peaks are more symmetrical than in the two-dimensional case. It should also be noticed that the total number of collisions has decreased considerably. In Fig. 22 the azimuthal and radial components of the relative impact velocity are plotted. When compared with the results in the two-dimensional case (Fig. 13), we see that the high values of collisional velocities have been removed. There is also an observable change in the relative shape of vw and vr histograms. Even though the azimuthal velocity component was more steeply concentrated on low values than was the radial one in the previous simulations, now this effect is more enhanced. In fact, now the impact velocities resemble more our earlier estimates of the magnitudes of the radial and azimuthal components at quadrature. The collisional change in semimajor axis is plotted in Fig. 23 as a function of the components of relative impact velocity. Now there is a big difference between the present and two-dimensional results, as should be expected because of the lower impact velocities. It is also

SIMULATIONS OF NEPTUNE’S RING ARCS

21

FIG. 20. Critical arguments (Cc , Cv , and Cp) plotted in the three-dimensional case for the particle centered at the CIR site.

remarkable that the collisional changes in eccentricity and longitude of pericenter of the particles have decreased (Fig. 24). The OVR apparently helps to make the arc more stable against interparticle collisions via several different physical mechanisms: reducing the number of collisions, removing the high-velocity collisions, and suppressing the azimuthal velocity component relatively more than the radial one by concentrating the collisions more effectively at quadrature. Naturally, part of the relative impact velocity is transferred to the vertical direction, but an impulse in the z direction does not affect the semimajor axis and eccentricity of a colliding particle. To what degree do all these mechanisms improve arc stability? In Fig. 25 the evolution of number of particles in the arc is compared between the present and two-dimensional simulations. The number of particles that will stay confined has doubled. In 8 years (6815Torb), 20–30% of the original particles remain in the arc. In the simulations with s 5 15.5 m and s 5 62 m the initial optical depths are the same. In the simulation with s 5 124 m the initial optical depth is four times higher (i.e., the number of parti-

cles is 25), which makes the comparison with the other simulations somewhat difficult. We have also done a threedimensional simulation with 25 particles of size s 5 248 m. In this simulation the collisions remove all but two particles within less than 2 years, after which these two particles never collide during the rest of the simulation (about 8 years) and thus stay confined. It would be interesting to know how the number of particles in the arc evolves during longer time scales in different simulations but the CPU time limitations do not allow us to do longer simulations. However, one simulation (the dotted line in Fig. 25) was extended up to 100 years, and is shown in Fig. 26. To save CPU time unconfined particles were removed from the simulation at two stages; by the end of the simulation we were simulating only 100 confined particles. If unconfined particles could collide with the arc particles, this procedure would be incorrect; however, because the unconfined particles are not resonantly forced by either the OLR or the OVR, these particles will have different eccentricities and inclinations and their periapses and nodes will precess freely. The probability of collision, or even the crossing of orbits within a

22

¨ NNINEN AND PORCO HA

FIG. 21. Phase distribution of colliding particles. The simulation with particle size s 5 15.5 m is used. In this simulation the inclinations are initially properly adjusted with respect to the OVR. Compare with Fig. 12.

distance equal to a particle radius, is probably totally insignificant. Figure 26 shows the evolution of the number of arc particles up to 100 years for this simulation. During the first 10 years the particle number has dropped to about 20%, but the strong decrease in the collision frequency slows down the evolution. During the next 90 years the number of particles drops down to roughly 10%. At the end of the simulation there are still 68 confined particles out of 800 particles initially. 8. EFFECT OF THE COEFFICIENT OF RESTITUTION

The elasticity of collisions, described by «, the coefficient of restitution, depends on vc . In general, the lower the coefficient of restitution, the higher the energy lost and the lower will be the equilibrium velocity dispersion in a system of colliding particles. In the simulations discussed so far we have used Eq. (5), which relates « to vc for ice, and we have chosen vc 5 0.063 cm/sec; however, we have also studied how arc stability depends on the elasticity of collisions by choosing vc 5 0.0016, 0.010, and 0.063 cm/ sec. The value vc Q 0.01 cm/sec is the value found for icy particles in laboratory measurements (Bridges et al. 1984).

In these simulations we have used 100 particles with size s 5 62 m. The particle inclinations are initially properly adjusted with respect to the vertical resonance as in the simulations of the previous section. In Fig. 27, we compare the evolution in the number of arc particles remaining within the corotation zone among simulations with different vc . First thing to notice is that the elasticity affects the efficiency of the collisional removal process, but not dramatically. The second thing is that the effectiveness of removal is neither a decreasing nor an increasing function of vc . Why? As discussed earlier, the fact that collisions occur at a preferential orbital phase (i.e., quadrature), where orbits cross, maintains the velocity dispersion at a minimum value so that the usual equilibrium between the velocity dispersion and the coefficient of restitution is never achieved. Though the amount of damping is still directly proportional to «, smaller values of « do not significantly alter the relative impact velocities as these depend on the kinematics near streamline crossing. Normally, softer impacts, i.e., impacts with smaller vc , should remove particles less effectively from the arcs. Our simulations show that when vc 5 0.0016 cm/sec this is indeed the case. However, it appears that the collision

SIMULATIONS OF NEPTUNE’S RING ARCS

23

FIG. 22. Velocity distribution of azimuthal and radial components of the relative impact velocity for the three dimensional case. The collisional data are collected from the simulation of 800 particles with particle size s 5 15.5 m.

frequency also depends on the elasticity of collisions: simulations with smaller vc now have higher impact frequency. For example, the number of collisions that occurred during the three simulations with vc 5 0.0016, 0.010, and 0.063 cm/sec, respectively, were 107, 91, and 79. Thus as collisions become more elastic (or particles become harder or vc gets higher), collisions more effectively remove particles from the arc, until the impact frequency drops to the point where this trend reverses. Furthermore, it appears that the value of vc determined experimentally in the laboratory, i.e., vc Q 0.01 cm/sec, leads to the most efficient removal. 9. DISCUSSION

The dynamics of ‘‘large’’ Neptune arc particles is governed by three different resonances due to the nearby satellite Galatea. The confinement mechanism itself requires only closely spaced corotation–inclination and outer Lindblad resonances. However, as we have demonstrated in this paper, a nearby outer vertical resonance, also due to Galatea, increases arc stability against collisional effects, so long as particle orbits have proper forced inclinations and orientations. A fourth resonance, mentioned by Foryta

and Sicardy (1996), is a form of parametric resonance, is considerably weaker than the OVR, and does not have any significant effect on the arc dynamics. The reasons for this are largely geometrical. First, the strong inclination gradient and particle nodal line locking generated by the OVR decrease the impact frequency. Also, the inclination gradient transfers some of the relative impact velocity into the z direction decreasing the collisional change in semimajor axis that is the most important reason for a removal of a particle from the arc. Furthermore, the probability of a collision occurring outside quadrature has been observed to be practically zero when the particles are on forced inclined orbits. As a result the number of surviving particles is significantly increased when particles are properly phased with respect to the OVR. There are still other significant effects due to the OVR that have not yet been discussed. After a particle is removed from the arc due to a sufficiently strong collision, the differential rotation will cause it to have an encounter with the arc after a few years. If all the particles were in one plane, the probability for a collision between the already removed particle and an arc particle would not be insignificant; however, the gradient in inclination and the fact

24

¨ NNINEN AND PORCO HA

FIG. 23. Collisional change in the (geometric) semimajor axis plotted as a function of azimuthal and radial components of the relative impact velocity for the three dimensional case.

that the nodes of the arc particles do not precess freely as do the nodes of the unconfined particles further reduces the probability for a high-impact-velocity collision. We have observed, as have Foryta and Sicardy (1996), that particles which are not in periodic orbits with respect to the satellite’s vertical perturbations may spontaneously, without experiencing any impacts, leave the CIR site. This occurs at the point in the libration cycle when the inclination becomes so small that the corotation–inclination resonance becomes weak and confinement becomes unlikely. In principle, collisions could bring about a series of inclination changes such that initially confined particles eventually leak from one corotation zone into another; however, we have not observed this in our simulations. Alternatively, the five Neptunian arcs may be the result of particle migration from one corotation zone to another during an early time when the arcs were tending toward stability. If the arcs are in fact the result of the disruption of a larger body, then many particles would certainly have been initially lost completely from the arc region, some would have been caught directly into corotation zones, and others may have been caught into a cycle of corotation zone migration until collisions either removed them also or set them on a proper orbit for stability.

In the two-dimensional simulations, we found that the evolution of particle number density within the arcs depends on the particle size: particle density decreases more rapidly when the colliding particles are bigger than the critical size, found to be between 30 and 60 m. No such size dependency was observed when the arc particles were smaller than the critical size. Apparently, when the particle size is bigger than the critical size, finite size effects will add a significant and extra component to the velocity dipersion. In the three-dimensional simulations in which the particles were placed on periodic orbits (or nearly so) with respect to all resonances, we did not observe similar finite size effects despite the fact that we used particles ranging from 15.5 to 248 m. Recall in the three-dimensional case that there is already a considerable z component to the relative impact velocity, and that a change in the z component of a particle’s velocity due to a collision will not affect the semimajor axis and therefore confinement. Consequently, we suspect that in the three-dimensional case, the critical particle size above which the velocity dispersion increases dramatically must be larger than even the largest particles we considered in our simulations. We found also that the process of particle removal from the arc depends mildly on the elasticity of interparticle

SIMULATIONS OF NEPTUNE’S RING ARCS

25

FIG. 24. Collisional change in the (geometric) eccentricity of colliding particles plotted as a function of azimuthal and radial components of the relative impact velocity for the three dimensional case. Also, the collisional change in the (geometric) pericenter is shown as a function of the azimuthal and radial impact velocity components.

collisions, though not in a monotonic way. The effect of the coefficient of restitution on the evolution of arc density was studied by using the Bridges et al. (1984) coefficient of restitution and varying the scaling factor vc. It was found that the more elastic the impact, the more effective is particle removal from the arc; however, it was also found that the collision frequency is higher when impacts are softer. The combined effect of these opposing mechanisms is that the impacts are most destructive of arc stability for ice particles having vc Q 0.01 cm/sec. Our numerical simulations in which the arc particles are placed on nearly closed orbits with respect to the OVR, CIR, and OLR demonstrate that a significant population of large particles (s . 10 m) can be confined at least for 100 years. Given the dramatic drop in collisional frequency at the very beginning of each of our simulations, and the fact that arc stability apparently relies on a proper geometric phasing of the orbits of particles remaining in the arcs,

it is plausible that large particle confinement may last for time scales much longer than 100 years. The problem of dust creation in the arcs is still one begging for solution. Presumably, particle impacts are the source of dust within the arcs, yet the collisional frequency among the large, equal-sized bodies in our simulations over long time scales is not great enough to replenish the dust reservoir and account for Voyager observations. We will investigate in the future the dynamics of arc dust particles, as well as the possibility that the dust-forming collisions are those between the large particles we have studied here and a putative population of smaller (i.e., millimeter- to centimeter-sized) ones. ACKNOWLEDGMENTS We thank Dr. P. Goldreich and Dr. H. Salo for discussions. This work was supported by NASA Grant NAGW-3411 and the Academy of Fin-

26

¨ NNINEN AND PORCO HA

FIG. 25. Number of particles confined in the arc plotted as a function of time for several collisional simulations to compare the effects of the outer vertical resonance on the stability of the arcs. The initial values for the simulations with same particle size are similar, except the fact that in one simulation the particles are in equatorial plane of the planet, having thus zero inclinations, and in the other simulation the forcing of the OVR has been taken into account when calculating the initial values.

FIG. 26. Evolution of number of particles confined in the arc followed 100 years in a simulation of initially 800 particles with s 5 15.5 m. This is three-dimensional simulation; i.e., the particle inclinations are initially set to the value forced by OVR using our empirical equation.

SIMULATIONS OF NEPTUNE’S RING ARCS

27

FIG. 27. Evolution of the number of arc particles compared between three simulations with particles of different elasticity. The Bridges et al. coefficient of restitution is used with different vc scaling factors. The scaling factor vc 5 0.01 cm/s refers to icy particles, and the others are used as comparison. The scaling factor vc 5 0.063 cm/s is used otherwise throughout this paper. In these comparison simulations we have used 100 particles with s 5 62 m. land. We also thank referees Dr. B. Sicardy and Dr. I. Mosqueira for valuable comments.

GOLDREICH, P., S. TREMAINE, AND N. BORDERIES 1986. Towards a theory for Neptune’s arc rings. Astron. J. 92, 490–494.

REFERENCES

HA¨ NNINEN, J., AND H. SALO 1992. Collisional simulations of satellite Lindblad resonances. Icarus 97, 228–247.

AARSETH, S. J. 1972. Direct integration methods of the N-body problem. In Gravitational N-Body Problem (M. Lecar, Ed.), pp. 373–387. Reidel, Dordrecht. AARSETH, S. J. 1985. Direct methods for N-body simulations. In Multiple Time Scales (J. E. Brackbill and B. I. Cohen, Eds.), pp. 377–418. Academic Press, Orlando, FL. BRIDGES, F. G., A. HATZES, AND D. N. C. LIN 1984. Structure, stability and evolution of Saturn’s rings. Nature 97, 333–335. COLWELL, J. E., AND L. W. ESPOSITO 1990. A model of dust production in the Neptune ring system. Geophys. Res. Lett. 17, 1741–1744. FORYTA, D. W., AND B. SICARDY 1996. The dynamics of the neptunian Adams ring’s arcs. Icarus 123, 129–167. FRANKLIN, F., M. LECAR, AND W. WIESEL 1984. Ring particle dynamics in resonances. In Planetary Rings (R. Greenberg and A. Brahic, Eds.), pp. 562–588. Univ. of Arizona Press, Tucson.

HA¨ NNINEN, J, AND H. SALO 1994. Collisional simulations of satellite Lindblad resonances: II. Formation of narrow ringlets. Icarus 108, 325–346. HORANYI, M., AND C. PORCO 1993. Where exactly are the arcs of Neptune. Icarus 106, 525–535. LISSAUER, J. 1985. Shepherding model for Neptune’s arc ring. Nature 318, 544–545. OWEN, W. M., R. M. VAUGHAN, AND S. P. SYNNOTT 1991. Orbits of the six new satellites of Neptune. Astron. J. 101, 1511–1515. PORCO, C. 1991. An explanation for Neptune’s ring arcs. Science 253, 995– 1001. PORCO, C. C., P. D. NICHOLSON, J. N. CUZZI, J. J. LISSAUER, AND L. W. ESPOSITO 1995. Neptune’s ring system. In Neptune and Triton (D. P. Cruikshank, Ed.), pp. 703–804. Univ. of Arizona Press, Tucson. SMITH, B. A., AND 64 COLLEAGUES 1989. Voyager 2 at Neptune: Imaging science results. Science 246, 1422.