MIGRATION OF A MOONLET IN A RING OF SOLID PARTICLES: THEORY AND APPLICATION TO SATURN'S PROPELLERS

The theory of disk–satellite interactions has for the most part been developed after Voyager's encounter with Saturn. The satellites that orbit beyond the outer edge of the rings perturb the dynamics of the particles composing the rings. This leads to an exchange of angular momentum between the rings and the satellites and to the formation of density waves in the rings. Using two different methods, Lin & Papaloizou (1979) and Goldreich & Tremaine (1979, 1980) have calculated the total torque exerted by a satellite on a disk interior (or exterior) to its orbit. This torque is called the one-sided Lindblad torque because it can be computed as the sum of the torques exerted at Lindblad resonances with the secondary body. In the lowest order, local approximation, the inner and outer torques are equal and opposite: the torque exerted by a satellite on a disk located inside its orbit is negative, with the same absolute value as the positive torque exerted on a disk located outside the orbit.

Reciprocally, a disk exerts a torque on the secondary body. When the strictly local approximation is relaxed, the inner and outer torques are not exactly equal and opposite (Ward 1986). Their sum, called the differential Lindblad torque, is generally negative. As a consequence, the orbital angular momentum of a body embedded in a disk decreases and so does its semimajor axis (on circular Keplerian orbits, the orbital angular momentum is proportional to the square root of the semimajor axis). This is planetary migration of type I (Ward 1997). So far, this phenomenon has been mainly studied in the frame of planets embedded in protoplanetary gaseous disks (see Papaloizou et al. 2007 for a review).

The Cassini spacecraft has been orbiting the Saturnian ring system since 2004, offering the possibility to observe the coupled evolution of the ring system and the satellites. Due to short orbital timescales (1 year is equivalent to about 700 orbits of the A ring), it may be possible to observe the exchange of angular momentum between the two systems. One of the most striking discoveries of the Cassini spacecraft is the observation of propeller-shaped features in the A ring (located between 122,000 and 137,000 km from Saturn), with a longitudinal extent of about 3 km (Tiscareno et al. 2006, 2008; Sremčević et al. 2007). They are most likely caused by the presence of moonlets about 100 m in size, embedded in the ring, and scattering ring particles (Spahn & Sremčević 2000). As they are embedded in the ring, these small bodies should exchange angular momentum with the ring and migrate (Crida et al. 2009). This migration could be detected by Cassini observations through the cumulative lag, or advance with time t of the orbital longitude ϕ induced by a small variation of the semimajor axis and the angular velocity Ω (δϕ = δΩ × δt), offering for the first time the possibility to directly confront the planetary migration theory with observations, and to give insights and constraints on the physical properties of the rings and the moonlets.

In this paper, we address the question of the theoretical migration rate of these propellers, using both numerical and analytical approaches. The theory is then compared with observations. In Section 2, we review the standard theory of type I migration; the differences between migration in protoplanetary disks and in Saturn's rings are explained, showing the need for a new calculation of the migration rate of embedded moonlets. This rate is given in Section 3 for a homogeneous, axisymmetric disk with a flat surface density profile, as a result of numerical computation in Section 3.1, and analytical calculation in Section 3.2. In Section 4, we consider the effect of density fluctuations in the rings, in particular the role of short-lived gravitating clumps, also called gravity wakes, which are known to be numerous in the A ring (Colwell et al. 2006; Hedman et al. 2007). We then conclude in Section 5 on what our model tells us on the properties of the rings, given the observed migration rates.

In protoplanetary gaseous disks, the perturbation caused by a terrestrial planet leads to the formation of a one-armed spiral density wave, leading the planet in the inner disk, and trailing behind the planet in the outer disk. This wave is pressure-supported and generally called the wake, but it has nothing to do with the gravity wakes mentioned above: the latter are local features, while the planet wake spirals through the whole disk. The planet wake carries angular momentum so that the angular momentum given by the planet to the disk is not deposited locally (e.g., Crida et al. 2006, their Appendix C). Therefore, the disk profile is hardly modified in this linear regime. Still, the negative torque exerted by the outer disk on the planet through the wake is larger in absolute value than the positive torque from the inner disk. Without going into the details (for which the reader is refereed to Ward 1997), the main reason the outer disk wins over the inner disk lies in pressure effects: from the dispersion equation of a pressure-supported wave, one finds that the location r_L,m where the wave with azimuthal mode number m, corresponding to the mth Lindblad resonance with the planet, is launched, is not exactly the location of the resonance given by Kepler's law. The shift is not symmetrical with respect to the planet position for inner and outer resonances, but favors the outer ones. As a consequence, the planet feels a negative total torque, called the differential Lindblad toque, and given by (Tanaka et al. 2002)

$$\begin{equation} T_{\rm diff} = -C q^2 \Sigma {r_p}^4 {\Omega _p}^2 h^{-2}\;, \end{equation} \tag{ 1 }$$

where the index p refers to the planet, r_p being the radius of its orbit and Ω_p its angular velocity, q is the planet to primary mass ratio, and Σ is the surface density of the disk in the neighborhood of the planetary orbit. Finally, h is the aspect ratio of the disk, being the ratio between its scale height and the distance to the central body, being proportional to the square root of the gas pressure. In a typical protoplanetary disk, h ≈ 0.05. The numerical coefficient C is given by C = 2.340 − 0.099ξ, where ξ is the index of the power law of the density profile: Σ ∝ r^−ξ.

This result is robust. In particular, the value of the negative torque is almost independent of the slope of the density profile ξ. This is due to the so-called pressure buffer: the resonances are shifted when the density gradient varies (Ward 1997). If the disk were pressureless, then the expression of C would be completely different. Also, the aspect ratio h in Equation (1) appears because r_L,m does not converge toward r_p when m tends to infinity, but toward r_p(1 ± 2h/3), due to pressure effects. To sum up, the gas pressure plays a fundamental role in type I migration.

It should be mentioned for completeness that, in addition to the differential Lindblad torque discussed above, the horseshoe drag—exerted on the planet by the gas on horseshoe orbits around the planetary orbit—plays a significant role in type I migration (see, e.g., Ward 1991; Masset 2001; Baruteau & Masset 2008; Kley & Crida 2008; Paardekooper & Papaloizou 2009; Paardekooper et al. 2010).

