Exploring Shear-free Ringlet Formation with Direct Simulations of Saturn's B Rings

Following Wisdom & Tremaine (1988; also see Perrine et al. 2011; Ballouz et al. 2017), we perform local simulations of ring particles in a comoving patch frame. The equations of motion are solved in the Hill approximation, a corotating local coordinate system (Hill 1878). The equations of motion for a test particle are given by

$$\begin{eqnarray}&&\ddot{x}={F}_{x}+3{{\rm{\Omega }}}^{2}x+2{\rm{\Omega }}\dot{y},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\ddot{y}={F}_{y}-2{\rm{\Omega }}\dot{x},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\ddot{z}={F}_{z}-{{\rm{\Omega }}}^{2}z,\end{eqnarray} \tag{ 3 }$$

where ${\boldsymbol{F}}=({F}_{x}$, F_y, F_z) is the acceleration due to particle self-gravity and contact forces, x, y, and z are the coordinates of the particle in the local coordinate system (whose origin is located at the center of the patch; positive x is in the radial direction away from Saturn, positive y is in the orbital direction, and positive z points out of the ring plane according to the right-hand rule), and the derivatives are with respect to time. These equations of motion allow us to use shear-periodic boundary conditions in the plane of the patch.

The collisional component of pkdgrav incorporates a soft-sphere discrete element method (SSDEM; Schwartz et al. 2012). In SSDEM, spherical particles are allowed to overlap slightly, which is meant to represent deformation at the point of contact. Particle contacts can last many time steps, with reaction forces dependent on the degree of overlap and contact history. By allowing particles to overlap, multicontact and frictional forces can be modeled. Our implementation of SSDEM uses a spring/dashpot mechanism for the contact forces. In this model, two overlapping particles will apply a normal contact force, ${{\boldsymbol{F}}}_{N},$ and a tangential stick-slip force, ${{\boldsymbol{F}}}_{T}.$ In the current implementation, contact forces are determined by restoring forces with optional damping and/or friction (static, rolling, and/or twisting):

$$\begin{eqnarray}{{\boldsymbol{F}}}_{N} & = & -{k}_{N}x\hat{{\boldsymbol{n}}}+{C}_{N}{{\boldsymbol{u}}}_{N},\\ {{\boldsymbol{F}}}_{T} & = & \min ({k}_{T}{{\boldsymbol{\delta }}}_{T}+{C}_{T}{{\boldsymbol{u}}}_{T},{\mu }_{S}| {{\boldsymbol{F}}}_{N}| {{\boldsymbol{\delta }}}_{T}/| {{\boldsymbol{\delta }}}_{T}| ),\end{eqnarray} \tag{ 4 }$$

where k_N and k_T are the normal and tangential spring constants, respectively. Here the min() function indicates taking the vector quantity of smallest magnitude. The damping parameters (C_N, C_T) are related to the conventional normal and tangential coefficients of restitution used in hard-sphere implementations, ε_n and ε_t (see Schwartz et al. 2012 for further details). x is the mutual overlap of the two particles in length units, ${{\boldsymbol{\delta }}}_{T}$ is the sliding displacement from the equilibrium contact point, and $\hat{{\boldsymbol{n}}}$ is the unit vector from the particle's center to its neighbor. ${{\boldsymbol{u}}}_{N}$ and ${{\boldsymbol{u}}}_{T}$ are the normal and tangential components of the total relative velocity between the two particles at the point of contact, respectively. The coefficient of static friction, μ_s, determines the maximum amount of tangential force that can be supported by the contact point.

Twisting resistance arises from the slip and friction at the contact region due to a difference in the rotation rate of the particles in a direction along the normal vector $\hat{{\boldsymbol{n}}}$. Rolling resistance can come about from a number of sources including slip and friction on the contact surface and shape effects. Both resistances are modeled as functions of the relative twisting or rotational angular velocity that can reach critical values based on a shape parameter, β, after which they maintain a maximum resisting torque (see Section 2.2 of Zhang et al. 2017). β represents a statistical measure of real particle shape, which has been found to be one of the primary physical mechanisms for rotational resistance that also plays an important role in modeling interparticle cohesion. See Zhang et al. (2017) for the full equations describing the torque contributions to the contact forces.

