Modeling Saturn's D68 clumps as a co-orbital satellite system

The D68 ringlet is the innermost feature in Saturn's rings. Four clumps that appeared in D68 around 2014 remained evenly spaced about 30 degrees apart and moved very slowly relative to each other from 2014 up until the last measurements were taken in 2017. D68's narrowness and the distribution of clumps could either indicate that we have a collection of source bodies in a co-orbital configuration or imply that an outside force confines the observed dust and any source bodies. In this paper we explore the possibility that these four clumps arose from four source bodies in a co-orbital configuration. We find that there are no solutions with four masses that produce the observed spacings. We therefore consider whether an unseen fifth co-orbital object could account for the discrepancies in the angular separations and approach a stable stationary configuration. We find a range of solutions for five co-orbital objects where their mass ratios depend on the assumed location of the fifth mass. Numerical simulations of five co-orbitals are highly sensitive to initial conditions, especially for the range of masses we would expect the D68 clumps to have. The fragility of our D68 co-orbital system model implies that there is probably some outside force confining the material in this ringlet.

A narrow ringlet referred to as D68 lies near the inner edge of Saturn's D ring, about 67,630 km from Saturn's center. From its discovery in Voyager images (Showalter 1996) through much of the Cassini mission, investigation of D68 focused on its radial profile and phase angle properties (Hedman et al. 2007). Later studies brought attention to its longitudinal brightness variations (Hedman et al. 2014). In 2014-15, four bright clumps formed and remained relatively evenly spaced with small longitudinal variations about mean separations of 26 • , 32 • , and 29 • (Hedman 2019). Hedman (2019) investigated these clumps in depth and designated them T (trailing), M (middle), L (leading), and LL (leading leading). The most likely explanation for the sudden increase in brightness in the four clump regions of the ringlet is that fine material was released by collisions into or among larger objects located near or within D68. These hypothetical larger objects are called source bodies, whose minimum sizes can be constrained by estimating the amount of material associated with each clump from phase-corrected normal equivalent area values, and whose maximum sizes can be constrained by the fact that they have not been observed directly. The range of masses that would correspond to these size constraints is 10 5 − 10 10 kg. The narrowness of the D68 ringlet and the distribution of clumps could either indicate that there is a collection of source bodies in a co-orbital configuration or imply that there is some outside force confining this material. In this paper we test the first idea by modeling the D68 clumps as a co-orbital satellite system.
The study of the dynamics of co-orbital systems is motivated by the many cases of co-orbital systems we find in our solar system. We are especially interested here in systems in which the co-orbitals have comparable masses. The best known of such systems are the horseshoe orbits of Janus and Epimetheus (Dermott & Murray 1981). Co-orbital asteroids have been suggested as the source of Venus's zodiacal dust ring (Pokorný & Kuchner 2019). Finally, the ring arcs in the Neptunian system have been proposed to be confined by either a corotation resonance with a moon on a separate orbit (Goldreich et al. 1986;Porco 1991;Salo & Hanninen 1998;Namouni & Porco 2002) or a co-orbital resonance with an undetected moon or even multiple moons shar-ing the same orbit (Lissauer 1985;Sicardy & Lissauer 1992;Renner et al. 2014).
In Section 2, we analyze potential stable configurations. In Section 3, we describe how we use numerical simulations to investigate these scenarios. In Section 4, we discuss some remarks for co-orbital systems as well as the possibilities for D68.

ANALYSIS OF POTENTIAL STABLE CONFIGURATIONS
Here we first review the theory of stable co-orbital objects, and we then apply the theory to the D68 clumps.

Theory
Salo & Yoder (1988) originally examined stationary configurations of equal-mass co-orbital satellites for small N (N ≤ 9) using a simple first-order theory, neglecting terms of the order (m/M ) 3/2 , where m and M are the masses of the satellite and the primary. A numerical search revealed three distinct types of stationary solutions, of which we are here concerned with only one, which Salo & Yoder (1988) label Type Ia: an equilibrium where all the N satellites are most concentrated on the same side of the common orbit. The case where N = 2 is the well known Trojan configuration, with an angular separation of 60 • . Type Ia configurations are stable for 2 ≤ N ≤ 8 but are not found for N ≥ 9 (Salo & Yoder 1988). This study, motivated by the D68 clumps, focuses on configurations with N = 4 and N = 5. Renner & Sicardy (2004) generalized the work of Salo & Yoder (1988) for similar but not necessarily equal masses, which is what we expect for the D68 clumps. When we define φ i as the longitude of satellite i and ξ i = ∆r i /r 0 as its relative radial excursion with respect to its average radius r 0 , the relevant equations of motion become (Renner & Sicardy 2004) where m j is the mass of satellite j and The derivative of f (φ) is For a co-orbital stationary configuration (Renner & Sicardy 2004), Equation 6 can be written in matrix form: Because the N × N matrix is antisymmetric and depends only on the longitudinal separations φ i between the bodies, for N ≥ 3 one can always find a set of relative masses that satisfies these equations for any given set of angular separations. This solution, however, might not be physical because one or more of the masses could be negative or zero.