In contrast to gaseous protoplanetary disks, pressure effects in Saturn's rings are not important. The aspect ratio h is on the order of 10⁻⁷. The spiral density waves that are observed in the A ring are gravity-supported, not pressure-supported. Thus, the standard theory of type I migration does not apply. In particular, the spiral planet wake does not appear. The interaction of the moonlet responsible for the propeller structure with the disk is observed to take place within a few hundred kilometers. Resonances with m ≳ 10³ are located within this distance and should play a significant role. However, in the standard type I migration the important resonances have m ∼ 1/h ∼ 10⁷. Therefore, Equation (1) cannot be directly applied to a moonlet in Saturn's rings. A new approach is needed, adapted to the two main characteristics of the problem: the fact that the rings are made of solid particles and the fact that the interaction is taking place very close to the moonlet.

In this section, we compute the interaction between a moonlet and a ring test particle. In this analysis, the gravity of the other ring particles is neglected. This leads to the torque exerted on the moonlet by an initially unperturbed, homogeneous ring. In Section 3.1, the computation is performed numerically. In Section 3.2, it is derived analytically. Both results are in agreement, and a corresponding migration rate for the moonlet is given and discussed in Section 3.3.

From now on, m denotes the mass of the moonlet (and not the order of a resonance). The moonlet, has a circular orbit of radius r_m around the central planet of mass M. The gravitational potential due to the moonlet is Ψ. The radial and azimuthal components of the equation of motion of a ring particle in two-dimensional cylindrical polar coordinates (r, ϕ) are

$$\begin{equation} \frac{d^2 r}{dt^2}-r\left(\frac{d\phi }{dt}\right)^2 = -\frac{\partial \Psi }{\partial r}-\frac{GM}{r^2} \end{equation} \tag{ 2 }$$

$$\begin{equation} {\rm and} \quad r\frac{d^2 \phi }{dt^2}+2\left(\frac{d r}{dt}\right)\left(\frac{d\phi }{dt}\right) = -\frac{1}{r}\frac{\partial \Psi }{\partial \phi }. \end{equation} \tag{ 3 }$$

The angular velocity of the moonlet is $\omega =\sqrt{GM/{r_m}^3}$. Let r₀ be the radius of the initially circular orbit of a test particle and $\Omega =\sqrt{GM/{r_0}^3}$ its angular velocity. We note b = r₀ − r_m is the impact parameter and $\hat{b}=b/r_H$ is the normalized impact parameter, where $r_H=r_m\left(\frac{m}{3M}\right)^{1/3}$ is the Hill radius of the moonlet.

3.1. Numerical Computation of the Ring–Moonlet Interaction

In this section, the numerical integration of the above equations of motion is performed in order to find the trajectories of ring particles in the presence of a perturbing moonlet in the frame corotating with the moonlet. To measure the tiny asymmetry between the inner and the outer part of the ring, the full equations are integrated, without linearization or simplification. A Bulirsch–Stoer algorithm (Press et al. 1992) is used, and a Taylor expansion is performed in the code when necessary to subtract accurately large numbers, in order to achieve machine double precision (10⁻¹⁶). We have checked that the Jacobi constant is conserved to this precision along the trajectories. Examples of obtained trajectories are given in Figure 1.

Zoom In Zoom Out Reset image size

It is well known within the framework of the restricted three-body problem that if $|\hat{b}|$ is small enough the test particle has a horseshoe-shaped orbit, while if $|\hat{b}|$ is larger than ∼2.5, the test particle is circulating, and scattered into an eccentric orbit. This can be seen in Figure 1. We perform many numerical integrations with various b in the case of a moonlet of mass m = 3 ×  10⁻¹²M, starting the particle at an azimuth |ϕ₀| = 3000 r_H/r_m = 0.3 (where the moonlet is at ϕ = 0). This angle is large enough so that at this location the influence of the moonlet is negligible and the orbital parameters of the test particle are not disturbed, as will be checked later. We find that the horseshoe regime occurs for $\hat{b}<1.8$ and the scattered regime occurs for $\hat{b}>2.5$. For $1.774<\hat{b}<2.503$, however, the trajectory approaches the center of the moonlet to within a distance smaller than 0.95 r_H. In that case, if one assumes the moonlet is a point mass, the test particle eventually leaves the Hill sphere, either on a horseshoe or a circulating trajectory, but the outcome changes several times with increasing $\hat{b}$. In the case we are concerned about here, the moonlet most likely almost fills its Roche lobe, and therefore we stop the integration of the trajectory as soon as the distance between the test particle and the moonlet is less than 0.95 r_H, assuming a collision.

The specific orbital angular momentum J = r²(dϕ/dt) of the test particles is computed along the trajectories. Angular momentum is exchanged during the close encounter with the moonlet. For $\hat{b}\geqslant 2.503$, the test particle is scattered onto an eccentric orbit of larger angular momentum than the initial one, which results in a gain in angular momentum. The variation of orbital angular momentum along the trajectory is shown in the bottom panel of Figure 2 for the case $\hat{b}=3$, where the top panel is the trajectory. The difference in angular momentum between the initial circular orbit at $\phi _0=0.3\,{\rm {sgn}}(b)$ and the end of the integration, when |ϕ| = 0.3 again, is noted ΔJ. In the figure, only the interval −0.05 < ϕ < 0.05 is displayed for convenience. Most of the exchange of angular momentum occurs when |ϕ| < 0.01.

Zoom In Zoom Out Reset image size

Figure 3 shows |ΔJ| (top thick curve) as a function of $\hat{b}$, in units of the specific angular momentum of the moonlet J_m = r²_mω. For $0<\hat{b}\leqslant 1.774$, the horseshoe trajectory corresponds to a U-turn toward the central planet, and to a loss of angular momentum for the test particle. More precisely, as for circular orbits J ∝ r^1/2, one expects for such a U-turn $\Delta J/J = \frac{1}{2}\frac{\Delta r}{r}=-b/r_m$; this is indeed the case for $\hat{b}<1.3$. In the case where the test particle collides with the moonlet, we assume that it gives all its orbital angular momentum to the moonlet: ΔJ = r²₀Ω − J_m, so that ΔJ/J ≈ b/(2r_m). This also appears in Figure 3. The opposite holds for b < 0.

Computing ΔJ as a function of b to numerical precision enables us to also compute the difference between the inner and outer disk: δJ(b) = ΔJ(b) + ΔJ(−b). This quantity is small with respect to ΔJ(b), but nonetheless well determined and converged in our simulations: ΔJ(b) + ΔJ(−b) is constant after the encounter to a precision better than 0.5% for all |ϕ|>0.02. This validates our choice of ϕ₀. In Figure 3, the bottom thick dashed curve shows δJ in the same scale as |ΔJ|. We see that δJ>0 for all b>0 and that δJ ≪ ΔJ, with

$$\begin{equation} \delta J/\Delta J \approx 5\times 10^{-4}\,\hat{b} \end{equation} \tag{ 4 }$$

for circulating trajectories, and

