Numerical Simulations of (10199) Chariklo's Rings with a Resonant Perturber

For the simple case of spherical central bodies and ring particles, the classical Roche limit is given by a_Roche = R(2δ)^1/3, where R is the radius of the primary and δ is the density ratio of the primary to the ring particles. Assuming a spherical equivalent radius for Chariklo of 129 km (Leiva et al. 2017) and δ = 1, the Roche limit at Chariklo is <100 km. In order for Chariklo's rings to reside within the Roche limit (a_Roche ⪆ 400 km), the density of the ring particles would have to be ≤0.2 the density of the nucleus. In comparison, the density of Saturn's ring particles is typically ≥0.5 the density of the planet, and these are likely the least dense ring particles of the giant-planet systems (e.g., Tiscareno et al. 2013).

A more detailed definition of the Roche radius is ${a}_{\mathrm{Roche}}={(3M/\gamma {\rho }^{{\prime} })}^{1/3}$, where M is the body's mass, ${\rho }^{{\prime} }$ is the density of the orbiting material, and the factor describing particle shape γ = 0.85 for classical calculations while γ = 1.6 has been preferred to represent the tidal destruction limit (e.g., Porco et al. 2007; Sicardy et al. 2019; Morgado et al. 2023). For M = 6 × 10¹⁸ kg (the solution for a Jacobi triaxial ellipsoid, Leiva et al. 2017), water-ice particles with ${\rho }^{{\prime} }=1\,{\rm{g}}\,{\mathrm{cm}}^{-3}$, and either value of γ, Chariklo's Roche limit is <300 km. Considering an ellipsoidal shape for the nucleus, nonspherical particles, and reasonable density ratios (δ < 5), other calculations for the Roche limit are also <400 km (Melita et al. 2017; Sicardy et al. 2019). The Roche limit can reach past the rings by assuming a mass of M = 8 × 10¹⁸ kg (the solution for a Maclaurin spheroid; Leiva et al. 2017), classical particle shape γ = 0.85, and particle density typical of the small, inner, Saturnian satellites, ${\rho }^{{\prime} }=0.4\,{\rm{g}}\,{\mathrm{cm}}^{-3}$, or lower (Thomas & Helfenstein 2020).

Morgado et al. (2023) pointed out that more elastic collisions than those at Saturn could prevent accretion in rings at distances well beyond the Roche limit; however, such impacts would cause viscous spreading. Until the physical characteristics of Chariklo's nucleus and ring particles are better known, the Roche limit will not be well constrained. Given the range of reasonable values, it is likely to be near or interior to the rings. Our explorations here thus focus on how Chariklo's rings can be maintained as well as why they do not accrete. The idea of a shepherd moon is bolstered by Saturn's F ring, which is outside the Roche limit, while perturbations from Prometheus and Pandora likely hinder reaccretion of material (e.g., Murray et al. 2005).

Our code stems from Lewis et al. (2011), which was originally developed for N-body simulations with geometry similar to that of Saturn's outer Encke gap edge with the satellite Pan. This code was used to demonstrate that a single satellite orbiting near a dense ring can cause the ring to collapse on a timescale of a few orbital periods (Lewis & Stewart 2005, 2009). Here, smaller-scale systems are modeled by scaling down the nucleus mass and the orbital distances and sizes of the rings. The simulations employ a linearized, closed-form solution to Hill's equations and utilize guiding-center coordinates (Stewart 1991). A few million particles are located in a local cell that slides downstream past a nearby satellite. The cell has open boundaries in the radial direction, and the azimuthal boundaries are semiperiodic. The boundary conditions are designed to preserve the gradient in the epicyclic phase that is induced by the passage of the satellite. The nucleus and satellite are both treated as point masses.

The original version of the code employs hard-sphere, dissipative collisions with a velocity-dependent coefficient of restitution (from Bridges et al. 1984), meaning that the collisions are treated as discrete events and handled in order by time. All ring particles are spheres. Particle spin and self-gravity are model options: when turned on, particle self-gravity is calculated using a tree method. The same tree is used for both gravity calculations and to find nearby particle pairs that need to be checked for collisions.

