© The Authors 2026

1 Introduction

Since 2013, narrow and dense rings have been discovered around four small objects of the outer Solar System: the centaur Chariklo (Braga-Ribas et al. 2014), the dwarf planet Haumea (Ortiz et al. 2017), the large trans-Neptunian object Quaoar (Morgado et al. 2023; Pereira et al. 2023), and possibly the centaur Chiron (Ortiz et al. 2023) (see also the reviews by Sicardy et al. (2024) and Sicardy et al. (2025a)). Current observations show that the rings of Chariklo and Haumea are close to the second-order 1/3 spin orbit resonance (SOR) with the central body, meaning that the ring particles complete one revolution, while the central body completes three rotations. The two rings of Quaoar are also near second-order resonance, the main ring being again near the 1/3 SOR, while the fainter ring orbits near the 5/7 SOR.

In Sicardy et al. (2025b) (Paper I hereafter), we studied the topological structure of the phase portraits corresponding to resonances of the first to fifth order. It appears from this work that among all these resonances, only those of order one (known as Lindblad resonances) and order two are expected to significantly perturb a collisional ring where eccentricities are damped by dissipative impacts. The reason for this is that only for these two kinds of resonances is the origin of the phase portrait (corresponds to zero eccentricity) not a stable (elliptical) fixed point.

This paper is the numerical counterpart of Paper I, first to validate the special status of the first- and second-order SORs mentioned above, and second to explore the ability of the 1/3 SOR to confine collisional rings. As mentioned in Paper I, only the first-order resonances force streamlines that do not self-intersect. As such, they can create nonsingular spiral waves that have been studied for decades in galactic and planetary ring contexts. In contrast, any periodic orbit forced by a resonance of order greater than one has at least one self-intersecting point, thus preventing in principle a regular response from a collisional disk. This is true in particular with the 1/3 SOR that forces streamlines with one self-intersecting point.

Here, we use numerical simulations to explore how the 1/3 SOR can confine material, in spite of the self-intersection problem. Compared to existing N-body simulation studies for small-body rings, our models include simultaneously all three ingredients necessary for a realistic modeling of SOR resonances: (i) a rotating nonaxisymmetric central body, (ii) particle physical collisions, and (iii) an azimuthally complete 3D ring. For example, Michikoshi & Kokubo (2017) performed large-N-body collisional simulations using an azimuthally complete ring, but assumed a spherical central body. Gupta et al. (2018) made collisional simulations with a nonspherical potential, but since the simulations were done in a local co-moving frame, they were limited to a case of axisymmetric nonrotating central potential. The two studies mentioned above thus missed SOR resonance effects.

On the other hand, Sickafoose & Lewis (2024) suggested that Chariklo’s rings could be stabilized by a Lindblad resonance with a nearby shepherd satellite. Their simulations used a modified local code, where the semiperiodic azimuthal boundary conditions were designed to correspond to the streamlines of the assigned first-order resonance. It is thus not obvious how to apply this method to 1/3 SOR or to higher order resonances in general. The mechanism itself, the collisional confinement of rings at first-order resonances, was already found in simulations of Hänninen & Salo (1994, 1995), which followed azimuthally complete rings. In the current study we demonstrate how a similar confinement takes place in rings perturbed by second-order resonances relevant for the observed systems.

This paper is organized as follows. In Section 2, we summarize our simulation method. In Section 3, we review the results of noncollisional test-particle simulations and their good agreement with the analytical calculations of Paper I. In Sections 4–7, collisional simulations are reported; after providing an illustration of the 1/3 SOR confinement, we establish a scaling of the simulation results to real systems and estimate the sizes of mass anomaly capable of confining ring material around small irregular bodies (Sect. 4), explore the normal modes excited inside confined ringlets (Sect. 5), and the outward migration of ring material through various resonances (Sect. 6). In Section 7, we demonstrate that the long-term stabilization of the ringlet is possible via a torque from a small additional external satellite. However, even in this case the location and confinement of the ringlet is determined by the 1/3 SOR. Finally, Section 8 shows that the confinement mechanism also works in self-gravitating rings. This paper describes in more detail the preliminary results obtained by Salo et al. (2021); Sicardy et al. (2021); Salo & Sicardy (2024) and Sicardy & Salo (2024).

2 Simulation method

Our numerical simulations treat a ring of inelastically colliding particles around a nonaxisymmetric central body rotating at angular speed Ω_B. Several models for the potential of the central body can be used (see details in Appendix A): (i) a homogeneous triaxial ellipsoid, and (ii) a spherical body with a dimensionless mass anomaly µ at its surface. These two models can also be combined so that (iii) the mass anomaly is attached to the triaxial ellipsoid. Finally, (iv) the ring particles can be perturbed by an additional satellite orbiting the central body.

Calculations extending over tens of thousands of central body revolutions are performed with IDL (Interactive Data Language) using a RK4 integrator in double precision. Impacts between ring particles are treated as soft-sphere collisions by including visco-elastic pressure forces between colliding, slightly penetrating particles. The treatment of collisions follows the schemes presented by Salo (1995), except that azimuthally complete rings are now simulated instead of local co-moving ring patches. In all our simulations, the adopted coefficient of restitution is ϵ_n = 0.1 (for more details, see Appendix B). Most of our simulations ignore mutual gravity between particles. However, Sect. 8 reports preliminary simulations with particle-particle calculation of ring self-gravity.

In this paper, a, e, L, and ϖ denote the osculating orbital semimajor axis, eccentricity, true longitude and longitude of pericenter of the particles, respectively, as calculated in the center-of-mass reference frame. The Appendix A also provides the definitions of the reference radius R_ref, the elongation ϵ_elon, and the oblateness f of a homogeneous triaxial ellipsoid, used later in this paper.

3 Results of test-particle integrations

3.1 Triaxial central body vs. mass anomaly

Spin-orbit resonances (SORs) occur near commensurabilities between Ω_B and the mean motion n of the ring particles (see Paper I). One kind of SOR is the corotation resonance that occurs for n = Ω_B, i.e., at the synchronous orbit with radius a_cor. This resonance is not studied here because the corresponding equilibrium points (similar to the Lagrange triangular points L₄ and L₅) are dynamically unstable or close to instability in the cases of Chariklo and Haumea (Sicardy et al. 2019). Moreover, they are also generally unstable against dissipative collisions since they correspond to local maxima of potential.

The other SORs occur for jκ = m(n − Ω_B), where κ is the particle epicyclic frequency, m is the azimuthal number of the resonance (positive or negative) and j > 0 is its order. These SORs excite the orbital eccentricities of the particles and, as shown in this paper, can confine rings in narrow regions. If the particles’ precession rate $\dot{\varpi }=n-\kappa$ is small compared to n, the above condition reads $$\frac{n}{{\text{Ω}}_{\text{B}}}\approx \frac{m}{m-j}.$$(1)

Following the notation of Paper I, m ≥ 1 corresponds to inner resonances (n > Ω_B), while m ≤ −1 corresponds to outer resonances (n < Ω_B). As mentioned in the Introduction, only resonances with j = 1 and j = 2 are able to excite the orbital eccentricities of particles colliding inelastically. Following Paper I, we define the quantity $$\overline{a}=a+a_0\left(\frac{m-j}{j}\right)e^2$$(2)

as the modified semimajor axis. It is a local expression of the Jacobi constant for a given m/(m − j) SOR, where a₀ is the semi-major axis of the circular orbit at exact resonance. The modified semimajor axis is essentially a measure of the semimajor axis a for e ≪ 1. The advantage of a over a is that it is on the average constant near a resonance, so that particles trapped into this resonance will tend to move vertically in the diagrams showing the eccentricity plotted versus $\overline{a}$.

The results given in this paper are applicable to any ring around any body with a given mass anomaly µ or triaxial shape (and thus with given reference radius R_ref, elongation ϵ_elon and oblateness f ). Unless otherwise explicitly stated, times will be expressed in terms of the number of rotations of the central body (corresponding to 2π time units) and the lengths (including the particle radius R) will be expressed in units of the corotation orbit radius a_cor. The mass anomaly, if present, is located at the distance R_ref. To relate these quantities to more physical cases, we can consider for instance the case of Chariklo, using the parameters given in Table 1. The cases of Haumea and Quaoar can also be considered, using the values provided in the Table 2 of Paper I.

Notes. Numerical values are the same as in Paper I, from Leiva et al. (2017); Morgado et al. (2021).

Examples of test-particle responses to various outer SORs are displayed in Fig. 1, in terms of the maximum eccentricity e_max achieved by particles initially on circular orbits. Orbits around a spherical central body with a mass anomaly µ = 10⁻⁴ (upper frame) and an ellipsoidal body with elongation ϵ_elon = 10⁻² and oblateness f = 0 (lower frame) were integrated for 10 000 revolutions of the central body.

In case of mass anomaly, the perturbation potential includes a full range of m-components, leading to a strong response at first-order resonances where n/Ω_B = 1/2, 2/3, 3/4…, corresponding to m = −1, −2, −3… with j = 1. Even with a mass anomaly as small as µ = 10⁻⁴, test particles reach orbital eccentricities as high as e_max = 0.01–0.05. Meanwhile, due to the π-symmetry of the ellipsoidal potential, only even m components are allowed. Consequently, the only first-order resonances present are those corresponding to n/Ω_B = 2/3, 4/5, …. Thus, the response at n/Ω_B = 1/2 for an elongated body, visible in Fig. 1, actually corresponds to a second-order resonance with m = −2, j = 2, and is thus noted 2/4. Also shown in this figure are the analytical responses calculated in Paper I, showing a good agreement with our numerical integrations.

In this paper, we focus on the 1/3 SOR, near which the main rings of Chariklo, Haumea, and Quaoar are observed (Sicardy et al. 2025a), noting that Quaoar’s fainter ring is close to another SOR with n/Ω_B ≈ 5/7 (Pereira et al. 2023). In the case of a mass anomaly, the 1/3 SOR corresponds to a second-order resonance with m = −1, j = 2. Conversely, for an ellipsoidal body and its associated π-symmetry, the 1/3 SOR corresponds to a fourth-order SOR with m = −2, j = 4, thus noted 2/6. Concerning the 5/7 SOR, it can be created only by a mass anomaly, thus corresponding to a second-order SOR with m = −5, j = 2.

With a mass anomaly µ = 10⁻⁴, no visible signature is seen in Fig. 1 at the 1/3 SOR location. Increasing µ to 10⁻² (inset in the upper panel of Fig. 1), a clear response to the 1/3 SOR is seen, in agreement with the theoretical expectation. On the other hand, for an elongated body, no response is noticeable near the 2/6 SOR, even with an elongation as high as ϵ_elon = 0.2, corresponding to the shape of Chariklo. This absence of response is as expected for a fourth-order resonance: as shown in Fig. 1 of Paper I based on the topology of phase portraits, only first-and second-order resonances excite eccentricities of particles initially on circular orbits.

Fig. 1 Comparison between the numerical and theoretical responses of test particles to various SORs. Numerical integrations followed the motion of 10 000 test particles initially distributed on circular orbits, during 10 000 revolutions of the central body. The maximum eccentricities emax reached by these particles are plotted in black as a function of ΩB/n = (a/acor)3/2, and are compared with the analytical estimates of Paper I (red or green curves). Upper panel: case of a spherical body with a mass anomaly µ = 10−4. Lower panel: case of a homogeneous triaxial ellipsoid with elongation ϵelon = 0.01 and oblateness f = 0. In the case of mass anomaly the strongest responses are at the outer first-order Lindblad resonances corresponding to commensurabilities n/ΩB = m/(m − 1), with m = −1, −2, −3… The inset in the upper panel is a zoomed-in image of the ΩB/n = 3 region, using a 100 times larger mass anomaly of µ = 0.01. The response to the second-order 1/3 resonance is now visible, and is compared to the analytical estimate in green. In the case of an elongated body, the π-symmetry of the perturbation imposes even values m = −2, −4, −6… for the first-order resonances. In this case the strongest second-order resonance (m = −2 and j = 2) is also visible, with the analytical response plotted in green. The inset in the lower panel is a zoomed-in image of the ΩB/n = 3 region. The fourth-order 2/6 resonance has no signature even when using a large value ϵelon = 0.20.

3.2 Scaling of first- and second-order resonances

From here on we concentrate on first- and second-order SOR resonances with a mass anomaly. We checked that our simulations correctly reproduce the expected scaling of the particle responses against µ. Figure 2 illustrates the responses of test particles near the first-order 2/3 and second-order 1/3 SORs for various µ’s, compared to the values of e_max displayed in Figs. 3 and 4 of Paper I. We denote by e_peak the largest possible value of e_max near a given resonance and W_res is the width the resonance given in table 1 of Paper I. For first-order resonance, W_res is close to the full width at half maximum (FWHM) of the e_max distribution, while for second-order resonance, it corresponds to the interval where e_max is nonzero. In the particular case of Chariklo, and following the methodology presented in Paper I to derive the resonance strengths, we obtain $$\begin{array}{l}{e}_{\text{peak}}\approx 0.93{\mu }^{1/3},W_{\text{res}}\approx 2.52{\mu }^{2/3}(\text{first}-\text{order}2/3\text{SOR}),\hfill \\ e_{\text{peak}}\approx 0.41{\mu }^{1/2},W_{\text{res}}\approx 0.26\mu (\text{second}-\text{order}1/3\text{SOR}).\hfill \end{array}$$(3)

The numerical factors are specific to the Chariklo case, while the µ-scaling depends only on the resonance order.

The left column of Fig. 2 shows that the general shape of the e_max distribution, the numerical values of e_peak and W_res and the µ-scaling are correctly reproduced in our integrations¹. The right column shows the time T_res required for a particle initially on a circular orbit to reach the maximum value e_max. Numerical integrations imply that at the resonance $$\begin{array}{ll}{T}_{res}\approx 0.25{\mu }^{-2/3}\hfill & (\text{first}-\text{order}2/3\text{SOR}),\hfill \\ T_{res}\approx 40{\mu }^{-1}\hfill & (\text{second}-\text{order}1/3\text{SOR}).\hfill \end{array}$$(4)

The above µ scalings confirm the scaling of timescales obtained in Paper I. Figure 3 shows a good agreement between numerical integrations and the analytical estimates for the linear growth rate in the case of first-order resonance (T_lin ∝ µ^−2/3; Eq. (32) in Paper I), and the exponential growth rate in the second-order resonance (T_exp ∝ µ⁻¹; Eq. (35) in Paper I).

For the typical threshold value µ = 10⁻³ that we later estimate for confining ring material (which corresponds to a ∼10 km mountain on Chariklo), we obtain at the first-order 2/3 resonance (orbital radius a_2/3 = 1.31) the values e_peak = 0.093, W_res = 0.0252 and T_res ≈ 600. For the second-order 1/3 resonance at a_1/3 = 2.08, we obtain e_peak = 0.013, W_res = 0.00026 and T_res ≈ 40 000. The long timescale for the excitation of eccentricities at the 1/3 resonance fully explains the absence of any signature in Fig. 1. For µ = 10⁻⁴, even if the maximum eccentricities would be about 0.004, the growth of resonance eccentricities would have required a 40-fold longer length of integration than depicted in Fig. 1. Moreover, the resonance width (≈3 × 10⁻⁵) would have been much too small to be discernible in this figure.

Because of the large strength and radial extent of the first-order resonances, it is worth estimating what their effect is at the distance of the second-order 1/3 SOR we are interested in. This is illustrated in Fig. 4, comparing the theoretical amplitudes e_max in the vicinity of 1/3 SOR at a/a_cor ≈ 2.08. For µ = 0.1, the eccentricities associated with the 2/3 SOR at the same distance are ≈0.04, or nearly 30% of the e_peak due to the 1/3 SOR. Moreover, the maximum eccentricities associated with first-order resonances are reached very rapidly, compared to the slow growth of second-order perturbations. Not surprisingly, collisional simulations performed for 1/3 resonance with large values of µ also show a clear m = 2 undulation in their initial evolution phase (see, e.g., the T = 500 frames in Fig. 8). However, the m = 2 undulation has no effect on the 1/3 resonance confinement observed in collisional simulations². In addition, although the values of e_peak associated with the first-order resonances decrease more slowly than that of the second-order 1/3 SOR (the ratio of the e_peak’s scale as µ^−1/6; see Eq. (3)), the local ratio of eccentricity amplitudes at the 1/3 location drops as µ^1/2. This stems from the fact that the widths of the first-order SOR’s shrink proportionally to µ^2/3.

