VISCOSITY IN PLANETARY RINGS WITH SPINNING SELF-GRAVITATING PARTICLES

Angular momentum transport is a key process that drives the large-scale structures of planetary rings. It determines the timescale of radial spreading, the behavior of wave structures, and the condition of gap formation in rings in the vicinity of embedded moonlets, and it is also related to the stability of dense rings (see, e.g., Goldreich & Tremaine 1978; Stewart et al. 1984; Schmidt et al. 2009). Furthermore, recent studies show that Saturn's small moons orbiting just outside the main rings likely formed by an accretion of radially spreading particles from the outer edge of the rings (Charnoz et al. 2010; Canup 2010); thus, the spreading rate of the ring is also related to the formation of these moons.

In planetary rings, angular momentum transport is caused by collisions and gravitational interactions between particles and satellites, and the rate of the transport is expressed in terms of viscosity. In most previous works on the viscosity of planetary rings, effects of the self-gravity of the rings were either neglected or only taken into account in an approximate manner (e.g., Goldreich & Tremaine 1978; Borderies et al. 1985; Araki & Tremaine 1986). Wisdom & Tremaine (1988) performed local simulation of planetary rings with nongravitating particles and studied their equilibrium properties, including viscosity. In addition to the local component of angular momentum transport due to particle radial random motion examined by Goldreich & Tremaine (1978), Wisdom & Tremaine (1988) also evaluated the nonlocal component due to collision. They found that the nonlocal component dominates the viscosity in dense rings, and they confirmed that the values of the viscosity agree well with the result of Araki & Tremaine (1986). They also examined effects of self-gravity in an approximate manner by enhancing the vertical frequency of particle orbital motion, and they found that the viscosity is enhanced by a factor of five. Richardson (1994) calculated ring viscosity using N-body simulation of self-gravitating rings, but did not examine the angular momentum transport due to gravitational interactions between particles.

On the other hand, Daisaka et al. (2001) studied viscosity in self-gravitating rings using local N-body simulation. In addition to the local and nonlocal components studied by Wisdom & Tremaine (1988), Daisaka et al. (2001) calculated the gravitational component of viscosity based on the formulation of Takeda & Ida (2001), who studied angular momentum transport in proto-lunar disks. Takeda & Ida (2001) and Daisaka et al. (2001) found that the viscosity is strongly enhanced when particle disks become gravitationally unstable and gravitational wakes are formed. From their numerical results, Daisaka et al. (2001) derived a formula of the viscosity in self-gravitating dense rings given by ν ∼ CG²Σ²/Ω³, where G is the gravitational constant, Σ and Ω are the surface density and angular velocity of the ring, respectively, and C is a correction factor that depends on the ratio of the Hill radius of the particles to their physical size (see Section 5). However, the range of parameters studied by Daisaka et al. (2001) was somewhat limited; for example, they did not examine viscosities in dense rings far from the planet, where temporary gravitational aggregates can form. Also, they assumed that particles are smooth spheres and did not examine the effects of surface friction of particles on the viscosity.

Spins of ring particles arise as a natural outcome of oblique impacts between particles with surface friction. Spin rate of ring particles is an important parameter in modeling thermal emission from Saturn's rings, and it has been studied in detail, using analytic calculation, three-body orbital integration, and N-body simulation (Salo 1987; Richardson 1994; Ohtsuki 2005, 2006a, 2006b; Ohtsuki & Toyama 2005; Morishima & Salo 2006). These studies show that the spin period of the particles is on the order of their orbital period in rings of equal-sized particles, while small particles spin faster than larger ones when size distribution is included. Considering energy balance between viscous gain and dissipation due to an inelastic collision with surface friction, Morishima & Salo (2006) examined the effects of particle surface friction on viscosity in rings consisting of nongravitating particles, but such effects on the viscosity in self-gravitating rings were not studied.

On the other hand, Tanaka et al. (2003) derived a new formulation for the calculation of the total viscosity (i.e., the sum of the local, nonlocal, and gravitational components) in particle disks. Under Hill's approximations in the three-body problem, they derived an expression of angular momentum flux in terms of the changes of orbital elements of constituent particles due to mutual collisions and gravitational interactions. In the case of planetary rings with low optical depth, they showed that viscosity can be evaluated by three-body orbital integration. They also derived a formula to calculate viscosity in dense rings using N-body simulation. In this case, viscosity in rings in a quasisteady state can be obtained from energy dissipation at inelastic collisions between particles (see Section 2). In the case of rings with the effect of vertical overall self-gravity, Salo et al. (2001) also found agreement between viscous heating and energy dissipation due to inelastic collision using N-body simulation. Tanaka et al. (2003) confirmed that viscosities in planetary rings calculated using their new formula agree with those obtained by the method used by Daisaka et al. (2001). This new method of viscosity calculation can also be used in the case of particles with surface friction, and the method is convenient when investigating detailed parameter dependence of the total viscosity of the rings.