$$\begin{equation*} \delta J/|\Delta J| \approx 1.17\times 10^{-4}\,\hat{b} \end{equation*}$$

for horseshoe orbits. In the following subsection, the empirically found Equation (4) is analytically derived and justified.

Zoom In Zoom Out Reset image size

3.2. Analytic Model for the Ring–Moonlet Interaction

In this section, we consider only circulating trajectories. Developing the exchange of angular momentum during an encounter with the moonlet ΔJ to second order, we can find the asymmetry δJ.

3.2.1. Solution for the Perturbed Moonlet Orbit

Let us start again from Equations (2) and (3). The ring particle is assumed to be on an unperturbed circular orbit of radius r₀ = r_m + b. It orbits with angular velocity $\Omega = \sqrt{GM/{r_0}^3}$ such that ϕ = Ωt, where without loss of generality we have defined the origin of time t = 0 to be when the particle is at ϕ = 0. Under the perturbation induced by Ψ, the particle moves to r = r₀ + x, and ϕ = Ωt + y/r_m, where x and y are assumed to be small. Linearizing Equations (2) and (3) about the circular orbit state, we obtain equations for x and y in the form

$$\begin{equation} \frac{d^2 x}{dt^2}-2\Omega \frac{d y}{dt} -3\Omega ^2 x={-}\left. \frac{\partial \Psi }{\partial r}\right|_0 \ \ \ {\rm and} \ \ \ \end{equation} \tag{ 5 }$$

$$\begin{equation} \frac{d^2 y}{dt^2}+2\Omega \frac{d x}{dt}={-}\left. \frac{1}{r_0}\frac{\partial \Psi }{\partial \phi }\right|_0. \end{equation} \tag{ 6 }$$

Here, the subscript 0 denotes evaluation on the unperturbed particle orbit.

We suppose that the perturbation is induced by a moonlet of mass m that is on a circular orbit of radius r_m and has an angular velocity ω. Then its azimuthal coordinate ϕ_m = ωt + ϕ_m,0, with ϕ_m,0 being a constant. The perturbing potential

$$\begin{equation} \Psi = \frac{Gm}{\sqrt{r_0^2+r_m^2-2r_mr_0\cos (\phi -\phi _m)}} \end{equation} \tag{ 7 }$$

becomes a function of time through substituting ϕ − ϕ_m = (Ω − ω)t − ϕ_m,0 therein.

Thus, we have

$$\begin{equation} \left. \frac{1}{r_0}\frac{\partial \Psi }{\partial \phi }\right|_0 \equiv \left. \frac{1}{r_0(\Omega -\omega)}\frac{\partial \Psi }{\partial t}\right|_0. \end{equation} \tag{ 8 }$$

Using this in Equation (6) and integrating with respect to time, we obtain

$$\begin{equation} \frac{d y}{dt}+2\Omega x = {-}\left. \frac{\Psi }{r_0(\Omega -\omega)}\right|_0, \end{equation} \tag{ 9 }$$

which when combined with Equation (5) gives an equation for x in the form

$$\begin{equation} \frac{d^2 x}{dt^2} + \Omega ^2x = -\left.\left(\frac{\partial \Psi }{\partial r} +\frac{2\Omega \Psi }{r_0(\Omega -\omega)}\right) \right|_0=S. \end{equation} \tag{ 10 }$$

3.2.2. Solution of the Linearized Equations

To solve Equation (10), we note that the perturbing potential Equation (7) evaluated on the unperturbed orbits is a periodic function of time with period 2π/|ω − Ω| = 2π/β. Thus, we should look for a periodic response. In order to do this, we have to introduce a small frictional term into Equation (10) to enable transients to decay and a net torque on the moonlet to be set up. When the frictional term is small, it is expected that the resulting torque should not depend on it (e.g., Goldreich & Tremaine 1980). Hence we add a frictional term γ(dx/dt) to the left-hand side of Equation (10), where γ/Ω is a small constant parameter so that it now reads

$$\begin{equation} \frac{d^2 x}{dt^2} + \gamma \frac{dx}{dt} +\Omega ^2 x=-\left.\left(\frac{\partial \Psi }{\partial r} +\frac{2\Omega \Psi }{r_0(\Omega -\omega)}\right) \right|_0=S. \end{equation} \tag{ 11 }$$

As the potential is periodic in time we can adopt a Fourier series of the form

$$\begin{equation} S = \sum _{n=-\infty }^{n=\infty }S_n\exp (in\beta t)\, \end{equation} \tag{ 12 }$$

where it is implicit that the real parts of such complex expressions are to be taken, and

$$\begin{equation} S_n = \frac{1}{2\pi }\int _0^{2\pi /\beta } S(t)\exp (-in\beta t)dt{.} \end{equation} \tag{ 13 }$$

The periodic solution of Equation (11) is now readily written down as

$$\begin{equation} x = \sum _{n=-\infty }^{n=\infty }\frac{S_n\exp (in\beta t)}{(\Omega ^2-n^2\beta ^2+i\gamma n \beta)}. \end{equation} \tag{ 14 }$$

We may write this in terms of a Green's function defined through

$$\begin{equation} G(\tau) =\frac{1}{2\pi } { \sum _{n=-\infty }^{n=\infty }\frac{\exp (in\beta \tau)}{(\Omega ^2-n^2\beta ^2+i\gamma n \beta)}}. \end{equation} \tag{ 15 }$$

Then the solution for x may be written as

$$\begin{equation} x =\beta \int _0^{2\pi /\beta } S(t-t^{\prime }) G(t^{\prime })dt^{\prime }. \end{equation} \tag{ 16 }$$

Note that as the orbit of a ring particle relative to the planet is periodic, the solution given by Equation (16) includes the effects of infinite numbers of repeating encounters. However, we wish to consider the case when dissipative effects, although weak, are strong enough to recircularize orbits between encounters in which case they will be independent of each other. This condition requires that γ/|ω − Ω| = γ/β ≫ 1. This is equivalent to requiring that the damping timescale be short compared to the relative orbital period between moonlet and ring particle. On account of the length scale of the encounters of interest being comparable to the Hill radius of the moonlet, this is much longer than the orbital period itself, so that we may adopt the ordering

$$\begin{equation} \gamma /|\omega -\Omega |=\gamma /\beta \gg 1 \gg \gamma /\Omega. \end{equation} \tag{ 17 }$$

In order to make use of the above ordering, we write down the form of the Green's function derived in the Appendix (see Equation (A7)) valid for 0 < t < 2π/β:

$$\begin{equation} G(\tau) = \frac{e^{-\gamma \tau /2}\sin (\omega _{\gamma }\tau)- e^{-\gamma \pi /\beta }\sin (\omega _{\gamma }(\tau -2\pi /\beta))}{\omega _{\gamma }\beta [1+e^{-2\gamma \pi /\beta }-2e^{-\gamma \pi /\beta }\cos (2\pi \omega _{\gamma }/\beta)]}\, \end{equation} \tag{ 18 }$$

where $\omega _{\gamma }=\sqrt{\Omega ^2-\gamma ^2/4}.$ The function is defined elsewhere through its periodicity with period 2π/β. Making use of the inequality Equation (17), we may replace the Green's function (Equation (18)) by the simple expression

$$\begin{equation} G(\tau) = \frac{\exp (-\gamma \tau /2)\sin (\omega _{\gamma }\tau)}{\omega _{\gamma }\beta }{.} \end{equation} \tag{ 19 }$$

Then the solution (Equation (16)) gives

$$\begin{equation} x = A(t) e^{-\gamma t/2}\sin (\omega _{\gamma } t) + B(t) e^{-\gamma t/2}\cos (\omega _{\gamma } t), \end{equation} \tag{ 20 }$$

where

$$\begin{equation} A(t)= \frac{1}{\omega _{\gamma }}\int ^t_{t-2\pi /\beta }S(t^{\prime }) \exp (\gamma t^{\prime }/2) \cos (\omega _{\gamma } t^{\prime })dt^{\prime } \end{equation} \tag{ 21 }$$

and

$$\begin{equation} B(t)= - \frac{1}{\omega _{\gamma }}\int ^t_{t-2\pi /\beta }S(t^{\prime }) \exp (\gamma t^{\prime }/2) \sin (\omega _{\gamma } t^{\prime })dt^{\prime }. \end{equation} \tag{ 22 }$$

To make use of the above expressions, we consider the situation when the ring particle has a close encounter with the moonlet at time t = 0, thus we take ϕ_m,0 = 0 (we note that a non-zero ϕ_m,0 can be dealt with by rotating the coordinate system and shifting the origin of time). The source term S is then expected to be highly peaked around t' = 0, and almost all of the contributions to the above integrals will occur for |t'| ≲ 2π/Ω. Furthermore, during this dynamical interaction, dissipation will be negligible. Thus, if we are interested in times after the main interaction, but before significant dissipation takes place, we may set γ = 0 and extend the limits of the integration to ±∞. However, in practice one may have to apply a cutoff to the potential at large distances from the moonlet in order to do that (see below). But this should not matter if the important interaction occurs when the moonlet and ring particle are close.

Then we simply have

$$\begin{equation} A = \frac{1}{\Omega }\int ^{\infty }_{-\infty }S(t^{\prime }) \cos (\Omega t^{\prime }) dt^{\prime } \end{equation} \tag{ 23 }$$

and

$$\begin{equation} B = - \frac{1}{\Omega }\int ^{\infty } _{-\infty }S(t^{\prime }) \sin (\Omega t^{\prime }) dt^{\prime }. \end{equation} \tag{ 24 }$$

Thus, A and B are constants representing epicyclic oscillation amplitudes induced after the close approach of the ring particle to the moonlet.

We remark that the approximations made in obtaining Equations (23) and (24) relate to how dissipation is treated. There has been no assumption that the particle trajectories are symmetric on opposite sides of the moonlet so that curvature effects remain fully incorporated during particle moonlet encounters. When dissipation is negligible during the encounter, an epicyclic oscillation is established immediately afterward. The assumption that dissipation circularizes orbits between encounters implies that we should consider approaching ring particles to be on circular orbits. The above discussion indicates that errors associated with this assumption are exponentially small.

3.2.3. Angular Momentum Transfer

For the setup considered here, symmetry considerations imply that S(t) is an even function of time (see below) so that B = 0. The generation of the epicyclic oscillation is associated with an angular momentum transfer between the moonlet and particle. To find this, we first note that

$$\begin{equation} \Delta J = \sqrt{GM}\Big(\sqrt{a_f(1-e^2)}-\sqrt{r_0}\Big), \end{equation} \tag{ 25 }$$

where a_f and e are the post-encounter semimajor axis and eccentricity of the particle, respectively. We also note that the Jacobi constant implies that the change of the particle orbital energy and angular momentum are related by ΔE = GM[1/(2r₀) − 1/(2a_f)] = ω ΔJ. This can be used to eliminate a_f in Equation (25) after which ΔJ may be found correct to second order in e ≡ A/r₀ with the result that ΔJ = Ω²A²/[2(ω − Ω)]. This in turn may be simply determined after evaluating A. Note that ΔJ < 0 for particles interior to the moonlet which have Ω>ω and conversely ΔJ>0 for particles orbiting exterior to the moonlet.

3.2.4. Development of the Perturbing Potential

We now consider

$$\begin{equation} S(t) =-\left.\left(\frac{\partial \Psi }{\partial r} +\frac{2\Omega \Psi }{r_0(\Omega -\omega)}\right) \right|_0. \end{equation} \tag{ 26 }$$

We begin by recalling that

$$\begin{equation*} \Psi = - \frac{Gm}{\sqrt{r_0^2+r_m^2-2r_0r_m\cos (\phi -\phi _m)}} \end{equation*}$$

$$\begin{equation} = - \frac{Gm}{\sqrt{r_0^2+r_m^2-2r_0r_m\cos (\beta t)}}. \end{equation} \tag{ 27 }$$

In order to evaluate the Fourier transform as specified by Equation (23), which was derived under the assumption that the interaction occurs only near the closest approach, we must truncate the potential at large |t|. As the encounter takes place over a time ≪1/β, this can be achieved by replacing cos(βt) in Equation (27) by 1 − β²t²/2. Note that a dimensionless estimate of the error involved is of order (β/ω)² ∼ (r_H/r_m)², where r_H is the Hill radius of the moonlet. This is small enough that the leading order asymmetry in the angular momentum transferred to orbits with the same impact parameter on either side of the disk can be estimated.

As the first stage in evaluating the Fourier transform of S specified in Equation (23) that gives the epicyclic amplitude, we evaluate

$$\begin{equation*} C = \frac{1}{\Omega }\int ^{\infty }_{-\infty } \Psi \cos (\Omega t)dt \end{equation*}$$

$$\begin{equation} = - \frac{1}{\Omega }\int ^{\infty }_{-\infty } \frac{Gm\cos (\phi) }{\frac{\beta }{\Omega }\sqrt{(r_0-r_m)^2\Omega ^2/\beta ^2+r_0r_m\phi ^2}}\frac{d\phi }{\Omega }{.} \end{equation} \tag{ 28 }$$

This can also be expressed as