Fig. 2 Responses of test particles to the first-order 2/3 and second-order 1/3 SORs. Left column: symbols represent the maximum eccentricities emax reached by particles initially on circular orbits for three values of the mass anomaly µ, near the 2/3 SOR (upper panel) and the 1/3 SOR (lower panel). The gray dashed curves are the analytical estimates of emax displayed in Figs. 3 and 4 of Paper I. The points are plotted against the normalized distance to the resonance, Δa/Wres = (ā − a0)/Wres, where a0 is the radius of the circular orbit at exact resonance, ā is the modified semimajor axis (Eq. (2)), and Wres is the width of the resonance (Eq. (3)). Right column: same, but with the timescales of Tres necessary to reach the maximum eccentricities emax. This figure shows that our numerical integrations reproduce satisfactorily the calculated distribution of emax, as well as the µ-scaling expected from Eqs. (3) and (4).

Fig. 3 Comparison of eccentricity evolution at the first-order 2/3 and second-order 1/3 SORs. The black curves follow the evolution of test particles released near exact resonance, while the red dashed lines display the theoretical prediction, i.e., a linear growth rate for the first-order SOR and an exponential growth rate for the second-order SOR.

Fig. 4 Comparison of the maximum eccentricity emax reached by particles at the second-order 1/3 SOR (blue) compared with the first-order 2/3 (red) and 1/2 (green) SORs. Three values µ = 0.1, 0.01, 0.001 are compared. For µ = 0.1, the theoretical emax associated with the 2/3 and 1/2 resonances at the location of 1/3 resonance are 30% and 8% of that due to 1/3 resonance, respectively. For other µ values the ratios at the 1/3 SOR location scale proportionally to µ1/2.

Fig. 5 Three cases showing the combined effects of collisions and the 1/3 SOR on the ring confinement. The left frames show the time evolution of the vertical angular momentum (Lz) distribution of the particles around the 1/3 SOR at Lz = 1.44. The magenta lines are the average value of Lz and the red lines show the RMS of the particle eccentricities; the full y-range of the frame corresponds to e = 0.1. The insets are polar plots of the particle positions, shown again in the right column in cartesian coordinates. For better viewing, in the cartesian projections the width of the ring around its center position at each azimuth and the deviations of this center position from the overall mean distance have both been expanded by a factor of five. We compare three cases: (i) A ring of collisionless test particles perturbed by a µ = 0.1 mass anomaly on an otherwise spherical central body (upper row); (ii) Colliding particles around a spherical central body without perturbation (µ = 0, middle row); and (iii) A ring of colliding particles with µ = 0.1 (lower row). All simulations use the same initial conditions with 30 000 particles placed initially in an annulus r = 2.02–2.14 straddling the 1/3 SOR at semimajor axis a = 2.08. In the collisional simulations the particle radius R = 10−3 corresponds to about 200 m for Chariklo’s ring particles and yields an initial geometric optical depth τ0 = 0.06. In the case of a mass anomaly, the perturbation is turned on linearly during the first 50 revolutions of the central body.

4 Results from collisional simulations

4.1 Illustrative example

We first demonstrate the crucial role of physical impacts on the response to resonances. This is illustrated in Fig. 5 using a large mass anomaly µ = 0.1. The particles with radii R = 10⁻³ are initially placed on a wide annulus that straddles the 1/3 resonance semimajor axis a_1/3 ≈ 2.08. Recalling that R is measured in units of the corotation radius a_cor, this corresponds from Table 1 to a physical radius of ≈200 m for Chariklo’s ring particles.

Figure 5 shows density plots of the angular momentum distribution L_z evolving with time, together with snapshots of particle positions at the end of the simulation, both in polar and cartesian systems. In the case of noncolliding test particles (upper panels of Fig. 5), the time evolution of the system near the resonance location exhibits a growth of eccentricities, with maximum amplitudes and timescales behaving as illustrated in Fig. 2 (see the peak in the eccentricity RMS at T ≈ 400). Most notably, there is no secular change of particle mean distances, so that there is very little change in the L_z distribution³. When plotting instantaneous particle positions, a gap develops at the resonance due to excited eccentricities, since the particles spend most of their time near the extremes of their orbital epicycles.

In the case of unperturbed but mutually colliding particles (middle panels of Fig. 5), the ring evolution shows a rapid collisional damping of inclinations and eccentricities, with a timescale of a few tens of impacts per particle, and a gradual viscous spreading taking place on longer timescales determined by the ring viscosity ν. This is seen as a widening of the ring and the L_z distribution in Fig. 5.

Finally, when both impacts and resonant perturbations are included (lower panels of Fig. 5), the behavior of the ring is drastically different. The simulation leads, after initial viscous spreading, to the excitation of eccentricities by the 1/3 SOR in parallel with the accumulation of particles around the resonance. At the same time the mean L_z of the ring jumps and the system settles to a narrow ringlet just outside the resonance. In the simulation of Fig. 5 this accumulation and confinement takes place at T ≈ 2500–3000. We discuss in Sect. 4.3 the conditions under which such behavior occurs. However, on longer timescales, not covered in Fig. 5, the ringlet slowly leaks away particles at its outer edge. We return to this process in Sect. 7.

4.2 The ring confinement at the 1/3 resonance

We now describe the various phases of the ring confinement process, using as an example a simulation with more realistic parameter values. The particles have a radius R = 1.25 × 10⁻⁴ (about 25 m for Chariklo’s ring particles) and are perturbed by a mass anomaly µ = 3 × 10⁻³. Due to smaller µ the timescale of the evolution is much longer than in the simulation of Fig. 5. For example, it takes about 150 000 central body revolutions before the resonance accumulation starts. The total span of the simulation, T ≈ 257 000 corresponds to about 206 years in case of Chariklo.

The nine upper panels of Fig. 6 display three phases of the ring evolution. Phase I corresponds to the damping of the initial eccentricities followed by a radial viscous spreading of the ring material; during Phase II, the 1/3 resonance excites the orbital eccentricities of the particles, forming a ringlet which gathers material from the background material. The eccentricity then reaches the peak value e_peak (Eq. (3)). A similar surge of eccentricities and its ensuing damping were visible in a simulation with µ = 0.1 (see the red line in the lower left panel of Fig. 5). This damping leads to a quasi-steady-state that constitutes the Phase III. A ringlet is now confined, in which the damping of eccentricities by collisions is balanced by the resonance excitation. The Fig. 7 synthesizes the three phases, in which 1600 snapshots of the system have been stacked in the (ā, e)-space.

It summarizes the history of the ring material, from the initial viscous spreading phase to the steady-state situation.

The three lower panels of Fig. 6 show the ring in polar coordinates at three selected times. At T = 163 120, a streamline excited by the 1/3 SOR has appeared, exhibiting a tendency of self-crossing as expected from the Fig. 6 of Paper I. At T = 169 840, this streamline is gathering particles from the unexcited background ring material. This gathering can be understood by the fact that the background particles near the point A move slower than the particles in the streamline. Consequently, they receive angular momentum during collisions, so that their semimajor axes increase. Conversely, the semimajor axes of background particles near point B decrease during collisions. This globally leads to the confinement of the particles near the resonance radius. At T = 257 280, all the background particles have been gathered into a nonintersecting streamline dominated by a m = 1 azimuthal mode.

The results of a complementary run using a larger mass anomaly µ = 0.1 are provided in Fig. C.2. It shows better the initial self-intersecting streamline forced by the 1/3 SOR, and its further transformation into a nonintersecting streamline dominated by a m = 1 azimuthal mode. In this case, two other ringlets form inside and outside the central ringlet, the outer ringlet being eventually swallowed by the central ringlet during the run.

4.3 Scaling to physical systems

We now address the applicability of simulation results to real systems. Firstly, we study the role of impact frequency in the transition between a test-particle system and a collisional ring. Then, we provide the criterion for which the resonant confinement is expected to win over the viscous spreading due to collisions. These are important factors since a fully realistic simulation of an azimuthally complete dense ring including collisions between presumably meter-sized particles would imply an unmanageable number of such particles⁴. In practice, a smaller optical depth and larger than real particles must be employed in simulations, calling for a scaling between the particle size R, mass anomaly µ, and optical depth τ in simulated and real physical systems.

In the nongravitating simulations with constant rebound coefficient ϵ_n, the collisional steady state of an unperturbed simulation system is determined by three parameters: the particle radius R, the dynamical optical depth τ, and the coefficient of restitution ϵ_n. The influence of perturbation depends on the mass anomaly µ and on the particular resonance(s) acting on the ring. For a nonaxisymmetric central body, the strength of the perturbation also depends on the oblateness f and elongation ϵ_elon. To a lesser extent, the initial radial width of ring W₀ also matters, as it determines the maximum number of particles that can be perturbed by the resonance. Ideally, W₀ should be large compared to both R and the resonance width W_res. On the other hand, the initial values of the vertical ring thickness and velocity dispersion are not critical, since the collisional damping of eccentricities and inclinations rapidly leads in the unperturbed case to a steady state with velocity dispersion c ∼ Rn, the precise proportionality factor depending on ϵ_n (Salo et al. 2018). The importance of τ follows from the fact that in an unperturbed 3D ring the impact frequency w_c is proportional to τn. The basic formula for viscosity is ν ≈ w_cλ², where the radial mean free path λ ≈ c/n at low τ rings. This implies $$v=k_{visc}\tau R^{2}n,$$(5)

where k_visc is a numerical factor on the order of unity. With our standard value ϵ_n = 0.1, we have k_visc ≈ 3.5.

When scaling the simulations to physical systems, two questions at least arise: (1) What the minimum frequency of impacts w_c is that makes the ring behave as a collisional ring, in contrast to a mere ensemble of independent test particles. (2) What the parameter regime is in which the timescale for the resonance build-up is shorter than the viscous spreading time due to collisions. These two questions are addressed in the next two subsections.

Fig. 6 Simulation with 40 000 ring particles of radius R = 1.25 × 10−4 (corresponding to 25 m for Chariklo’s ring particles) perturbed by a mass anomaly µ = 3 × 10−3. The time T is given in units of Chariklo’s rotation period (about 7 h, Table 1), so that the maximum time shown here (T = 256 280) corresponds to about 206 years. Upper nine panels: density maps of the particles in the modified semimajor-eccentricity (ā, e) space (see Eq. (2)). The gray spiky curve is the value of emax shown in the lower left panel of Fig. 2. Three phases of the ring evolution are displayed. Phase I corresponds to a rapid damping of the initial eccentricities accompanied by a radial spreading caused by collisions; During Phase II, the 1/3 SOR excites the orbital eccentricities of the particles up to the predicted peak value epeak (Eq. (3)), while confining concomitantly the material near the resonance (ā ≈ 2.08) over a timescale Tres (Eq. (4)). Finally, during Phase III, the eccentricities damp due to dissipative collisions and a quasi-steady-state is reached, where the eccentricity damping is balanced by the resonance excitation. Lower three panels. Density maps of the particles in polar coordinates (L, r) space at three selected times shown in the upper panels. The white dotted lines indicate the location of the 1/3 SOR. At T = 163 120 a streamline forced by the 1/3 SOR has appeared, with a tendency of self-crossing around L = 260 deg. It is gathering material from the background unexcited particles. At T = 169 840, the accumulation process is going on, where a well-formed ringlet with various azimuthal modes collects more background material. At T = 257 280, all the ring material has been accumulated onto the ringlet, which is now dominated by a m = 1 azimuthal mode, with the presence of two kinks. The colors in the plots indicate the density of particles (in arbitrary units, density increasing from blue to red). Five movies generated from this simulation are available online. Movie1 illustrates the evolution in the (ā, e) space during Phase I, between T=0 and T=99 840 (corresponding to 0 and 79.77 years for Chariklo, respectively), using snapshots stored every 240 Chariklo rotations (about 70 days). Movie2 and Movie3 show the same for Phase II (between T=100 000 and T=199 840) and Phase III (between T=200 000 and T=257 280), respectively. The third movie also plots the position of an individual particle (white dot), showing its capture into the ring and its subsequent motion inside the cloud of points. Movie4 and Movie5 display the same snapshots in the longitude-radius space corotating with Chariklo, for Phases II and III, respectively.

Fig. 7 Upper panel: overview of the three phases shown in the upper panels of Fig. 6. The plot shows a stack of 1600 snapshots of the system from T = 120 000 to T = 200 000 (96 to 160 years respectively in the Chariklo case) with time steps ∆T = 50 Chariklo’s revolutions (about 15 days). The three phases I, II, and III described in the text are indicated by the arrows. Lower panel: phase III quasi-stationary stage obtained at the end of our run (we note the change in radial scale compared with the upper panel). A total of 780 snapshots taken from T = 256 500 to T = 257 280 (204.94 to 205.57 years) with time steps ∆T = 1 Chariklo’s revolution (about 7 h) have been stacked. A quasi steady state is reached, where the eccentricity damping due to collisions is balanced by the resonance excitation. On the long term (not reached here), a leakage of particles toward larger semimajor axes would be observed, as illustrated in Fig. 17, where a larger mass anomaly µ = 0.03 is used.

4.3.1 Influence of impact frequency

For nongravitating rings, the impact frequency w_c ≈ 3nτ, corresponding to about 20τ impacts/particle/ring orbital period. A condition sometimes quoted for a collisional ring is to have at least one impact/particle/ring orbital period. This condition, w_c > n/(2π), would correspond to τ ≳ 0.05. However, since the resonance timescales at 1/3 SOR are much longer than the ring orbital period (see Eq. (4)), a much smaller impact frequency can in practice be expected to modify the test-particle evolution.

We explore the transition from a noncollisional to a collisional ring by conducting simulations with successively larger values of initial optical depths τ₀. Fig. 8 shows snapshots of the ring near the 1/3 SOR, using µ = 0.1 and τ₀ = 0.00025–0.06. The corresponding unperturbed impact frequency w_c increases from about 0.01 to about one impact/particle/orbit, indicating on average one-hundred to one ring orbital periods between impacts (marked as T_c in the second column). For the second-order 1/3 SOR, T_res ∼ 400 central body rotation periods (Eq. (4)), i.e., about 130 ring orbital periods. This timescale refers to the time it takes to reach the maximum eccentricity starting from circular orbits; the e-folding time T_exp of eccentricity growth is about ten times shorter.

Thus, for the smallest τ₀ explored in Fig. 8, the impact timescale T_c is roughly equal to the resonance timescale T_res. Consequently, the effect of impacts is weak: most of the particles initially close to the resonance are able to cross radially the ring without colliding. A gap similar to that in the upper panel of Fig. 5 (collisionless test particles) opens at the resonance. This gap gets slowly filled with time, as inter-particle impacts lead to increased eccentricities throughout the ring. No sign of particle accumulation in the resonance region is seen during the time span of the simulation. In contrast, the ring appears very diffuse.

When τ increases to τ = 0.001 and τ = 0.002 (corresponding to the L_z maps colored green in the last column of Fig. 8), both the opening of the initial gap and the rapid growth of eccentricities take place. The system eventually forms a nearly circular ringlet just outside of the 1/3 resonance, surrounded by a population of hot ring particles not trapped by the resonance. We note a weak undulation of the ringlet with azimuthal number m = 2, caused by the distant 2/3 resonance which is still quite strong at the 1/3 SOR distance (see Fig. 4).

For τ ≥ 0.004, corresponding to T_c ≲ 10, the behavior is reminiscent of that of Fig. 5: the whole system goes through the initial accumulation, resonance excitation and confinement phases (L_z maps colored in red in the last column). The eventual leakage of particles outside the resonance radius is obvious for the τ ≳ 0.01 runs, the dispersal of the ringlet getting faster with larger τ.