In this work, we examine the viscosity of planetary rings consisting of spinning, self-gravitating particles. Using an N-body simulation code similar to the one developed by and used in Daisaka & Ida (1999) and Daisaka et al. (2001) and on the basis of the above viscosity calculation formula derived by Tanaka et al. (2003), we calculate ring viscosity for a wide range of parameters, including cases of particles with surface friction. We also examine viscosities in dense rings far from the central planet, where temporary gravitational aggregates can form. In Section 2, we describe our numerical methods. Numerical results of our N-body simulation are presented in Sections 3 and 4. First, in Section 3, we describe the effects of surface friction on viscosity in rings with low optical depth. Then, we show results for dense rings of spinning self-gravitating particles in Section 4. In Section 5, we derive semianalytic expressions of ring viscosity based on our numerical results. Our conclusions are summarized in Section 6.

We adopt the method of local N-body simulation (e.g., Wisdom & Tremaine 1988; Richardson 1994; Salo 1995; Daisaka & Ida 1999; Ohtsuki & Emori 2000; Daisaka et al. 2001) and use a code based on Daisaka & Ida (1999) and Daisaka et al. (2001). We erect a rotating Cartesian coordinate system with the origin at the center of the square simulation cell that moves on a circular orbit with a semimajor axis a₀ at the Keplerian angular velocity Ω = (GM_c/a³₀)^1/2 (M_c is the mass of the central planet). The x-axis points radially outward, the y-axis points in the direction of the orbital motion, and the z-axis is normal to the equatorial plane. In this case, the equations of motion of particle i are written as

$$\begin{eqnarray} &&\ddot{{x}}_{i} = 2 \dot{{y}}_{i}\Omega + 3x_{i}\Omega ^{2} + \sum ^N_{j \ne i}\frac{ G m_{j} (x_{j} - x_{i})}{r_{ij}^{3}} \nonumber \\ &&\ddot{{y}}_{i} = -2 \dot{{x}}_{i}\Omega + \sum ^N_{j \ne i}\frac{ G m_{j}(y_{j} - y_{i})}{r_{ij}^{3}} \\ &&\ddot{{z}}_{i} = - {z}_{i}\Omega ^{2} + \sum ^N_{j \ne i}\frac{ G m_{j} (z_{j} - z_{i})}{r_{ij}^{3}} \nonumber, \end{eqnarray} \tag{ 1 }$$

where r_ij = [(x_i − x_j)² + (y_i − y_j)² + (z_i − z_j)²]^1/2, N is the particle number, m_j is the mass of particle j, and the last terms on the right-hand side of Equation (1) denote gravitational forces between particles. We directly calculate the gravitational forces between particles using GRAPE-7, which is a special-purpose hardware for calculating gravitational forces, and orbits of particles are integrated with the second-order leap-frog method. In this work, we assume that all particles have identical sizes of 1 m radius.

When collisions between particles are detected, velocity changes are calculated based on the hard-sphere model, including surface friction (e.g., Salo 1987; Richardson 1994; Ohtsuki & Toyama 2005). Velocity changes are described in terms of the normal and tangential restitution coefficients, ε_n and ε_t, where 0 ⩽ ε_n ⩽ 1 and −1 ⩽ ε_t ⩽ 1. (Perfectly smooth spheres have ε_t = 1.) Let the normal and the tangential components of the relative velocity (v) of colliding particles to the tangent plane be v_n and v_t, respectively, and the normal and the tangential components of the relative velocity of the two contacting points at the time of impact be u_n and u_t, respectively. Then, the normal and the tangential components of the relative velocity of the two contacting points after impact are given as

$$\begin{equation} \begin{array}{l} \mbox{\boldmath $u$}_\mathrm{n}^{\prime } = \mbox{\boldmath $v$}_\mathrm{n}^{\prime } = -\varepsilon _\mathrm{n} \mbox{\boldmath $v$}_\mathrm{n}, \\ \mbox{\boldmath $u$}_\mathrm{t}^{\prime } = \varepsilon _\mathrm{t} \mbox{\boldmath $u$}_\mathrm{t}. \end{array} \end{equation} \tag{ 2 }$$

Using these relations, we calculate changes in velocities and rotation rates of the colliding particles (Appendix A). In this study, we assume that ε_t is constant and examine the dependence of viscosity on ε_t. As for the normal restitution coefficient, in addition to cases with constant ε_n's, we also examine cases with the velocity-dependent coefficient based on laboratory impact experiments (Bridges et al. 1984), which is given as

$$\begin{eqnarray} &&\varepsilon _{\mathrm{n}}(v)= {\min} \big\lbrace 0.32 (v/v_{c} ) ^{-0.234}, \quad 1\big\rbrace \end{eqnarray} \tag{ 3 }$$

where v_c = 1 cm s⁻¹.

We calculate ring viscosities based on the formulation derived by Tanaka et al. (2003). Deriving an expression for angular momentum flux in terms of changes in particle orbital elements and assuming a quasisteady state in a local ring region, Tanaka et al. (2003) showed that the averaged viscosity in the ring region can be obtained by calculating energy dissipation due to inelastic collisions using N-body simulation. In the case of N-body simulation using a square simulation cell of width L, the viscosity averaged over the cell and a time interval T is written as (Tanaka et al. 2003)