Results
We first considered the observed configuration with four masses separated by angles of 29 • , 32 • , and 26 • because these are the observed separations (Hedman 2019).
These separations are closer than the expected separations for an equal-mass situation: 41.498 • , 37.356 • , and 41.498 • (Salo & Yoder 1988). We therefore solved the above equations for arbitrary masses, using Gaussian elimination, which involved re-ordering the rows, and found that the solution contains a mass that is calculated as either zero or a small negative number on the order of 10 −16 − 10 −19 , depending on the order in which the rows are solved (most likely a numerical issue involving the limit of double precision numbers): when we normalize the relative masses such that their sum is 1. Thus, there does not exist a physical solution for the stationary configuration with four objects that would produce four comparably bright clumps.
There are two possible ways that the clumps could still reflect a collection of co-orbital source bodies: the four source bodies could have been librating around the equilibrium location or there could be another massive body in the system that did not produce a visible clump. It is certainly possible for there to be only four non-stationary clumps and for this to be a transient phenomenon. In fact, Hedman (2019) identified slow changes in the clumps' azimuthal separations over time that could be evidence for libration. It is unlikely, however, for the clumps to be on the edge of a libration cycle, due to how azimuthally compact the whole configuration is. The most compact state of a configuration of three is in the symmetric mode when the outer bodies are at their closest approach to the middle body. A similar symmetric mode in a system of four bodies would require the outer two bodies to converge at a faster rate than the middle two bodies. The observed drift rates, however, show the opposite trend, with the middle two clumps drifting at a faster rate than the outer clumps (Hedman 2019). If, however, the dust around four source bodies was stirred up by an object that passed nearby, it is certainly possible this object could have missed other source bodies in the D68 ringlet. We therefore consider whether there could be an unseen fifth object, whose mass could account for the angular separations we observe between the four known clumps. We explore the approximately 270-degree span of longitudes ahead of Clump LL and behind Clump T. Using the same equations of motion (Renner & Sicardy 2004), we find physically realistic solutions in two regions, one centered 33 • ahead of Clump LL and one centered 32 • behind Clump T. These regions each span about 22 • in longitude and are mapped out in Figure 1. The relative masses of the clumps that correspond to these solutions are plotted in Figure 2. The horizontal axis shows the longitude of Object 5 in the same longitude reference system used by Hedman (2019). The left-hand side of the split horizontal axis corresponds to a configuration in which Object 5 is trailing the other D68 clumps; the right-hand side corresponds to a configuration in which Object 5 is leading the other D68 clumps. In more compact configurations (when Object 5 is near longitudes -80 and 55), the middle and outer masses are greater than the second and fourth masses. In less compact configurations (when Object 5 is near longitudes -95 and 75), Object 5's mass would be more than double that of any other mass, and the clump farthest away from Object 5 becomes the least massive while the other three would require comparable masses.
Figure 2. This plot shows the relative masses of the five co-orbitals for each possible configuration. Compact configurations are characterized by the more massive bodies in positions 1, 3, and 5. Extended configurations, by contrast, have the most massive object on one end, with the other masses tending to decrease with increasing distance.