To conclude, all the simulations of Fig. 5 where T_c ≲ 0.1T_res ≈ T_exp show an evolution similar to what is displayed in Figs. 6 and 7. Taking into account the scaling between T_res and µ leads to the following empirical condition for a ring near the 1/3 SOR to be collisional (instead of an ensemble of test particles): $$\tau \gtrsim (0.01-0.04)\mu .$$(6)

Here the larger limit corresponds to the formation of confined eccentric ringlet while the lower limit corresponds to T_c where the first signs of collective behavior appear.

4.3.2 Competition between resonance accumulation and viscous spreading

The eccentricity growth forced by a resonance takes place on timescales of ∼T_res. The growing epicyclic excursions induce collisions between particles in and close to the resonance zone. Due to the dissipative nature of the impacts, the particles originating further away from the resonance tend to remain closer to resonance after impacts as their eccentricities are damped. In order to accumulate particles at the resonance, the viscous removal of particles from resonance must not be too fast. To get a condition for net accumulation, we follow the heuristic arguments presented in Franklin et al. (1980).

The displacements of particles diffusing away from the resonance zone are given by $$\text{Δ}{r}_{visc}\approx \sqrt{vt},$$(7)

while the epicyclic excursions due to growing eccentricities transport particles to the resonance zone from distances $$\text{Δ}{r}_{res}=e_{res}(t)a_0$$(8)

These particles preferentially collide with particles near the resonance zone and end up with mean distances close to the resonance. In order to obtain a net accumulation at the resonance, we must have ∆r_visc ≲ ∆r_res during the whole time range t ≲ T_res. This implies $\sqrt{v{T}_{\text{res}}}\lesssim e_{\max }{a}_0$, leading to the condition $$v\lesssim {\left(e_{\max }{a}_0\right)}^2/T_{\text{res}}=k_{\text{res}}{\left(e_{\text{peak}}{a}_0\right)}^2/T_{\text{res}}.$$(9)

Here, the prefactor k_res << 1 takes into account the fact that the average e_max amplitudes near the resonance are less than the peak value e_peak, and also that the timescales to reach e_max are generally longer than T_res (see the right panel of Fig. 2).

We now plug in Eq. (5) for the viscosity, while using the formulae (3) and (4) for e_peak and T_res at the 1/3 SOR. This provides a simple condition for the parameter regime in which the resonance confinement should be possible, $$k{\mu }^2\gtrsim \tau R^2,$$(10)

where k ∼ 0.05k_res/k_visc. We note that this estimate suggests that the threshold µ for confinement is expected to scale as ${\mu }_{\text{thresh}}\propto \sqrt{v}$. In order to check the τ, R and µ dependence of the accumulation threshold and to estimate the value of k, we turn to numerical simulations.

Fig. 8 Transition from noncollisional to collisional rings. This is explored by simulations with various initial optical depths τ0 from 0.00025 to 0.06, all using a mass anomaly µ = 0.1. Snapshots of the particle distributions in polar coordinates are shown at various times, using the radial range 1.88–2.28 centered around the resonant semimajor axis a1/3 ≈ 2.08. The label Tc in the second column indicates the average time interval between impacts, in units of the particle orbital periods. The rightmost column shows the time evolution of the vertical angular momentum Lz distribution with time up to 30 000 central body revolutions. The lowermost simulation (τ0 = 0.06) is the same as the one shown in the last row of Fig. 5, where it is displayed up to T = 5000.

4.3.3 Empirical criterion for ringlet formation

Figure 9 displays a grid of 1/3 resonance simulations all with the same initial optical depth τ₀ = 0.015, using values of µ between 0.003 and 0.1, and particle radii from 0.000125 to 0.004 (25 m to 800 m when scaled to Chariklo’s ring system). This figure illustrates the competition between resonance accumulation and viscous dispersal, covering for each µ a range of runs both leading and not leading to confinement. In particular, this survey confirms that for a fixed τ, and as predicted by Eq. (10), the minimum µ required for the resonance confinement increases linearly with R. Thus, a much smaller size of mass anomaly can be expected to lead to resonance confinement when the particle size and ring viscosity approach more realistic small values. The black unit-slope line in the right panel delineates the boundary between the two regimes. Except for the largest value µ = 0.1, the boundary is reasonably well approximated by this unit-slope line, in agreement with Eq. (10). The boundary provides an estimated value k ∼ 4 × 10⁻⁵, so that Eq. (10) for 1/3 SOR can be re-expressed as $$\mu \gtrsim 160\sqrt{\tau }R,$$(11)

where we recall that the particle radius R is measured in units of the corotation radius a_cor. For Chariklo system this yields $$\mu \gtrsim {10}^{-3}\sqrt{\tau }{R}_{\text{phys}},$$(12)

where R_phys is the particle radius in meters. As found in test-particle integrations, for µ = 0.1 the eccentricity amplitudes are about 30% weaker than predicted by the theoretical formula for e_peak. It also takes longer to reach the maximum eccentricity than predicted by the fitted T_res based on smaller µ’s. Together these effects weaken the perturbation by roughly a factor of 2, explaining the reduced threshold value of R for µ = 0.1 compared to the black line limit.

Additional tests were performed to check the τ dependence of the accumulation boundary: according to Eq. (11) accumulation requires τ ≲ τ_lim ≈ 4 × 10⁻⁵(µ/R)². In case of Fig. 8 with µ = 0.1, R = 10⁻³, the implied τ_lim ≈ 0.4, consistent with the fact that accumulation was seen even in the case of the largest τ₀ = 0.06. Appendix C.3 reports similar surveys using R = 10⁻³ with µ = 0.05 and 0.03, in which case the predicted τ_lim ≈ 0.1 and 0.035, respectively. Simulations covering a range of τ₀’s confirm the expected trend of reduced τ_lim ∝ µ², though the limiting values observed in simulations are roughly 30% smaller than estimated, again most likely due to the large µ.

In addition to the value τ = 0.015 used in the simulations, an extrapolated curve for τ = 1 is shown as a red line in the right panel of Fig. 9. It assumes that the linear relation ν ∝ τ holds all values of τ, which is not strictly true when τ starts to approach unity (see Salo et al. 2018). In any case, the red curve indicates that confinement should take place in dense Chariklo type rings provided that mass anomaly exceeds about 0.001, assuming typical one-meter ring particles. We note that this estimate includes only collisions (see Sect. 8 for a refined estimate in case gravitational viscosity is taken into account).

Fig. 9 Transition between accumulation and dispersal, a test of Eq. (10). Left panel: same as the left column of Fig. 5, but exploring a grid of simulations near the 1/3 SOR. The particle radii increase from left to right in the range R = 1.25 × 10−4 − 4 × 10−3, corresponding to R = 25-0800 m in the case of Chariklo’s rings, while the mass anomaly increases upward from 0.003 to 0.1. The number of particles and the initial width of the rings were chosen in order to provide the same initial optical depth τ0 = 0.015 for all simulations. The vertical axes span values of Lz between 1.35 and 1.55 for µ = 0.1 and 0.03, and between 1.4 and 1.5 for µ = 0.01 and 0.003. The small insets covering the radial range 1.9–2.4 show snapshots of the ring in polar coordinates at the end of the run; the label indicates the duration of the run in Chariklo rotation periods. Increasing the particle size (and thus the viscosity) for a given perturbation strength prevents the resonance accumulation. Right panel: filled and open symbols respectively distinguish between simulations leading to resonance confinement and dispersal. The black line indicates the accumulation threshold for the simulations displayed in the left panel, following Eq. (10); the dashed portion of the line indicates the region where the linear relation fails since the resonance excitation at large µ ≳ 0.03 becomes somewhat weaker than predicted by the analytical formulas (see footnote 1). The shaded region bounded by the red line extrapolates the accumulation region to a denser ring with τ = 1.

5 Normal modes

In the case of forced Lindblad resonances ( j = 1 in Eq. (1)), the response of collisional rings is relatively straightforward: dissipative impacts force the particles to follow m-lobed non-intersecting streamlines that are close to the resonant periodic orbits (see Fig. C.4). In the more general case of a m/(m − j) SOR, the periodic orbits have |m|( j − 1) self-intersecting points (Sicardy 2020). Thus, in the case of the 1/3 SOR (m = −1, j = 2), there would be one intersecting point (see Fig. 6 in Paper I). In practice, collisions prevent this crossing of orbits from persisting (see the T = 3000–8000 frames of Fig. C.2). An interesting question is the resulting flow configuration the system settles to.

According to our simulations, the typical overall outcome is an excitation of a dominant |m| = 1 mode. For example, among the simulations shown in Fig. 9, only in one case a |m| = 2 ringlet was initially formed. Even in this case, when the simulation was continued further, the |m| = 1 mode eventually took over (Fig. 10). However, in addition to the |m| = 1 mode, additional Fourier modes are always present (see, e.g., the T = 30 000 frames in Fig. 8). The most striking example is the µ = 0.003 simulation of Fig. 6 which settles to a complicated azimuthal profile with two prominent kinks. As demonstrated below, such a ring response can be interpreted as the 1/3 SOR excitation being transferred to a superposition of several free outer Lindblad modes. We note that in this simulation the smallest particle size and perturbation are used, whereupon it can be assumed to mimic closest the possible behavior of real rings. In what follows we analyze in more detail the ringlet formed in this simulation.

Figure 11 shows the radius versus longitude profile of the µ = 0.003 simulation, at T = 250 000. In the uppermost frame the fitted |m| = 1 Fourier-component is superposed on the profile, while the lower frames show the residual profile after subtracting consecutive components |m| = 1, 2, 3. Clearly, both |m| = 2 and |m| = 3 modes are significant in addition to the |m| = 1 mode, whereas subtracting the |m| = 4 mode would not have much effect on the residual profile. In order to analyze the propagation of the modes, a shape model, $$r_m(L,t)=A_m\cos \left[m\left(L-{\text{Ω}}_{m}t\right)+{\varphi }_m\right]$$(13)

was fitted to each Fourier component over the time range ∆T = 200. Here r_m(L, t) is the radius of the streamline versus the true longitude L at time t, A_m is the amplitude of the mode and ϕ_m its phase, Ω_m is the pattern speed and ω_m = |m|Ω_m is the frequency. The results of this fit are collected to Fig. 12. The peaks of the periodogram for |m| = 2, 3, … all fall on a linear relation ω_m/Ω_B = (1 + |m|)/3, while for |m| = 1, ω_m ≈ 0.

This is a quite remarkable result as it shows that the system’s response is a superposition of free outer (since Ω_m > n) Lindblad modes corresponding to the condition κ = m(n − Ω_m), with m < −1. In such a mode, a particle executes exactly |m| radial excursions while completing one revolution in a frame rotating with the pattern speed Ω_m. The m = −1 mode corresponds to ${\text{Ω}}_m=\dot{\varpi }$, i.e., the locked precession mode of an essentially Keplerian ellipse. The other modes with m ≠ 1 correspond to $$\frac{{\text{Ω}}_m}{n}\approx \frac{m-1}{m},$$(14) $${\omega }_m\approx |m-1|n$$(15)

where the approximation stems from the fact that $\dot{\varpi }\ll n$. Taking into account that n = Ω_B/3 and m < 0, this is exactly the same relation as observed in the simulation.

Figure 12 indicates that the Fourier-modes with |m| ≥ 4 have small amplitudes, compared to |m| = 1, 2, 3. However, their amplitudes decay quite slowly, and most importantly, do not stem from noise but convey a true signal, corresponding to the kink, best seen in the lowermost residual profile of Fig. 11.

This correspondence is illustrated in Fig. 13, comparing the simulated ringlet to a toy model where |m| = 4–20 Fourier modes have been superposed. We assume here that the amplitudes drop as A_m ∼ |m|^−3/2 and that the phases are the same for all modes, roughly corresponding to the properties of the simulated modes (the π phase shift between ϕ_m of even and odd modes seen in last frame of 12 only moves the phase of the kink feature, not affecting its shape or propagation; on the other hand, a nonalignment of ϕ_m’s would destroy the kink feature). The result is a propagating kink, moving with the orbital speed of the ring, closely resembling what is seen in the simulation.

Since ω_m is a multiple of n, each mode is invariant through a time translation of T_orb, the orbital period of the particle. In the present case (a ring near the 1/3 SOR), this means that the ring recovers its initial shape after three rotations of the central body. We have made use of this in Fig. 14, illustrating the ring evolution over one full ring period (three central body rotations), using for each time stacks of 20 particle snapshots separated by ∆T = 3. In addition to the azimuthal profiles, two sets of cartesian projections are shown: in the middle frames the deviations of the ringlet mean distance are exaggerated, while in the right the same is done with width variations. As seen, the shape and width of the ringlet varies in a complicated cycle, however repeating regularly over time.

Also indicated in Fig. 14, are the regions where local shear reversal is taking place, defined as the locations where the non-diagonal component of the velocity dispersion tensor, T_rt =< ∆v_r∆v_t >, has negative values. Here ∆v_r and ∆v_t are the particle velocities with respect to the mean flow at their position, and brackets indicate a mean over a local region. In order to measure T_rt we used a method quite similar to that in Hänninen & Salo (1992). We divided the particles into twenty streamlines according to their Jacobi energies, and calculated the mean radial and tangential flow velocities along each such streamline, to finally obtain the T_rt along the streamline. Stacking of several snapshots was essential to reduce the noise in this process.

In an unperturbed case with Keplerian shear, T_rt > 0 throughout the ring, indicating that collisions on average transfer angular momentum outward; for example, impacts by inner particles approaching with positive relative v_r have on average positive v_t with respect to outer particle, providing a positive impulse. This situation corresponds to viscous spreading of unperturbed rings. Negative values on the other hand associate with inward transport (Borderies et al. 1983). The flux of angular momentum (per unit length of streamline) relates to T_rt by F(a, ϕ) = Σ(a, ϕ)aT_rt(a, ϕ), and the flux integrated over the streamline gives the angular momentum luminosity L_H(a). The situation where L_H(a) < 0 corresponds to the maintenance of sharp edges due to flux reversal (Borderies et al. 1983) or single-sided shepherding (Goldreich et al. 1995). The numerical simulations of Hänninen & Salo (1994) confirmed that such a mechanism can maintain sharp edges at first-order Lindblad resonances.

Compared to the Lindblad-resonance case (see Fig. C.4), where the flow pattern, and the regions of positive and negative T_rt are fixed in a frame rotating with the satellite, the current situation is more complicated, due to superposition of several modes making the pattern truly variable in time. Because of this the detailed analysis of the angular momentum transport is left for a later study.

Fig. 10 Transition from a |m| = 2 ringlet stage to a |m| = 1 ringlet stage in the µ = 0.03, R = 5 · 10−4 (corresponds to 100 m in the case of Chariklo) simulation of Fig. 9. In all other cases the ringlet shape is dominated by the |m| = 1 response once it has formed. Also in this case, the |m| = 1 response eventually replaces the |m| = 2 shape. During this transition at T ≈ 35 000, the torque exerted on the nearly circular ringlet temporarily vanishes (see the flat portion of the mean Lz curve in the upper frame).

Fig. 11 Fourier fits to the ring azimuthal profile in the µ = 0.003, R = 25 m simulation at T = 250 000 (same simulation as in Fig. 6). The uppermost frame shows r(L) profile, with a |m| = 1 fit superposed. The second frame shows the residual profile after subtracting the |m| = 1 component; the green curve is the |m| = 2 fit to the residual profile. The next frames repeat the procedure for |m| = 3 and |m| = 4.

Fig. 12 Proper Lindblad modes appearing in a ring confined at the 1/3 SOR. The outputs of the run shown in Fig. 6 have been used to detect normal modes in the confined rings between times 250 000 and 250 200, with steps ∆T = 0.05. The Lomb normalized periodogram power spectra of the modes described by Eq. (13) are displayed in the upper six panels as a function of ωm/ΩB for |m| = 1 to 6. Consistent with Eq. (15), the maximum power (red squares) is reached at the frequency ωm = (1 + |m|)ΩB/3. This is illustrated in the lower left panel. The lower middle panel shows the rapid decrease in the amplitude Am of the modes with |m|, while the right panel shows the distribution of phases ϕm.