$$\begin{eqnarray} \nu &=&\frac{4}{9\Omega ^{2} \Sigma TL^{2}}\sum _{j}^{{\rm inside}}\sum ^{{\rm all}}_{\displaystyle \scriptstyle k \atop \displaystyle \scriptstyle b_{k}>b_{j}}\sum _{{\rm all} \hspace{2.84526pt} {\rm collisions}}\Delta E_{\mathrm{col},jk}, \end{eqnarray} \tag{ 4 }$$

where Σ is the averaged surface density and ΔE_{col, jk}(> 0) is an amount of energy dissipated at a collision between particle j and particle k. ΔE_{col, jk} can be calculated, for example, from the change of the orbital elements for the relative motion of particles due to a collision (Appendix B). The first summation is taken over particle j whose guiding center is inside the simulation cell, while the second summation is over particle k, which have a semimajor axis of the guiding center (b) larger than that of particle j. The above expression can also be obtained by assuming a balance between energy loss due to inelastic collisions and the viscous gain due to the shear motion (Tanaka et al. 2003; see also Salo et al. 2001).

In the case of nongravitating rings or self-gravitating rings with low optical depth, N ∼ 10³ or an even smaller number of particles is sufficient to examine such equilibrium properties, e.g., velocity dispersion. However, in the case of self-gravitating dense rings, a sufficiently large simulation cell (thus, a sufficiently large number of particles) should be used to accurately obtain various dynamical quantities (Salo 1995). Following Salo (1995), we set the size of the simulation cell so that L ⩾ 4λ_cr, where λ_cr = 4π²GΣ/Ω² is the critical wavelength for axisymmetric gravitational instability. We performed simulations for sufficiently long time periods to achieve a quasisteady state; these were typically at least 20–50T_K (T_K is the Keplerian orbital period) for dense rings (10T_K for τ = 1 and r_h = 1, where r_h is defined by Equation (5)), and still longer periods in the case of dilute rings. The time average for calculating viscosities from Equation (4) is carried out after a quasisteady state is achieved.

In this section, we first show results for nongravitating rings. Then we describe, in detail, results for self-gravitating rings with a low optical depth. Numerical results for self-gravitating dense rings will be presented in Section 4.

Figure 1 shows numerical results for the viscosity in the case of nongravitating rings as a function of optical depth. First, to confirm the validity of our code for N-body simulation and viscosity calculation, we compare our results with Wisdom & Tremaine (1988) by performing simulations for the same parameters, i.e., smooth particles (ε_t = 1) with ε_n = 0.5. In Figure 1(a), the crosses represent results from Wisdom & Tremaine (1988), and the line shows our results. We find excellent agreement between the two results. Figure 1(b) shows our numerical results for the velocity-dependent normal restitution coefficient given by Equation (3) with ε_t = 1, 0.5, and 0. We find that the dependence of viscosity on the value of ε_t is different between the low and high optical depth cases, and this difference is related to the behavior of particle random velocity. In the case of low optical depth, random velocity is somewhat reduced with the increasing strength of surface friction (Salo 1987; Ohtsuki & Toyama 2005), and the reduced radial excursion of particles lowers the rate of angular momentum transport and viscosity (Goldreich & Tremaine 1978). On the other hand, in the case of high optical depth, the dependence of random velocity on τ or ε_t becomes rather weak (Wisdom & Tremaine 1988). In this case, viscosity becomes larger with decreasing ε_t, because viscous heating increases to balance the additional energy dissipation due to surface friction.

Zoom In Zoom Out Reset image size

Next, we examine the case of self-gravitating rings with a low optical depth. The relative importance of self-gravity to mutual collision between particles depends on the relative size of the Hill radius of the particles to their physical size (Ohtsuki 1993, 1999; Salo 1995; Daisaka et al. 2001). In the case of rings of equal-sized particles (mass m, radius R, internal density ρ) at a distance a₀ from the planet, the ratio of the mutual Hill radius R_H(=  a₀(2m/3M_c)^1/3) of colliding particles to the sum of their radii is written as (Daisaka et al. 2001)

$$\begin{eqnarray} r_\mathrm{h} & =& R_\mathrm{H}/(2R)\nonumber\\ & = & {\displaystyle 0.82 \left(\frac{M_\mathrm{c}}{5.7 \times 10^{29} \, {\rm g}}\right)^{-1/3}\! \left(\frac{\rho }{0.9 \, {\rm g} \, {\rm cm}^{-3}}\right)^{1/3}\! \left(\frac{a_0}{10^{10} \, {\rm cm}}\right)}. \nonumber\\ \end{eqnarray} \tag{ 5 }$$

