Issue 
A&A
Volume 524, December 2010



Article Number  A22  
Number of page(s)  13  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201015177  
Published online  22 November 2010 
Stochastic orbital migration of small bodies in Saturn’s rings^{⋆}
University of Cambridge, Department of Applied Mathematics and
Theoretical Physics, Centre for Mathematical Sciences, Wilberforce Road, Cambridge
CB3 0WA,
UK
email: hr260@cam.ac.uk
Received:
8
June
2010
Accepted:
11
August
2010
Many small moonlets that create propeller structures have been found in Saturn’s rings by the Cassini spacecraft. We study the dynamical evolution of such 20−50 m sized bodies, which are embedded in Saturn’s rings. We estimate the importance of various interaction processes with the ring particles on the moonlet’s eccentricity and semimajor axis analytically. For low ring surface densities, the main effects on the evolution of the eccentricity and the semimajor axis are found to be caused by collisions and the gravitational interaction with particles in the vicinity of the moonlet. For high surface densities, the gravitational interaction with selfgravity wakes becomes important.
We also perform realistic threedimensional, collisional Nbody simulations with up to a quarter of a million particles. A new set of pseudo shear periodic boundary conditions is used, which reduces the computational costs by an order of magnitude compared to previous studies. Our analytic estimates are confirmed to within a factor of two.
On short timescales the evolution is always dominated by stochastic effects caused by collisions and gravitational interaction with selfgravitating ring particles. These result in a random walk of the moonlet’s semimajor axis. The eccentricity of the moonlet quickly reaches an equilibrium value owing to collisional damping. The average change in semimajor axis of the moonlet after 100 orbital periods is 10100 m. This translates to an offset in the azimuthal direction of several hundred kilometres. We expect that such a shift is easily observable.
Key words: celestial mechanics / planets and satellites: rings / planets and satellites: dynamical evolution and stability / planetdisk interactions
Two movies are only available in electronic form at http://www.aanda.org
© ESO, 2010
1. Introduction
Small bodies called moonlets, which are embedded in Saturn’s rings, can create propellershaped structures by disturbing the rings. These have been predicted both analytically and numerically (Spahn & Sremčević 2000; Sremčević et al. 2002; Seiß et al. 2005). Only recently have they been observed by the Cassini spacecraft in both the A and B rings (Tiscareno et al. 2006, 2008).
These 20−100 m sized bodies can migrate within the rings, similar to protoplanets that migrate in a protostellar disc. Depending on the disc properties and the moonlet size, this can happen in either a smooth or in a stochastic (randomwalk) fashion. We refer to those migration regimes as types I and IV, respectively, in analogy to the terminology in discplanet interactions. Crida et al. (2010) show that, in the case of a laminar disc, there is a type I regime that might be important on very long timescales. This migration is qualitatively different to the migration in a pressure supported gas disc. However, the migration of moonlets in the A ring is generally dominated by type IV migration, at least on short timescales.
In this paper, we study the type IV migration regime in detail. The full mutual interaction of the ring particles with the moonlet and its consequent induced motion are considered both analytically and numerically.
We first review the basic equations governing the moonlet and the ring particles in a shearing box approximation in Sect. 2. Then, in Sect. 3, we estimate the eccentricity damping timescale due to ring particles colliding with the moonlet and ring particles on horseshoe orbits, as well as the effect of particles on circulating orbits. In Sect. 4 we estimate the excitation of the moonlet eccentricity caused by stochastic particle collisions and gravitational interactions with ring particles. This enables us to derive an analytic estimate of the equilibrium eccentricity. In Sect. 5 we discuss and evaluate processes, such as collisions and gravitational interactions with ring particles and selfgravitating clumps, which lead to a random walk in the semimajor axis of the moonlet.
We describe our numerical code and the initial conditions used in Sect. 6. We perform realistic three dimensional Nbody simulations of the ring system and the moonlet, taking into account a moonlet with finite size, a size distribution of ring particles, selfgravity, and collisions. The results are presented in Sect. 7. All analytic estimates are confirmed both in terms of qualitative trends and quantitatively to within a factor of about two in all the simulations that we performed. We also discuss the longterm evolution of the longitude of the moonlet and its observability, before summarising our results in Sect. 8.
2. Basic equations governing the moonlet and ring particles
Fig. 1 Trajectories of ring particles in a frame centred on the moonlet. Particles accumulate near the moonlet and fill its Roche lobe. Particles on trajectories labelled (a) are on circulating orbits. Particles on trajectories labelled (b) collide with other particles in the moonlet’s vicinity. Particles on trajectories labelled (c) collide with the moonlet directly. Particles on trajectories labelled (d) are on horseshoe orbits. Particles on trajectories labelled (e) leave the vicinity of the moonlet. 
We adopt a local righthanded Cartesian coordinate system with its origin in a circular Keplerian orbit with semimajor axis a and rotating uniformly with angular velocity Ω. This orbit coincides with that of the moonlet when it is assumed to be unperturbed by ring particles. The x axis coincides with the line joining the central object of mass M_{p} and the origin. The unit vector in the x direction e_{x}, points away from the central object. The unit vector in the y direction e_{y}, points in the direction of rotation and the unit vector in the z direction, e_{z}, points in the vertical direction being normal to the disc midplane.
In general, we consider a ring particle of mass m_{1} interacting with the moonlet that has a much higher mass m_{2}. A sketch of three possible types of particle trajectory in the vicinity of the moonlet is shown in Fig. 1. These correspond to three distinct regimes, a) denoting circulating orbits, b) and c) denoting orbits that result in a collision with particles in the vicinity of the moonlet and directly with the moonlet, respectively, and d) denoting horseshoe orbits. All these types of trajectories occur for particles that are initially in circular orbits, both interior and exterior to the moonlet, when at large distances from it.
Approximating the gravitational force due to the central object by its firstorder Taylor expansion about the origin leads to Hill’s equation, governing the motion of a particle of mass m_{1} of the form (1)where r = (x,y,z) is position of a particle with mass m_{1}, and Ψ is the gravitational potential acting on the particle due to other objects such as the moonlet.
The square amplitude of the epicyclic motion ℰ^{2} can be defined through (2)Note that neither the eccentricity e nor the semimajor axis a are formally defined in the local coordinate system. However, in the absence of any interaction with other masses (Ψ = 0), ℰ is conserved, and up to first order in the eccentricity we may make the identification ℰ = ea. We recall the classical definition of the eccentricity e in a coordinate system centred on the central object (3)Here, s is the position vector (a,0,0), and w is the velocity vector of the particle relative to the mean shear associated with circular Keplerian motion as viewed in the local coordinate system plus the circular Keplerian velocity corresponding to the orbital frequency Ω of the origin. Thus, (4)All eccentricities considered here are very small ( ~ 10^{8}−10^{7}) so that the difference between the quantities defined through use of Eqs. (2) and (3) is negligible, allowing us to adopt ℰ as a measure of the eccentricity throughout this paper.
Let us define another quantity, , that is also conserved for noninteracting particle motion, Ψ = 0, which is the x coordinate of the centre of the epicyclic motion and is given by (5)We identify a change in as a change in the semimajor axis a of the particle, again, under the assumption that the eccentricity is small.
2.1. Two interacting particles
We now consider the motion of two gravitationally interacting particles with position vectors r_{1} = (x_{1},y_{1},z_{1}) and r_{2} = (x_{2},y_{2},z_{2}). Their corresponding masses are m_{1} and m_{2}. The governing equations of motion are where the interaction gravitational potential is Ψ_{12} = −Gm_{1}m_{2}/r_{1} − r_{2}. The position vector of the centre of mass of the two particles is given by (8)The vector also obeys Eq. (1) with Ψ = 0, which applies to an isolated particle, because the interaction potential does not affect the motion of the centre of mass. We also find it useful to define the vector (9)Then, consistently with our earlier definition of ℰ, we have ℰ_{i} = ℰ_{i} . The amplitude of the epicyclic motion of the centre of mass is given by (10)Here φ_{12} is the angle between E_{1} and E_{2}. It is important to note that is conserved even if the two particles approach each other and become bound. This is true as long as frictional forces are internal to the two particle system and do not affect the centre of mass motion.
3. Effects leading to damping of the eccentricity of the moonlet
We begin by estimating the moonlet eccentricity damping rate associated with ring particles that either collide directly with the moonlet, particles in its vicinity, or only interact gravitationally with the moonlet.
3.1. The contribution due to particles impacting the moonlet
Particles hitting the moonlet in an eccentric orbit exchange momentum with it. Let us assume that a ring particle, identified with m_{1} has zero epicyclic amplitude, so that ℰ_{1} = 0 far away from the moonlet. The moonlet is identified with m_{2} and has an initial epicyclic amplitude ℰ_{2}. The epicyclic amplitude of the centre of mass is therefore (11)where we have assumed that m_{1} ≪ m_{2}.
The moonlet is assumed to be in a steady state where there is no net accretion of ring particles. Therefore, for every particle that either collides directly with the moonlet or nearby particles bound to it (and so itself becomes temporarily bound to it), one particle must also escape from the moonlet. This happens primarily through slow leakage from locations close to the L_{1} and L_{2} points such that most particles escape from the moonlet with almost zero velocity (as viewed from the centre of mass frame) and so do not change its orbital eccentricity. However, after an impacting particle becomes bound to the moonlet, conservation of the epicyclic amplitude associated with the centre of mass motion, together with Eq. (11), imply that each impacting particle will reduce the eccentricity by a factor 1 − m_{1}/m_{2}.
It is now an easy task to estimate the eccentricity damping timescale by determining the number of particle collisions per time unit with the moonlet or particles bound to it. To do that, a smooth window function W_{b + c}(b) is used, being unity for impact parameters b that always result in an impact with the moonlet or particles nearby that are bound to it, and being zero for impact parameters that never result in an impact.
The number of particles impacting the moonlet per time unit, dN/dt, is obtained by integrating over the impact parameter with the result that (12)where Σ is the surface mass density. Allowing b to be negative enables impacts from both sides of the moonlet to be taken into account.
Therefore, after using Eq. (11), we find that the rate of change of the moonlet’s eccentricity e_{2}, or equivalently of its amplitude of epicyclic motion ℰ_{2}, is given by (13)where τ_{e,collisions} defines the circularisation time arising from collisions with the moonlet. The natural unit for b is the Hill radius of the moonlet, r_{H} = (m_{2}/(3M_{p}))^{1/3}a, so that the dimensional scaling for τ_{e,collisions} is given by (14)which we find to also apply to all the processes for modifying the moonlet’s eccentricity discussed below. If we assume that W_{b + c}(b) can be approximated by a box function, being unity in the interval [1.5r_{H},2.5r_{H}] and zero elsewhere, we get (15)
3.2. Eccentricity damping due to the interaction of the moonlet with particles on horseshoe orbits
The eccentricity of the moonlet manifests itself in a small oscillation of the moonlet about the origin. Primarily ring particles on horseshoe orbits will respond to that oscillation and damp it. This is because only those particles on horseshoe orbits spend enough time in the vicinity of the moonlet, i.e. many epicyclic periods.
In Appendix A, we calculate the amplitude of epicyclic motion ℰ_{1f} (or equivalently the eccentricity e_{1f}) that is induced in a single ring particle in a horseshoe orbit undergoing a close approach to a moonlet that is assumed to be in an eccentric orbit. In order to calculate the circularisation time, we have to consider all relevant impact parameters. We begin by noting that each particle encounter with the moonlet is conservative and is such that for each particle, the Jacobi constant, applicable when the moonlet is in a circular orbit, is increased by an amount by the action of the perturbing force, associated with the eccentricity of the moonlet, as the particle passes by. Because the Jacobi constant, or energy in the rotating frame, for the moonlet and the particle together is conserved, the square of the epicyclic amplitude associated with the moonlet alone changes by Accordingly the change in the amplitude of epicyclic motion of the moonlet ℰ_{2}, consequent on inducing the amplitude of epicyclic motion ℰ_{1f} in the horseshoe particle, is given by (16)Note that this is different compared to Eq. (11). Here, we are dealing with a secondorder effect. First, the eccentric moonlet excites eccentricity in a ring particle. Second, because the total epicyclic motion is conserved, the epicyclic motion of the moonlet is reduced.
Integrating over the impact parameters associated with ring particles and taking particles streaming by the moonlet from both directions into account by allowing negative impact parameters gives (17)where τ_{e,horseshoe} is the circularisation time and, as above, we have inserted a window function, which is unity on impact parameters that lead to horseshoe orbits, otherwise being zero. Using ℰ_{1f} given by Eq. (A.20), we obtain (18)For a sharp cutoff of W_{d}(b) at b_{m} = 1.5r_{H} we find the value of the integral in the above to be 2.84. Thus, in this case we get (19)However, note that this value is sensitive to the value of b_{m} adopted. For b_{m} = 1.25r_{H},τ_{e,horseshoe} is a factor of 4.25 larger.
3.3. The effects of circulating particles
The effect of the response of circulating particles to the gravitational perturbation of the moonlet on the moonlet’s eccentricity can be estimated from the work of Goldreich & Tremaine (1980, 1982). These authors consider a ring separated from a satellite such that coorbital effects are not considered. Thus, their expressions may be applied to estimate effects of circulating particles. However, we exclude their corotation torques because they are determined by the ring edges and are absent in a local model with uniform azimuthallyaveraged surface density. Equivalently, one may simply assume that the corotation torques are saturated. When this is done only Lindblad torques act on the moonlet. These tend to excite the moonlet’s eccentricity rather than damp it.
We replace the ring mass in Eq. (70) of Goldreich & Tremaine (1982) by the integral over the impact parameter and switch to our notation to obtain (20)where W_{a}(b) is the appropriate window function for circulating particles. Assuming a sharp cutoff of W_{a}(  b  ) at b_{m}, we can evaluate the integral in Eq. (20) to get (21)where we have adopted b_{m} = 2.5r_{H}, consistent with the simulation results presented below. We see that, although τ_{e,circ} < 0 corresponds to growth rather than damping of the eccentricity, it scales in the same manner as the circularisation times in Sects. 3.1 and 3.2 scale (see Eq. (14)).
This timescale and the timescale associated with particles on horseshoe orbits, τ_{e,horseshoe}, are significantly larger than the timescale associated with collisions, given in Eq. (15). We can therefore ignore those effects for most of the discussion in this paper.
4. Processes leading to the excitation of the eccentricity of the moonlet
In Sect. 3 we assumed that the moonlet had a small eccentricity but neglected the initial eccentricity of the impacting ring particles. When this is included, collisions and gravitational interactions of ring particles with the moonlet may also excite its eccentricity.
4.1. Collisional eccentricity excitation
To see this, assume that the moonlet initially has zero eccentricity, or equivalently no epicyclic motion, but the ring particles do. As above we consider the conservation of the epicyclic motion of the centre of mass in order to connect the amplitude of the final epicyclic motion of the combined moonlet and ring particle to the initial amplitude of the ring particle’s epicyclic motion. This gives the epicyclic amplitude of the centre of mass after one impact as (22)Assuming that successive collisions are uncorrelated and occur stochastically with the mean time between consecutive encounters being τ_{ce}, the evolution of is governed by the equation (23)where can be expressed in terms of the surface density and an impact window function W_{b + c}(b) (see Eq. (12)). The quantity is the mean square value of ℰ_{1} for the ring particles. Using Eq. (2), this may be written in terms of the mean squares of the components of the velocity dispersion relative to the background shear, in the form (24)We find ⟨ e_{1} ⟩ ~ 10^{6} in numerical simulations for almost all ring parameters. This value is mainly determined by the coefficient of restitution (see e.g. Fig. 4 in Morishima & Salo 2006).
4.2. Stochastic excitation by circulating particles
Ring particles that are on circulating orbits, as given by path (a) in Fig. 1, exchange energy and angular momentum with the moonlet and therefore change its eccentricity. Goldreich & Tremaine (1982) calculated the change in the eccentricity of a ring particle due to a satellite. We are interested in the change of the eccentricity of the moonlet induced by a ring of particles and therefore swap the satellite mass with the mass of a ring particle. Thus, rewriting their results (Eq. (64) in Goldreich & Tremaine 1982) in our local notation without reference to the semimajor axis, we have (25)Supposing that the moonlet has very small eccentricity, it will receive stochastic impulses that cause its eccentricity to undergo a random walk that will result in increasing linearly with time, so that we may write (26)where d(1/t_{b}) is the mean encounter rate with particles that have an impact parameter in the interval (b,b + db) and W_{a}(b) is the window function describing the band of particles in circulating orbits. Setting d(1/t_{b}) = 3ΣΩ  b  db/(2m_{1}), we obtain (27)For high surface densities, an additional effect can excite the eccentricity when selfgravity wakes occur. The Toomre parameter Q is defined as (28)where is the velocity dispersion of the ring particles^{1}. When the surface density is high enough for Q to approach unity, transient clumps will form in the rings. The typical length scale of these structures is given by the critical Toomre wave length λ_{T} = 2π^{2}GΣΩ^{2} (Daisaka et al. 2001), which can be used to estimate a typical mass of the clump: (29)Whenever strong clumping occurs, m_{T} should be used in the above calculation instead of the mass of a single particle m_{1}: (30)where is the appropriate window function. For typical parameters used in Sect. 7, this transition occurs at Σ ~ 200 kg/m^{2}.
4.3. Equilibrium eccentricity
Putting together the results from this and the previous section, we can estimate an equilibrium eccentricity of the moonlet. Let us assume that the eccentricity e_{2}, or the amplitude of the epicyclic motion ℰ_{2}, evolves under the influence of excitation and damping forces as (31)Both the excitation terms and the decay timescales are constant within our approximations so that Eq. (31) describes a stable equilibrium. The equilibrium eccentricity is then found by setting the above equation equal to zero and solving for ℰ_{2}.
To make quantitative estimates we need to specify the window functions W_{a}(b), W_{b + c}(b), and W_{d}(b) that determine in which impact parameter bands the particles are in (see Fig. 1). To compare our analytic estimates to the numerical simulations presented below, we measure the window functions numerically. Alternatively, one could simply use sharp cutoffs at some multiple of the Hill radius. We already made use of this approximation as a simple estimate in the previous sections. The results may vary slightly, but not significantly.
However, as the window functions are dimensionless, it is possible to obtain the dimensional scaling of ℰ_{2} by adopting the length scale applicable to the impact parameter to be the Hill radius r_{H} and simply assume that the window functions are of order unity. As already indicated above, all of the circularisation times scale as τ_{e} ∝ Ωr_{H}G^{1}Σ^{1}, or equivalently ∝ . Assuming the ring particles have zero velocity dispersion, the eccentricity excitation is then due to circulating clumps or particles, and the scaling of due to this cause is given by Eqs. (27) and (30) by (32)where m_{i} is either m_{1} or m_{T}. We may then find the scaling of the equilibrium value of ℰ_{2} from consideration of Eqs. (31) and (13) as (33)This means that the expected kinetic energy in the noncircular motion of the moonlet is ∝ which can be viewed as stating that the noncircular moonlet motion scales in equipartition with the mass m_{i} moving with speed r_{H}Ω. This speed applies when the dispersion velocity associated with these masses is zero, indicating that the shear across a Hill radius replaces the dispersion velocity in that limit.
In the opposite limit, when the dispersion velocity exceeds the shear across a Hill radius and the dominant source of eccentricity excitation is collisions, Eq. (31) gives (34)so that the moonlet is in equipartition with the ring particles. Results for the two limiting cases can be combined to give an expression for the amplitude of the epicyclic motion excited in the moonlet of the form (35)where C_{i} is a constant of order unity. This indicates the transition between the shear dominated and the velocity dispersion dominated limits.
5. Processes leading to a random walk in the semimajor axis of the moonlet
We have established estimates for the equilibrium eccentricity of the moonlet in the previous section. Here, we estimate the random walk of the semimajor axis of the moonlet. In contrast to the case of the eccentricity, there is no tendency for the semimajor axis to relax to any particular value, so that there are no damping terms, and the deviation of the semimajor axis from its value at time t = 0 grows on average as for large t. Depending on the surface density, there are different effects that dominate the contributions to the random walk of the moonlet. For low surface densities, collisions and horseshoe orbits are most important. For high surface densities, the random walk is dominated by the stochastic gravitational force of circulating particles and clumps. In this section, we estimate the strength of each effect.
5.1. Random walk due to collisions
Let us assume, without loss of generality, that the guiding centre of the epicyclic motion of a ring particle is displaced from the orbit of the moonlet by in the inertial frame, whereas in the corotating frame the moonlet is initially located at the origin with . When the impacting particle becomes bound to the moonlet, the guiding centre of the epicyclic motion of the centre of mass of the combined object, as viewed in the inertial frame, is then displaced from the original moonlet orbit by (36)This is the analogue of Eq. (22) for the evolution of the semimajor axis. For an ensemble of particles with impact parameter b, the average centre of epicyclic motion is . Thus we can write the evolution of due to consecutive encounters as which should be compared to the corresponding expression for the eccentricity in Eq. (23).
5.2. Random walk due to stochastic forces from circulating particles and clumps
Particles on a circulating trajectory (see path (a) in Fig. 1) that come close the moonlet will exert a stochastic gravitational force. The largest contributions occur for particles within a few Hill radii. Thus, we can crudely estimate the magnitude of the specific gravitational force (acceleration) due to a single particle as (39)When selfgravity is important, selfgravity wakes (or clumps) have to be taken into account, as Sect. 4.2. Then a rough estimate of the strongest specific gravitational forces due to self gravitating clumps is (40)where 2r_{H} has been replaced by 2r_{H} + λ_{T}. This allows for λ_{T} being significantly larger than r_{H}, in which case the approximate distance of the clump to the moonlet is λ_{T}.
Following the formalism of Rein & Papaloizou (2009), we define a diffusion coefficient as D = 2τ_{c} ⟨ f^{2} ⟩ , where f is an acceleration, τ_{c} is the correlation time, and the angle brackets denote a mean value. We here take the correlation time associated with these forces to be the orbital period, 2πΩ^{1}. This is the natural dynamical timescale of the systems and has been found to be a reasonable assumption from an analysis of the simulations described below. By considering the rate of change of the energy of the moonlet motion, we may estimate the random walk in due to circulating particles and selfgravitating clumps as follows (Rein & Papaloizou 2009): This is only a crude estimate of the random walk undergone by the moonlet. In reality several additional effects might also play a role. For example, circulating particles and clumps are clearly correlated, the gravitational wakes have a large extent in the azimuthal direction, and particles that spend a long time in the vicinity of the moonlet have more complex trajectories. Nevertheless, we find that the above estimates are correct up to a factor 2 for all the simulations that we performed (see below).
5.3. Random walk due to particles in horseshoe orbits
Finally, let us calculate the random walk induced by particles on horseshoe orbits. Particles undergoing horseshoe turns on opposite sides of the planet produce changes of opposite sign. Encounters with the moonlet are stochastic and therefore the semimajor axis will undergo a random walk. A single particle with impact parameter b will change the semimajor axis of the moonlet by (45)Analogous to the analysis in Sect. 5.1, the time evolution of is then governed by This equation is identical to Eq. (37) except for a factor 4, as particles with impact parameter b will leave the vicinity of the moonlet at − b.
6. Numerical simulations
Initial simulation parameters.
We perform realistic threedimensional simulations of Saturn’s rings with an embedded moonlet and verify the analytic estimates presented above. The nomenclature and parameters used for the simulations are listed in Table 1. The first column lists the name of the simulation. The second column gives the surface density of the ring. The third gives the moonlet radius. The fourth column gives the timestep. The fifth and sixth columns give the lengths of the main box as measured in the xyplane and the number of particles used, respectively.
6.1. Methods
Fig. 2 Shearing box, simulating a small patch of a planetary ring system. The particles that leave an auxiliary box in the azimuthal (y) direction reenter the same box on the other side and get copied into the main box at the corresponding location. All auxiliary boxes are equivalent and there are no curvature terms present in the shearing sheet. 
The gravitational forces are calculated with a BarnesHut tree (Barnes & Hut 1986). The numerical scheme is similar to that used by Rein et al. (2010). Additionally, we implemented a symplectic integrator for Hill’s equations (Quinn et al. 2010). This turned out to be beneficial for accurate energy conservation at almost no additional cost when the eccentricity of the moonlet ( ~ 10^{8}) was small and integrations were undertaken over many hundreds of dynamical timescales.
To speed up the calculations further, we run two coupled simulations in parallel. The first adopts the main box that incorporates a moonlet. The second adopts an auxiliary box which initially has the same number of particles but which does not contain a moonlet. This is taken to be representative of the unperturbed background ring. We describe this setup as enabling us to adopt pseudo shear periodic boundary conditions. We consider the main box together with eight equivalent auxiliary boxes to be stacked as illustrated in Fig. 2. All auxiliary boxes are identical copies whose centres are shifted according to the background shear as in a normal shearing sheet. On account of this motion, auxiliary boxes are removed when their centres are shifted in azimuth by more than 1.5L_{y} from the centre of the main box and then reinserted in the same orbit on the opposite side of the main box so that the domain under consideration is prevented from shearing out. If a particle in the main box crosses one of its boundaries, it is discarded. If a particle in the auxiliary box crosses one of its boundaries, it is reinserted on the other side of this box, according to normal shear periodic boundary conditions. But in addition it is also copied into the corresponding location in the main box. We describe this procedure as applying pseudo shear periodic boundary conditions.
In a similar calculation, Lewis & Stewart (2009) use a very long box (about 10 times longer than the boxes used in this paper) to ensure that particles are completely randomised between encounters with the moonlet. We are not interested in the long wavelength response that is created by the moonlet. Effects that are most important for the moonlet’s dynamical evolution are found to happen within a few Hill radii. Using the pseudo shear periodic boundary conditions, we ensure that incoming particles are uncorrelated and do not contain prior information about the perturber. This setup speeds up our calculations by more than an order of magnitude.
The gravity acting on a particle in the main box, which also contains the moonlet, is calculated by summing over the particles in the main box and all auxiliary boxes. The gravity acting on a particle in the auxiliary box is calculated the standard way, by using ghost boxes that are identical copies of the auxiliary box. Tests have indicated that our procedure does not introduce unwanted fluctuations in the gravitational forces.
The moonlet is allowed to move freely in the main box. However, to prevent it leaving the computational domain, as soon as the moonlet has left the innermost part (defined as extending one eighth of the box size), all particles are shifted together with the box boundaries, such that the moonlet is returned to the centre of the box. This is possible because the shearing box approximation is invariant with respect to translations in the xy plane (see Eq. (1)).
Collisions between particles are resolved using the instantaneous collision model and a velocity dependent coefficient of restitution given by Bridges et al. (1984): (48)where v_{ ∥ } is the impact speed projected on the axis between the two particles. The already existing tree structure is reused to search for nearest neighbours (Rein et al. 2010).
6.2. Initial conditions and tests
Throughout this paper, we use a distribution of particle sizes, r_{1}, ranging from 1 m to 5 m with a slope of q = −3. Thus The density of both the ring particles and the moonlet is taken to be 4.0kg/m^{2}. This is at the lower end of what has been assumed reasonable for Saturn’s A ring (Lewis & Stewart 2009). The moonlet radius is taken to be either 50m or 25m. We found that using a larger ring particle and moonlet density only leads to more particles being bound to the moonlet. This effectively increases the mass of the moonlet (or equivalently its Hill radius) and can therefore be easily scaled to the formalism presented here. Simulating a gravitational aggregate of this kind is computationally very expensive, as many more collisions have to be resolved each timestep.
Fig. 3 Snapshots of the particle distribution and the moonlet (blue) in simulation EQ40050DTW. The wake is much clearer than in Fig. 4 as the box size is twice as large. 
The initial velocity dispersion of the particles is set to 1mm/s for the x and y components and 0.4mm/s for the z component. The moonlet is placed at a semimajor axis of a = 130000km, corresponding to an orbital period of P = 13.3h. Initially, the moonlet is placed in the centre of the shearing box on a circular orbit.
We run simulations with a variety of moonlet sizes, surface densities, and box sizes. The dimensions of each box, as viewed in the (x,y) plane, are specified in Table 1. Because of the small dispersion velocities, the vertical motion is automatically strongly confined to the midplane. After a few orbits, the simulations reach an equilibrium state in which the velocity dispersion of particles does not change anymore.
We ran several tests to ensure that our results are converged. Simulation EQ40050DT uses a ten times longer timestep than simulation EQ40050. Simulations EQ5050DTW, EQ20050DTW, and EQ40050DTW use a box that is twice the size of that used in simulation EQ40050. Therefore, four times more particles are used. No differences in the equilibrium state of the ring particles or in the equilibrium state of the moonlet have been observed in any of those test cases.
A snapshot of the particle distribution in simulation EQ40050DTW is shown in Fig. 3. Snapshots of simulations EQ20025, EQ20050, EQ40025, and EQ40050, which use the smaller box size, are shown in Fig. 4. In all these cases, the length of the wake that is created by the moonlet is much longer than the box size. Nevertheless, the pseudo shear periodic boundary conditions allow us to simulate the evolution of the moonlet accurately.
Fig. 4 Snapshots of the particle distribution and the moonlet (blue) for different surface densities after 25 orbits. The wake created by the moonlet is much longer than the size of the box and therefore hardly visible in these images. 
Fig. 5 Impact bands for simulation EQ20050 (left) and EQ20025 (right). The vertical lines indicate sharp cutoffs at 1.5r_{H} and 2.5r_{H}. 
7. Results
7.1. Impact bands
The impact band window functions W_{a}(b), W_{b + c}(b), and W_{d}(b) are necessary for estimating the damping timescales and the strength of the excitation. To measure those, each particle that enters the box is labelled with its impact parameter. The possible outcomes are plotted in Fig. 1: (a) being a circulating particle which leaves the box on the opposite side, (b) representing a collision with other ring particles close to the moonlet where the maximum distance from the centre of the moonlet has been taken to be twice the moonlet radius, (c) being a direct collision with the moonlet, and finally (d) showing a horseshoe orbit in which the particle leaves the box on the same side as it entered.
The impact band (b) is considered in addition to the impact band (c) because some ring particles will collide with other ring particles that are (temporarily) bound to the moonlet. Thus, these collisions take part in transferring the energy and momentum to the moonlet. Actually, if the moonlet is simply a rubble pile of ring particles as suggested by Porco et al. (2007), then there might be no solid moonlet core and all collisions are in the impact band (b).
Figure 5 shows the impact bands, normalised to the Hill radius of the moonlet r_{H} for simulations EQ20050 and EQ20025^{2}. For comparison, we also plot sharp cutoffs, at 1.5r_{H} and 2.5r_{H}. Our results show no sharp discontinuity because of the velocity dispersion in the ring particles which is ~ 5mm/s. This corresponds to an epicyclic motion of ~ 40m or equivalently ~ 0.8r_{H} and ~ 1.6r_{H} (for r_{2} = 50m and r_{2} = 25m, respectively), which explains the cutoff width of approximately one Hill radii. The impact bands are almost independent of the surface density and only depend on the Hill radii and the mean epicyclic amplitude of the ring particles, or, more precisely, on the ratio thereof.
7.2. Eccentricity damping timescale
To measure the eccentricity damping timescale, we first let the ring particles and the moonlet reach an equilibrium and integrate them for 200 orbits. We then change the velocity of the moonlet and the ring particles within 2r_{H}. The new velocity corresponds to an eccentricity of 6 × 10^{7}, which is well above the equilibrium value. We then measure the decay timescale τ_{e} by fitting a function of the form (49)The results are given in Table 2. The second column lists the simulation results. The third and fourth columns list the analytic estimates of collisional and horseshoe damping timescales, respectively. They are in good agreement with the estimated damping timescales from Sects. 3.1 and 3.2, showing clearly that the most important damping process, in every simulation considered here, is indeed through collisional damping as predicted by comparing Eq. (15) with Eqs. (19) and (21).
7.3. Mean moonlet eccentricity
The mean eccentricity of the moonlet is measured in all simulations after several orbits when the ring particles and the moonlet have reached an equilibrium state. To compare this value with the estimates from Sect. 2, we set Eq. (31) equal to zero and use the analytic damping timescale listed in Table 2. The analytic estimates of the equilibrium eccentricity are calculated for each excitation mechanism separately to disentangle their effects. They are listed in Table 3. The second column lists the simulation results. The third, fourth, and fifth columns list the analytic estimates of the equilibrium eccentricity assuming a single excitation mechanism. The last column lists the analytic estimate of the equilibrium eccentricity summing over all excitation mechanisms. The sixth column lists the analytic estimate for the mean eccentricity using the sum of all excitation mechanisms.
For all simulations, the estimates are correct within a factor of about 2. For low surface densities, the excitation is dominated by individual particle collisions. For higher surface densities, it is dominated by the excitation due to circulating self gravitating clumps (selfgravity wakes). The estimates and their trends are surprisingly accurate, as we have ignored several effects (see below).
7.4. Random walk in semimajor axis
The random walk of the semimajor axis a (or equivalently the centre of epicyclic motion in the shearing sheet) of the moonlet in the numerical simulations are measured and compared to the analytic estimates presented in Sect. 5. We ran one simulation per parameter set. To get a statistically meaningful expression for the average random walk after a given time , we average all pairs of and for which t − t′ = Δt. In other words, we assume the system satisfies the Ergodic hypothesis.
then grows like and we can fit a simple square root function. This allows us to accurately measure the average growth in after Δt = 100 orbits by running one simulation for a long time and averaging over time, rather than running many simulations and performing an ensemble average. The measured values are given in the second column of Table 4. The second column lists the simulation results. The third to sixth columns list the analytic estimates of the change in semimajor axis assuming a single excitation mechanism. The last column lists the total estimated change in semimajor axis summing over all excitation mechanisms. These values correspond to the analytic expressions in Eqs. (46), (37), (41), and (44), respectively.
For low surface densities, the evolution of the random walk is dominated by collisions and the effect of particles on horseshoe orbits. For high surface densities, the main effect comes from the stochastic gravitational force of circulating clumps (selfgravity wakes).
7.5. Longterm evolution and observability
Eccentricity damping timescale τ_{e} of the moonlet in units of the orbital period.
Equilibrium eccentricity ⟨ e_{eq} ⟩ of the moonlet.
Change in semimajor axis of the moonlet after 100 orbits, .
The change in semimajor axis is relatively small at a few tens of metres after 100 orbits ( = 50 days). The change in longitude that corresponds to this is, however, is much larger. We therefore extend the discussion in Rein & Papaloizou (2009), who derive the time evolution of orbital parameters undergoing a random walk, to the time evolution of the mean longitude. This is of special interest in the current situation, because in observations of Saturn’s rings, it is much easier to measure the shift in the longitude of an embedded moonlet than the change in its semimajor axis.
The time derivative of the mean motion is given by (50)where F_{θ} is the timedependent stochastic force in the θ (or in shearing sheet coordinates, y) direction, and we have assumed a nearly circular orbit, e ≪ 1. The mean longitude is thus given by the double integral The root mean square of the difference compared to the unperturbed orbit (n_{0}) is where g(t) is the autocorrelation function of F_{θ}, which has been assumed to be exponentially decaying with an autocorrelation time τ, thus g(t) = exp(−  t  /τ). In the limit t ≫ τ, the above equation simplifies to (57)On a shorter timescale, the other terms in Eq. (56) are important, leading to a linear and oscillatory behaviour. For uncorrelated noise (e.g. collisions), Eq. (57) is replaced by (58)where ⟨ Δv ⟩ is the average velocity change per impulsive event and τ is taken to be the average time between consecutive events.
Instead of a stochastic migration, let us also assume a migration in a laminar disc with constant migration rate τ_{a}, defined by (59)As above, we can calculate the longitude as a function of time as the following double integral and thus (62)Note that Eq. (62) is an exact solution, whereas Eq. (57) is a statistical quantity, describing the root mean square value in an ensemble average. In the stochastic and laminar case, (Δλ)^{2} grows like ~ t^{3} and ~ t^{4}, respectively. This behaviour allows a stochastic and laminar migration to be distinguished by observing Δλ over an extended period of time.
Fig. 6 Offset in the azimuthal distance due to the random walk in the semimajor axis a measured in simulations EQ5025, EQ5050, EQ20025, EQ20050, EQ40025 and EQ40050. Individual moonlets may show linear, constant or oscillatory growth. On average, the azimuthal offset grows like t^{3/2}. 
8. Conclusions
In this paper, we have discussed the dynamical response of an embedded moonlet in Saturn’s rings to interactions with ring particles both analytically and by the use of realistic threedimensional many particle simulations. Both the moonlet and the ring particle density were taken to be 4.0kg/m^{2}. Moonlets of radius 25m and 50m were considered. Particle sizes ranging between 1m and 5m were adopted.
We analytically estimated the eccentricity damping timescale of the moonlet due to both collisions with ring particles and the response of ring particles to gravitational perturbations by the moonlet. We found the effects due to the response of particles on horseshoe and circulating orbits to be negligible. On the other hand, the stochastic impulses applied to the moonlet by circulating particles were found to cause the square of the eccentricity to grow linearly with time as did collisions with particles with non zero velocity dispersion. A balance between excitation and damping processes then leads to an equilibrium moonlet eccentricity.
We also estimated the magnitude of the random walk in the semimajor axis of the moonlet induced by collisions with individual ring particles and the gravitational interaction with particles and selfgravity wakes. There is no tendency for the semimajor axis to relax to any particular value, so that there are no damping terms. The deviation of the semimajor axis from its value at time t = 0 grows on average as for large t.
From our simulations we find that the evolution of the eccentricity is indeed dominated by collisions with ring particles. For high surface densities (more than 200 kg/m^{2}), the effect of selfgravity wakes also becomes important, leading to an increase in the mean steady state eccentricity of the moonlet. When the particle velocity dispersion is larger than Ωr_{H}, we obtain approximate energy equipartition between the moonlet and ring particles as far as epicyclic motion is concerned.
Similarly, the random walk in the semimajor axis was found to be dominated by collisions for low surface densities and by selfgravity wakes for high ones. In addition we have shown that on average, the difference in longitude Δλ of a stochastically forced moonlet grows with time like t^{3/2} for large t.
The distance travelled within 100 orbits (50 days at a distance of 130 000 km) is, depending on the precise parameters, between 10 m and 100 m. This translates to a shift in longitude of several hundred kilometres. The shift in λ is not necessarily monotonic on short timescales (see Fig. 6). We expect that such a shift should be easily detectable by the Cassini spacecraft. And indeed, the recently published paper by Tiscareno et al. (2010) provides such an observation. These long term observations show that at least one propeller undergoes sustained nonKeplerian orbital motion.
This should be compared with Fig. 6 in this paper, however, the moonlet size differs by a factor four. It can be seen that the observed shift by several hundred kilometres in the azimuthal direction over a year can be easily explained by the mechanisms described in this paper. Future observations and a statistical analysis of the existing observations can further constrain the stochastic nature of the nonKeplerian motion.
There are several effects that have not been included in this study. The analytic discussions in Sects. 3, 4, and 5 do not model the motion of particles correctly when their trajectories take them very close to the moonlet. Material that is temporarily or permanently bound to the moonlet is ignored. The density of both the ring particles and the moonlet are at the lower end of what is assumed to be reasonable (Lewis & Stewart 2009). This has the advantage that for our computational setup only a few particles (with a total mass of less than 10% of the moonlet) are bound to the moonlet at any given time. In a future study, we will extend the discussion presented in this paper to a wider variety of particle sizes, moonlet sizes and densities.
Online material
∑ _{i = a,b,c,d}W_{i} ≳ 1 because all particles are shifted from time to time when the moonlet is to far away from the origin (see Sect. 6.1). In this process, particles with the same initial impact parameter might become associated with the more than one impact band.
Acknowledgments
We thank Aurélien Crida for reading an early draft of this paper. Hanno Rein was supported by an Isaac Newton Studentship and St John’s College Cambridge. Simulations were performed on the computational facilities of the Astrophysical Fluid Dynamics Group and on Darwin, the Cambridge University HPC facility.
References
 Barnes, J., & Hut, P. 1986, Nature, 324, 446 [NASA ADS] [CrossRef] [Google Scholar]
 Bridges, F. G., Hatzes, A., & Lin, D. N. C. 1984, Nature, 309, 333 [NASA ADS] [CrossRef] [Google Scholar]
 Crida, A., Papaloizou, J. C. B., Rein, H., Charnoz, S., & Salmon, J. 2010, AJ, 140, 944 [NASA ADS] [CrossRef] [Google Scholar]
 Daisaka, H., Tanaka, H., & Ida, S. 2001, Icarus, 154, 296 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Tremaine, S. 1980, ApJ, 241, 425 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Goldreich, P., & Tremaine, S. 1982, ARA&A, 20, 249 [NASA ADS] [CrossRef] [Google Scholar]
 Lewis, M. C., & Stewart, G. R. 2009, Icarus, 199, 387 [NASA ADS] [CrossRef] [Google Scholar]
 Morishima, R., & Salo, H. 2006, Icarus, 181, 272 [NASA ADS] [CrossRef] [Google Scholar]
 Porco, C. C., Thomas, P. C., Weiss, J. W., & Richardson, D. C. 2007, Science, 318, 1602 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Quinn, T., Perrine, R. P., Richardson, D. C., & Barnes, R. 2010, AJ, 139, 803 [NASA ADS] [CrossRef] [Google Scholar]
 Rein H., & Papaloizou, J. C. B. 2009, A&A, 497, 595 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rein, H., Lesur, G., & Leinhardt, Z. M. 2010, A&A, 511, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Seiß, M., Spahn, F., Sremčević, M., & Salo, H. 2005, Geophys. Res. Lett., 32, 11205 [Google Scholar]
 Spahn F., & Sremčević, M. 2000, A&A, 358, 368 [NASA ADS] [Google Scholar]
 Sremčević, M.,Spahn, F., & Duschl, W. J. 2002, MNRAS, 337, 1139 [Google Scholar]
 Tiscareno, M. S., Burns, J. A., Hedman, M. M., et al. 2006, Nature, 440, 648 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Tiscareno, M. S., Burns, J. A., Hedman, M. M., & Porco, C. C. 2008, AJ, 135, 1083 [Google Scholar]
 Tiscareno, M. S., Burns, J. A., Sremčević, M., et al. 2010, ApJ, 718, L92 [Google Scholar]