Fig. 13 Left column: time evolution of the ringlet at the end of the µ = 0.003 simulation. The dominant |m| = 1, 2, 3 modes have been subtracted in order to highlight the kink feature. Right column: toy model for the formation of a kink as a superposition of |m| = 4–20 modes. Each mode has the same phase ϕm while the amplitudes obey Am ∝ 1/|m|3/2. Profiles calculated with Eq. (13) are shown over one full ring period.

Fig. 14 Local shear reversal in 1/3 SOR ringlet. The evolution of the µ = 0.003, R = 25 m simulation is shown over one ring orbital period (or three central body revolutions). The left frames display the r(L) profiles, with radial range 2.05–2.12. The ring patches with negative angular momentum flux Trt are indicated in light blue, while the white bullet indicates the location of a tracer particle. The middle frames show cartesian projections, where the ring deviations from its mean distance (dashed curve) have been exaggerated by a factor of 40 at each azimuth, and the open square marks the location of the mass anomaly. Similarly, in the right frames (again a cartesian projection) the ring width variations have been exaggerated by a factor of 50, using the same color convention as in the left frames. For the construction of all the figures, 20 snapshots of the ring separated by ∆T = 3, 6, …57 have been stacked.

6 Migration of material from inner regions

So far, we have focused our attention to the confinement of ring material initially near the 1/3 SOR. However, the various rings observed so far around small objects were probably formed over a broad range of radii around the central body, and in very different contexts (Sicardy et al. 2025a). For instance, Haumea’s ring may have formed during the spewing of material associated with a spin up phase (Noviello et al. 2022), while Chariklo’s rings may originate from a cometary activity that launched material from its surface (Sicardy et al. 2025a).

In general, the resonances rapidly clear the corotation region, pulling the ring material inside the synchronous orbit down to the surface of the body, while pushing away the material outside the synchronous orbit. Using a toy model with Stokes-like friction, Sicardy et al. (2019) showed that this clearing may occur over decadal scales under the effect of a Chariklo elongation ϵ_elon = 0.2, and in a few centuries if a mass anomaly of µ ∼ 10⁻³ is present.

We have followed the migration of a colliding ring initially placed near the second-order 5/7 SOR. The Fig. 15 shows the time evolution of the particle angular momenta L_z. The crossing of each first- or second-order SOR by the material results in a jump in L_z, superimposed to a general positive drift of L_z, i.e., a positive torque exerted by the mass anomaly on the disk. The Fig 16 shows the evolution of the same system in the (ā, e) space. This figure confirms the expectation that each first-and second-order SOR excites the orbital eccentricities as predicted by Figs. 2 and 7, while collisions damp the eccentricities when the ring material evolves between resonances. For the case µ = 0.003, Fig. 16 shows that the characteristic timescale for the ring migration is some 10⁵ Chariklo’s rotations, corresponding to some centuries, confirming the results obtained by the toy model of Sicardy et al. (2019).

The elongations of bodies like Chariklo, Haumea, or Quaoar cause stronger resonances than a mass anomaly (see Figs. 7–9 of Paper I), and thus repel even more rapidly the ring material, with larger eccentricities. In that context, the weaker second-order 1/3 resonance may correspond to a protected zone where a ring may be confined.

Fig. 15 Time evolution of the vertical angular momentum Lz of particles that cross various resonances. The simulation contains 7500 particles with radius R = 0.001 (corresponding to 200 m in case of Chariklo), perturbed by a mass anomaly µ = 0.003. The ring starts near the 5/7 SOR and is pushed outward by a positive torque which takes it through the 2/3, 3/5, and 1/2 SORs. The red curve is the mean value of Lz. The green curve shows the mean eccentricity of the particles.

7 Effect of an outer satellite

Even though the 1/3 SOR is efficient in confining a ring on short timescales, our simulations show a leakage of material outward on the long term. This is shown in the upper rows of Figs. 17 and 18, which display the evolution in a simulation with µ = 0.03 up to T = 200 000. An eventual dispersal of the ringlet is unavoidable as the continuous torque exerted on the ringlet by the mass anomaly implies a gradual increase in its mean L_z, even after the rapid jump in L_z associated with the resonance passage is over (see Fig. 19).

The drift of the ringlet out of the resonance is not as rapid as the value of dL_z/dt would indicate, since most of the angular momentum is carried out by particles pushed out from the resonance, while the core of the ringlet stays confined. However, even the ringlet core would eventually erode out. Estimated from the profiles in Fig. 18, this might take about 10⁶–10⁷ central body revolutions in the µ = 0.03 simulation, corresponding to about 800–8000 years in the Chariklo case. Without an active confinement, the particles pushed outside the resonance experience continuous viscous spreading.

In order to avoid this leakage, a small satellite may be placed outside the ring, as illustrated in the lower row of Fig. 17. This hypothetical satellite exerts a negative torque on the ring material through Inner Lindblad Resonances (ILR) m/(m − 1), where m > 0. If acting alone, such a satellite excites spiral density enhancements in its inner Lindblad resonances, and the associated torque pushes the ring inward (lower left frame in Fig. 17 and the bottom row in Fig. 18). If combined with a mass anomaly, the satellite may stop the outward leakage of particles from the ringlet, and lead to a steady state where the outward torque by 1/3 SOR is balanced by the satellite induced negative torque. The torque associated with a m/(m − 1) resonance is classically given by (see, e.g., Sicardy et al. 2019) $${\text{Γ}}_m=\text{sign}\left(n_{\text{s}}-n\right)3{\pi }^2\left|m{(m-1)}^3\right|GM{\text{Σ}}_0{ϵ}^{\prime 2},$$(16)

where n_s is the satellite mean motion, Σ₀ is the ring surface density, and ϵ′ is a dimensionless coefficient defined in Eq. (17) of Paper I that measures the strength of the m/(m − 1) ILR. The coefficient ϵ′ is proportional to µ_s, the mass of the satellite relative to the mass of the primary (not to be confused with the mass anomaly µ of the central body).

The satellite can halt the outward leakage of the ring if |Γ_m| is larger than the viscous torque ${\text{Γ}}_v=3\pi n_0{a}_0^2\nu {\text{Σ}}_0$. Using Eq. (5), the condition |Γ_m| ≳ Γ_ν provides an order-of-magnitude estimation of the strength ϵ′ (and thus on the satellite mass µ_s) necessary to balance the ring outward diffusion, $${ϵ}^{\prime }\gtrsim \sqrt{\frac{k_{\text{visc}}\tau }{\pi m{(m-1)}^3}}\left(\frac{R}{a_0}\right),$$(17)

where a₀ ≈ 2.08 for the 1/3 SOR. Using the methodology of Paper I, it can be shown that ϵ′ ≈ 0.6/m for values of m larger than a few times unity. From our simulations with a rebound coefficient ϵ_n = 0.1, we obtain k_visc ≈ 3.5, knowing that larger values of k_visc would increase that factor a bit without changing its order of magnitude. Using these values, the Equation (17) can be re-expressed as $${\mu }_{\text{s}}\gtrsim 0.8\sqrt{\frac{\left|m\right|\tau }{|m-1{|}^3}}R$$(18)

when restricted to a 1/3 SOR ringlet.

We apply this equation to the run shown in Fig. 17 where m = 8, and R = 5×10⁻⁴. Inserting the peak optical depth at the ringlet core, τ = 0.15, would require µ_s ≳ 3 × 10⁻⁵. However, since the ringlet itself is confined by the 1/3 SOR, it is more relevant to use the τ of the tail in this estimate, τ ≈ 0.01, which indicates µ_s ≳ 7 × 10⁻⁶. This is of similar order of magnitude with the result of Fig. 17, where a satellite with mass µ_s = 2 × 10⁻⁶ prevents the viscous spreading of the ring. We note that the m-number has no special role in preventing the leaking: depending on mass and distance of the exterior satellite, a balance could be achieved at an ILR with a different m.

Fig. 16 Eccentricities e of particles vs. semimajor axis a. Unlike in Fig. 7, here we do not use the quantity ā because its definition depends on the particular resonance considered (Eq. (2)). However, because e remains small, this choice does not change the interpretation of this figure taken from the run shown in Fig. 15. We note that at time T = 21 600, the narrow second-order 3/5 SOR is able to significantly perturb the ring, in spite of the fact that the first-order 2/3 SOR is still expected to exert a nonnegligible perturbation, amounting to 10% of the eccentricity forced by the 3/5 SOR. At time T = 112 000, we see the damping effect of collisions as the ring evolves between the 3/5 and 1/2 SORs. The run is stopped at T=249 720 Chariklo’s rotations, corresponding to about 200 years. We note that third-order resonances, for example the 4/7 SOR at a=1.45, have no effect on the migration. Two movies generated from snapshots stored every 360 Chariklo rotations (about 105 days) are available online: Movie6 displays the time evolution of eccentricity versus semi-major axis and Movie7 shows the ring migration in a cartesian frame rotating with Chariklo.

8 Self-gravity

The simulations presented so far have not included mutual gravity between ring particles, but have concentrated on the collisional confinement at the resonances. However, self-gravity is expected to have significant influence on the ring dynamics. In low density rings the gravitational scattering in binary encounters enhances the steady-state velocity dispersion and thereby the ring viscosity, while for larger densities the continuously forming and dissolving self-gravity wakes dominate the dynamics (Salo 1992, 1995). The wakes increase the viscosity significantly, both via gravitational torques and due to increased velocity dispersion associated with wake motions (Daisaka et al. 2001). Finally, for sufficiently large planetocentric distances, the particles start to collect to gravity-bound aggregates. Although a fully realistic (see below) treatment of self-gravity would require much larger N than achievable in our current simulations, we here briefly address how self-gravity might affect the ring evolution, and in particular whether the confinement at 1/3 SOR could still work.

Fig. 17 Upper row: long-term evolution of a 1/3 SOR ringlet perturbed by a mass anomaly µ = 0.03, using 30 000 particles with radius R = 5 × 10−4 (corresponds to 100 m for Chariklo rings). Although the ringlet maintains a well-defined core, a slow outward leakage of particles appears, forming a faint tail. For better viewing, the width of the ring has been expanded by a factor of 10. Lower left panel: effect of a satellite located at 2.3 ≈ 1.1a1/3, with a mass µs = 2 × 10−6 relative to the central body. Without the 1/3 SOR, the ring spreads inward under the effects of spiral density perturbations forced by Lindblad resonances, shown here at T=100 000. Lower right panel: when both the satellite and the 1/3 SOR are present, the outward leakage of particles is prevented by the 8/7 Lindblad resonance with the satellite. Two movies generated from simulation snapshots stored every 360 Chariklo rotations (about 105 days) are available online. Movie8 displays the evolution of eccentricity versus semi-major axis in the simulation without an external satellite: after T ∼ 60 000 (about 50 years for Chariklo), the particles start to leak outwards due to a small torque associated with the 1/3 SOR. Movie9 shows the same for the simulation with the moonlet: the leakage observed when only the 1/3 SOR is acting on the ring is prevented, and the ring ends up being confined between the 1/3 SOR and the 8/7 inner Lindblad resonance with the satellite.

8.1 Scaling with the r_h parameter

The importance of self-gravity is conveniently described by the dimensionless parameter r_h (Ohtsuki 1993; Salo et al. 2018), that compares the size of particle pair’s gravitational Hill radius to their physical sizes. For a pair of identical particles at distance a from a spherical central body with radius R_B, $$r_{\text{h}}=\frac{R_{\text{H}}}{2R}={\left(\frac{\rho }{12{\rho }_{\text{B}}}\right)}^{\frac{1}{3}}\left(\frac{a}{R_{\text{B}}}\right),$$(19)

where ρ and ρ_B are the bulk densities of the particles and central body, respectively. Here R_H = (2M_p/3M_B)^1/3a is the radius of the Hill-sphere of two particles with masses M_p, inside which the pair’s mutual gravity dominates over the tidal pull from the central body at distance a. When r_h decreases, the particle pair extends more and more out from its Hill-sphere: r_h = 0 corresponds to the nongravitating case, while for r_h = 1, the net attraction between two synchronously rotating, radially aligned ring particles in contact equals the disruptive tidal force. The classical Roche limit, a/R_B ≤ 2.456 (ρ_B/ρ)^1/3 for the tidal disruption of a fluid body corresponds to r_h ≤ 1.072.

Both the influence of binary encounters and self-gravity wakes can be written in terms of r_h, in addition to τ and nR which alone were sufficient to describe the state of a nongravitating ring for a given coefficient of restitution. The velocity dispersion maintained by encounters is comparable to the two-body escape velocity, $v_{esc}=\sqrt{2G{M}_p/R}$. In terms of r_h, this velocity dispersion is $$\frac{c_{\text{esc}}}{nR}=4.9r_{\text{h}}^{3/2}$$(20)

For r_h ≳ 0.7, c_esc exceeds the typical dispersion maintained by impacts alone, c_imp ≈ (2 − 3)nR. At small τ when self-gravity wakes are weak, the viscosity is enhanced due to encounters roughly by a factor (c_esc/c_imp)² compared to the value given in Eq. (5). For r_h ∼ 1, this corresponds to a factor of ∼4 increase in viscosity. A same enhancement would be obtained by using a particle size increased by a factor of ∼2 in a nongravitating simulation.

Similarly, the Toomre critical wavelength and velocity dispersion can be written in terms of r_h as (see, e.g., Salo et al. 2018) $$\frac{{\lambda }_{\text{cr}}}{R}=48\pi \tau r_{\text{h}}^3,$$(21) $$\frac{c_{\text{cr}}}{nR}=12.8\tau r_{\text{h}}^3.$$(22)

The wake structure, with a typical radial spacing ∼λ_cr starts to emerge whenever the radial velocity dispersion maintained by impacts drops below (2 − 3)c_cr, corresponding to $\tau r_{\text{h}}{}^3\gtrsim 0.1$. Strong wakes imply substantially increased viscosity, both due to gravitational torques exerted by the wakes and due to increased random motions. The standard formula for the gravitational viscosity reads (Daisaka et al. 2001) $$v_{\text{grav}}\approx 190r_{\text{h}}{}^{11}{\tau }^{2}n{R}^2,$$(23)

and including the wake motions, this leads to the total ν_tot ≈ 2ν_grav. Comparing to Eq. (5), this implies an enhanced viscosity by a factor $\sim 100\tau r_{\text{h}}{}^{11}$ over a nongravitational system (see also Fig. B1 in Salo & Mondino-Llermanos 2025).

For r_h approaching unity, the wake structure becomes increasingly clumpy, and for r_h ≳ 1.1–1.2 the wakes degrade to semipermanent particle aggregates (Salo 1995; Karjalainen & Salo 2004). Accretion takes place also in low density rings via pairwise accumulation of particles (see Fig. 16.24 in Salo et al. (2018) for an illustration of various parameter domains dominated by impacts, encounters, self-gravity wakes, and gravitational accretion).

Fig. 18 Comparison of the semimajor axis distributions in the three simulations of Figs. 17 and 19. The upper frame is the run with mass anomaly, and the colors indicate profiles at various times T = 20 000 to 200 000. The dashed gray line is the initial distribution, and the arrow marks the 1/3 SOR location. We note the growing tail of the distribution in profiles for T ≳ 100 000. In the middle frame, the simulation includes both the mass anomaly and a satellite. The run extends to T = 150 000 revolutions, with no leakage or spreading of the ringlet. The bottom frame shows the simulation with a satellite alone, with arrows marking its inner Lindblad resonances. We note the inward spreading of the ring and the gradual accumulation of material to the m = 5 resonance, and the disappearance of the m = 9 peak.

Fig. 19 Evolution of mean angular momentum of the 1/3 SOR ring perturbed by a µ = 0.03 mass anomaly (black curve), by a µs = 2 × 10−6 satellite (green), or simultaneously by both (red). The snapshots of these simulations at T = 100 000 are displayed in Fig. 17.

8.2 Self-gravitating simulations

8.2.1 Accretion boundary