Figure 2 shows the plots of r_h as a function of the distance from Saturn for several values of ρ. When ρ = 0.9 g cm⁻³, r_h = 0.5, 0.71, and 1 (see, e.g., Figure 3) correspond to a/10¹⁰ cm = 0.61, 0.87, and 1.22, respectively. However, if ρ = 0.5 g cm⁻³, the above three values of r_h correspond to a/10¹⁰ cm = 0.74, 1.05, and 1.48, respectively. Although the actual physical properties of ring particles are poorly known, comparisons between the observations and dynamical studies suggest that it is likely they are underdense (e.g., Schmidt et al. 2009).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 3 shows the dependence of viscosity and random velocity on ε_n for low optical depth rings (τ = 0.01) with four different values of r_h, including the nongravitating case (r_h = 0). Figure 3(b) confirms the dependence of random velocity on ε_n and r_h (or the distance from the planet for a given density) obtained by previous studies (e.g., Salo 1995; Ohtsuki 1999). In the case of no or weak gravity, random velocity is determined by a mutual collision and the equilibrium random velocity increases with increasing ε_n, while the ε_n-dependence becomes weak in the case of larger r_h, where gravitational encounters between particles become important. More importantly, the random velocity significantly increases with increasing r_h, i.e., the increasing effect of gravitational encounters. The ε_n-dependence of viscosity shown in Figure 3(a) is similar to, but somewhat weaker than, that of random velocity. However, we find that the r_h-dependence of the viscosity is similar to that of random velocity, indicating that viscosity in self-gravitating rings with low optical depth can be roughly expressed in terms of random velocity, as in the case of the dilute nongravitating rings discussed above.

Figure 4 shows the dependence of viscosity and random velocity on ε_t in self-gravitating rings with low optical depth. We also show the results for the nongravitating case (r_h = 0) for comparison. The curves for the random velocity in the case of weak gravity (r_h = 0.5) are very close to those for the nongravitating case. In these cases, energy dissipation due to surface friction is most effective at ε_t = 0; thus, the random velocity becomes the smallest at ε_t ≃ 0 (Morishima & Salo 2006). Then, the reduced radial excursion of particles results in the reduction of viscosity (Figure 4(a)), as mentioned above. The behavior in the case of r_h = 0.71 is similar, although the ε_t-dependence of the viscosity seems rather weak. In the case of r_h = 1, on the other hand, gravitational encounters between particles play a major role in determining random velocity. As a result, the random velocity becomes almost independent of ε_t. In this case, impact velocity also becomes nearly independent of ε_t, and energy dissipation at the collision is the largest when surface friction is the strongest, i.e., at ε_t ≃ 0. Therefore, the viscosity in the case of r_h = 1 becomes larger when ε_t ≃ 0 by a factor of about two compared with the case of smooth particles (ε_t = 1). Here, we emphasize that both viscosity and random velocity increase with increasing r_h and, if we neglect the above rather weak ε_t-dependence, viscosity in self-gravitating rings with low optical depth seems to be expressed in terms of random velocity (Section 5).

Zoom In Zoom Out Reset image size

Next, we examine the case of dense, self-gravitating rings. First, we show results for smooth particles to compare with cases with surface friction shown later. The solid line with marks in Figure 5(a) shows the plots of the viscosity obtained by our N-body simulation as a function of optical depth (r_h = 0.82, ε_n = 0.5, ε_t = 1). In the case of low optical depth, viscosity is proportional to the optical depth, because angular momentum is transferred via successive binary collisions and gravitational encounters between particles, and the frequency of interactions is proportional to the particle surface number density. The coefficient of proportionality in the case of self-gravitating rings with low optical depth can be obtained from results of three-body orbital integration (Tanaka et al. 2003). In this case, equilibrium random velocity for given parameters (i.e., r_h, ε_n, and ε_t) is first calculated by solving the velocity evolution equation with the velocity stirring rates obtained from the three-body orbital integration (Ohtsuki 1999). Then, we evaluate the coefficient for the viscosity for this equilibrium velocity by also using the results of the three-body orbital integration, and the result is shown in Figure 5(a) with the dotted line. We find that the results obtained by the two different methods perfectly agree with each other in the case of low optical depth. On the other hand, the results of N-body simulation deviate from the three-body results for τ ≳ 0.2. This is because the viscosity is dominated by the angular momentum transfer related to the gravitational torques exerted by the wake structures formed in such dense, self-gravitating rings (Figure 6; Daisaka et al. 2001).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition to the smooth particles case, we show the results of our N-body simulation for self-gravitating rings of particles with surface friction in Figure 5(b). We find that the overall behavior in the case with surface friction is similar to the case with smooth particles. Figure 5(b) shows the numerical results for four combinations of ε_n (0.5 or ε_n(v) given by Equation (3)) and ε_t(1, 0.5, or 0). Figure 5(c) is the same as Figure 5(b), but the viscosity values are scaled by those in the case of ε_n = 0.5 and ε_t = 1, to facilitate a comparison among the four cases. We note that the viscosity values change systematically when the restitution coefficients are varied. In dense rings where gravitational wakes are formed, the balance between gravitational heating by wakes and inelastic dissipation leads to a quasisteady state with Q ≃ 1–2, where Q is the Toomre parameter, and the dependence of the particle random velocity (thus their impact velocity) on elastic properties is rather weak (Salo 1995). In this case, the average energy dissipation rate is larger when the effect of the surface friction is stronger, resulting in somewhat larger viscosity values (Figures 5(b) and (c)). We also note that the slope of the curve for the case with ε_n(v) in Figure 5(b) becomes somewhat steeper at τ ≳ 0.2. This is because the increased impact velocities at a larger optical depth results in the increase in energy dissipation and viscosity in the case with the velocity-dependent restitution coefficient. Figure 5(c) shows that varying restitution coefficients can change the viscosity in dense self-gravitating rings by a factor of about two.