Previous studies (Iwashita & Oda 1998; Mohamed & Gutierrez 2010) have shown that the inclusion of rotational resistances into a contact model allows numerical simulations to better match laboratory experiments. Our previous study (Ballouz et al. 2017) showed that rolling friction is an important micro-scale feature that can control whether kilometer-scale self-gravity wakes or vertically coherent viscous overstability structure develops in the ring. In that study, we used a rotational resistance implementation (Schwartz et al. 2012) that is suited for the modeling of the dynamic flow of granular material. For this work, we anticipate that the development of shear-free zones in our simulations may be better modeled with a rotational resistance implementation developed in Zhang et al. (2017) that is better suited for grains in a quasi-static state. This is because the new method uses a spring-dashpot model for rolling resistance that keeps track of the original point of contact and is able to account for static rolling friction, which is more appropriate for dense particle systems. There are some models that can accurately capture the rotational resistance properties of grains in a mixture of the two states (dynamic flow and quasi-static aggregates, Ai et al. 2011); however, a full-scale comparison of rotational resistance models is beyond the scope of this work.

In the quasi-static-state implementation we use here, two rotational resistances are modeled: twisting and rolling resistance (parameterized through the friction coefficients μ_R and μ_T). The SSDEM implementation of pkdgrav has been validated through comparison with laboratory experiments (Schwartz et al. 2014) and has been used for various solar system applications: size-sorting on asteroids (Maurel et al. 2017), avalanche dynamics (Yu et al. 2014), the stability of asteroid rotations (Zhang et al. 2017, 2018), and collisions between rubble piles at low speeds (Ballouz et al. 2014, 2015).

In this paper, we introduce the use of an interparticle cohesive force between particles in a ring-patch simulation due to Zhang et al. (2018), who applied it to the study of critical spin limits of asteroids. The cohesive forces between ring particles may arise from van der Waals forces due to molecular or atomic polarization effects between micro-sized particles (Scheeres et al. 2010). Zhang et al. (2018) implemented a cohesion model in pkdgrav based on hypothetical contact between microscopic grains characterized by the shape parameter β (introduced in the end of Section 2.1). In this model, β represents a statistical measure of the area where the interstitial grains are in contact. The effective contact area, A_eff, between two grains is given by

$$\begin{eqnarray}&&{A}_{\mathrm{eff}}=4{(\beta R)}^{2},\end{eqnarray} \tag{ 5 }$$

where R is the effective particle radius, R = R_pR_n/(R_p + R_n), and R_p and R_n are the radii of the particle and its neighbor, respectively. The cohesive contact force, ${{\boldsymbol{F}}}_{C},$ between the two particles is given by

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{C}={{cA}}_{\mathrm{eff}}\hat{{\boldsymbol{n}}},\end{eqnarray} \tag{ 6 }$$

where c is the interparticle cohesive tensile strength. The value of c reflects the physical properties of the interstitial microscopic grains: porosity, size distribution, surface energy, etc. The two particles only feel a cohesive force when they are in contact.

The value of the cohesive tensile strength for planetary ring particles is not well constrained in the literature. As a result, we tested a wide range of cohesion values in order to see its effect. For this study, after experimenting with cohesion coefficient values (c) from 0 to 10 Pa, we settled on the range 5 × 10⁻² to 7 Pa, as smaller values were indistinguishable from using zero and larger values resulted in coherent structures on the order of the patch size, which invalidate the sliding patch model constraints. Note that this range is consistent with the values of a few to 100 dyn found by Hatzes et al. (1991) in their ice ball experiments.

For the simulations presented in this paper, we set ${k}_{t}=\tfrac{2}{7}$ k_n in Equation (4) to keep the damped harmonic frequencies of the normal and tangential springs in phase—see Schwartz et al. (2012). An appropriate value of k_n is given by

$$\begin{eqnarray}&&{k}_{n}\sim m{\left(\displaystyle \frac{{v}_{\max }}{{x}_{\max }}\right)}^{2},\end{eqnarray} \tag{ 7 }$$

where m is the mass of a typical energetic particle, v_max is the maximum expected relative particle speed in the simulation, and x_max is the maximum allowed particle overlap, which we set to be ∼1% of the smallest particle radius. The timestep is set so that a typical collision (i.e., a full oscillation of the normal spring) takes about 30 steps. The typical collision speeds, v_coll, in a shearing patch frame are given by the Keplerian orbital frequency, Ω, and the typical radius, R_particle, of the ring particles:

$$\begin{eqnarray}&&{v}_{\mathrm{coll}}\sim {\rm{\Omega }}{R}_{\mathrm{particle}}.\end{eqnarray} \tag{ 8 }$$