Inserting the values of Table 1 into Eq. (19) implies that for the Chariklo ring system $$r_{\text{h}}=0.697{\left(\frac{\rho }{900\text{kg}{\text{m}}^{-3}}\right)}^{1/3}\frac{a}{a_{\text{cr}}}=1.45{\left(\frac{\rho }{900\text{kg}{\text{m}}^{-3}}\right)}^{1/3}\frac{a}{a_{1/3}},$$(24)

so that the 1/3 SOR distance corresponds to r_h ≳ 1 whenever ρ ≳ 300 kg m⁻³. Thus, even for relatively under-dense particles with ρ = 200–300 kg m⁻³, the C1R ring region should be strongly affected by self-gravity.

Figure 20 compares two simulations starting with a wide initial ring around the 1/3 SOR, spanning the range a = 1.55–2.80; for the adopted bulk density of particles (ρ = 300 kg m⁻³) this corresponds to r_h = 0.75–1.35. In the first simulation (upper row) a spherical central body is used, while the second (lower row) uses µ = 0.3. A large particle radius R = 0.005 (corresponding to 1 km in Chariklo’s case) is adopted, yielding an initial optical depth τ₀ = 0.11, which should be sufficiently large for self-gravity wakes to become discernible. A large τ also speeds up the formation of gravitational aggregates at large distances. Indeed, in both simulations particle aggregates start to form rapidly, within first few ring orbital periods, in the region r_h ≳ 1.15 (marked by the solid line on the right axis of the T = 20 frames). Inside this distance self-gravity wakes form, inclined by ∼20^◦ with respect to the tangential direction and with radial spacing ∼0.1, consistent with Eq. (21). There are also a few elongated streaks visible at the border zone between the wake and accretion regions: these are formed by particle aggregates recently destroyed by tidal forces or due impacts. Eventually, the aggregates manage to merge, and in the end of both runs at T = 150, about 95% of the particles beyond r_h = 1.15 have collected into a single aggregate. A significant viscous spreading is revealed by the inward motion of the ring inner edge in the µ = 0 simulation.

For the simulation with a mass anomaly, a large µ was chosen, in order to make the resonance perturbations significant in spite of the large R and the very short duration of the run. Indeed, the particles are gradually pushed outward due to resonance torques: in the frame at T = 20, a weak m = 2 undulation is visible, related to the 1/2 SOR at a = 1.59. Conversely, the large amplitude m = 1 pattern at the T = 150 snapshot is related to 1/3 SOR. However, the ring viscosity is much too large to allow an efficient resonance confinement to take place. We also note that self-gravity wakes and particle clumps are totally absent in the perturbed region at the T = 150 frame (lower right panel of Fig. 20).

8.2.2 Resonance confinement

The condition derived in Section 4.3 for the resonance confinement implies $\sqrt{v}$. To check whether the resonance confinement is possible in the presence of self-gravity in spite of the viscosity enhancement, we added particle gravity to a non-gravitating simulation that has already formed a confined ringlet. We selected as a starting point the run shown in the upper-left corner of Fig. 9, corresponding to µ = 0.1 and R = 0.0005 (100 m in Chariklo’s case), For this value of µ, even a 4-fold increase in R (i.e., a factor of 16 in viscosity) still permits a confinement in the nongravitating case. We note that this radius R is ten times smaller than in the examples of Fig. 20. Three simulations with ρ = 300, 450, and 900 kg m⁻³, corresponding to r_h = 1.0, 1.15, and 1.44, respectively, started from the ringlet stage at T = 10 000. The time evolution of particle distributions with various ρ’s is followed in Fig. 21. The first and second columns display polar snapshots at T = 50 (about 17 ring revolution periods after the inclusion of self-gravity), and at T = 400 (the end of the self-gravitating run), while the right column shows the time evolution of the ringlet mean width (uppermost frame) and the semimajor axis distributions.

With ρ = 300 kg m⁻³ the ringlet remains well-confined with sharp inner and outer boundaries. Nevertheless, its width has become about two-fold compared to the nongravitating case. This increase in the width takes place rapidly during the first few ring revolutions, accompanied by a similar increase in the velocity dispersion. During later evolution, the width of the ring as well as the dispersion of the semimajor axis is practically constant. With larger bulk densities a rapid formation of particle aggregates takes place. At T = 50 with ρ = 450 kg m⁻³, the ringlet has developed several azimuthal gaps, connected to the largest individual aggregates which scatter particles out from the ringlet core. For ρ = 900 kg m⁻³, the whole ringlet has become quite fuzzy. Nevertheless, after the initial evolution the ringlet spreading seems to stop (ρ = 450 kg m⁻³) or even become reversed (ρ = 900 kg m⁻³).

Fig. 20 Self-gravitating simulations with R = 0.005 (corresponding to a physical radius of 1 km in the Chariklo case), assuming particles with bulk density ρ = 300 kg m−3. The initial distribution extends to r = 1.55–2.8, which corresponds to the range rh = 0.75–1.35. With N = 24 000 particles the initial τ0 = 0.11. The upper row shows two snapshots of the simulation with no mass anomaly, while the lower row shows a simulation with a µ = 0.3 mass anomaly. The location of the 1/3 SOR (a = 2.08) corresponds to rh ≈ 1.0. In the T = 20 frames, the tick marks on the right vertical axis indicate the radii where rh = 1.15. The red and green solid lines in the µ = 0.3 frames indicate the locations of 1/2 (left) and 1/3 (right) SOR. The insets in the T = 150 frames zoom in on the largest aggregate formed during the run, covering a 0.4 × 0.4 region with a correct aspect ration. The slice of the aggregate through the z = 0 plane is shown.

8.2.3 Cyclic formation and destruction of aggregates

Based on Fig. 21 it seems plausible that a ringlet can remain confined at the 1/3 SOR for values ρ ≤ 900 kg m⁻³. However, due to the short duration of the simulation, it is unclear whether an actual steady-state has been achieved in the runs with ρ = 450 kg m⁻³ or 900 kg m⁻³. Therefore, another simulation with ρ = 450 kg m⁻³ was conducted, now extended over 20-times longer duration but reducing the particle number by one half (N = 15 000) to speed-up the calculation. Fig. 22 shows the evolution of ringlet width and velocity dispersion⁵, and also quantifies the amount of particles in aggregates, separately in the small (10 < N_group < 100; gray color) and large groups (N_group > 100; blue), and in the largest aggregate (red).

The initial evolution in the simulation is practically similar to the corresponding run in Fig. 21: the velocity dispersion (c_z ∼ Rn) is initially less than the escape velocity of individual particles (c_esc ∼ 6Rn) and therefore small particle groups are rapidly formed. However, the gravitational stirring by the growing groups heats the system, soon preventing further pairwise accretion. The already formed groups gradually grow through merging (and occasionally destroy each other in impacts), until at T ∼ 200, the system is dominated by one large aggregate, with mass on the order of 4% of the total ringlet mass (i.e., involving ∼600 particles), exceeding the total mass in other aggregates. At T ∼ 300 this large aggregate looses half of its mass by tidal leakage and at T = 500 it breaks completely. When this happens, the system starts to cool down due to collisional dissipation as there is no more gravitational stirring by the aggregates.

Interestingly, the aggregate formation/destruction appears to be a self-regulating process: once the system has cooled down sufficiently (c_z ∼ 5Rn, at T ≈ 1750), the accretion starts again, going through the formation phase of increasingly large particle groups, until at T ≈ 2300 the last large aggregate is again tidally destroyed, and the cooling phase starts again. This cycle repeats throughout the simulation, though with somewhat diminished variations. The alternating cycle is also visible in the snapshots of the system, the ringlet appearance varying from sharper to fuzzier, depending of whether aggregates are present or not (see the insets in Fig. 22: the individual aggregates are too small to be discernible in the plots, except indirectly in the snapshot for T = 6300, displaying the system just after a very large particle aggregate with 2000 particles has been tidally shredded, the debris being still visible).

The simulation of Fig. 22 suggests that quite an interesting type of behavior is possible near the accretion boundary at r_h ∼ 1.1–1.2. To check that the observed tidal disruption of large aggregates is not a numerical artifact, but indeed due to the perturbation, an additional simulation was conducted. It started from the above simulation at T = 6250, shortly before the breakup of the big aggregate, and used a time step reduced to one-half (1/1200 of the ring orbital period). Also in this case the aggregate broke up, although a bit later than in the original run. In both runs the break-up happened via gradual stretching in the radial direction. On the other hand, another experiment starting from T = 6250, where the single large aggregate was isolated and all other particles removed, led to a stable aggregate (surviving to the end of the run lasting to T = 7000). This suggest that the tidal break-up is induced by impacts from particles perturbed by the resonance, rather than some secular instability due to errors accumulated during the orbital integration of the aggregate.

Fig. 21 Self-gravitating simulations of N=30 000 particles with different bulk densities. Simulations continue from the ringlet stage (at T = 10 000) of the nongravitating simulation of Fig. 9 with µ = 0.1, R = 0.0005 (100 m in Chariklo’s case): the snapshot of the initial state is shown in the upper left frame. The other frames display self-gravitating simulations with ρ= 300, 450, and 900 kg m−3, corresponding to rh = 1.0, 1.15, and 1.44, respectively at the distance of the 1/3 SOR. Times are counted from the beginning of the self-gravitating simulation. The upper right frame displays the evolution of the ringlet mean width, $W=\sqrt{12\langle {\left(r-r_{\text{fit}}\right)}^2\rangle }$, where rfit denotes the fitted ringlet central line at the azimuth of the particle and the averaging is over all particles (for a uniform particle distribution this formula with the factor of 12 recovers the full width of the ringlet). The remaining frames on the right show the semimajor axis distribution at different instants for the three self-gravitating runs.

9 Discussion

Our main result is that a mass anomaly µ leads to the confinement of ring material near the 1/3 SOR, provided that the strength of the resonance overcomes the spreading effect of collisions (see Eq. (10) and Fig. 9). In analogy with Saturn’s dense rings, we may assume that the dynamics of the dense rings observed so far is dominated by large particles with radius on the order of a meter (Cuzzi et al. 2018). Using for instance a corotation radius a_cor ∼ 200 km in the Chariklo case (Table 1), we obtain R ∼ 5 × 10⁻⁶ in Eq. (11). Considering that τ_Q1R ∼ 1, this gives a threshold value µ ≳ 10⁻³ that permits ring confinement. If the mass anomaly µ is caused by a mountain of height h at the surface of a body with radius R_ref, then µ ∼ 0.5(h/R_ref)³. Using R_ref = 115 km for Chariklo, then µ ∼ 10⁻³ corresponds to h ∼ 10 km, a plausible value when compared with other small objects of the Solar System.

Eq. (11) shows that the threshold value µ is proportional to R, which is the particle size normalized to a_cor. Consequently, for a given particle size in meters, µ scales like $a_{\text{cor}}^{-1}$. Since Haumea’s and Quaoar’s corotation radii are respectively six and ten times greater than for Chariklo (Table 2 of Paper I), the threshold value of µ is on the order of 10⁻⁴ for these two bodies.

It remains to be seen if these objects can sustain such mass anomalies. They can be compared with various small bodies of the Solar System that exhibit large topographic features⁶. For instance Ceres (R_ref ∼ 470 km) has mountains with h/R_ref ∼ 8.5 × 10⁻³, while Vesta (R_ref ∼ 270 km) has features with h/R_ref ∼ 8.4 × 10⁻². This corresponds to µ = 3 × 10⁻⁷ and µ = 3 × 10⁻⁴, respectively. The large Trans-Neptunian Object (307261) 2002 MS₄ (R_ref ∼ 400 km) possesses a large depression with depth ∼45 km and an extension of some 320 km, corresponding to µ of a few times −10⁻² (Rommel et al. 2023).

As noted in Section 3.1, while Chariklo’s, Haumea’s and Quaoar’s main rings orbit close to the 1/3 SOR, Quaoar’s fainter ring Q2R is close to the second-order 5/7 SOR. This is the starting point of the simulation shown in Fig. 16. However, this run uses Chariklo’s parameters, and the ring is rapidly evacuated from this resonance due to the large effects of the nearby 3/4 and 2/3 SORs. In the case of Quaoar, the 5/7 SOR is better isolated from the 3/4 and 2/3 SORs (Fig. 9 of Paper I), and future numerical simulations could test the ability of the 5/7 SOR to confine a ring. This said, we note in Fig. 16 that the material is temporarily trapped in another second-order SOR corresponding to n/Ω_B = 3/5, showing that the 1/3 resonance is not the only one that can confine a ring.

Our simulations show that the ring confinement is accompanied by the excitations of various free Lindblad modes with azimuthal wave numbers m (Fig. 12). This was observed most strikingly in our best simulation with µ = 0.003 and R = 25 m, where the resonance excitation was eventually divided between normal modes with several values of m (Fig. 12). In simulations with larger µ’s and R′, signs of similar excitations were seen, but not all normal modes were present. The presence of such modes could be revealed using multi-chord stellar occultations by providing accurate orbital radii and ring widths at various longitudes. Unfortunately, even for the best-observed occultations by Chariklo’s system, the azimuthal sampling of the rings is still too coarse, Chariklo’s center cannot be pinned down accurately, and the ring precession rate $\dot{\varpi }$ is poorly constrained. This prevents any accurate analysis of these modes. Hints for the presence of m = −1 and m = −2 Lindblad modes in Chariklo’s main ring C1R has been reported by Gomes-Junior et al. (2025). Although encouraging, these detections still require confirmation through more detailed observations.

As discussed in Section 7 and shown in Fig. 17, the 1/3 SOR confinement by a mass anomaly alone is not sufficient to explain the long-term presence of a ring. A small satellite is required to balance the outward viscous spreading of the confined ringlet. Eq. (18) shows that the mass of this satellite is proportional to the size of the ring particles R. Assuming meter-sized particles, this yields R = 5 × 10⁻⁶ in the Chariklo case, which is 100 times smaller that the radius used for the run shown in Fig. 17. Eq. (18) then indicates that a satellite with mass µ_s ∼ 10⁻⁷ can stabilize the outward spreading of a 1/3 SOR ringlet. This value is comparable to the masses obtained by Sickafoose & Lewis (2024) – several times 10⁻⁶ – estimated from their numerical simulations of Chariklo’s rings, taking particle radii of a few meters confined by various inner and outer Lindblad resonances.

The value µ_s = 10⁻⁷ corresponds to an icy moonlet as small as 0.5 km in radius in the case of Chariklo. The presence of such small satellites can be expected as Chariklo’s rings are close to the classical Roche limit, allowing the mixture of satellites and rings in the 1/3 SOR region. The same is true for Haumea’s and Quaoar’s rings (Hedman 2023). We note from Eq. (18) that again ${\mu }_{\text{s}}\propto a_{\text{cor}}^{-1}$. Thus, for Haumea and Quaoar, the threshold values of µ_s are ∼1.5 × 10⁻⁸ and ∼10⁻⁸, corresponding to icy moonlets of radii on the order of 2 km and 1 km, respectively.

The numerical estimates of the various threshold values for µ and µ_s are based on nongravitating simulations. However, for the range of plausible bulk densities of ring particles, ρ = 300–900 kg m⁻³, self-gravity can be expected to have a significant effect in Chariklo’s rings. In terms of the Hill dimensionless parameter r_h (see Section 8), this range of densities corresponds to r_h = 1.0–1.44. This implies that when applying Eq. (11), the viscosity enhancement due to gravitational encounters (τ ≲ 0.1) and/or self-gravity wakes (τ ≳ 0.1) are both important. The enhancement due to gravitational encounters was estimated to be on the order of 5–10, while wakes can lead to even 100-fold viscosity when r_h ∼ 1. In terms of threshold value for µ which is proportional to square-root of viscosity, such enhanced viscosities would require roughly 3 and 10 times larger minimum size for the mass anomaly in order to obtain confinement, or similarly, a smaller maximum particle size for a given value of µ.