Figure 7 shows the τ-dependence of viscosity for various values of r_h, including the nongravitating case (r_h = 0; dashed line), with surface friction and the velocity-dependent ε_n. As in the case shown in Figure 5, viscosity in dilute rings is proportional to the optical depth, while it significantly increases at large τ in the case with self-gravity. We confirmed that Q takes on values smaller than 2–3 when viscosities begin to deviate from the linear τ-dependence. The deviation begins at τ ∼ 0.1 in the case of strong gravity (r_h = 1). On the other hand, in the case of weak gravity (r_h = 0.5), the deviation is very small and the viscosity values are close to those for the nongravitating case, even at τ ≃ 1. In this case of weak self-gravity, the tendency of wake formation is not strong, even at τ ≃ 1 (Figure 8). The wake structures for τ = 1 become notable with increasing r_h; even temporary aggregates can form in the case of r_h = 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The τ-dependence of viscosity, often expressed in terms of β ≡ dln ν/dln τ, is related to the condition for the onset of viscous overstability in dense rings; Salo et al. (2001) found that β > β_cr ∼ 1 as a rough condition for the onset of overstability. In the high-τ regime where viscosity is enhanced by the effect of self-gravity, the slope of the curves shown in Figure 7 seems to change with increasing optical depth. For example, in the case of r_h = 0.82, viscosity increases more rapidly than the ν∝τ² dependence (i.e., β > 2) when the effect of gravitational wakes becomes important at τ ≃ 0.2–0.3, then β ≃ 2 at τ ≃ 0.5–1. Similar behavior can be seen in the cases of r_h = 0.71 and 1. Figure 14.8 (upper right panel) of Schmidt et al. (2009), drawn by Heikki Salo, shows the change of β in the case of r_h = 0.82 in more detail, where β increases from 1 to about 3 as τ increases from 0.1 to 0.5, then decreases to about 2 as τ further increases to ∼1. In the range of parameters shown here and studied in other works (e.g., Daisaka et al. 2001; Schmidt et al. 2009), we find that 1 ≲ β ≲ 3 in self-gravitating dense rings, satisfying the condition for the onset of overstability. Simulations of dense rings for much longer periods with a larger simulation cell are necessary for a more detailed study of viscous overstability (Daisaka et al. 2001; Salo et al. 2001).

On the basis of our numerical results presented above, here we derive a semianalytic expression for viscosity in rings of equal-sized particles. As we discussed above, ring viscosity depends on various parameters. We have shown that varying restitution coefficients can change ring viscosity for a given τ and r_h by a factor of about two. On the other hand, ring viscosity changes by orders of magnitude when τ, r_h, or both are varied. In the following, we will use our numerical results in the case of ε_n = ε_n(v) and ε_t = 0.5 and derive semianalytic expressions that approximately reproduce the dependence of these results on τ and r_h.

First, in the case of low optical depth, viscosity is determined by particle velocity dispersion c, as shown by Goldreich & Tremaine (1978) for nongravitating rings. We also discussed this in Section 3 for the case of self-gravitating rings. Viscosity in this case can be written as

$$\begin{eqnarray} &&\nu _{{\rm low}} = \alpha c^{2} \tau / \Omega, \end{eqnarray} \tag{ 6 }$$

where α is a constant. In optically thin rings, the velocity dispersion c can be written as (e.g., Salo 1995; Ohtsuki 1999)

$$\begin{eqnarray} &&c \sim \max \big\lbrace 2R\Omega, \; \sqrt{\strut 2Gm/R} \big\rbrace, \end{eqnarray} \tag{ 7 }$$

depending on either inelastic collisions or gravitational encounters dominate particle velocity evolution. Substituting Equation (7) into Equation (6), viscosity in low optical depth rings can be expressed in terms of optical depth, distance from the planet, and particle internal density, as

$$\begin{eqnarray} \frac{\nu _{{\rm low}}}{R^{2}\Omega } &\sim& \max \bigg\lbrace 4, \; 15 \times \left(\frac{\rho }{0.9 \, {\rm g} \, {\rm cm}^{-3}}\right) \left(\frac{a}{10^{10} \, {\rm cm}}\right)^{3}\nonumber\\ &&\times \left(\frac{M_{\mathrm{c}}}{5.7 \times 10^{29} \, {\rm g}} \right)^{-1} \bigg\rbrace \alpha \tau. \end{eqnarray} \tag{ 8 }$$

We find that α ∼ 0.35 from the comparison with our numerical results, which is roughly consistent with the nongravitating case (Goldreich & Tremaine 1978).