In order to be conservative with our estimation of typical collision speeds, we set v_max to be an order of magnitude larger than v_coll. For a ring patch at a Keplerian orbit 100,000 km from Saturn (B Ring) and ring particles with R_particle = 1 m, we use k_n = 4.2 × 10³ kg s⁻², and a timestep of 0.074 s.

The friction properties of ring particles can influence the types of structures that are able to develop in opaque ring patches (Ballouz et al. 2017); however, they are poorly constrained. Here, we use a set of parameters, μ_s = 1.0, μ_R = 1.05, μ_T = 1.3, β = 0.5, as nominal values of the friction coefficients in order to mimic the interaction of sand particles of medium hardness that have an angle of repose of ∼30° (Jiang et al. 2015; Zhang et al. 2017).

We simulate patches located in the B ring of dynamical optical depths of 0.8, 1.4, and 1.8 and surface densities between 530 and 2400 kg m⁻². The dynamical optical depth τ_dyn of a ring is defined as

$$\begin{eqnarray}&&{\tau }_{\mathrm{dyn}}=\int \pi {R}_{\mathrm{particle}}^{2}n({R}_{\mathrm{particle}}){\rm{d}}{R}_{\mathrm{particle}},\end{eqnarray} \tag{ 9 }$$

where R_particle is the particle radius and n(R_particle) is the surface number density of ring particles whose sizes are between R_particle and R_particle + dR_particle. For equal-size particles, this reduces to ${\tau }_{\mathrm{dyn}}=N\pi {R}_{\mathrm{particle}}^{2}/S$, where N is the number of particles and S is the area of the patch. We consider the middle optical depth value to be nominal in our simulations as it is typical of the B ring, and it corresponds to a number of particles that can be feasibly simulated in a timely fashion using our code.

The values of the patch surface density in our simulations cover the range found in Saturn's B ring (400–1400 kg m⁻²; Hedman & Nicholson 2016). In order to ensure that we allow sufficient space for the formation of large-scale structures, we simulate a patch that is at least 4 λ_cr in the radial direction and 4 λ_c in the azimuthal direction, where λ_c is the Toomre wavelength, defined as:

$$\begin{eqnarray}&&{\lambda }_{c}=\displaystyle \frac{4{\pi }^{2}G\sigma }{{{\rm{\Omega }}}^{2}},\end{eqnarray} \tag{ 10 }$$

where G is the gravitational constant and σ is the surface density. Previous simulations (Salo 1992; Daisaka & Ida 1999) have shown that the typical spacing of wakes is close to λ_cr.

For each simulation, we choose a single combination of particle density and particle radius (monodisperse distribution) that corresponds to a constant optical depth. This is done in order to ensure λ_c is kept constant across all simulations that use the same dynamical optical depth, allowing us to draw more accurate comparisons of simulation outcomes across different particle density and cohesion cases. The main parameters for the simulations are summarized in Table 1.

Typical run times were between 0.3 and 40 kSU, where one SU is equal to 1 hr of wall clock time on a single core of a CPU, and 1000 SU = 1 kSU.

In order to determine the effectiveness of cohesive forces to clumping particles together, we developed a procedure to measure the sizes of particle aggregates in the patch. While aggregates may easily be distinguished by the human eye, naiive computational procedures are often fooled by tenuous connections or overlaps that lead to overestimates of aggregate sizes. Our approach divides the patch into grid cells, calculates the projected particle number in each cell (particle number density), extends search branches from the highest-density cell, and checks along each branch to determine the aggregate semimajor axis, semiminor axis, ellipticity, and orientation. Since self-gravity tends to keep particles in the midplane, most of the aggregates formed can be analyzed as two-dimensional structures in the x–y plane. As a result, it is not necessary to analyze the aggregate size in 3D, which saves computational time. We use 20 radial branches for each search. For expediency and accuracy, the branches' extent is limited to one-fourth of the patch from the highest particle density cell. We found this to be sufficient since aggregates are smaller than one-fourth of the patch size; furthermore, coherent structures of the order of the patch size invalidate the assumptions of the sliding patch model. However, the clump-finding algorithm adopted here only uses the distribution of surface number density of particles and does not check if the identified aggregate is actually forming a shear-free structure or not.