The self-gravitating simulations of Section 8 indicate that resonantly confined ringlets are surprisingly resilient against gravitational accretion, whereas unperturbed rings experience rapid accretion once r_h ≳ 1.15. In our self-gravitating experiments the starting point was taken from a confined ringlet formed in a nongravitating simulation: after switching-on self-gravity the ringlet settled into a more diffuse, dynamically hotter state, while still remaining confined. During the initial evolution, when the velocity dispersion was still small, aggregates rapidly formed, but were eventually destroyed via tidal stretching. Near the accretion limit r_h = 1.15 a quite interesting cyclic behavior was seen, with the system alternating between states where aggregates were abundant or nearly absent.

The low τ simulations of Section 8 covered only the influence of gravitational encounters and particle accretion, as the current particle number is far too small for realistic modeling of self-gravity wakes: Fig. 20 was merely an illustration of the regimes where wakes/accretion take place, using a very large R. Indeed, a realistic modeling of a dense narrow ring with self-gravity wake structure would require R << λ_cr << W, whereas in the nongravitating simulations it is sufficient to have R << W, where λ_cr is the Toomre critical wavelength and W is the ring width. This extra intermediate length scale implies at least a factor of 10 smaller particles, increasing N by factor of at least 100 for a fixed τ (in fact, near Roche limit, λ_cr/R ∼ 100, indicating still larger separation of scales). Moreover, the optical depth should be on the order of 0.1–1, further increasing N. Thus at least 10⁶–10⁷ particles would be needed, instead of the few tens of thousands employed in the current simulations⁷.

Fig. 22 Self-gravitating simulations with µ = 0.1, R = 0.0005 (100 m in Chariklo’s case), assuming particles with bulk density ρ = 450 kg m−3. The initial particle distribution is similar to that in Fig. 21, except that τ has been reduced by a factor of 2 by selecting only every other particle from the confined ringlet in the nongravitating simulation at T = 10 000. The upper frame displays the evolution of the ringlet mean width, together with snapshot corresponding to various width maxima and minima; the snapshot at T = 6300 displays the system just after tidal breakup of the largest aggregate. The lower frame displays the number of particles in the largest aggregate (red) and the total number of particles in groups with more than 100 (blue) and 10 particles (gray; the total number of particles N = 15 000). Also shown as a green curve is the vertical dispersion cz (labels on the right).

10 Conclusions

In connection with the analytical calculations of Paper I, we investigated the effects of spin-orbit resonances (SORs) on a collisional ring surrounding an irregular body, using fully 3D collisional simulations. Our main conclusion is that the 1/3 SOR between the ring and the central body is efficient in confining ring material, while a small outer satellite may prevent the long-term outward spreading of material. More specifically, our main results are the following:

1. Our simulations confirm the theoretical expectation that only first- and second-order SORs effectively disturb a collisional disk;

2. We can reproduce the peak eccentricity and the resonance width expected from our analytical calculations (Figs. 2, 6 and 16). This validates the use of the (ā, e) space to study the ring evolution (Fig. 6);

3. The dense mesh of first- and second-order resonances around Chariklo, Haumea, and Quaoar excite large eccentricities and rapidly clear up a colliding disk up to the 1/2 resonance (Figs. 15–16). In that context, the outermost second-order 1/3 SOR forced by a mass anomaly is a quieter place, immune from the interactions with the inner SORs (Fig. 1);

4. The 1/3 SOR first forces a self-intersecting streamline (Figs. 6 and C.2) that gathers the background material of the disk into a ringlet. The self-intersection problem is avoided by the fact that the eccentricity excitation caused by the 1/3 SOR is transferred into superposed non-self-intersecting free Lindblad modes, which eventually leads to the ring confinement (Figs. 6 and 12);

5. We obtain a condition for the confinement of a ringlet at the 1/3 SOR. It expresses the fact that the collisional viscous spreading is counteracted by the 1/3 SOR confinement effect. This condition, 4 × 10⁻⁵µ² ≳ τR², relates the optical depth τ, the radius R of the particles normalized to the radius of the synchronous orbit, and the value µ of the mass anomaly (see Eq. (11) and Fig. 9). Application to Chariklo shows that a mass anomaly µ ≳ 10⁻³ can confine material at the 1/3 SOR, which corresponds to typical topographic features (mountains or craters) on the order of 10 km;

6. A long-term outward viscous spreading of the 1/3 SOR ringlet is observed in our simulations. We find that this leakage can be prevented through Lindblad resonances raised by an outer moonlet of mass µ_s ≳ 10⁻⁷ (Fig. 17) relative to Chariklo, corresponding to a subkilometer moonlet. However, it is important to note that even in this case the ring confinement is still ensured by the 1/3 SOR;

7. These results can be applied to Haumea’s and Quaoar’s rings by noting that the threshold values of µ and µ_s mentioned above scale as the inverse of the radius of the synchronous orbit a_cor. Therefore, the threshold values of µ necessary to confine a 1/3 SOR ring are ∼10⁻⁴ for these two bodies. The threshold value of the satellite mass that can prevent the viscous spreading of such a ring is on the order of 10⁻⁸;

8. Our preliminary experiments with self-gravity indicate that the confinement mechanisms also work for self-gravitating particles, although due to the enhanced viscosity, the threshold $\mu \propto \sqrt{v}$ is increased. At low τ the enhancement in ν is due to gravitational encounters and amounts to less than a factor of 10 compared to nongravitating values. In large τ, the gravitational wakes may enhance the viscosity even by a factor of 100. The implied increase in threshold µ is thus by a factor of 3–10.

Our results provide encouraging evidence that second-order resonances can confine ring material around a nonaxisymmetric object. However, several pieces of the puzzle are still missing to fully explain this confinement.

The initial material confinement observed at the 1/3 SOR (see the lowermost panels of Fig. 6) can easily be understood by the fact that the velocity field in the forming ringlet is reversed when compared to the background Keplerian motion: at its outermost radial position, a ring particle moves slower than the background particles on circular orbits, and vice versa at its innermost position. Things get more complicated when the ring becomes narrow. The confinement is then reminiscent of the single-sided shepherding seen in the simulations of Hänninen & Salo (1994, 1995) that described ring confinement at Lindblad resonances forced by a satellite. This result was explained analytically by the negative angular momentum luminosity forced by the satellite (Goldreich et al. 1995). Unfortunately, this analytical model is not readily applicable to the 1/3 SOR, due to the impossibility of avoiding self-intersection for a pure 1/3 SOR streamline. Thus, the next important step is to explain how the 1/3 SOR excitation is transferred to free Lindblad modes, and how these modes lead to the reversal of angular momentum flux.

While analytical expressions for the torques exerted by Lindblad resonances have been derived decades ago, there is no such expressions for second-order resonances, while such torque is observed in our simulation (see, e.g., the passage through the 3/5 SOR in Fig. 15 that show a clear jump in angular momentum. The calculation of this torque would be important in order to assess the long-term stability of a 1/3 SOR ringlet. At this point, and considering the difficulty to describe how free Lindblad modes are excited at the 1/3 SOR, it is not clear if this analytical expression is obtainable.

Another important topic to study further is the effect of self-gravity on the ring confinement. Our small-N simulations have addressed mainly the influence of gravitational encounters, while a much larger N would be needed to study realistic dense rings with self-gravity wakes. The gravitational accretion of particles in perturbed versus unperturbed rings would also merit a dedicated study. In particular, the cyclic behavior of perturbed rings near the accretion boundary, potentially relevant for example to the time-dependent structure of Saturn’s F-ring, should be verified by independent simulations.

Our simulations have focused on the case of the 1/3 SOR using Chariklo’s parameters. It is now important to have simulations more specifically focused on Haumea and on the 1/3 and 5/7 SORs around Quaoar.

Independent of rings around nonaxisymmetric bodies, it is also important to extend this work to the confinement of material at second-order resonances caused by forming planets in circumstellar disks or in the proto-Solar System.

Data availability

The movies associated to Figs. 6, 16 & 17 are available at https://www.aanda.org

Acknowledgements

This work has been supported by the French ANR project Roche, number ANR-23-CE49-0012.

References

Appendix A Potentials acting on the ring

A.1 Triaxial ellipsoid

The analytical expansions and numerical values for the potential of homogeneous triaxial ellipsoids are given in the Appendix B of Paper I. We denote the total mass of the body by M and the semiaxes of the ellipsoid by A, B, C, with A > B > C. The reference radius R_ref, elongation ϵ_elon and oblateness f of the body are defined by⁸ $$\frac{3}{R_{\text{ref}}^2}=\frac{1}{A^2}+\frac{1}{B^2}+\frac{1}{C^2},$$(A.1) $${ϵ}_{\text{elon}}=\frac{A^2-B^2}{2R_{\text{ref}}^2},$$(A.2) $$f=\frac{A^2+B^2-2C^2}{4R_{\text{ref}}^2}.$$(A.3)

Outside the ellipsoidal body, but near to its equatorial plane (|z| << r) we may approximate its potential in co-rotating cylindrical coordinates as $$\begin{array}{ll}U(r,\theta ,z)=\hfill & -\frac{GM}{r}\left[1+\frac{f}{5}{\left(\frac{R_{\text{ref}}}{r}\right)}^2+{\displaystyle \sum }_{p=1}^{\infty }{\left(\frac{R_{\text{ref}}}{r}\right)}^{2p}{S}_p{ϵ}_{\text{elon}}^p\cos (2p\theta )\right]\hfill \\ \hfill & -\frac{1}{2}\frac{GM}{r^3}\left[1+\frac{9f}{5}{\left(\frac{R_{\text{ref}}}{r}\right)}^2\right]z^2.\hfill \end{array}$$(A.4)

In this formula the origin is fixed to the center of the body and the angle θ is measured from the longest semi-axis A. The method to calculate the coefficients S_p is given in the Appendix B of Paper I, and in practice we include terms up to p = 10.

We assume that the central body rotates counterclockwise around its shortest axis (z-axis), so that at time t the longest axis makes an angle Ω_Bt with respect to the x-axis of the inertial system (initially the central body long axis is aligned with the x-axis). At each force evaluation step, the nonrotating cartesian coordinates of particles are transformed to the system aligned with the body ellipsoid, by rotation with the angle Ω_Bt around the z-axis. The accelerations components in this system are calculated from −∇U and then transformed back to the nonrotating system by rotation with the angle −Ω_Bt.

A.2 Mass anomaly

In the Appendix A of Paper I, the potential caused by a point-like mass anomaly on the central body was expanded in rotating cylindrical coordinates centered onto the body itself, the most convenient approach to investigate the structure of the various resonances at play.

In numerical simulations, it is simpler and computationally faster to write the forces directly in the cartesian center-of-mass system. Let M stand for the total mass of the central body (including the mass anomaly), and define µ as the dimensionless fractional mass contained in the mass anomaly, located at the distance R_ref from the center of the body along its longest semi-axis. In the center-of-mass system and at time t, the center of the body and the mass anomaly have respective coordinates r₀ = (x₀, y₀, 0) and r₁ = (x₁, y₁, 0), where $$\begin{array}{l}{x}_0=-\frac{\mu }{1+\mu }{R}_{\text{ref}}\cos \left({\text{Ω}}_{\text{B}}t\right),y_0=-\frac{\mu }{1+\mu }{R}_{\text{ref}}\sin \left({\text{Ω}}_{\text{B}}t\right),\hfill \\ x_1=\frac{1}{1+\mu }{R}_{\text{ref}}\cos \left({\text{Ω}}_{\text{B}}t\right),y_1=\frac{1}{1+\mu }{R}_{\text{ref}}\sin \left({\text{Ω}}_{\text{B}}t\right).\hfill \end{array}$$(A.5)

The acceleration felt by a ring particle i is $$\begin{array}{ll}{\ddot{r}}_i\hfill & =-GM(1-\mu )\frac{r_{i0}}{r_{i0}^3}-GM\mu \frac{r_{i1}}{r_{i1}^3}\hfill \\ \hfill & =-{\nabla }_i\left[-\frac{GM(1-\mu )}{\left|r_{i0}\right|}-\frac{GM\mu }{\left|r_{i1}\right|}\right],\hfill \end{array}$$(A.6)

where the abbreviation r_ij = r_i − r_j is used.

In case iii) where the mass anomaly is combined with a tri-axial central body shape, the first point-mass potential term in Eq. A.6 is replaced by that calculated from Eq. A.4.

A.3 Jacobi constant

Since the central body potential U is time-independent in the frame rotating with constant angular velocity ${\text{Ω}}_{\text{B}}\widehat{z}$, the equations of the motion conserve the Jacobi energy, $$E_{\text{J}}=E_{\text{kin}}+U+{\text{Ω}}_{\text{B}}{L}_z$$(A.7) $$=\frac{1}{2}\left({\dot{x}}^2+{\dot{y}}^2+{\dot{z}}^2\right)+U(x,y,z)-{\text{Ω}}_{\text{B}}(x\dot{y}-y\dot{x}),$$(A.8)

where E_kin and L_z denote the kinetic energy and angular momentum z-component in the center-of mass inertial frame. This is a useful quantity that can be used for checking the accuracy of the numerical integrations, in case impacts are ignored.

A.4 Satellites

In the case when an additional satellite with mass M_s is included, the integrations are also performed in the system centered at the mass-center of the system formed by the central body and its mass anomaly. In this system, the acceleration of the satellite located at r_s is obtained from Eq. A.6, except that the total mass M_s + M replaces M (see, e.g., Valtonen & Karttunen 2006, p. 31): $${\ddot{r}}_s=-G\left(M+M_{\text{s}}\right)(1-\mu )\frac{r_{\text{s}0}}{r_{\text{s}0}{}^3}-G\left(M+M_{\text{s}}\right)\mu \frac{r_{sl}}{r_{sl}{}^3}.$$(A.9)

For the ring particles, the perturbing acceleration (added to Eq. A.6) due to the satellite includes both the direct and indirect parts: $$\begin{array}{ll}\text{Δ}{\ddot{r}}_i=\hfill & -G{M}_{\text{s}}\frac{r_{is}}{r_{is}^3}\hfill \\ \hfill & -G{M}_{\text{s}}(1-\mu )\frac{r_{\text{s}0}}{r_{\text{s}0}^3}-G{M}_{\text{s}}\mu \frac{r_{\text{s}1}}{r_{\text{s}1}^3}.\hfill \end{array}$$(A.10)

Appendix B Treatment of impacts

For particle impacts we employ the visco-elastic "shock absorber" technique, initially introduced by Salo (1995) for self-gravitating local simulations pertaining to Saturn’s rings. This method involves integrating the motion of particle pairs through each collision event. During the impact the particle pairs experience a repulsive force that is directly proportional to the extent of their mutual overlap. Additionally, they encounter a viscous force proportional to the perpendicular component of their relative velocity. The combination of these force components ensures that the colliding particles rebound with a decreased post-collisional relative speed. The two parameters contained in this linear force model can be written in terms of the desired coefficient of restitution ϵ_n = |v_n′|/|v_n|, which specifies the ratio of post and pre-collisional velocity differences, and the duration of the impact T_dur (see Salo et al. 2018). In all our current simulations we employ ϵ_n = 0.1.

The advantage of such a soft-particle force-method, compared to the more-standard treatment of impacts as instantaneous velocity changes, is its ability to handle multiple simultaneous impacts and gravitational sticking of particles, without leading to artificial overlaps of particles. In addition, the technique allows for realistic inclusion of ring self-gravity, even when self-gravity is strong enough to lead to gravitationally bound particle groups (see Salo 1995; Mondino-Llermanos & Salo 2022). Recently, Torii et al. (2024, 2025) adopted the same method in their global simulations of embedded satellites in Saturn’s rings.

The downside of the soft-particle method is that impacts need to be resolved with very small time steps, which increases the CPU consumption. To alleviate this problem, we use two tricks. The first is to scale-up the impact duration to be longer than the actual physical impact duration (a fraction of second), but still keep T_dur short compared to orbital period T_orb: tests reported in Salo et al. (2018) indicate that impact durations on the order of T_dur = 0.0025T_orb yield results that are indistinguishable from treating impacts as instantaneous. Similar result was found by Torii et al. (2025). For too long impact duration, the steady-state velocity dispersion following from the balance between viscous gain and collisional dissipation would be modified. For example, for the Chariklo ringlets, T_dur = 0.0025T_orb corresponds to ∼100 seconds. Another trick to speed up the calculations is to use a separate treatment of those particles which are currently colliding, compared to those which are not. The short time steps ∆t_small << T_dur, dictated by the duration of impact, are only used for the integration of currently colliding pairs, while the other particles are integrated with longer time steps ∆t << T_orb, whose length is determined by the need to integrate accurately the orbital motion.