Next, we examine the case of high optical depth. From a dimensional analysis (see also Lynden-Bell & Kalnajs 1972; Ward & Cameron 1978), Daisaka et al. (2001) argued that the viscosity in dense rings is expected to be proportional to that of a self-gravitating disk,

$$\begin{eqnarray} &&\nu _{\scriptscriptstyle {\rm GD}} = G^{2} \Sigma ^{2} / \Omega ^{3}, \end{eqnarray} \tag{ 9 }$$

with a correction factor C(r_h), which depends on the relative importance of inelastic collisions to the self-gravity of a ring. In terms of τ and r_h, Equation (9) can be rewritten as (Daisaka et al. 2001)

$$\begin{eqnarray} &&\nu _{\scriptscriptstyle {\rm GD}}/ (R^{2}\Omega ) = (144/ \pi ^{2}) r_{\mathrm{h}}^{6} \tau ^{2}. \end{eqnarray} \tag{ 10 }$$

Daisaka et al. (2001) found C(r_h) = 26r⁵_h from a comparison with their numerical results.

Although our results of N-body simulations show that the dependence of the viscosity on optical depth is somewhat stronger than that given by Equation (9) (i.e., dln ν/dln τ > 2) when viscosity begins to deviate from the ν∝τ relation with increasing τ (Section 4), the τ-dependence of our numerical results for the high-τ cases is roughly consistent with Equation (9). Therefore, following Daisaka et al. (2001), we evaluate the correction factor $C = \nu / \nu _{\scriptscriptstyle {\rm GD}}$ using our numerical results (Figure 9(a)). We find that a correction factor given by

$$\begin{eqnarray} &&C_{\scriptstyle {\rm wake}}(r_{\mathrm{h}} ) \simeq 53 r_{\mathrm{h}}^{5} \end{eqnarray} \tag{ 11 }$$

seems to well reproduce our numerical results for dense rings with gravitational wakes (i.e., r_h ∼ 0.7–0.8 and τ ≳ 0.5). The fifth-power dependence of the correction factor on r_h is consistent with Daisaka et al. (2001), but the numerical coefficient is a factor of about two larger than Daisaka et al. (2001). This is due to the difference in restitution coefficients adopted in these simulations. We used our numerical results for the case with ε_n = ε_n(v) and ε_t = 0.5, while Daisaka et al. (2001) derived their correction factor from their simulations with smooth particles with normal restitution coefficients independent of impact velocity. In fact, we confirmed that our numerical results as well as those obtained by Schmidt et al. (2009) for ε_n = 0.5 and ε_t = 1 are consistent with Daisaka et al. (2001; Figure 9(b)). As we discussed in Section 4, including surface friction as well as using the velocity-dependent normal restitution coefficient somewhat increases the viscosity values, which explains the above difference in the correction factor. Thus, it should be noted that our semianalytic expressions shown below include uncertainty of a factor of 2–3 due to uncertainty in elastic properties of ring particles.

Zoom In Zoom Out Reset image size

Daisaka et al. (2001) argued that the r_h-dependence of the correction factor is related to "compressibility" of particle disks due to self-gravity. When r_h is significantly smaller than unity (e.g., r_h ∼ 0.5), collisions suppress the tendency of formation of gravitational wakes; thus, the enhancement of viscosity due to self-gravity is insignificant. This compressibility discussed in Daisaka et al. (2001) can be quantitatively demonstrated by calculating the physical optical depth τ_phys using a spatial distribution of particles obtained by N-body simulation; here, τ_phys is defined so that exp (− τ_phys) is the average fraction of skewers that penetrate the ring without piercing any particles (Figure 10; Wisdom & Tremaine 1988; Robbins et al. 2010). In the case of large τ and r_h, we confirm significant reduction in τ_phys compared with the dynamical optical depth τ(= NπR²/L²); thus, there is a stronger tendency of particle clumping and wake formation, which results in a significant increase in viscosity.

Zoom In Zoom Out Reset image size

Figure 11(a) shows the plots of our analytic results obtained by $\nu = \max \lbrace \nu _{{\rm low}}, C_{\scriptstyle {\rm wake}} \nu _{\scriptscriptstyle {\rm GD}} \rbrace$ ($C_{\scriptstyle {\rm wake}}$ is given by Equation (11)) together with numerical results of our N-body simulation. We confirm excellent agreement in both low and high optical depth regimes for r_h = 0.5–0.82, as expected. On the other hand, the above analytic expression significantly overestimates the numerical results for r_h = 1, as can be predicted from the plots of the correction factor in Figure 9(a). In the cases of r_h ≃ 1, the relatively large size of the Hill sphere compared with the particle physical size allows the formation of temporary aggregates rather than gravitational wakes (Salo 1995), as we mentioned above (Figure 8). Therefore, the above disagreement between the analytic result for the viscosity that is valid in the case of wake formation and the numerical results for r_h ≃ 1 is not surprising; viscosity in rings with temporary aggregate formation, which likely corresponds to the outer A ring of Saturn, should be treated separately from the case of rings with gravitational wakes. Our numerical results show that the τ²-dependence suggested by Equation (9) also seems to be approximately valid in the case of r_h  =  1 and τ  ≃  1, but the correction factor should be separately determined. Since our numerical results for r_h ≃ 1 suggest that the dependence of the correction factor on τ is rather weak (Figure 9(a)), here we adopt