Branches are drawn in 18° intervals for symmetry. Cells are sampled along each branch until the projected number of particles in the cell drops below a threshold fraction of the maximum central projected particle number. If a branch passes partly between two cells, the average particle number in each neighboring cell is used. The threshold fraction can be adjusted for different cases since low-cohesion simulations are more likely to form wakes instead of aggregates. For example, if a smaller threshold is used in low-cohesion analysis, the aggregate semimajor axis is likely to be identified as the whole vertical wake structure. Figure 1 shows results of the analysis using a threshold of 0.1 (left) compared to a threshold of 0.5 (right) for a cohesion value of 6.6 × 10⁻³ Pa. Because the particle per cell value decreases rather slowly along the gravity wake structures, instead of converging on an aggregate shape, the search continues along the wake. Ultimately, a threshold fraction value of 0.5 was chosen through trial and error based on visual inspection. The semimajor axis of the structure is then calculated as half of the maximum extent achieved along two oppositely oriented branches.

Zoom In Zoom Out Reset image size

In order to ensure the algorithm is not sensitive to patch size, we performed the analysis for two different patch sizes, as shown in Figure 2. Three different snapshots taken over three orbits after equilibrium is established were used to perform the average. Error bars show the 1σ standard deviation. Three different cohesive strengths were analyzed for each patch size. For consistency, the three simulation time steps were chosen to be the same for both patch sizes. Since the results are fairly similar in each case for the two patch sizes, we conclude that the smaller patch size is sufficient for characterizing the results of these simulations.

Zoom In Zoom Out Reset image size

We find that particles in Saturn's opaque rings have a nonintuitive relationship between their cohesive strengths and the maximum structure size that they can form. It was expected that ring particles were more likely to form aggregates as interparticle cohesion increases; however, we find that large values of cohesion frustrates aggregate growth, and self-gravity wake structure is destroyed. Furthermore, we find that this process tends to cause a sharp increase in the z velocity dispersion of the ring.

From visual comparisons of the early stages of a low-cohesion simulation and a high-cohesion simulation, as shown in Figure 7, we find a possible explanation. We observe that ring particles with high cohesion form large aggregates extending in the radial direction, and these aggregates are then rotated by the tidal shear, collide, and break apart, injecting a massive amount of energy into the system. This may explain why the system has a "hotter" equilibrium state where the z velocity dispersion is high and no large aggregates are able to persist.

Zoom In Zoom Out Reset image size

The evolution of the total kinetic energy per unit mass for low and high cohesion values further supports this hypothesis (Figure 8). This quantity is calculated by summing up the translational kinetic energy ($\tfrac{1}{2}{{mv}}^{2}$) and the rotational kinetic energy ($\tfrac{1}{2}I{\omega }^{2}$) from each particle and dividing by the total mass of the system. Here the moment of inertia $I=\tfrac{2}{5}{{mR}}^{2}$ for a spherical particle of mass m and radius R. The analysis shows a systematic difference between the high-cohesion run and the low-cohesion run at an early stage. With high cohesion, the system gained a huge amount of kinetic energy, presumably through tidal shearing as explained in the previous paragraph, near the start of the simulation. As the large aggregates are broken down by collisions, the energy injection slows down and the energy dissipation from collisions rapidly increases. This causes a rapid decrease in the total kinetic energy until an equilibrium state is reached, where the energy injection from tidal shear is balanced by the energy dissipation. Therefore, the system reaches a more uniform "hotter" equilibrium, compared to the low-cohesion case, which gradually gains energy from the tidal shearing before approaching a slightly lower equilibrium value compared to the high-cohesion case.

Zoom In Zoom Out Reset image size

Visual inspection of the snapshots of the early stages of a simulation with high restitution coefficients compared to that of low restitution coefficients with the same cohesive strength reveals a negative relation between the radial structure size and the restitution coefficients (see Figure 9), since lower restitution coefficients correspond to a higher energy lost when particles collide, which potentially increases the structure size. Larger structures can carry more energy, and this could explain why simulations with low restitution coefficients result in an earlier transition from wakes to structureless form.

Zoom In Zoom Out Reset image size

The hypothesis that bigger aggregates form with higher energy loss can also be confirmed by analyzing the particle speed distribution. We plotted the magnitude of the velocities (v) for all particles (see Figure 10) using $v=\sqrt{{v}_{x}^{2}+{v}_{{y}_{\mathrm{noShear}}}^{2}+{v}_{z}^{2}}$, where v_x, ${v}_{{y}_{\mathrm{noShear}}}$, v_z are the particle velocity components in the x, y, z direction after subtracting the shear speed in the y direction. The result shows larger aggregates forming for the high-energy-loss runs at a cohesive strength value of 7 Pa. Although the boundary condition might have played a role in breaking up the particles due to their large radial extent, there are significant clumps undergoing high rotation well inside the patch. The same test also shows that low-cohesion runs at the same moment in time have no particles that have speeds greater than 0.5 cm s⁻¹ or any large aggregates forming.