$$\begin{equation} C = -\frac{2 GmK_0(\xi _0)}{\Omega \beta \sqrt{r_0r_m}}, \end{equation} \tag{ 29 }$$

where $\xi _0 = (\Omega |r_0-r_m|)/(\beta \sqrt{r_0r_m})$ and K_j denotes the modified Bessel function of the second kind of order j.

3.2.5. Total Angular Momentum Exchange

We may now use the above expression together with Equation (26) to evaluate the epicyclic amplitude Equation (23) (noting that the radial derivative is with respect to r₀ with other quantities held fixed) so obtaining

$$\begin{eqnarray} A &=& -\frac{2 Gm}{\Omega \beta r_0\sqrt{r_0r_m}} \times \left(K_0(\xi _0)\left[\frac{1}{2}-\frac{2\Omega }{(\Omega -\omega)}\right]\right.\nonumber \\ &&\left. +\ K_1(\xi _0)\xi _0\left[\frac{1}{2} +\frac{r_m}{(r_0-r_m)}\right]\right). \end{eqnarray} \tag{ 30 }$$

The associated angular momentum exchanged is then given by

$$\begin{equation*} \Delta J = \frac{2 (Gm)^2}{r_0^3r_m(\omega -\Omega)^3} \times \left(K_0(\xi _0)\left[\frac{1}{2}-\frac{2\Omega }{(\Omega -\omega)}\right]\right. \end{equation*}$$

$$\begin{equation} \hfill \left. + K_1(\xi _0)\xi _0\left[\frac{1}{2} +\frac{r_m}{(r_0-r_m)}\right]\right)^2. \end{equation} \tag{ 31 }$$

In a strictly local approximation under which the inner and outer sides are symmetric, the contributions from orbits equidistant from the moonlet would cancel, leaving the net result to be determined by the surface density profile. However, although we have assumed the interactions are local, we did not assume symmetry between the exterior and interior orbits. Accordingly, we evaluate the difference in the magnitude of ΔJ evaluated from orbits equidistant from the moonlet: r₀ = r_m ±  b. The leading order contribution to ΔJ is symmetric in b. The lowest order contribution is antisymmetric and accordingly leads to cancellation between the two sides. We make use of the expansions ξ₀ =  2/3 − b/(2r_m) + O((b/r_m)²) and 2Ω/(Ω − ω) = −4r_m/(3b)(1 − b/(4r_m)) + O(b/r_m) together with standard properties of Bessel functions to write

$$\begin{equation} \Delta J = \frac{64 (Gm)^2r_m}{243\omega ^3b^5} \left(2K_0(2/3)+K_1(2/3)\right)^2 \left(1+\alpha \frac{b}{r_m} \right), \end{equation} \tag{ 32 }$$

where