$$\begin{eqnarray} &&C_{\scriptstyle {\rm agg}} \simeq 30 \end{eqnarray} \tag{ 12 }$$

as a simple remedy for the correction factor in the case of temporary aggregate formation. Figure 11(b) shows a comparison between the results obtained from the revised semianalytic expression with our results of N-body simulation, where we find excellent agreement for the entire parameter regime examined here.

Zoom In Zoom Out Reset image size

In summary, our numerical results can be well approximated by the following expression:

$$\begin{eqnarray} &&\nu = \max \left\lbrace \nu _{{\rm low}}, C \nu _{\scriptscriptstyle {\rm GD}} \right\rbrace, \end{eqnarray} \tag{ 13 }$$

where ν_low and $\nu _{\scriptscriptstyle {\rm GD}}$ are given by Equations (8) and (9), respectively, and the correction factor is C = min {C_wake(r_h), C_agg}, where C_wake(r_h) and C_agg are given by Equations (11) and (12), respectively. Again, it should be noted that viscosity values obtained by the above expression include an uncertainty factor of about 2–3 due to uncertainty in the elastic properties of particles.

In this work, using local N-body simulation, we examined viscosity in self-gravitating rings for a wide range of parameters. We calculated ring viscosity based on the formulation developed by Tanaka et al. (2003), in which viscosity is evaluated by calculating the average energy dissipation rate due to inelastic collisions in N-body simulation. We investigated the dependence of viscosity on various parameters in detail, including the effects of particle surface friction.

In the case of self-gravitating rings with low optical depth, we found that viscosity is determined by particle random velocity, as in the case of nongravitating rings with low optical depth. The effects of surface friction on viscosity were found to depend on the importance of the mutual gravity of the particles relative to collisions, which can be expressed in terms of the ratio of the mutual Hill radius of colliding particles to their physical size (r_h). When the effects of gravitational encounters are unimportant (r_h ∼ 0.5) and particle random velocity is determined by inelastic collisions, the inclusion of surface friction slightly reduces both random velocity and viscosity. On the other hand, when r_h ∼ 1 and gravitational encounters do play a major role in particle velocity evolution, the inclusion of surface friction results in the increase of viscosity, since the viscous heating increases to balance the increased energy dissipation at collisions due to surface friction. We found that varying restitution coefficients can change the viscosity in dilute rings by a factor of about two.

On the other hand, viscosity is significantly enhanced in dense self-gravitating rings where gravitational wakes are formed (Daisaka et al. 2001). We found that viscosity in dense rings with a given τ and r_h is larger when collisions are more dissipative, for example, by the inclusion of surface friction, and the enhancement of viscosity due to this effect can be a factor of about two. Following Daisaka et al. (2001), we assumed that the dependence of viscosity on the optical depth can be approximated by that for a self-gravitating disk (Equation (9)), so we evaluated a correction factor using our numerical results of N-body simulation. We found that the expression for the correction factor of the dense rings with gravitational wakes is similar to the one obtained by Daisaka et al. (2001); the difference of a factor two in the numerical coefficient can be explained by the difference in restitution coefficients that were adopted. On the other hand, in the case of the outer rings where temporary aggregates form due to self-gravity, the above expression overestimates the viscosity values; thus, we derived a revised expression for such a case. We confirmed that the semianalytic expressions we obtained reproduce our numerical results well for the entire range of parameters examined here.

The expression we obtained for the viscosity in the case with low optical depth is consistent with previous theoretical works for nongravitating rings (e.g., Goldreich & Tremaine 1978). However, the expressions for the viscosity (or the expression for the correction factor) for dense rings were empirically obtained by comparing numerical results. Although the strong r_h-dependence of the viscosity in the case of the rings with gravitational wakes has been explained by ring compressibility (i.e., the degree of density contrast created as a result of the formation of gravitational wakes), more detailed theoretical modeling of viscosity in self-gravitating collisional disks would be useful for a better understanding of the dynamics in dense planetary rings.

This work was supported by the Japanese Society for Promotion of Science and NASA's Outer Planets Research Program and Planetary Geology and Geophysics Program. We thank Ryuji Morishima, Heikki Salo, Jürgen Schmidt, and Hidekazu Tanaka for their discussion and useful comments. A portion of the numerical simulations was performed using the GRAPE-7 system at the Center for Computational Astrophysics of the National Astronomical Observatory of Japan.

When two particles, i and j, collide with each other, the changes of their velocities and rotation rates are calculated as follows. We assume that ring particles are spheres, and we let the mass, radius, and spin angular velocity vector of particle j be given by m_j, R_j, and ω_j, respectively. Suppose that two particles, i and j, collide with each other with relative velocity of the centers of the two bodies in the rotating coordinate system, v, and that their relative position at impact is given by r. Let v_n and v_t denote the normal and the tangential component of v to the tangent plane, respectively. In this case, the normal and the tangential components of the relative velocity of the two contacting points at the time of impact can be written as (e.g., Araki & Tremaine 1986; Richardson 1994; Ohtsuki 2006a, 2006b)