Adjustable model parameters include (i) cell size and aspect ratio; (ii) initial distribution of particles; (iii) nucleus mass; (iv) satellite mass and location; (v) ring particle minimum and maximum size, size-distribution power law, compositional density, and (vi) self-gravity. Initial simulations were conducted in a square cell ∼28 × 28 km, large enough to span both C1R and C2R radially; however, the rings do not interact with each other dynamically and can thus be studied separately. For the remainder of the simulations, the size of the cell was selected to be ∼7 × 7 km in order to span the anticipated width of only the larger ring, C1R. All simulations use an initial condition with the surface densities of particles having Gaussian distribution(s) in the radial direction and uniform distribution(s) in the azimuthal direction. The radial distribution is centered at the observed ring location(s) from Braga-Ribas et al. (2014). We note that these simulations are all for perturbed circular rings.

For the rings to be present for long periods, they need to be confined yet not accrete into a single body. We explore confinement and accretion with separate simulations. To look at confinement by a single shepherd, we ran a suite of simulations with parameters listed as numbers 1 through 9 in Table 2: parameters for the accretion simulations are shown as simulation number 10. These parameters were chosen to begin understanding the effect of satellite location, satellite mass, and the characteristics of the ring material (through particle sizes, particle densities, and model ring optical depth, τ_m ).

In these simulations, ring material is confined as a result of the gravitational interaction between the ring particles and the satellite. A collisionally damped wake is generated by the satellite in the ring, which reduces ring width. Physically, a strong satellite wake causes negative flux of angular momentum across the ring, resulting in particles moving to regions of higher density. In the reference frame rotating with the satellite, the forced eccentricity induced by the passage of the moon appears as a standing wave moving away from it. This wave is sheared due to differential rotation, causing regions of compression and rarefaction indicative of wakes (Showalter et al. 1986). The diffusion coefficient becomes negative, which we refer to as "negative diffusion," when the orbits of particles with nearby semimajor axes cross. The wake damping is an effective confinement mechanism over a range of properties; however, it requires appropriate particle sizes and optical depths for particle collisions to be dynamically significant.

Model data for a full ring are created by joining together values from the cells at different times over the course of one synodic period. Once the system has reached equilibrium, defined by ring characteristics having constant values, the orbital elements of particles at a given location are determined by their position relative to the perturber, and any synodic period of data is a valid snapshot of the full ring. The resonant particles reach a state of relative equilibrium fairly quickly; however, breaks appear in the unequilibrated material outside of resonance, as the perturbations from one satellite pass to the next do not yet align. As a result, the model ring plots in this work are continuous for the material in resonance.

2.2.1. Soft-sphere Model

2.2.2. Simulation Times

To study accretion, we wrote a data-analysis tool specifically to track the formation and evolution of moonlets. We define a moonlet as a group of particles of a sufficient size where all the particles are in contact and gravitationally bound. A particle is gravitationally bound to a group if the gravitational potential energy is larger than the relative kinetic energy. To make moonlet identification computationally tractable for simulations with millions of particles, we start by binning the particles into a grid where the grid size is several particle diameters across. We then search for groups in order from the grid cell with the most particles down to a minimum threshold. Each particle that is found to be bound is tagged as being in that particular group. The search moves outward in rings of cells around the central one, and the particles in each new cell are tested in order based on distance from the center of mass when that cell is started. If any particles are encountered that are already part of another group, they are ignored. Groups that do not meet a minimum size requirement are also ignored.

The processing of the first time step uses only the grid to pick a location at which to start looking for groups. After the first step, for each group found in the previous step, the "core" particle of the group is defined to be the particle whose center is closest to the center of mass of the group. The search for groups begins with each of those "core" particles. This method makes the code more capable of tracking the same groups, and it is particularly important when moonlets run back into any dense ring wakes. We continue tracking cores for several time steps even if the surrounding groups do not reach the threshold, further improving our ability to track moonlets moving through any wakes and to track separate moonlets after a nonmerging collision. Checking whether a new particle in a moonlet overlaps existing ones can also be computationally expensive. So, once a certain number of particles are in a group, we create a small spatial grid with a cell size equal to a particle diameter. This way, the search for overlapping particles only has to include particles in adjacent cells.

The bulk of our simulations were focused on exploring the confinement of ring material by single-sided shepherding. We simulated a two-ring system with different satellite locations, and then a one-ring system with different satellite masses and model ring optical depths. Chariklo was assumed to have a mass of 7 × 10¹⁸ kg (the middle of the range from Leiva et al. 2017).