$$\begin{equation} \alpha = \frac{3}{4}+ \frac{(6K_1(2/3)+3K_0(2/3))}{(4K_0(2/3)+2K_1(2/3)}=2.46. \end{equation} \tag{ 33 }$$

The first-order term of Equation (32) was already given by Goldreich & Tremaine (1980). It is plotted as a straight green dashed line in Figure 3. Our expansion to second order enables us to go further, and to give the expression of the magnitude of the asymmetry between the two sides of the disk:

$$\begin{equation} \frac{\delta J}{\Delta J} = 2 \alpha |b|/r_m = 4.92|b|/r_m{.} \end{equation} \tag{ 34 }$$

It is such that for an orbit with a given impact parameter, the angular momentum exchanged in the outer disk is the larger one.

In the case studied numerically, we had r_H = 10⁻⁴, so that $|b|/r_m=10^{-4}\hat{b}$. Then, Equation (34) remarkably agrees with the numerical fit (Equation (4)). The light blue dot-dashed curve in Figure 3 displays $4.92\times 10^{-4}\,\hat{b}\,\Delta J$.

In the context of the above, we note that approximations made in obtaining Equation (31) such as effectively starting and truncating the interaction at some finite though large distance from the moonlet could conceivably lead to changes comparable to those given by Equation (34). However, such changes are again approximately symmetric for trajectories on both sides of the moonlet and thus approximately cancel so we do not expect such effects to significantly alter Equation (34).

3.3. Migration Rate and Discussion

If the surface density of ring particles is Σ, the total rate of angular momentum transferred to the moonlet is

$$\begin{equation} \frac{dJ}{dt}= -\int \!\!\int _{\rm disk}\Sigma \ \Delta J\, \frac{|\omega -\Omega |}{2\pi }dr\,rd\phi, \end{equation} \tag{ 35 }$$

where the integral is taken over the disk. The particles exterior to the moonlet contribute negatively while those interior to the moonlet contribute positively. The cumulative torque exerted by the moonlet on the region of the ring located within a distance b to its orbit then reads

$$\begin{equation} T_c(b) = \int _{-b}^b \Sigma (r_m+b^{\prime }) (\Delta J(b^{\prime })) |\omega -\Omega | db^{\prime }. \end{equation} \tag{ 36 }$$

The normalized cumulative torque

$$\begin{equation*} T_c(b)\big/\big[(m/M)^{4/3}\big(\Sigma /Mr_m^{\,-2}\big)\big] \end{equation*}$$

is plotted in Figure 4. The proportionality to Σ is obvious; that T_c ∝ (m/M)^4/3 is numerically verified for 3 × 10⁻¹⁵ ⩽ m/M ⩽ 3 × 10⁻⁹, and has been already found analytically by Ward (1991) for the horseshoe drag in a similar context.

Zoom In Zoom Out Reset image size

Most of the total torque comes from scattered, circulating particles, in particular the ones with the smallest impact parameter $\hat{b}\approx 2.5$. This makes the total torque sensitive to the physical size of the moonlet (taken as 0.95 r_H here), as some particles colliding with the moonlet could be circulating if it were smaller.

The role of the horseshoe drag appears to be non-negligible, amounting to ∼4.1(Σ/Mr⁻²_m)(m/M)^4/3 Mr²_mω². The expression of Ward (1991) for the torque arising from material executing horseshoe turns, called the horseshoe drag, is for a Keplerian disk with a flat density profile:

$$\begin{equation} T_{\rm HS} = \frac{9}{8}\Sigma {w}^4\omega ^2, \end{equation} \tag{ 37 }$$

where w is the half-width of the horseshoe region. In our case, w = 1.774 r_H, which gives T_HS = 2.6(Σ/Mr⁻²_m)(m/M)^4/3Mr²_mω². The agreement is good because Ward's analysis is based only on geometrical effects and angular momentum variation in a Keplerian disk, without any pressure effect. Therefore, it also applies in Saturn's ring. We remark that taking w = 2r_H in Equation (37) gives a perfect match with what we find numerically for the total horseshoe drag.

In conclusion, from Figure 4, the total torque felt by a moonlet of mass m on a circular orbit of radius r_m around a planet of mass M can be written as

$$\begin{equation} T = - 17.8\,\left(\frac{\Sigma }{Mr_m^{\,-2}}\right)\left(\frac{m}{M}\right)^{4/3} Mr_m^{\ 2}\omega ^2{.} \end{equation} \tag{ 38 }$$

Note that to get the same dependency of the type I torque in the parameters of the system, one has to assume h ∝ r_H/r_m in Equation (1); however, in a protoplanetary disk, h is fixed and independent of the mass of the secondary body, so that this proportionality would not be justified.

The torque is related to the migration speed through T = 0.5 m r_mΩ(dr_m/dt). Hence we deduce that

$$\begin{equation} \frac{dr_m}{dt} = - 35.6\,\frac{\Sigma r_m^{\ 2}}{M}\left(\frac{m}{M}\right)^{1/3}\ r_m\Omega {.} \end{equation} \tag{ 39 }$$

Here the migration rate is proportional to the mass of the moonlet to the power 1/3, in contrast to standard type I migration where dr_m/dt ∝ m. A numerical application to the case of an m = 10⁻¹⁸M_Saturn = 5.68 × 10⁸ kg moonlet in orbit in the A ring of density Σ = 400 kg m⁻² at r_m= 130,000 km from Saturn gives dr_m/dt = −0.23 m yr⁻¹. Increasing the mass by two orders of magnitude to correspond to a radius of ∼200 m speeds up the migration rate by a factor of only ∼4.5 to ∼−1 m yr⁻¹.

After time t, a migrating propeller will be shifted longitudinally with respect to a corresponding non-migrating one by a distance r_m δϕ = 3Ω |dr_m/dt| t²/4. For the above parameters, this gives 713 [t/(1 yr)]² m. A shift of this magnitude is potentially detectable on a timescale of a year to a few years (Porco et al. 2004). Actually, migration of propellers has already been detected: during one time period of nearly a year, a particular propeller has been seen moving outward at a rate of ∼110 m yr⁻¹, and during a later similar time period, the same propeller has been seen moving inward at a rate of ∼40 m yr⁻¹ (Tiscareno et al. 2010).

These observations are not compatible with the above theory. But we recall that the process of migration for a moonlet described above assumed a smooth particle disk with constant surface density. Here we note that there are features and mechanisms that might produce a significantly faster migration rate, possibly in both directions inward and outward, and non-constant in time. One can first think of a radial density gradient: as there is no pressure buffer here, this would directly affect the balance between the torques from the inner and outer parts of the ring. This would also affect the torque from the horseshoe region, which could turn positive. However, if the migration is governed by the gradient of some quantity, it seems likely that the moonlet would have approached an extremum in that quantity, and thus should have attained a migration rate comparable to that estimated for a constant surface density.

Another possibility resulting in the moonlet migrating faster than what the previous calculation indicates, and possibly outward, is a runaway migration in a planetesimal disk (see Ida et al. 2000 and Levison et al. 2007, for a review), similar to the type III migration of planets in protoplanetary disks (Masset & Papaloizou 2003). In this regime, the migration of the moonlet in the disk leads to a positive feedback on its migration rate, because of the material of the inner (resp. outer) disk making horseshoe U-turns to the outer (resp. inner) disk. This speeds up the migration, possibly leading to a runaway. However, this inevitably leads to an asymmetry in the horseshoe region, while the propeller structures observed are rather symmetrical.

Finally, the A ring of Saturn is not homogeneous. It is close to gravitational instability, which should lead to the formation of gravity wakes and density fluctuations. The effect of these density fluctuations on the moonlet is studied in the following section.

The analytic calculations and numerical simulations in the previous chapters assume an inflow of particles on circular orbits only perturbed by the nearby moonlet. However, we know that Saturn's A ring is marginally gravitationally stable (Daisaka et al. 2001). The Toomre Q parameter (Toomre 1964), which is a measure of the importance of self-gravity, is expected to be on the order of 2–7, indicating that the ring particles' mutual gravity is indeed a strong effect. It leads to the regular formation and dispersion of gravity wakes, which are local density enhancements elongated in parallel directions by the Keplerian shear. Those overdensities give rise to stochastic forces which act on the embedded moonlet.

A very similar effect is expected to occur in protoplanetary disks. These disks are thought to be turbulent due to their magnetorotational instability (Balbus & Hawley 1991). The turbulent fluctuations create overdensities which interact gravitationally with embedded small-mass planets. The stochastic forces make the planet undergo a random walk. An analytic model of this random walk has been derived by Rein & Papaloizou (2009). In the following, we apply this model to moonlets embedded in Saturn's rings. To do that, we need to get an estimate of the amplitude of the stochastic forces.

Here, we only calculate a first estimate of these forces. A comprehensive parameter space survey of all possible ring parameters is beyond the scope of this work and discussed in a separate paper (Rein & Papaloizou 2010).

4.1. Numerical Calculations

We perform three-dimensional simulations of ring particles, in a shearing box, similar to Salo (1995). The simulations are done in a local cube of size H with shear periodic boundary conditions, and the origin of the box is fixed at a semimajor axis of a= 130,000 km. A BH tree code (Barnes & Hut 1986) is used to calculate the self-gravity between ring particles and resolve inelastic collisions. Collisions between particles are resolved using the instantaneous collision model and a velocity-dependent coefficient of restitution given by (Bridges et al. 1984)

$$\begin{equation} \epsilon (v) = \mathrm{min}\left\lbrace 0.34\times \left(\frac{v}{1\,\mathrm{cm\;s}^{-1}} \right)^{-0.234}, 1\right\rbrace, \end{equation} \tag{ 40 }$$

where v is the impact speed projected on the vector joining the centers of the two particles. The code is described in more detail in Rein et al. (2010).

The size of ring particles is not well constrained. Therefore, and to be able to scale to different locations in Saturn's rings, we perform multiple simulations. For a given simulation, all the particles are assumed to be spherical and have the same size (or radius), which varies from simulation to simulation from 0.52 to 13 m. The simulation parameters are listed in Table 1. The nomenclature and physical parameters are, for easy comparison, the same as in Lewis & Stewart (2009), as our simulations are similar to theirs.

Table 1. Simulation Parameters

Name	ra (m)	ρp (g cm−3)	Σp (kg m−2)	H (m)	N
L2	13	0.5	885	5000	4,808
S1	0.52	0.7	49.7	1000	120,548
S3	1.3	0.7	246	1000	38,188

Notes. The first column identifies the simulation, following the convention of Lewis & Stewart (2009). The second and third columns give the size (or radius) of the particles and their density. The fourth and fifth columns give the surface density and the size of the computational domain, respectively. The last column lists the number of particles.

Download table as:  ASCIITypeset image

The moonlet is not taken into account in the simulations. We measure the specific gravitational force ${\hat{\bf f}}$ (or acceleration) felt by a passive test particle sitting at the origin. We calculate the force in two different ways in order to avoid the singularity at the origin and to account for the physical size of the moonlet. In the first case, we use a cutoff at the moonlet's radius d and exclude all particles within that radius from the force calculation. In the second case, we use a smoothed gravitational force per unit mass in the form

$$\begin{equation} {\hat{\bf f}} = -\frac{Gm_{\rm part}}{|{\hat{\bf r}}|^2+d^2} {\hat{\bf r}}, \end{equation} \tag{ 41 }$$

where ${\hat{\bf r}}$ is the vector linking the origin to the particle and m_part is the mass of the particle. The smoothing length d is set equal to the moonlet's size. In a self-consistent simulation, one should include the moonlet with its real physical size. However, this goes beyond the scope of this paper and will be considered in future work (Rein & Papaloizou 2010). Our purpose here is to estimate the underlying stochastic fluctuations in the migration rate that occur independently of the moonlet. This procedure is reasonable as long as the moonlet is in a steady state, namely if it does not accumulate or lose a large amount of mass over one orbit. In all our simulations we assume a moonlet size of d = 200 m.

4.2. Results

We measure the amplitude and the correlation time of the stochastic forces in all simulations. The results are listed in Table 2. We also plot the time evolution of the azimuthal (y) force component in Figure 5. Whereas the correlation time in all simulations is almost the same, the amplitude of the stochastic fluctuations varies by almost a factor of 10³. The forces in the vertical direction are negligible and are not presented here. The diffusion coefficient, being a measure of the strength of stochastic forces, is defined as D = 2〈f²〉τ (see Rein & Papaloizou 2009), where 〈f²〉^1/2 and τ are the root mean square value and the approximate correlation time, respectively, of the specific stochastic force in one direction.

Zoom In Zoom Out Reset image size

Table 2. Simulation Results

Simulation		Correlation Time (s)		Diffusion Coefficient (m2 s−3)		Q
		τx	τy	Dx	Dy
L2	Cutoff	5000	7000	9.61 × 10−12	14.97 × 10−12	4.1
	Smooth	6000	9000	5.34 × 10−12	9.11 × 10−12
S1	Cutoff	2000	4000	1.92 × 10−17	2.26 × 10−17	7.2
	Smooth	3000	8000	1.59 × 10−17	1.69 × 10−17
S3	Cutoff	2000	6000	7.50 × 10−16	30.24 × 10−16	2.5
	Smooth	4000	10000	8.87 × 10−16	20.48 × 10−16

Notes. The first column gives the name of the simulation as defined in Table 1 and the method used in the computation of the force. The second and third columns give the correlation time in the x (radial) and y (azimuthal) directions, respectively. The fourth and fifth columns list the diffusion coefficients. The sixth column is the Toomre Q parameter as measured in the simulation.

Download table as:  ASCIITypeset image

The change in the semimajor axis a due to the effect of stochastic forces with diffusion coefficient D after time t is given by

$$\begin{eqnarray} \Delta a &=& \frac{2}{\omega } \sqrt{Dt}, \end{eqnarray} \tag{ 42 }$$

where ω is the mean motion of the moonlet (Rein & Papaloizou 2009). Note that in this regime the acceleration of the moonlet does not depend on its mass (as can be seen in Equation (41)). Therefore, the migration rate and the diffusion coefficient are independent of the mass of the moonlet; this might be an observational indication for this migration regime.

Assuming an initial semimajor axis of a = 130,000 km and D ∼ 10⁻¹⁷ m² s⁻³ as found in simulation S1, one can calculate the expected difference in the semimajor axis after one orbit due to stochastic forces which turns out to be Δa = 0.01 m. For simulation S3, assuming D ∼ 10⁻¹⁵ m² s⁻³, one finds Δa = 0.11 m. For the simulation L2, assuming D ∼ 10⁻¹¹ m² s⁻³, one finds Δa = 10.5 m. These translate to random walks with standard deviation given as a function of time by $\Delta a = 0.27 \sqrt{t/(1\;{\rm yr})}$ m, $\Delta a = 2.7 \sqrt{t/(1\;{\rm yr})}$ m, and $\Delta a = 270 \sqrt{t/(1\;{\rm yr})}$ m, respectively.

The results might depend on a variety of properties of the ring particles such as the coefficient of restitution, internal density, and size distribution. However, our results can be seen as an upper limit on the strength of the stochastic forces.

We use the coefficient of restitution given by Equation (40). Additional simulations have shown that reducing the coefficient of restitution makes the system gravitationally unstable forming persistent aggregates. One could circumvent this fate by decreasing the internal density of particles, so that they are well within their Roche limit. And indeed, the internal densities used in our simulations are rather high (compare to Porco et al. 2008). However, because the simulations are close to forming persistent clumps, Q is close to unity and the gravitational wakes that occur are as strong as they can possibly get. Thus, the forces are as large as they can possibly get in a steady-state ring of a fixed surface density.

4.3. Discussion

In this paper, we have calculated the differential torque exerted on a moonlet by the outer and the inner disk with a smooth, flat surface density profile. We performed both an accurate numerical integration and a second-order analytical calculation. These approaches were found to be in excellent agreement where their domains of validity overlap. The migration rate found in this case is proportional to the mass of the moonlet to the power 1/3. It is about −1 m yr⁻¹ for a 200 m radius moonlet in the A ring. This is way too low to explain the observed migration of the propellers in Saturn's rings. Nonetheless, density fluctuations in the rings, due to their proximity to gravitational instability, can lead to stochastic torques on a moonlet, which may dominate on the timescale of the Cassini mission, to an extent that depends mainly on the surface density of the rings. These stochastic torques may account for the observations.

The possibility that the migration of propellers is induced by stochastic processes rather than by a regular type-I-like migration is therefore very exciting: this may help to infer the local surface density and therefore the size of the ring particles. Indeed, both quantities are linked through the optical depth, which is observationally well constrained. Our estimate of Σ ∼ 700 kg m⁻² is equivalent to ∼10 m size particles in the A ring for a single-sized population. This is in good agreement with previous estimates of the surface density from the structure of the gravity waves (Tiscareno et al. 2007, and references therein). For the A ring, radio occultations of the Voyager spacecraft (Zebker et al. 1985) find that the radius r_a of particles follows: 0.1 m (assumed) <r_a < 11 m, with a power index ∼−3. For the 28 Sgr occultation (French & Nicholson 2000), a range 1 m <r_a < 20 m is found, with a power index between −2.7 and −3, and an effective size (the single average size accounting for the fluctuations in photon count) being about 7 m (Cuzzi et al. 2009).

Fortunately, the Cassini mission has been extended to 2017. In the meantime, numerous observations of the propellers will hopefully give a clear picture of their orbital evolution for a period of 10 years, representing about 10,000 orbits. This will allow us to test the hypothesis presented in this paper, and in particular to check whether the propellers really are in stochastic migration, whereas first results seem to favor the stochastic hypothesis.

Among the questions that still need to be addressed is why all the propellers seem to be gathered in the A ring, in places apparently devoid of density waves. Indeed the propellers seem gathered in three narrow radial bands of only about 1000 km width (Tiscareno et al. 2008). This is especially surprising since the A ring is densely populated by numerous density waves launched by the nearby small moons (Atlas, Prometheus, Pandora, Janus, and Epimetheus). Is there a systematic mechanism that would eject the propellers away from density waves? Or does the present location of propellers just reflect the initial location of the parent body, assuming that the propeller population comprises the fragments resulting from the destruction of an ancient moon orbiting within the rings?

If the moonlets really undergo stochastic migration, then Equation (42) may strongly constrain the age of the propellers, which cannot be larger than the time needed to diffuse over Δa>1000 km. Unfortunately, as long as D is unknown Equation (42) does not provide any useful information. However, considering that Δa is proportional to the square root of the time, and assuming that a moonlet migrates about 100 m in 1 year, one finds that Δa = 1000 km for t = 100 million years. Assuming |Δa| = 40 m in 1 year, we find that it requires ∼625 million years to diffuse over 1000 km. Note that for 200 m radius moonlets, the spreading due to stochastic migration equates to the contraction of the semimajor axes occurring as a result of smooth, regular migration after ∼10⁴ yr, and then Δa ≈ 10 km. Thus, it would take about a million years to migrate through Δa = 1000 km in this case, largely through the action of the non-stochastic, regular migration process, if that can be assumed to operate smoothly and simultaneously with the stochastic migration process. A 100 m radius moonlet would migrate through only 500 km in the same period, so that the smooth, regular migration process spreads a population of moonlets of various sizes over 1000 km in one to two million years. These could be the times since the catastrophic disruption of a small moon orbiting at 130,000 km from Saturn that was broken into smaller moonlets by a meteoritic impact. On the other hand, an estimate of the lifetime of a Pan size moon (∼14 km in radius) against today's cometary flux is provided by Dones et al. (2009) and gives a range between 100 Myr and 16 Gyr, depending on the size distribution of impactors. Therefore, the recent occurrence of such an event, about 4–4.5 billion years after solar system formation, is possible. Note also that an age of about 100 Myr is coherent with some estimates of Saturn's ring age despite the lack of fully satisfactory explanation for their origin (see Charnoz et al. 2009a for a review).

We see that the question of a propeller's migration is inextricably linked to the issue of the origin of Saturn's moons embedded in the rings, which is still a mystery. Porco et al. (2007) and Charnoz et al. (2007) have jointly proposed that small moons embedded in the rings could be aggregates of material on an initial shard denser than ice. When destroyed by meteoritic bombardment, these could release dense chunks of material that could explain the origin of the propellers. However, the origin of Saturn's ring system is still a matter of debate (Harris 1984; Charnoz et al. 2009b). Knowledge of the age of the propellers could provide important constraints on the age of the main ring system and its embedded moons, as there are strong indications that these could have about the same age, provided these moonlets hide a dense shard (Charnoz et al. 2007; Porco et al. 2007). Understanding the migration rate of the propellers is therefore an important piece of this puzzle.

We thank J. Burns for stimulating discussions and the organizers of the "Dynamics of Discs and Planets" workshop at the Isaac Newton Institute in Cambridge where these took place, as well as M. Tiscareno for providing us with migration rates. Hanno Rein was supported by an Isaac Newton Studentship, STFC, and St John's College, Cambridge. We also thank an anonymous referee for providing a speedy and constructive report.

Here we evaluate the Green's function defined by Equation (15) as

$$\begin{equation} G(\tau) = \frac{1}{2\pi } { \sum _{n=-\infty }^{n=\infty }\frac{\exp (in\beta \tau)}{(\Omega ^2-n^2\beta ^2+i\gamma n \beta)}}. \end{equation} \tag{ A1 }$$

To perform the summation, we use the general result that if for a general periodic function

$$\begin{equation} g(\tau) = \sum _{n=-\infty }^{n=\infty }b(n)\exp (in\beta t), \end{equation} \tag{ A2 }$$

with period 2π/β, and b(n) being defined as an integrable function, we set

$$\begin{equation} F(\tau) = \frac{1}{2\pi }\int ^{\infty }_{-\infty }b(n)\exp (in\beta \tau)dn, \end{equation} \tag{ A3 }$$

then

$$\begin{equation} g(\tau) =2\pi \sum _{n=-\infty }^{n=\infty }F(\tau +2\pi n/\beta). \end{equation} \tag{ A4 }$$

We set

$$\begin{equation} b(n)= \frac{1}{2\pi (\Omega ^2-n^2\beta ^2+i\gamma n \beta)}. \end{equation} \tag{ A5 }$$

Then the integral in Equation (A3) defining F(τ) is readily performed by contour integration with the result that for t>0

$$\begin{equation} F(\tau) = \frac{\exp (-\gamma \tau /2)\sin (\omega _{\gamma }\tau)}{2\pi \omega _{\gamma }\beta }, \end{equation} \tag{ A6 }$$

otherwise F(τ) = 0. Here $\omega _{\gamma }=\sqrt{\Omega ^2-\gamma ^2/4}.$ Using the above to evaluate the sum in Equation (A4) as a geometric progression yields g(τ) ≡ G(τ) for 0 < τ < 2π/β as

$$\begin{eqnarray} &&\hspace{-1pt}G(\tau)\nonumber\\ && = \frac{\exp (-\gamma \tau /2)\sin (\omega _{\gamma }\tau)- \exp (-\gamma \pi /\beta)\sin (\omega _{\gamma }(\tau -2\pi /\beta))}{\omega _{\gamma }\beta \left[1+\exp (-2\gamma \pi /\beta)-2\exp (-\gamma \pi /\beta)\cos (2\pi \omega _{\gamma }/\beta)\right]},\nonumber\\ \end{eqnarray} \tag{ A7 }$$

the function is determined elsewhere by its periodicity with period 2π/β.