NUMERICAL INVESTIGATIONS OF THE CONFIGURATIONS' STABILITY
To further examine the dynamics of a co-orbital system at the semi-major axis of the D68 ringlet and to investigate stability limits, we numerically simulated the motion of point masses at the longitudes of the clump peaks, adding in a fifth point mass at one of the locations permitted by the methods found in Renner & Sicardy (2004). For orbital simulations, we used the hybrid symplectic/Bulirsch-Stoer algorithm in the Mer-cury6 package (Chambers 1999). Our orbital simulations considered Saturn as the central mass and included terms up to J 6 in its gravitational field. The constants used for these simulations were taken from Jacobson et al. (2006) and Archinal et al. (2018), converted to the units used in Mercury6, and are found in Table 1. We used a time step of 0.02 days, which for D68 corresponds to about one tenth of an orbit.
For the sake of simplicity, we focused on one specific stable solution with the corresponding angular separation of Object 5 in order to do numerical simulations, though other configurations were also investigated, both on the leading and trailing sides, to ensure that our conclusions are general. We focus on a configuration with Object 5 ahead of Clump LL by 33 • , as specified in Table 2.
We explored perturbations to this configuration in semi-major axis and longitude, modifying the initial semi-major axis or longitude for some of the bodies. We  also varied their absolute mass, while keeping their relative masses constant, as calculated above (Renner & Sicardy 2004). Although the highest mass range we expect for the clump source bodies is 10 9 −10 10 kg because they have not been observed directly (Hedman 2019), we also consider much more massive configurations because these evolve more quickly and in this way clarify how these systems respond to perturbations. Thus we consider three different situations: one with extreme masses of 10 20 − 10 21 kg (i.e., similar to Enceladus, Tethys, and Dione), one with with masses of 10 13 −10 14 kg (i.e., similar to Polydeuces, Pallene, and Daphnis), and one with masses of 10 8 − 10 9 kg, close to that expected for the D68 source bodies.
In each simulation, we plot the longitudinal evolution of the bodies with respect to a reference longitude, which is calculated for each timestep as where λ i is the mean longitude of body i. This equation works well when the longitudinal oscillations are small. This type of plot gives a quick sense of stability and of orbital evolution.
We verify that the stationary points found using the method of Renner & Sicardy (2004) are indeed stable by placing objects there and finding they do not evolve in 1000-year simulations with high masses (10 20 − 10 21 kg; see Figure 3). Here we do not explore perturbations in initial longitude or semi-major axis for the high-mass case because the larger masses complicate scalings to the real system. We consider two types of perturbation, longitudinal and radial, in the medium-mass case, 10 13 −10 14 kg. The objects are massive enough that it is easier to demonstrate both stable libration and more chaotic mutual encounters. First, we consider a longitudinal shift in which the system begins in a more compact configuration, and we find that the masses oscillate stably around the solution (see Figure 4).
Second, we consider radial perturbations in which we modify the initial semi-major axis. We define a critical semi-major axis separation ∆a crit which separates small oscillatory motion like that shown in Figure 5 from the sort of motion shown in Figure 6. We explored through simulations the allowable perturbations to semi-major axis using the mode in which Clump LL is given a positive ∆a and Clump M is given a negative ∆a, just as in Figure 5. We found that, for these relative masses in this specific perturbation mode, the critical semi-major axis separation's relation to absolute mass is best represented as  m are small enough that when two of the bodies approach each other, they exchange energy and angular momentum in such a way as to begin receding from each other, similar to the periodic orbital swap of Janus and Epimetheus. Perturbations of ∆a = 200 m are too much for a stable configuration, which results in bodies looping around to approach the other side of the co-orbital system and eventual spreading into multiple orbits via gravitational interactions with the other bodies. Figure 5. With sufficiently small initial perturbations to semi-major axes (50 m for Daphnis-scale co-orbitals), the bodies oscillate around a stable solution, which is indicated by the dashed lines.
To apply our numerical simulations to the D68 clumps, however, we must also consider the dynamics Figure 6. With large enough initial perturbations to semimajor axes (200 m for Daphnis-scale co-orbitals), the system becomes unstable when some of the bodies encounter each other. Figure 7. With realistic masses, even semi-major axis perturbations of one meter result in system instability. Although low-mass co-orbital systems are fragile, stability could be achieved with the help of an external force.
in a low-mass case, 10 8 − 10 9 kg. For the low-mass case, ∆a crit = 27 cm, which is confirmed by Figure 7. With only 1-m perturbations (a 2-m separation in semi-major axis), the point masses drift by each other, with the three closest approaches between the centers of any two bodies as 266 m, 365 m, and 400 m. 1 With a density of 0.5 g/cm 3 , spherical bodies of these masses would have radii ranging from 49 to 89 m. Thus, although these closest approaches would not be collisions, they would still be close enough gravitational encounters to provide significant perturbations, in a range of 4 to 6 Hill radii for the largest mass. We consider such a system to be fragile.