3.1.1. Satellite Location

When modeling both rings, simulations 1–4 in Table 2, we used model optical depths of 0.4 and 0.06 for the inner and outer rings, respectively. The satellite was assumed to have a mass of 3 × 10¹³ kg and location from 320 to 470 km. Ring-particle sizes were 3–10 m in radius, following a standard power-law distribution ∼r^−q with fixed index q = 3 (following the power law determined for Saturn's rings; see, e.g., Cuzzi et al. 2009; Eckert et al. 2021).

Figures 1 and 2 show the results from these simulations, aimed at exploring the primary effects of a satellite in a Chariklo-like ring system. When there was no satellite in the system, the ring material dispersed as expected. When there was a satellite but it was not in a resonance with the initial ring locations, the material also dispersed. The inner ring width as a function of time was not linear, as in the case of no satellite, because the material was effected by the satellite; however, the net result remained that the material was unconfined. In both cases, the expansion of ring material was slow, as shown by the m-scale increases in the widths of the inner ring over multiple orbits on the right side of Figure 1.

Zoom In Zoom Out Reset image size

A satellite in resonance with the rings acted to confine the particles to km-level widths over a relatively short timescale. Examples are shown in Figure 2 for a satellite in both an inner and outer mean-motion resonance (MMR) with the rings. In these cases, the inner ring was on a first-order MMR with the satellite and another resonance was near the inner part of the outer ring. The effects of the satellite perturbation are apparent in both rings as azimuthal and radial variations in the distributions of particles. We note that the effect of a satellite in an interior resonance can be considered dynamically similar to that of one in an exterior resonance, or that of a non-axisymmetric nucleus that is in a rotational resonance with the rings (as considered by Sicardy et al. 2019).

Zoom In Zoom Out Reset image size

3.1.2. Satellite Mass and Ring Optical Depth

To explore the effect of the satellite mass, as well as a broader range of ring-particle characteristics, we carried out simulations 5–9 in Table 2 for the inner ring only. For these simulations, we fixed the satellite location at the 4:5 MMR. This resonance was chosen based on the initial simulations, because it confined ring material at the approximate locations of C1R and C2R within a reasonable timescale. Placing the satellite in other resonances will change the strength and number of perturbations, resulting in different ring characteristics and timelines for confinement (e.g., Figure 2). We explored a range of satellite masses from 7.5 × 10¹² to 6 × 10¹³ kg (∼10⁻⁷–10⁻⁶ the mass of Chariklo) and a range of ring model optical depths from 0.02 to 0.4. Ring-particle sizes were 1–10 m in radius, following a power-law distribution with fixed index q = 3. The particle density was 0.5 g cm⁻³.

Figure 3 shows the ring guiding-center full width at half maximum (FWHM) over time of a subset of these simulations, for three different satellite masses and a ring of two different model optical depths. Setting an initial condition of τ_m = 0.4 in these simulations produced clumping of material, as well as an increased collision rate, and it slowed the computations to the point of being unreasonable. Simulations 1–4 did not demonstrate this, because particle self-gravity was not included, but here we took the more physically realistic approach of increasing particle density and reducing the initial model optical depths. See Section 4 for further discussion.

Zoom In Zoom Out Reset image size

From Figure 3, satellites with the lowest mass (1.5 × 10¹³ kg) maintained but did not confine the ring to the observed widths. A satellite with a slightly higher mass contracted the ring to the few-km widths observed for C1R. However, as optical depth increased, the mass of the satellite required for confinement also increased. Satellites with higher mass contracted the ring even further, to the scale of a single km.

The width of the simulated ring from the middle right panel in Figure 3 is shown in Figure 4. This plot demonstrates that the actual ring width varied over azimuth. The azimuthal width variations in this simulation were approximately two km.

Zoom In Zoom Out Reset image size

In Figure 5, we have combined the results from simulations 5–9 to roughly map ring confinement as a function of satellite mass and initial optical depth of the inner ring. Material is defined to be confined with a visual inspection of the density of guiding centers over the course of the simulation. Satellite mass and ring optical depth were positively correlated, in an approximately linear relationship, in terms of confining material to the few-km width observed for Chariklo's inner ring. In addition, tighter confinement occurs for simulations that have lower initial optical depths (and thus lower initial collisional rates).

Zoom In Zoom Out Reset image size

We next explored the impact of the shepherd-moon wakes on accretion. The confinement mechanism from single-sided shepherding produces wakes with a significantly higher density than the general ring material, where particles can easily form clumps or gravitationally bound moonlets. Lewis & Stewart (2005) found that material leaving the satellite wakes in simulations of the Encke gap formed extended gravity wakes. While this mechanism could enhance accretion, the dynamics of the satellite wakes produces regions of rarefaction in addition to the compression. This occurs when the perturbed streamlines separate. Our simulations addressed two hypotheses relating to accretion in the Chariklo system: (i) particles that are not sufficiently bound together will shear out, and (ii) collisional interactions between the wakes and the moonlets will cause more erosion than additional accretion.

We ran two simulations with parameters that favored accretion: one used the standard hard-sphere model and the other used the soft-sphere model, following Section 2.2.1. The satellite was 3 × 10¹³ kg in the 4:5 MMR, the initial ring optical depth was 0.1, and the particle density was 0.5 g cm⁻³ (Table 2). To promote accretion, these simulations considered a smaller range of particle sizes (1–1.1 m in radius with q = 3) and Chariklo had a smaller mass (5 × 10¹⁸ kg). Using the size from Morgado et al. (2021), this mass corresponds to a density of 0.62 g cm⁻³ for Chariklo. These parameters give a scaled Hill radius of 1.27, which is significantly larger than the value of 1.1 found to lead to accretion in previous ring studies (see Section 3.2 of Michikoshi & Kokubo (2017), for references). In other words, the ring was well outside the classical Roche limit in these simulations.

The accretion-testing simulations were run for just over 30 orbits, long enough for the resonant satellite wakes to form and for self-gravity to begin the process of accreting particles. We found and tracked moonlets from the time of their formation to when they fell apart or the end of the simulation, following the method in Section 2.3. Results from the simulations at three specific time steps are shown in Figure 6. Figures 7–9 contain close-ups of a selection of the moonlets detected during these time steps.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 7 shows the state of the hard-sphere simulation during the 32nd orbit, in which the density wakes were easily visible and there were two moonlets. Moonlets formed as clumps of material fell out of the passing satellite wake, where densities were elevated; however, the material was only tenuously bound and sheared out with even small perturbations. The top moonlet lasted the longest of any in this simulation: no moonlets were stable longer than half an orbital period.

Figure 8 shows the state of the soft-sphere simulation early in the 27th orbit. As with the hard-sphere simulation, material came out of the satellite wake in clumps. Unlike the hard-sphere simulation, the highest-density clumps formed moonlets that were stable for many orbits. Figure 8 highlights nine of the more than 40 moonlets present at this time, most of which display tails of sheared material.

Zoom In Zoom Out Reset image size

Figure 9 shows the soft-sphere simulation near the beginning of the 31st orbit. At this point, a reasonably large fraction of the ring material was contained in just over 30 moonlets. The satellite wake was not as defined because the effective optical depth had been dramatically reduced. The state of the ring was similar to the oligarchic growth phase seen in planetary-formation simulations. The moonlets also had crossing orbits, so they could collide. Such collisions were rare and more energetic than earlier in the simulation, although the collisions at this time were generally not destructive.

Zoom In Zoom Out Reset image size

A plot of the number of moonlets detected along with the numbers of particles in those moonlets as a function of time for both the hard- and soft-sphere simulations is shown in Figure 10. To remove any "false positives" in the wake peaks, only clumps of particles that survived for at least one-third of an orbit were counted. In both simulations, the number of moonlets varied throughout an orbital period because the satellite wake created and destroyed clumps of material. No moonlets survived long-term in the hard-sphere simulation. Indeed, when the criteria to rule out false positives was increased to one-half of an orbit, there was never more than one moonlet present at any time in the hard-sphere simulation. The oscillations were created by the fact that moonlets formed as material dropped out of the wake peak, then fell apart because of various collisions with other ring material over the course of the next orbit.

Zoom In Zoom Out Reset image size

In the soft-sphere simulation, the number of moonlets increased drastically to nearly 50 before leveling out to roughly 30 stable moonlets at the end of the simulation. These simulations included roughly 2.5 million particles with a narrow size distribution. In roughly 30 orbits, nearly a third of the simulation material was entrained into approximately 30 larger moonlets. The final moonlet count leveled off when the wakes no longer had sufficient material to form new moonlets and the mean intercollision time for the moonlets was large, much like in the oligarchic growth phase seen in planetary simulations.

The satellite that produced the ring widths shown in Figure 4 would have a radius of ∼3 km, assuming a reasonable density of 1.5 g cm⁻³. Although Chariklo is not currently known to have any satellites, the existence of a few-km-sized moon at Chariklo is not unrealistic. The angular size of a 2 km body at Chariklo's current distance is <0.1 mas, with the entire ring system extending ∼75 mas. Direct imaging, for example with the Hubble Space Telescope, has not been able to resolve the rings nor anything else closer than 1000 km from Chariklo's center (e.g., Bérard et al. 2017). In addition, a few binary Centaurs have been detected (e.g., Noll et al. 2006, 2008) and 50%–60% of one of the Centaur's parent TNO populations, the Plutinos, are thought to consist of tight binary systems (Di Sisto et al. 2010; Thirouin & Sheppard 2018).

An important factor in the simulations and the observations is ring optical depth. Optical depth is related to both particle sizes and surface number densities: a simplistic example is that higher optical depth can be caused by a few large particles or many small particles (e.g., van de Hulst 1957). The particle size distribution is thus a key parameter for ring modeling, which is currently unknown for Chariklo.

Assuming that the apse alignment is maintained by the ring's self-gravity, Pan & Wu (2016) estimated a typical particle size of a few meters. On the other hand, Michikoshi & Kokubo (2017) argued that the Centaur rings are observable only given ring particles smaller than a meter or with the existence of shepherding satellites, otherwise the viscous spreading occurs on too short a timescale (100 yr). Analysis of dual-wavelength data (visible and red) from a stellar occultation by Morgado et al. (2021) showed no dependence of ring opacity with observational wavelength, suggesting that the particle sizes are larger than a few microns.

The simulations presented here have a minimum particle radius of 1 m. The size distributions of planetary rings are typically well-described by a power law with a differential slope that is characterized by a steep increase in the number of particles at smaller sizes. For this reason, it is impractical from a computational standpoint for ring models to consider particles at scales smaller than roughly a meter for a cell this size, artificially limiting the optical depth of systems that can be effectively simulated (or ignoring particle self-gravity, as in simulations 1–4 in Table 2). However, stellar-occultation observations are sensitive to small particles, down to submicron sizes for visible-wavelength data. Therefore, caution must be taken when comparing modeled and observed optical depths.

The optical depths reported for our simulations are geometric optical depths of the fractional area coverage, calculated by

$$\begin{eqnarray}&&{\tau }_{m}=\displaystyle \frac{\displaystyle \sum _{i}\pi {r}_{i}^{2}}{A},\end{eqnarray} \tag{ 1 }$$

where r_i is the radius of the ith particle and A is the area of the simulation cell in the ring plane. Observed optical depths are calculated from stellar-occultation data. Stellar occultations measure the transmission, T, of light when a star is blocked by ring material. The optical depth along the line of sight, τ₀, is calculated by $T={e}^{-{\tau }_{0}}$, and the normal optical depth, τ_N, is calculated given a ring opening angle, B, by ${\tau }_{{\rm{N}}}={\tau }_{0}\sin | B|$ (for a polylayer ring, e.g., Elliot et al. 1984). Reported optical depths for Chariklo's rings, the proposed rings at Quaoar, and the recent recalculation of the characteristics of the features at Chiron have included a factor of two in this calculation, ${\tau }_{{\rm{N}}}=\tfrac{{\tau }_{0}}{2}\sin | B|$ (Braga-Ribas et al. 2014; Bérard et al. 2017; Braga-Ribas et al. 2023; Morgado et al. 2023; Pereira et al. 2023). This factor is based on inferred physical properties of the rings, assuming that the rings are thin and the occultation optical depth is larger, by an extinction efficiency factor of two, than the area physically filled, similar to what is observed for the Uranian rings (e.g., Cuzzi 1985; Roques et al. 1987). Published values for normal optical depths for material surrounding Chiron were based on the standard calculation (Sickafoose et al. 2020, 2023). It is important to take this factor of two into consideration when comparing published ring characteristics and determining model parameters.

Here, we have used the normal optical depths for Chariklo given by Braga-Ribas et al. (2014) while not employing self-gravity, or we have used significantly lower model optical depths with reasonable particle densities. A factor-of-two increase from the published optical depth for C1R would be too dense for our simulations to proceed reasonably, noting again that occultations can detect much smaller particles sizes than those in our simulations—and thus model and observed optical depths need not be equivalent.

The model presented here could apply to other small-body ring systems. The Centaur Chiron is similar in size to Chariklo (spherical estimate of 218 ± 20 km and ellipsoidal with semi-axes of 114, 98, and 62 km; Lellouch et al. 2017; Braga-Ribas et al. 2023) with a similar rotational period (5.92 hr; Marcialis & Buratti 1993), and its orbit ranges from 8.5 to 19 au. Stellar occultations in the early 1990s revealed an asymmetric dust coma at Chiron as well as narrow jets (Elliot et al. 1995; Bus et al. 1996). A combination of stellar occultations, rotational light curves, and long-term photometric and spectroscopic variations have been used to support the idea of a two-ring system with additional surrounding material (Ortiz et al. 2015; Ruprecht et al. 2015; Sickafoose et al. 2020). The proposed Chiron rings are at 300 and 309 km, are 2.5–4.5 km in width, and vary azimuthally in width by 1.5 km (inner ring) and 0.5 km (outer ring) (Sickafoose et al. 2020). However, the most recent occultation observations suggest that the situation is not as simple as a stable, two-ring system (Ortiz et al. 2023; Sickafoose et al. 2023). Further observations are needed to better characterize material around Chiron, but the demonstrated applicability of our simulations for Chariklo suggests that they could also be relevant for any stable rings at Chiron.

Haumea was the second small body with a confirmed ring system (Ortiz et al. 2017). Unlike Chariklo and Chiron, Haumea is a large TNO with even faster rotation and known satellites: assuming a triaxial ellipsoid, the size is estimated to be 2322 × 1704 × 1026 km, its orbit ranges from 35 to 52 au, and the rotational period is 3.9 hr (Rabinowitz et al. 2006; Ortiz et al. 2017). Haumea has a collisional history, which resulted in the formation of its two satellites and more than a dozen smaller TNOs known as Haumea-family objects (e.g., Brown et al. 2007; Ragozzine & Brown 2009; Schlichting & Sari 2009; Leinhardt et al. 2010). A single, coplanar ring was detected at Haumea, with a width of ∼70 km at a radius of ∼2287 km, which is inside the Roche limit (Ortiz et al. 2017). In size, body shape, and ring width, the Haumea system is different from Chariklo and Chiron. Although the known satellites could suggest that Haumea is a good candidate for our simulations, they orbit far enough away from the ring that they do not affect it dynamically. Our simulations also specifically maintain narrow rings a few km in width, and they are thus likely to be less applicable to Haumea than to the Centaur systems.

Most recently, one or two possible rings were reported around Quaoar (Morgado et al. 2023; Pereira et al. 2023). Quaoar is a ∼1110 km diameter TNO, orbiting from 42 to 45 au, with a rotational period of either 8.84 hr (single-peaked) or 17.68 hr (double-peaked) (Trujillo & Brown 2003; Brown & Trujillo 2004; Braga-Ribas et al. 2013; Morgado et al. 2023). It is known to have one satellite, Weywot (e.g., Fraser & Brown 2010). Quaoar does not have a homogeneous ring: the majority of Quaoar stellar-occultation light curves have not detected any material or only during either ingress or egress (Morgado et al. 2023; Pereira et al. 2023). There is a wide range of reported normal optical depths, with only one significant detection having a low-error measurement where τ_N ⪆ 0.1 (Morgado et al. 2023; Pereira et al. 2023). The proposed rings are located at 2520 and 4057 km, with possibly all—and certainly the more distant—material being well beyond the Roche limit (Pereira et al. 2023). Quaoar is known to have a moon, Weywot, for which the more distant ring is proposed to be in a 6:1 MMR, considering the double-peaked period (Morgado et al. 2023; Pereira et al. 2023). If stable rings were confirmed at Quaoar, this could be the most applicable system to model using our simulations.