$$\begin{equation} \begin{array}{l} \mbox{\boldmath $u$}_\mathrm{n} = \mbox{\boldmath $v$}_\mathrm{n} = (\mbox{\boldmath $v$} \cdot \mbox{\boldmath $\lambda $}) \mbox{\boldmath $\lambda $}, \\ \ \mbox{\boldmath $u$}_\mathrm{t} = \mbox{\boldmath $v$}_\mathrm{t} + \mbox{\boldmath $\Omega $}\times \mbox{\boldmath $r$} + \mbox{\boldmath $\lambda $}\times (R_i \mbox{\boldmath $\omega $}_i + R_j \mbox{\boldmath $\omega $}_j), \end{array} \end{equation} \tag{ A1 }$$

where λ ≡ r/|r| is the unit vector pointing from the center of particle i to that of particle j, and Ω ≡ (0, 0, Ω) is the Keplerian angular velocity vector of the coordinate system. In the above equation, we used the relation r = (R_i + R_j)λ at impact. In terms of the normal and the tangential restitution coefficients, ε_n and ε_t, the normal and the tangential components of the relative velocity of the two contacting points after impact are given by Equation (2). From Equation (2) and the conservation of linear and angular momenta, the change of the normal and the tangential components of the relative velocity of the centers of the two spheres are given as

$$\begin{equation} \begin{array}{l}\Delta \mbox{\boldmath $v$}_\mathrm{n} = -(1+\varepsilon _\mathrm{n}) \mbox{\boldmath $v$}_\mathrm{n}, \\ \ \Delta \mbox{\boldmath $v$}_\mathrm{t} ={\displaystyle -\frac{{\cal K}}{{\cal K}+1}(1-\varepsilon _\mathrm{t}) \lbrace \mbox{\boldmath $v$}_\mathrm{t} + R_\mathrm{p} \mbox{\boldmath $\lambda $}\times (\mbox{\boldmath $\omega $}- \mbox{\boldmath $\Omega $})\rbrace, } \end{array} \end{equation} \tag{ A2 }$$

where R_p ≡ R_i + R_j, and

$$\begin{equation} \mbox{\boldmath $\omega $} \equiv (R_i \mbox{\boldmath $\omega $}_i + R_j \mbox{\boldmath $\omega $}_j)/(R_i + R_j). \end{equation} \tag{ A3 }$$

In the above, ${\cal K}$ is the coefficient of the moment of inertia, e.g., $I_j = {\cal K}m_jR_j^{2}$. In the present work, we set ${\cal K}=2/5$, assuming a homogeneous sphere. The change of spin angular velocity of particle j due to the collision is written as

$$\begin{equation} \Delta \mbox{\boldmath $\omega $}_j = {\displaystyle \frac{\mu (1-\varepsilon _\mathrm{t})}{({\cal K}+1)m_jR_j} \mbox{\boldmath $\lambda $}\times \lbrace \mbox{\boldmath $v$}_\mathrm{t} + R_\mathrm{p} \mbox{\boldmath $\lambda $}\times (\mbox{\boldmath $\omega $}- \mbox{\boldmath $\Omega $}) \rbrace,} \end{equation} \tag{ A4 }$$

where μ is the reduced mass.

When a local ring region is in a quasisteady state, its average viscosity can be calculated by summing up the amount of energy dissipated at collisions that occur in the region during a certain time interval (Equation (4)). Following Tanaka et al. (2003, their Equation (43)), we calculate the amount of energy dissipated at a collision between particles j and k as

$$\begin{equation} \Delta E_{\mathrm{col}, jk} = - \frac{1}{2} \mu h^{2} a_{0}^{2} \Omega ^{2} \left(\Delta \tilde{e}^{2}_{\mathrm{col}, jk} + \Delta \tilde{i}^{2}_{\mathrm{col}, jk} - \frac{3}{4} \Delta \tilde{b}^{2}_{\mathrm{col}, jk} \right). \end{equation} \tag{ B1 }$$

In the above, μ = (1/m_j + 1/m_k)⁻¹ (m_j is the mass of particle j, and m_j = m in the equal-sized case studied in the present work); h = (2m/3M_c)^1/3; and $\Delta \tilde{e}^{2}_{\mathrm{col}, jk}$, $\Delta \tilde{i}^{2}_{\mathrm{col}, jk}$ and $\Delta \tilde{b}^{2}_{\mathrm{col}, jk}$ are the changes in the square of the scaled relative orbital elements $\tilde{e}$, $\tilde{i},$ and $\tilde{b}$ due to the collision, where $\tilde{e}$ and $\tilde{i}$ are the eccentricity and the inclination of the relative motion scaled by h (see Ohtsuki 1999), and $\tilde{b}$ is the difference in the semimajor axes of the two particles scaled by ha₀. Note that the effects of time variation of the total self-gravitational energy, random kinetic energy, or rotational energy of particles do not appear in Equation (B1) because we assume the region is in a quasisteady state.