Appendix A: Response calculation of particles on horseshoe orbits
A.1. Interaction potential
Due to some finite eccentricity, the moonlet undergoes a small oscillation about the origin. Its Cartesian coordinates then become (X,Y), with these considered small in magnitude. The components of the equation of motion for a ring particle with Cartesian coordinates (x,y) ≡ r_{1} are where the interaction gravitational potential due to the moonlet is (A.3)with the cylindrical coordinates of the particle and moonlet being (r,φ) and (R,Φ), respectively. This may be expanded correct to first order in R/r in the form (A.4)The moonlet undergoes small amplitude epicyclic oscillations such that X = ℰ_{2}cos(Ωt + ϵ),Y = −2ℰ_{2}sin(Ωt + ϵ) where e is its small eccentricity and ϵ is an arbitrary phase. Then Ψ_{1,2} may be written as (A.5)where (A.6)gives the part of the lowest order interaction potential associated with the eccentricity of the moonlet.
We here view the interaction of a ring particle with the moonlet as involving two components. The first term operates when ℰ_{2} = 0 due to the first term on the righthand side of Eq. (A.4) and results in standard horseshoe orbits for ring particles induced by a moonlet in circular orbit. The second term, Eq. (A.6), perturbs this motion when ℰ_{2} is small. We now consider the response of a ring particle undergoing horseshoe motion to this perturbation. In doing so, we make the approximation that the variation in the leading order potential Eq. (A.4) due to the induced ring particle perturbations may be neglected. This is justified by the fact that the response induces epicyclic oscillations of the particle, which are governed by the dominant central potential.
A.2. Response calculation
Setting x → x + ξ_{x},y → y + ξ_{y}, where ξ is the small response displacement induced by Eq. (A.6), and linearising Eqs. (A.1), we obtain the following equations for the components of ξFrom these we find a single equation for ξ_{x} in the form (A.9)When performing the time integral on the righthand side of the above, because we are concerned with a potentially resonant epicyclic response, we only retain the oscillating part.
The solution to Eq. (A.9), which is such that ξ vanishes in the distant past (t = −∞) when the particle is far from the moonlet, may be written as (A.10)where After the particle has had its closest approach to the moonlet and moves to a large distance from it, it will have an epicyclic oscillation with amplitude and phase determined by In evaluating the above we note that, although F vanishes when the particle is distant from the moonlet at large  t  , it also oscillates with the epicyclic angular frequency Ω, which results in a definite non zero contribution. This is the action of the coorbital resonance. It is amplified by the fact that the encounter of the particle with the moonlet in general occurs on a horseshoe libration timescale that is much longer than Ω^{1}. We now consider this unperturbed motion of the ring particles.
A.3. Unperturbed horseshoe motion
The equations governing the unperturbed horseshoe motion are Eqs. (A.1) with (A.15)We assume that this motion is such that x varies on a timescale that is much longer than Ω^{1} so that we may approximate the first of these equations as x = −2/(3Ω)(dy/dt). Consistent with this, we also neglect x in comparison to y in Eq. (A.15). From the second equation we then find that (A.16)which has a first integral that may be written (A.17)where as before b is the impact parameter, or the constant value of  x  at large distances from the moonlet. The value of y for which the horseshoe turns is then given by y = y_{0} = 24Gm_{2}/(9Ω^{2}b^{2}). At this point we point out that we are free to choose the origin of time t = 0 such that it coincides with the closest approach at which y = y_{0}. Then in the approximation we have adopted, the horseshoe motion is such that y is a symmetric function of t.
A.4. Evaluation of the induced epicyclic amplitude
Using Eqs. (A.6), (A.9), and (A.17), we may evaluate the epicyclic amplitudes from Eq. (A.13). In particular we find (A.18)When evaluating Eq. (A.13), consistent with our assumption that the epicyclic oscillations are fast, we average over an orbital period assuming that other quantities in the integrands are fixed, and use Eq. (A.17) to express the integrals with respect to t as integrals with respect to y. We then find α_{∞} = −A_{‘∞}sinϵ and
β_{∞} = −A_{‘∞}cosϵ where (A.19)with r^{2} = 8Gm_{2}/(3Ω^{2}y_{0})(1 − y_{0}/y) + y^{2}. From this we may write (A.20)where ℰ_{1f} is the final epicyclic motion of the ring particle and the dimensionless integral ℐ is given by (A.21)We notice that the dimensionless quantity with r_{H} being the Hill radius of the moonlet. It is related to the impact parameter b by η = 2^{6}(b/r_{H})^{6}. Thus for an impact parameter amounting to a few Hill radii, in a very approximate sense, η is of order unity, y_{0} is of order r_{H}, and the induced eccentricity ℰ_{1f} is of order ℰ_{2}.
All Tables
Eccentricity damping timescale τ_{e} of the moonlet in units of the orbital period.
All Figures
Fig. 1 Trajectories of ring particles in a frame centred on the moonlet. Particles accumulate near the moonlet and fill its Roche lobe. Particles on trajectories labelled (a) are on circulating orbits. Particles on trajectories labelled (b) collide with other particles in the moonlet’s vicinity. Particles on trajectories labelled (c) collide with the moonlet directly. Particles on trajectories labelled (d) are on horseshoe orbits. Particles on trajectories labelled (e) leave the vicinity of the moonlet. 

In the text 
Fig. 2 Shearing box, simulating a small patch of a planetary ring system. The particles that leave an auxiliary box in the azimuthal (y) direction reenter the same box on the other side and get copied into the main box at the corresponding location. All auxiliary boxes are equivalent and there are no curvature terms present in the shearing sheet. 

In the text 
Fig. 3 Snapshots of the particle distribution and the moonlet (blue) in simulation EQ40050DTW. The wake is much clearer than in Fig. 4 as the box size is twice as large. 

In the text 
Fig. 4 Snapshots of the particle distribution and the moonlet (blue) for different surface densities after 25 orbits. The wake created by the moonlet is much longer than the size of the box and therefore hardly visible in these images. 

In the text 
Fig. 5 Impact bands for simulation EQ20050 (left) and EQ20025 (right). The vertical lines indicate sharp cutoffs at 1.5r_{H} and 2.5r_{H}. 

In the text 
Fig. 6 Offset in the azimuthal distance due to the random walk in the semimajor axis a measured in simulations EQ5025, EQ5050, EQ20025, EQ20050, EQ40025 and EQ40050. Individual moonlets may show linear, constant or oscillatory growth. On average, the azimuthal offset grows like t^{3/2}. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.