Typically, we choose the constants of the visco-elastic force model in such a way that T_dur ∼ 0.0015T_orb, and employ an integration time step ∆t_small = 0.1T_dur for particles currently colliding. With these values, the maximum penetration in impacts typically stays below 1% of particle radius.⁹ For orbital integration, we use $\text{Δ}t=\frac{1}{300}{T}_{\text{orb}}$ unless otherwise indicated. In the self-gravitating simulation described in Section 8, two-times smaller T_dur, ∆t, and ∆t_small were used.

In practice, it is important to keep the number of pairs requiring small time steps ∆t_small as low as possible. To achieve this we construct in the beginning of each orbital integration step a list of particle pairs that may collide during the next large time step ∆t. This list is constructed in two stages: we first make an initial list of all particle pairs whose mutual distances are less than a given threshold distance R_lim1. Search of these close pairs makes use of the fact that the system is a narrow annulus and orders the particles according to their azimuthal coordinate. In the second stage, we utilize particle velocities and use a first-order Taylor-expansion to estimate the minimum separation each pair included to the initial list can achieve during the orbital integration step ∆t. If this minimum separation is smaller than a threshold separation R_lim2, then the pair is accepted to a refined list of potentially colliding pairs.

In the calculation of forces on particles in RK4, all particles which are members of the refined list of pairs (potential impactors), are integrated with the small time step ∆t_small. Not all of them are currently experiencing an impact: at each small time step, we check whether the pair actually overlaps and if so, then impact forces are added. On the other hand, those particles which are not in the list of potential impactors, are integrated with the large time step ∆t. At low optical depth systems we concentrate on, only a small fraction of particles are impacting at a given instant of time, and therefore such a splitting of particles to impacting and non-impacting particles gives easily an order of magnitude speed-up. In practice, it is important to keep the list of potential impactors as short as possible, while still making sure that all actual impacts are found. This is achieved by dynamically adjusting both R_lim1 and R_lim2. The first threshold is based on the velocity dispersion of the system, and the second one is based on what has been, during the last few integrations steps, the maximum pre-step separation actually leading to an impact.

In our self-gravitating simulations, the additional forces between particles were calculated after every ∆t, together with the time derivatives of the forces. Due to small N in the current simulations, a particle-particle force calculation method was used, including all particle pairs within a certain limiting distance R_grav. This limiting distance was chosen to be on the order of 100 particle diameters and at least 10 times bigger than the largest aggregate forming in the simulation. In practice the calculation of forces was done together with finding the nearby impact pairs. Since R_grav >> R_lim1, the CPU-time consumption in finding the pairs was substantially larger in self-gravitating simulations compared to nongravitating simulations, in particular if large aggregates formed. During the small integration time steps, a linear Taylor approximation of forces was used: in Karjalainen & Salo (2004) it was found that such an approximation improves significantly the accuracy of the gravity calculation. In particular, it removed the secular rotational instability of the aggregates seen in case constant forces were used during the step ∆t.

Appendix C Additional simulations

In this section, we present complementary simulations to test some of the results obtained in the main text.

C.1 Effect of nearby first-order 2/3 resonance on 1/3 SOR

Figure 4 indicated that the theoretical eccentricity amplitudes excited by the first-order 2/3 resonance are very large at the vicinity of the 1/3 SOR, even 30% for µ = 0.1. To demonstrate that this resonance overlap does not affect the confinement seen at collisional 1/3 resonance simulations, we have conducted additional simulations, where only the m = 2 and m = 1 terms of the mass anomaly potential are taken into account. As seen in Fig. C.1, the m = 2 term is not able to lead to any secular change in L_z distribution at the 1/3 SOR. On the other hand, using just the m = 1 leads to essentially similar evolution as using the full potential.

C.2 Wide initial distribution

Figure 6 shows the confinement of a ring near the 1/3 SOR using a small mass anomaly µ = 0.003. In order to better show the formation of the initial self-intersecting streamline forced by the resonance, we have performed a simulation with a large mass anomaly µ = 0.1 (see Fig. C.2). It clearly shows the formation of a self-intersecting streamline forced by the 1/3 SOR and its transformation into a single ringlet. Due to large initial width of the particle distribution, too additional ringlets accumulate in the same resonance site, the outer of them eventually merging with the central feature.

Fig. C.1 Four collisional simulations of 1/3 SOR with µ = 0.1, R = 0.001 (corresponds to 200 m for Chariklo’s rings). Upper row, left: Standard simulation using full mass anomaly potential. Upper row, right: Using an expansion of the potential, including m=1 and m=2 terms. Lower row: Using the m = 1 (left) or m = 2 (right) terms of the potential expansion.

C.3 Additional tests of confinement criterion

The Fig. 9 explored the effects of various parameters of the runs (particle radius R and mass anomaly µ) in order to test Eq. 10. The Fig. C.3 complements this testing by examining more specifically the effect of τ on the ring confinement for a given particle radius of R = 0.001.

According to Eq. 10, τ_lim is proportional to µ² for a fixed R. To test this behavior, two additional series of simulations were performed with varying τ₀, with µ = 0.05 and 0.03 (Fig. C.3). According to Eq. 10 the expected τ_lim ≈ 0.1 and 0.035, respectively, but based on simulations the actual limits are about 0.06 and 0.025, respectively. We note that in the case of µ = 0.05, all the simulations with τ₀ up to 0.2 eventually lead to accumulation at the resonance zone. However, this does not contradict the above threshold criterion, since the optical depth at the resonance zone drops during the simulation due to viscous spreading (see the caption of Fig. C.3).

C.4 Shear reversal in first-order resonance

Figure C.4 shows examples of the response of collisional ring to first-order 2/1 and 2/3 satellite resonances, in a frame rotating with the satellite. In case of 2/1, the m = 2 response is oriented nearly perpendicular to the instantaneous direction of the satellite, while for 2/3 the orientation is nearly aligned with the satellite. Due to impacts, there is a slight offset in the response (seen best in the right frames, emphasizing the ring width variations), leading to negative (positive) torque exerted on the ring by the outer (inner) satellite. Also shown in the right frames are the regions where local shear is reversed: blue regions indicate where T_rt =< ∆v_r∆v_t > is negative, ∆v_r and ∆v_t denoting the particles’ radial and tangential velocities compared to the mean flow velocity at its location.

Fig. C.2 Time evolution of a ring perturbed by a large mass anomaly µ = 0.1, to be compared with Fig. 6. The ring is initially centered around 1/3 resonance distance a0 with a full width of 0.6. This is slightly larger than the theoretical maximum of epicyclic excursions, a0epeak ≈ 0.27 and over a factor of 10 larger than the width of the resonance zone, Wres ≈ 0.026. Particle size R = 0.001 corresponds to 200 m for Chariklo’s rings. The three upper panels show the self-intersecting streamline forced by the 1/3 SOR. This self-crossing disappears around T = 6 000 to form a streamline dominated by a m = 1 azimuthal mode. In contrast to the case of narrower initial rings (Fig. 6), two supplementary rings form both interior and exterior to the main ring at T = 6 000–8 000. However, the outer ringlet is rapidly swallowed by the main ringlet (at T ≈ 10 000), while the inner ringlet survives until the end of the simulation. The last panel on the lower right (red color palette) shows the time evolution of the Lz distribution until the end of the simulation at T = 30 000.

Fig. C.3 Effect of optical depth on the resonance accumulation. The time evolution of the Lz distribution is shown with various initial optical depths τ0 and a fixed particle radius R = 0.001 (corresponds to 200 m for Chariklo’s rings). Upper row: Mass anomaly of µ = 0.05; lower row: µ = 0.03. Because R is fixed, different values of τ0 correspond to different numbers of particles. The orange curve shows the evolution of the optical depth at the resonance zone, calculated from the number of particles with |a − a0| < Wres (axis labels on the right), while the green curve indicates the excess radial velocity dispersion at the same region (full scale corresponds to a five-fold dispersion compared to unperturbed dispersion). The dashed lines indicate the estimated τlim: for τres < τlim the excitation of resonance perturbations (indicated by the growing oscillations of cres) starts, either initially in the case τ0 ≲ tlim, or after a sufficient drop in τres caused by the initial viscous diffusion. The latter cases are indicated by green arrows and provide an estimate of τlim. Based on the simulations τlim ≈ 0.06 and 0.025 for µ = 0.05 and 0.03, respectively.

Fig. C.4 Forced response of a collisional ring to first-order Lindblad resonances: 2/1 inner resonance (m = 2, j = 1) in the upper row and 2/3 outer resonance (m = −2, j = 1) in the lower row. In both cases the ringlet is centered at unit distance, and is composed of 7 500 particles with radius 0.0005; the satellite has a mass µs = 2 × 10−4 relative to the spherical central body. In the frames 50 snapshots have been stacked in a coordinate system corotating with the satellite (the satellite location is shown by the white square). In the left frame, the particle deviations from the ringlet mean distance have been exaggerated by a factor of 10. In the right frames the ringlet width has been exaggerated by a factor of 20. In the right frames the blue (red) contours indicate regions where Trt is negative (positive). The line indicates the direction where the ring width is the smallest. Local shear reversal, corresponding to negative Trt, is seen in the ring regions where particles are approaching the width minimum.

All Tables

All Figures

Fig. 1 Comparison between the numerical and theoretical responses of test particles to various SORs. Numerical integrations followed the motion of 10 000 test particles initially distributed on circular orbits, during 10 000 revolutions of the central body. The maximum eccentricities emax reached by these particles are plotted in black as a function of ΩB/n = (a/acor)3/2, and are compared with the analytical estimates of Paper I (red or green curves). Upper panel: case of a spherical body with a mass anomaly µ = 10−4. Lower panel: case of a homogeneous triaxial ellipsoid with elongation ϵelon = 0.01 and oblateness f = 0. In the case of mass anomaly the strongest responses are at the outer first-order Lindblad resonances corresponding to commensurabilities n/ΩB = m/(m − 1), with m = −1, −2, −3… The inset in the upper panel is a zoomed-in image of the ΩB/n = 3 region, using a 100 times larger mass anomaly of µ = 0.01. The response to the second-order 1/3 resonance is now visible, and is compared to the analytical estimate in green. In the case of an elongated body, the π-symmetry of the perturbation imposes even values m = −2, −4, −6… for the first-order resonances. In this case the strongest second-order resonance (m = −2 and j = 2) is also visible, with the analytical response plotted in green. The inset in the lower panel is a zoomed-in image of the ΩB/n = 3 region. The fourth-order 2/6 resonance has no signature even when using a large value ϵelon = 0.20.

In the text

Fig. 2 Responses of test particles to the first-order 2/3 and second-order 1/3 SORs. Left column: symbols represent the maximum eccentricities emax reached by particles initially on circular orbits for three values of the mass anomaly µ, near the 2/3 SOR (upper panel) and the 1/3 SOR (lower panel). The gray dashed curves are the analytical estimates of emax displayed in Figs. 3 and 4 of Paper I. The points are plotted against the normalized distance to the resonance, Δa/Wres = (ā − a0)/Wres, where a0 is the radius of the circular orbit at exact resonance, ā is the modified semimajor axis (Eq. (2)), and Wres is the width of the resonance (Eq. (3)). Right column: same, but with the timescales of Tres necessary to reach the maximum eccentricities emax. This figure shows that our numerical integrations reproduce satisfactorily the calculated distribution of emax, as well as the µ-scaling expected from Eqs. (3) and (4).

In the text

	Fig. 3 Comparison of eccentricity evolution at the first-order 2/3 and second-order 1/3 SORs. The black curves follow the evolution of test particles released near exact resonance, while the red dashed lines display the theoretical prediction, i.e., a linear growth rate for the first-order SOR and an exponential growth rate for the second-order SOR.
In the text

Fig. 4 Comparison of the maximum eccentricity emax reached by particles at the second-order 1/3 SOR (blue) compared with the first-order 2/3 (red) and 1/2 (green) SORs. Three values µ = 0.1, 0.01, 0.001 are compared. For µ = 0.1, the theoretical emax associated with the 2/3 and 1/2 resonances at the location of 1/3 resonance are 30% and 8% of that due to 1/3 resonance, respectively. For other µ values the ratios at the 1/3 SOR location scale proportionally to µ1/2.

In the text

Fig. 5 Three cases showing the combined effects of collisions and the 1/3 SOR on the ring confinement. The left frames show the time evolution of the vertical angular momentum (Lz) distribution of the particles around the 1/3 SOR at Lz = 1.44. The magenta lines are the average value of Lz and the red lines show the RMS of the particle eccentricities; the full y-range of the frame corresponds to e = 0.1. The insets are polar plots of the particle positions, shown again in the right column in cartesian coordinates. For better viewing, in the cartesian projections the width of the ring around its center position at each azimuth and the deviations of this center position from the overall mean distance have both been expanded by a factor of five. We compare three cases: (i) A ring of collisionless test particles perturbed by a µ = 0.1 mass anomaly on an otherwise spherical central body (upper row); (ii) Colliding particles around a spherical central body without perturbation (µ = 0, middle row); and (iii) A ring of colliding particles with µ = 0.1 (lower row). All simulations use the same initial conditions with 30 000 particles placed initially in an annulus r = 2.02–2.14 straddling the 1/3 SOR at semimajor axis a = 2.08. In the collisional simulations the particle radius R = 10−3 corresponds to about 200 m for Chariklo’s ring particles and yields an initial geometric optical depth τ0 = 0.06. In the case of a mass anomaly, the perturbation is turned on linearly during the first 50 revolutions of the central body.

In the text

Fig. 6 Simulation with 40 000 ring particles of radius R = 1.25 × 10−4 (corresponding to 25 m for Chariklo’s ring particles) perturbed by a mass anomaly µ = 3 × 10−3. The time T is given in units of Chariklo’s rotation period (about 7 h, Table 1), so that the maximum time shown here (T = 256 280) corresponds to about 206 years. Upper nine panels: density maps of the particles in the modified semimajor-eccentricity (ā, e) space (see Eq. (2)). The gray spiky curve is the value of emax shown in the lower left panel of Fig. 2. Three phases of the ring evolution are displayed. Phase I corresponds to a rapid damping of the initial eccentricities accompanied by a radial spreading caused by collisions; During Phase II, the 1/3 SOR excites the orbital eccentricities of the particles up to the predicted peak value epeak (Eq. (3)), while confining concomitantly the material near the resonance (ā ≈ 2.08) over a timescale Tres (Eq. (4)). Finally, during Phase III, the eccentricities damp due to dissipative collisions and a quasi-steady-state is reached, where the eccentricity damping is balanced by the resonance excitation. Lower three panels. Density maps of the particles in polar coordinates (L, r) space at three selected times shown in the upper panels. The white dotted lines indicate the location of the 1/3 SOR. At T = 163 120 a streamline forced by the 1/3 SOR has appeared, with a tendency of self-crossing around L = 260 deg. It is gathering material from the background unexcited particles. At T = 169 840, the accumulation process is going on, where a well-formed ringlet with various azimuthal modes collects more background material. At T = 257 280, all the ring material has been accumulated onto the ringlet, which is now dominated by a m = 1 azimuthal mode, with the presence of two kinks. The colors in the plots indicate the density of particles (in arbitrary units, density increasing from blue to red). Five movies generated from this simulation are available online. Movie1 illustrates the evolution in the (ā, e) space during Phase I, between T=0 and T=99 840 (corresponding to 0 and 79.77 years for Chariklo, respectively), using snapshots stored every 240 Chariklo rotations (about 70 days). Movie2 and Movie3 show the same for Phase II (between T=100 000 and T=199 840) and Phase III (between T=200 000 and T=257 280), respectively. The third movie also plots the position of an individual particle (white dot), showing its capture into the ring and its subsequent motion inside the cloud of points. Movie4 and Movie5 display the same snapshots in the longitude-radius space corotating with Chariklo, for Phases II and III, respectively.