DISCUSSION AND IMPLICATIONS
To emphasize how fragile the system is, we can estimate the impulse required to perturb a moonlet's semimajor axis by 1 m, similar to what has been done in Hedman & Bridges (2020). For nearly circular orbits, the standard orbital perturbation equations can describe the rate of change of semi-major axis over time as (Burns 1976;Hedman 2018) where the mean motion n = GM/a 3 1751.7 • /day, F p is the azimuthal component of the perturbing force, F G = GM m/a 2 = n 2 am is the gravitational force on the moonlet from the planet's center, M is the planet's mass, and m is the moonlet's mass. The moonlet will thus undergo a semi-major axis change δa upon receiving an azimuthal impulse With the range of masses we use for our low-mass case, the impulse required to perturb a moonlet's semi-major axis by 1 m ranges from 4.2 × 10 4 kg m/s to 2.7 × 10 5 kg m/s. For any collision between a moonlet and interplanetary debris, the impact speed would be comparable to the D68 orbit speed v = na 24 km/s. Dividing the range of impulses by this orbit speed, we get a range of masses roughly from 2 kg to 10 kg for the interplanetary impactor. Assuming a density of 0.5 g/cm 3 , the piece of interplanetary debris would need to be 0.1 m to 0.2 m in radius. The estimated cumulative influx rate Φ for debris of this size is around 10 −17 /m 2 /s (Tiscareno et al. 2013). Thus the rough timescale t 1/ΦA on which we can expect such a collision, using the crosssectional area A for moonlets with radii of 49 m to 89 m, corresponds to a range from 130,000 years to 420,000 years, but this is the impact timescale for just one of the objects. Because an impact into any of the objects can break the system, we can adjust the system timescale to about 40,000 years by adding their cross-sectional areas together. We therefore cannot expect a co-orbital configuration at D68 to last longer than a few tens of thousands of years. For this reason we call the system fragile and find it unlikely that a co-orbital system could explain the orbital evolution of the clumps or the ringlet. Table 3. Predicted corotation resonance locations, the rM14 column corresponding to our predictions based on the pattern periods reported in Marley (2014), the rM19 column corresponding to our predictions based on the pattern speeds reported in Mankovich et al. (2019)  Consequently, we look for other resonances that could drive the orbital evolution of the clumps or the ringlet. It is unlikely that a corotation resonance with any satellite is responsible for the clumping of material into ring arcs. A 30 • separation between clumps would be the result of a 12-fold pattern at the D68 semi-major axis. A 12-fold pattern could be caused by a 13:12 corotation resonance with an external perturber or an 11:12 corotation resonance with an internal perturber. A 13:12 corotation resonance with an external perturber would require a perturber at a semi-major axis of 71,300 km, which is not as far out as D72, the structure closest to D68. An 11:12 corotation resonance with an internal perturber would require a perturber at a semi-major axis of 63,800 km, which is a few thousand km away from Saturn's equatorial radius (60,268 km). There is no evidence for any moons or ringlets in these regions. Moreover, no results came from a numerical search for corotation resonances up to fourth-order between D68 and Janus, Mimas, Enceladus, Tethys, Dione, Rhea, or Titan.
It is possible that a resonance of some sort with Saturn itself could be responsible for the D68 clumping. The outer Lindblad resonance of Saturn's = 5, m = 3 oscillation mode is located in the D68 region, reported first at 67,625 km ± 550 km (Marley & Porco 1993) and more recently at roughly 67,550 km (Marley 2014). Although Lindblad resonances do not confine material, each such resonance can be associated with a corotation resonance, which can confine material. To locate these corotation resonances, we computed the radii at which the mean motion (using the second-order equation from Renner & Sicardy 2006) matches the pattern speeds associated with the modes reported in Marley (2014) and Mankovich et al. (2019). The modes that produce corotation resonances near D68 are listed in Table 3. Because the pattern speed is dependent on Saturn's structure, any of these modes could possibly be responsible for providing a corotation resonance to confine D68 material. Mode splitting or mixing could also be involved (Fuller 2014), allowing the locations of these resonances to fall at slightly different radii than what we can compute. For a corotation resonance of one of these modes to explain the D68 clumping, it would require a planetbased angle that moves at a speed comparable to D68's mean motion. Although the set of angular separations among the clumps favors a 12-fold pattern, it is also possible for them to be confined to within the libration longitude of one or a few stable points of a lower m mode, and then be spaced out within that external potential. Perhaps there is a set of co-orbital moonlets that are trapped together and librating within a larger potential, similar to Renner et al. (2014)'s model of the Neptune ring arcs. Radial oscillations of ± 10 km have been observed for the D68 ringlet with an estimated period of 14-15 years (Hedman et al. 2014), though the clumps are drifting more slowly than the rest of the ringlet (Hedman 2019). These radial oscillations could be evidence for that libration.
In conclusion, we have tested and ruled out long-term stable co-orbital configurations as an explanation for the spacing of the D68 clumps. We therefore predict that either the clumping is a transient phenomenon, or that an external mechanism is trapping the clumps in this region.