Zoom In Zoom Out Reset image size

We also analyzed the radial velocity of particles (x direction). For rotating clumps, we will expect a smooth transition of x velocity from small to large as the particle distance from the center increases, which is what we see for the clumps identified in Figure 11 by the red circles. These two tests further support our hypothesis that shear-induced rotations of large clumps inject energy at high cohesion values.

Zoom In Zoom Out Reset image size

The same test for velocity in the x direction is performed on a larger patch size (10 × 10 λ_c) in order to ensure the boundary condition is not the main factor in breaking up large aggregates. Figure 12 shows the aggregates in the same color scale shown in Figure 10. The red solid line carved out one of the aggregates. Compare to that of a smaller patch size, the difference in aggregate size is visually invisible.

Zoom In Zoom Out Reset image size

If our explanation is correct, shear-free rings around Saturn might not be possible since any small perturbation could start a large aggregate rotating for sufficiently cohesive particles, resulting in high collision speeds and destruction of the aggregate.

While we have demonstrated that cohesion limits the growth of subkilometer structures (∼100 m) in the B-ring, the analysis of Tremaine (2003) had suggested that annuli of "solid" rings up to 100 km in size may form through interparticle cohesion if the tensile strength of the material was large enough. Here, we analyze the results of our simulation in light of these previous theoretical provisions.

By considering the shear rate of a solid annulus in Saturn's ring, Tremaine (2003) analytically demonstrated that the maximum tensile stress, C, that could be withstood by a rotating slab of width, Δr, and bulk density, ρ, that is located in the B ring (i.e., 10⁵ km from Saturn) can be given by

$$\begin{eqnarray}&&C=4\times {10}^{4}\left(\displaystyle \frac{\rho }{0.3\,{\rm{g}}\,{\mathrm{cm}}^{-3}}\right){\left(\displaystyle \frac{{10}^{10}\mathrm{cm}}{r}\right)}^{3}{\left(\displaystyle \frac{{\rm{\Delta }}r}{100\mathrm{km}}\right)}^{2}\,\mathrm{Pa}.\end{eqnarray} \tag{ 12 }$$

We can compare the results of our simulations to these theoretical expectations by equating C to the expected bulk tensile strengths that we simulated. Here, we note that the bulk tensile strength differs from the interparticle cohesive strength, c. It is difficult to determine a precise relationship between these two values as the bulk tensile strength of an aggregate will depend on the average number of grain interactions, grain properties such as shape and sizes, and the geometrical structure of the aggregate. Nevertheless, we estimate this relationship between c and C from simulations of spherical aggregate spin destabilization by Zhang et al. (2018), who demonstrate the following empirical relationship using the same cohesion model and friction parameters as those of this study,

$$\begin{eqnarray}&&c\simeq 120\times {C}^{0.75}.\end{eqnarray} \tag{ 13 }$$

By equating Equations (12) and (13), we evaluate the necessary interparticle cohesive strength for a given Δr, for a structure bulk density of 0.3 g cm⁻³ (Figure 13).

Zoom In Zoom Out Reset image size

We find that for the range of interparticle cohesion we explored here (black dashed lines), the theoretical maximum structure size ranges from approximately a few to 100 m in size. However, the analysis does not account for the formation of self-gravity wakes that create large-scale structure without the aid of cohesion. The typical sizes of these wakes are set by the ring surface density, given by Equation (10). If self-gravity is considered, as was done in our simulations, then the cohesion acts to create structures that are larger than wake structure, which will require tensile strengths that exceeded the limit set by Equation (12) (solid blue line) for stability. Therefore, it is likely that the combination of self-gravity and cohesion that leads to the unstable growth of structure. Nevertheless, there may be some conditions in which cohesion may build stable large-scale structures. The green dashed–dotted curves in Figure 13 represent two hypothetical scenarios for the growth of the structure strength with size for the case of a ring that starts with self-gravity wakes that are a few 100 m in size. For Case A, increases in structure size lead to growth in structure strengths that exceed the estimated growth given by Equation (12), while Case B represents a shallower growth in structure strength in size. If real ring particles can create structures that behave like Case A, then "solid" ring annuli should exist above a certain size (intersection between blue and green curve). For structure strength growth as represented by Case B, a solid ring will always be unstable against shear deformation, regardless of size. Measuring the dependence of structure strength with size and interparticle cohesive strength is beyond the scope of this study; however, this analysis provides a baseline for future theoretical investigations for large-scale structure stability through grain cohesion.