In the text

Fig. 7 Upper panel: overview of the three phases shown in the upper panels of Fig. 6. The plot shows a stack of 1600 snapshots of the system from T = 120 000 to T = 200 000 (96 to 160 years respectively in the Chariklo case) with time steps ∆T = 50 Chariklo’s revolutions (about 15 days). The three phases I, II, and III described in the text are indicated by the arrows. Lower panel: phase III quasi-stationary stage obtained at the end of our run (we note the change in radial scale compared with the upper panel). A total of 780 snapshots taken from T = 256 500 to T = 257 280 (204.94 to 205.57 years) with time steps ∆T = 1 Chariklo’s revolution (about 7 h) have been stacked. A quasi steady state is reached, where the eccentricity damping due to collisions is balanced by the resonance excitation. On the long term (not reached here), a leakage of particles toward larger semimajor axes would be observed, as illustrated in Fig. 17, where a larger mass anomaly µ = 0.03 is used.

In the text

Fig. 8 Transition from noncollisional to collisional rings. This is explored by simulations with various initial optical depths τ0 from 0.00025 to 0.06, all using a mass anomaly µ = 0.1. Snapshots of the particle distributions in polar coordinates are shown at various times, using the radial range 1.88–2.28 centered around the resonant semimajor axis a1/3 ≈ 2.08. The label Tc in the second column indicates the average time interval between impacts, in units of the particle orbital periods. The rightmost column shows the time evolution of the vertical angular momentum Lz distribution with time up to 30 000 central body revolutions. The lowermost simulation (τ0 = 0.06) is the same as the one shown in the last row of Fig. 5, where it is displayed up to T = 5000.

In the text

Fig. 9 Transition between accumulation and dispersal, a test of Eq. (10). Left panel: same as the left column of Fig. 5, but exploring a grid of simulations near the 1/3 SOR. The particle radii increase from left to right in the range R = 1.25 × 10−4 − 4 × 10−3, corresponding to R = 25-0800 m in the case of Chariklo’s rings, while the mass anomaly increases upward from 0.003 to 0.1. The number of particles and the initial width of the rings were chosen in order to provide the same initial optical depth τ0 = 0.015 for all simulations. The vertical axes span values of Lz between 1.35 and 1.55 for µ = 0.1 and 0.03, and between 1.4 and 1.5 for µ = 0.01 and 0.003. The small insets covering the radial range 1.9–2.4 show snapshots of the ring in polar coordinates at the end of the run; the label indicates the duration of the run in Chariklo rotation periods. Increasing the particle size (and thus the viscosity) for a given perturbation strength prevents the resonance accumulation. Right panel: filled and open symbols respectively distinguish between simulations leading to resonance confinement and dispersal. The black line indicates the accumulation threshold for the simulations displayed in the left panel, following Eq. (10); the dashed portion of the line indicates the region where the linear relation fails since the resonance excitation at large µ ≳ 0.03 becomes somewhat weaker than predicted by the analytical formulas (see footnote 1). The shaded region bounded by the red line extrapolates the accumulation region to a denser ring with τ = 1.

In the text

Fig. 10 Transition from a |m| = 2 ringlet stage to a |m| = 1 ringlet stage in the µ = 0.03, R = 5 · 10−4 (corresponds to 100 m in the case of Chariklo) simulation of Fig. 9. In all other cases the ringlet shape is dominated by the |m| = 1 response once it has formed. Also in this case, the |m| = 1 response eventually replaces the |m| = 2 shape. During this transition at T ≈ 35 000, the torque exerted on the nearly circular ringlet temporarily vanishes (see the flat portion of the mean Lz curve in the upper frame).

In the text

Fig. 11 Fourier fits to the ring azimuthal profile in the µ = 0.003, R = 25 m simulation at T = 250 000 (same simulation as in Fig. 6). The uppermost frame shows r(L) profile, with a |m| = 1 fit superposed. The second frame shows the residual profile after subtracting the |m| = 1 component; the green curve is the |m| = 2 fit to the residual profile. The next frames repeat the procedure for |m| = 3 and |m| = 4.

In the text

Fig. 12 Proper Lindblad modes appearing in a ring confined at the 1/3 SOR. The outputs of the run shown in Fig. 6 have been used to detect normal modes in the confined rings between times 250 000 and 250 200, with steps ∆T = 0.05. The Lomb normalized periodogram power spectra of the modes described by Eq. (13) are displayed in the upper six panels as a function of ωm/ΩB for |m| = 1 to 6. Consistent with Eq. (15), the maximum power (red squares) is reached at the frequency ωm = (1 + |m|)ΩB/3. This is illustrated in the lower left panel. The lower middle panel shows the rapid decrease in the amplitude Am of the modes with |m|, while the right panel shows the distribution of phases ϕm.

In the text

Fig. 13 Left column: time evolution of the ringlet at the end of the µ = 0.003 simulation. The dominant |m| = 1, 2, 3 modes have been subtracted in order to highlight the kink feature. Right column: toy model for the formation of a kink as a superposition of |m| = 4–20 modes. Each mode has the same phase ϕm while the amplitudes obey Am ∝ 1/|m|3/2. Profiles calculated with Eq. (13) are shown over one full ring period.

In the text

Fig. 14 Local shear reversal in 1/3 SOR ringlet. The evolution of the µ = 0.003, R = 25 m simulation is shown over one ring orbital period (or three central body revolutions). The left frames display the r(L) profiles, with radial range 2.05–2.12. The ring patches with negative angular momentum flux Trt are indicated in light blue, while the white bullet indicates the location of a tracer particle. The middle frames show cartesian projections, where the ring deviations from its mean distance (dashed curve) have been exaggerated by a factor of 40 at each azimuth, and the open square marks the location of the mass anomaly. Similarly, in the right frames (again a cartesian projection) the ring width variations have been exaggerated by a factor of 50, using the same color convention as in the left frames. For the construction of all the figures, 20 snapshots of the ring separated by ∆T = 3, 6, …57 have been stacked.

In the text

Fig. 15 Time evolution of the vertical angular momentum Lz of particles that cross various resonances. The simulation contains 7500 particles with radius R = 0.001 (corresponding to 200 m in case of Chariklo), perturbed by a mass anomaly µ = 0.003. The ring starts near the 5/7 SOR and is pushed outward by a positive torque which takes it through the 2/3, 3/5, and 1/2 SORs. The red curve is the mean value of Lz. The green curve shows the mean eccentricity of the particles.

In the text

Fig. 16 Eccentricities e of particles vs. semimajor axis a. Unlike in Fig. 7, here we do not use the quantity ā because its definition depends on the particular resonance considered (Eq. (2)). However, because e remains small, this choice does not change the interpretation of this figure taken from the run shown in Fig. 15. We note that at time T = 21 600, the narrow second-order 3/5 SOR is able to significantly perturb the ring, in spite of the fact that the first-order 2/3 SOR is still expected to exert a nonnegligible perturbation, amounting to 10% of the eccentricity forced by the 3/5 SOR. At time T = 112 000, we see the damping effect of collisions as the ring evolves between the 3/5 and 1/2 SORs. The run is stopped at T=249 720 Chariklo’s rotations, corresponding to about 200 years. We note that third-order resonances, for example the 4/7 SOR at a=1.45, have no effect on the migration. Two movies generated from snapshots stored every 360 Chariklo rotations (about 105 days) are available online: Movie6 displays the time evolution of eccentricity versus semi-major axis and Movie7 shows the ring migration in a cartesian frame rotating with Chariklo.

In the text

Fig. 17 Upper row: long-term evolution of a 1/3 SOR ringlet perturbed by a mass anomaly µ = 0.03, using 30 000 particles with radius R = 5 × 10−4 (corresponds to 100 m for Chariklo rings). Although the ringlet maintains a well-defined core, a slow outward leakage of particles appears, forming a faint tail. For better viewing, the width of the ring has been expanded by a factor of 10. Lower left panel: effect of a satellite located at 2.3 ≈ 1.1a1/3, with a mass µs = 2 × 10−6 relative to the central body. Without the 1/3 SOR, the ring spreads inward under the effects of spiral density perturbations forced by Lindblad resonances, shown here at T=100 000. Lower right panel: when both the satellite and the 1/3 SOR are present, the outward leakage of particles is prevented by the 8/7 Lindblad resonance with the satellite. Two movies generated from simulation snapshots stored every 360 Chariklo rotations (about 105 days) are available online. Movie8 displays the evolution of eccentricity versus semi-major axis in the simulation without an external satellite: after T ∼ 60 000 (about 50 years for Chariklo), the particles start to leak outwards due to a small torque associated with the 1/3 SOR. Movie9 shows the same for the simulation with the moonlet: the leakage observed when only the 1/3 SOR is acting on the ring is prevented, and the ring ends up being confined between the 1/3 SOR and the 8/7 inner Lindblad resonance with the satellite.

In the text

Fig. 18 Comparison of the semimajor axis distributions in the three simulations of Figs. 17 and 19. The upper frame is the run with mass anomaly, and the colors indicate profiles at various times T = 20 000 to 200 000. The dashed gray line is the initial distribution, and the arrow marks the 1/3 SOR location. We note the growing tail of the distribution in profiles for T ≳ 100 000. In the middle frame, the simulation includes both the mass anomaly and a satellite. The run extends to T = 150 000 revolutions, with no leakage or spreading of the ringlet. The bottom frame shows the simulation with a satellite alone, with arrows marking its inner Lindblad resonances. We note the inward spreading of the ring and the gradual accumulation of material to the m = 5 resonance, and the disappearance of the m = 9 peak.

In the text

	Fig. 19 Evolution of mean angular momentum of the 1/3 SOR ring perturbed by a µ = 0.03 mass anomaly (black curve), by a µs = 2 × 10−6 satellite (green), or simultaneously by both (red). The snapshots of these simulations at T = 100 000 are displayed in Fig. 17.
In the text

Fig. 20 Self-gravitating simulations with R = 0.005 (corresponding to a physical radius of 1 km in the Chariklo case), assuming particles with bulk density ρ = 300 kg m−3. The initial distribution extends to r = 1.55–2.8, which corresponds to the range rh = 0.75–1.35. With N = 24 000 particles the initial τ0 = 0.11. The upper row shows two snapshots of the simulation with no mass anomaly, while the lower row shows a simulation with a µ = 0.3 mass anomaly. The location of the 1/3 SOR (a = 2.08) corresponds to rh ≈ 1.0. In the T = 20 frames, the tick marks on the right vertical axis indicate the radii where rh = 1.15. The red and green solid lines in the µ = 0.3 frames indicate the locations of 1/2 (left) and 1/3 (right) SOR. The insets in the T = 150 frames zoom in on the largest aggregate formed during the run, covering a 0.4 × 0.4 region with a correct aspect ration. The slice of the aggregate through the z = 0 plane is shown.

In the text

Fig. 21 Self-gravitating simulations of N=30 000 particles with different bulk densities. Simulations continue from the ringlet stage (at T = 10 000) of the nongravitating simulation of Fig. 9 with µ = 0.1, R = 0.0005 (100 m in Chariklo’s case): the snapshot of the initial state is shown in the upper left frame. The other frames display self-gravitating simulations with ρ= 300, 450, and 900 kg m−3, corresponding to rh = 1.0, 1.15, and 1.44, respectively at the distance of the 1/3 SOR. Times are counted from the beginning of the self-gravitating simulation. The upper right frame displays the evolution of the ringlet mean width, $W=\sqrt{12\langle {\left(r-r_{\text{fit}}\right)}^2\rangle }$, where rfit denotes the fitted ringlet central line at the azimuth of the particle and the averaging is over all particles (for a uniform particle distribution this formula with the factor of 12 recovers the full width of the ringlet). The remaining frames on the right show the semimajor axis distribution at different instants for the three self-gravitating runs.

In the text

Fig. 22 Self-gravitating simulations with µ = 0.1, R = 0.0005 (100 m in Chariklo’s case), assuming particles with bulk density ρ = 450 kg m−3. The initial particle distribution is similar to that in Fig. 21, except that τ has been reduced by a factor of 2 by selecting only every other particle from the confined ringlet in the nongravitating simulation at T = 10 000. The upper frame displays the evolution of the ringlet mean width, together with snapshot corresponding to various width maxima and minima; the snapshot at T = 6300 displays the system just after tidal breakup of the largest aggregate. The lower frame displays the number of particles in the largest aggregate (red) and the total number of particles in groups with more than 100 (blue) and 10 particles (gray; the total number of particles N = 15 000). Also shown as a green curve is the vertical dispersion cz (labels on the right).

In the text

	Fig. C.1 Four collisional simulations of 1/3 SOR with µ = 0.1, R = 0.001 (corresponds to 200 m for Chariklo’s rings). Upper row, left: Standard simulation using full mass anomaly potential. Upper row, right: Using an expansion of the potential, including m=1 and m=2 terms. Lower row: Using the m = 1 (left) or m = 2 (right) terms of the potential expansion.
In the text

Fig. C.2 Time evolution of a ring perturbed by a large mass anomaly µ = 0.1, to be compared with Fig. 6. The ring is initially centered around 1/3 resonance distance a0 with a full width of 0.6. This is slightly larger than the theoretical maximum of epicyclic excursions, a0epeak ≈ 0.27 and over a factor of 10 larger than the width of the resonance zone, Wres ≈ 0.026. Particle size R = 0.001 corresponds to 200 m for Chariklo’s rings. The three upper panels show the self-intersecting streamline forced by the 1/3 SOR. This self-crossing disappears around T = 6 000 to form a streamline dominated by a m = 1 azimuthal mode. In contrast to the case of narrower initial rings (Fig. 6), two supplementary rings form both interior and exterior to the main ring at T = 6 000–8 000. However, the outer ringlet is rapidly swallowed by the main ringlet (at T ≈ 10 000), while the inner ringlet survives until the end of the simulation. The last panel on the lower right (red color palette) shows the time evolution of the Lz distribution until the end of the simulation at T = 30 000.

In the text

Fig. C.3 Effect of optical depth on the resonance accumulation. The time evolution of the Lz distribution is shown with various initial optical depths τ0 and a fixed particle radius R = 0.001 (corresponds to 200 m for Chariklo’s rings). Upper row: Mass anomaly of µ = 0.05; lower row: µ = 0.03. Because R is fixed, different values of τ0 correspond to different numbers of particles. The orange curve shows the evolution of the optical depth at the resonance zone, calculated from the number of particles with |a − a0| < Wres (axis labels on the right), while the green curve indicates the excess radial velocity dispersion at the same region (full scale corresponds to a five-fold dispersion compared to unperturbed dispersion). The dashed lines indicate the estimated τlim: for τres < τlim the excitation of resonance perturbations (indicated by the growing oscillations of cres) starts, either initially in the case τ0 ≲ tlim, or after a sufficient drop in τres caused by the initial viscous diffusion. The latter cases are indicated by green arrows and provide an estimate of τlim. Based on the simulations τlim ≈ 0.06 and 0.025 for µ = 0.05 and 0.03, respectively.

In the text

Fig. C.4 Forced response of a collisional ring to first-order Lindblad resonances: 2/1 inner resonance (m = 2, j = 1) in the upper row and 2/3 outer resonance (m = −2, j = 1) in the lower row. In both cases the ringlet is centered at unit distance, and is composed of 7 500 particles with radius 0.0005; the satellite has a mass µs = 2 × 10−4 relative to the spherical central body. In the frames 50 snapshots have been stacked in a coordinate system corotating with the satellite (the satellite location is shown by the white square). In the left frame, the particle deviations from the ringlet mean distance have been exaggerated by a factor of 10. In the right frames the ringlet width has been exaggerated by a factor of 20. In the right frames the blue (red) contours indicate regions where Trt is negative (positive). The line indicates the direction where the ring width is the smallest. Local shear reversal, corresponding to negative Trt, is seen in the ring regions where particles are approaching the width minimum.

In the text