Turbulence in Extrasolar Planetary Systems Implies that Mean Motion Resonances are Rare
aa r X i v : . [ a s t r o - ph ] M a y TURBULENCE IN EXTRASOLAR PLANETARY SYSTEMSIMPLIES THAT MEAN MOTION RESONANCES ARE RARE
Fred C. Adams , , Gregory Laughlin , and Anthony M. Bloch , Michigan Center for Theoretical PhysicsPhysics Department, University of Michigan, Ann Arbor, MI 48109 Astronomy Department, University of Michigan, Ann Arbor, MI 48109 Lick Observatory, University of California, Santa Cruz, CA 95064 Department of Mathematics, University of Michigan, Ann Arbor, MI 48109
ABSTRACT
This paper considers the effects of turbulence on mean motion resonances in ex-trasolar planetary systems and predicts that systems rarely survive in a resonant con-figuration. A growing number of systems are reported to be in resonance, which isthought to arise from the planet migration process. If planets are brought together andmoved inward through torques produced by circumstellar disks, then disk turbulencecan act to prevent planets from staying in a resonant configuration. This paper studiesthis process through numerical simulations and via analytic model equations, whereboth approaches include stochastic forcing terms due to turbulence. We explore howthe amplitude and forcing time intervals of the turbulence affect the maintenance ofmean motion resonances. If turbulence is common in circumstellar disks during theepoch of planet migration, with the amplitudes indicated by current MHD simulations,then planetary systems that remain deep in mean motion resonance are predicted tobe rare. More specifically, the fraction of resonant systems that survive over a typicaldisk lifetime of ∼ P bound ≈ C/N orb1 / , where the dimensionless parameter C ∼ −
50 and N orb is thenumber of orbits for which turbulence is active. Subject headings:
MHD — planetary systems — planetary systems: formation — plan-ets and satellites: formation — turbulence
1. INTRODUCTION
An increasing number of the observed extrasolar planetary systems have been discovered tocontain multiple planets. A growing subset of these multiple planet systems have period ratios 2 –close to the ratios of small integers and hence are candidates for being in mean motion resonance(e.g., Mayor et al. 2001, Marcy et al. 2002, Butler et al. 2006). In true resonant configurations, theorbital frequencies not only display integer ratios, but, in addition, the relevant resonant arguments(angular variables composed of the osculating orbital elements) display constrained oscillatory timeevolution (for further discussion, see Murray & Dermott 1999; hereafter MD99; see also Peale 1976).As a result, these resonant configurations represent rather special dynamical states of the planetarysystem and their existence (if/when observed) places interesting constraints on the formation andlong term evolution of these systems. More specifically, as we show in this paper, planetary systemsare easily knocked out of mean motion resonance. One source of perturbations acting on the planetsis turbulent fluctuations in the circumstellar disks that originally formed the planets. Since thesedisks are also necessary ingredients in the current picture of planetary migration (e.g., Papaloizou& Terquem 2006), and since coupled migration seems necessary to put planets into mean motionresonance (e.g., Lee & Peale 2002; Lee 2004; Moorhead & Adams 2005; Beaug´e et al. 2006), theeffects of disk turbulence on the resonances is important to understand.This paper analyzes the effects of turbulence on planets in or near mean motion resonanceusing both a semi-analytic treatment ( §
2) and direct numerical integrations ( §
2. SEMI-ANALYTIC MODEL EQUATIONS2.1. Pendulum Model for the Resonance
For the sake of definiteness, we consider the case in which a larger planet on an external orbitperturbs a smaller planet on an interior orbit. The larger body is then assumed to have orbitalelements that do not change with time (or at least vary much less than those of the smaller innerplanet). We are thus making the most extreme approximation and keeping only the leading orderterms in order to obtain an analytic description. Following MD99, the equation of motion for thelibration angle ϕ for two planets in resonance reduces to that of a pendulum, i.e., d ϕdt + ω sin ϕ = 0 , (1)where the natural oscillation frequency ω is given by ω = − j C r Ω e | j | . (2)Here, Ω and e are the mean motion and eccentricity of the inner planet. The integers j and j depend on the type of resonance being considered. The parameter C r , which is taken to be constant,is given by C r = µ Ω αf d ( α ) , (3)where µ = m out /M ∗ and where m out is the mass of the outer (perturbing) planet. The quantity αf d ( α ) results from the expansion of the disturbing function (MD99), where the parameter α is theratio of the semimajor axes of the two planets, i.e., α ≡ a /a . This approximation neglects termsof order O ( µ ). Finally, we note that more complicated analytic models for mean motion resonancescan be derived (e.g., Holman & Murray 1996; Quillen 2006), but they are qualitatively similar tothe pendulum equation considered here.For much of this analysis, we are interested in the 2:1 mean motion resonance since observedextrasolar planetary systems are (sometimes) found in or near such a configuration, and since theseresonances are generally the strongest. For this case, the integers j = − j = −
1. For the 2:1resonance, αf d ( α ) ≈ − /
4, where α is the ratio of semimajor axes of the two planets. With thesespecifications, the natural oscillation frequency ω of the libration angle is given approximately by ω ≈ µ Ω e . (4)In typical cases, where planet masses are 1000 times smaller than the stellar masses, and theeccentricity e ∼ .
1, we find that ω / Ω ∼ − . As a result, the period of oscillation for thelibration angle is typically ∼
100 orbits.For completeness, we note that for first order resonances, the oscillation frequency ω oftenhas additional contributing terms (MD99). The order of the additional terms differs from that ofequation (2) by a factor of order O ( µ/e ), where e is the eccentricity of the inner planet. Under 4 –some circumstances of interest, µ ∼ − and e ∼ .
1, so this correction can be of order unity. Weare thus considering only the simplest form of the pendulum equation in this analysis, in order togain a general understanding of the effects of turbulence on resonant systems. In principle, however,care must be taken to include all of the relevant terms.
Given the pendulum model for the resonance (from MD99), the next step is to include theeffects of additional perturbations produced by turbulence in the circumstellar disk material thatremains after the planets have been formed. Since this disk is generally thought to induce planetarymigration, some disk material is likely to be present, for a typical time scale of 1 – 10 Myr. Undermany circumstances, these disks are expected to be turbulent (Balbus & Hawley 1991) due to themagneto-rotational instability (MRI). Although MRI can be shut down if the ionization rate isnot high enough (Gammie 1996), star forming regions can experience enhanced cosmic ray fluxes,and hence enhanced ionization rates, due to cosmic rays from supernovae becoming trapped inthe magnetic fields of molecular clouds (Fatuzzo et al. 2006). In this work, we consider turbulentfluctuations to be present and consider our previous model (Laughlin et al. 2004; hereafter LSA)to provide the parameters of the perturbations. Note that the LSA model describes turbulenceis disks with no gaps, applicable to low-mass embedded planets, whereas this paper primarilyconsiders giant planets that will clear gaps. As a result, we must modify the LSA formalism toinclude this effect (see below).In principle, the turbulent fluctuations in the disk can affect either planet, i.e., the outerperturbing planet (in this analysis) or the inner perturbed planet. If the outer planet is subjectto turbulent fluctuations, they have the effect of shaking the effective base of the pendulum. As aresult, they modify the effective potential energy of the pendulum, and the stochastic term entersinto the equation of motion as follows: d ϕdt + [ ω + ξ ( t )] sin ϕ = 0 , (5)where ξ ( t ) is the stochastic forcing term. This is the form expected when the circumstellar disksurrounds both planets and primarily influences the outer one. In contrast, if the disk (the torquesproduced by turbulence in the disk) primarily influences the inner planet, then the fluctuationsprovide an additional acceleration term and the equation of motion takes the form d ϕdt + ω sin ϕ = − e ξ ( t ) . (6)For the sake of simplicity, we will primarily consider the first case (eq. [5]). Notice also that bothtypes of fluctuations can be included.The inclusion of turbulent fluctuations thus provides a stochastic forming term in the pendulumequation of motion (1). The forcing term ξ ( t ) can be considered as a series of “kicks” produced 5 –by short-term forcing due to turbulence. To completely specify this forcing effect, one needs todetermine the spectrum of fluctuation amplitudes, the timing of the fluctuations, and any possiblecorrelations in the fluctuation timing. This latter property determines the so-called “color” of thenoise (fluctuations) and can have (perhaps) surprising effects on the long-term evolution of thestochastic pendulum and hence the resonance angles of the planetary system (Mallick & Marcq2004, hereafter MM04).The formulation for turbulence developed in LSA shows that the turbulent fluctuations arereloaded into the disk on roughly an orbital time scale. In order to consider the turbulent torquesto be fully independent, we must use the longest relevant orbital period, which corresponds to thatof the outer disk edge (in the calculations of LSA). This period at the outer disk edge is abouttwice that of the planet being forced. Since the outer (perturbing) planet is being forced by theturbulence in the present context, and since the outer planet is in a 2:1 mean motion resonancewith the inner planet, this turbulent forcing period is approximately four times the period of theinner planet.Next we need to account for the fluctuations in the equation of motion. To the same order ofapproximation used to derive the pendulum equation for the resonance, the time rate of change ofthe mean motion is given approximately by d Ω dt ≈ −
13 Ω (cid:18) J dJdt (cid:19) , (7)where J is the angular momentum of the inner planetary orbit. Here we consider the turbulentfluctuations to provide discrete changes in the angular momentum on time scales comparable tothe orbital period P D at the disk edge, where (as discussed above) we expect P D ≈ P orb = 8 π/ Ω.Thus, the time scale ∆ t for independent stochastic perturbations to act is given by ∆ t ≈ π/ Ω.Since d ϕ/dt ∼ d Ω /dt (see equation [8.39] from the derivation of MD99), the forcing term ξ k ∝ Ω( ˙
J /J ) ∼ Ω[(∆ J ) k /J ] δ ([ t ] − ∆ t ). As a result, we can write the stochastic differential equation inthe form d ϕdt + (cid:20) ω −
13 Ω (cid:18) ∆ JJ (cid:19) k δ ([ t ] − ∆ t ) (cid:21) sin ϕ = 0 . (8)In this equation, [(∆ J ) /J ] k represents the fractional change in the angular momentum of the planetfor a given realization of the turbulence. The subscript ‘ k ’ indicates that this increment of angularmomentum is delivered at discrete time intervals, which can be counted and labeled with an integer.To account for this timing, the Dirac delta function is periodic with period ∆ t , so that the quantity[ t ] is the time measured mod-∆ t .The amplitude [(∆ J ) /J ] k is calculated in LSA for disks without gaps (see also Nelson &Papaloizou 2003, hereafter NP03, and the discussion below). The resulting torques produce angularmomentum perturbations with zero mean and well-defined RMS amplitudes of order [(∆ J ) /J ] rms = j rms ∼ . R to account for the absence of the full complement of disk material near the planet. In practice,we find that Γ R ≈ .
1, where this value can be estimated in multiple ways. First, we note thatnumerical simulations of gap formation in turbulent disks (NP03) show that the gaps are notcompletely cleared out; instead, the surface density at the gap minimum is reduced by a factorof ∼
20 from the unperturbed level (see Fig. 1 of NP03). Thus we expect the reduction factorto lie in the range 0 . ≤ Γ R ≤
1. Next we note that material from both the outer part of thegap (where the reduction in surface density is less severe) and outside the gap (where there is noreduction) provide some contribution to the torques. Indeed, the dominant contribution of torquesis expected to come from intermediate length scales (see Johnson et al. 2006), specifically scaleslarger than the disk scale height, but not too much larger (given that gravitational forces decreasewith distance). To model this effect, we start with the formalism of LSA and remove disk materialcorresponding to a quadratic gap profile with zero surface density at the center (e.g., Goldreich &Sari 2003) and then calculate the reduction factor Γ R as a function of gap width. Let w g denote thehalf-gap width, i.e., the location on either side of the planet where the quadratic gap profile joinsonto the unperturbed surface density distribution. The planet is assumed have semi-major axis a p , which also defines the location of the gap center, where the surface density is taken to vanish.Then we find that Γ R ≈ w g /a p = 1, 1.5, and 2, respectively.Although the exact value of the reduction factor will depend on the details of the gap shape, weuse Γ R = 0.1 as a representative value for this work. We also note that these systems have twoplanets, so that each planet can clear disk material, but each planet can also experience torques;we assume here that these two competing effects cancel.For completeness, we note that the surface density structure in the local neighborhood of amassive planet can be strongly perturbed, with well-defined wakes and even a circumplanetary diskwithin the gap (e.g., see Fig. 10 of NP03). Although global transport seems to be unaffected bythese distortions, their effects on stochastic torques remain unclear and should be considered infuture work.Since the forcing amplitude [(∆ J ) /J ] k plays an important role, it is useful to have an heuristicunderstanding of the expected value. For a circumstellar disk with surface density Σ, the torqueexerted on a planet can be written in terms of the benchmark value T D = 2 πG Σ rm p . The torquedue to turbulence is expected to be a fraction f T of this physical scale, where numerical simulationssuggest that f T ≈ .
05 (e.g., Johnson et al. 2006; Nelson 2005). This torque acts over a timeinterval to produce a change in angular momentum. As argued above, this calculation uses atime interval of four orbital periods, long enough so that the turbulent torques can be taken to beindependent. The torque thus acts on a timescale 4 P orb = 8 π ( r /GM ∗ ) / and the net change inangular momentum (∆ J ) = 4 P orb f T T D . Next, we include the reduction factor Γ R due to the diskgap (where we expect Γ R ∼ . J of the planet is givenby J = m p [ GM ∗ r (1 − e )] / . We can absorb the eccentricity dependence into the definition of f T .Putting these results together, we can write the fractional change in angular momentum over one 7 –time interval (4 P orb ) in the form (cid:18) ∆ JJ (cid:19) k = f T Γ R π Σ r M ∗ ≈ − (cid:18) Σ1000g cm − (cid:19) , (9)where the second equality is scaled for a 1 AU orbit and a typical surface density at that radius.Note that we expect the surface density profile to have a nearly power-law form Σ ∼ r − p , where p ∼ /
2, so the quantity Σ r increases slowly with radius. Since giant planets form at larger radiallocation and migrate inwards, the torques will be somewhat larger during the earlier phases ofmigration. With the mean motion resonance modeled by a pendulum equation, and with the character-istics of the turbulent fluctuations specified, we now find the time evolution of the system. If wedefine a dimensionless time variable τ ≡ ω t , (10)where ω is the oscillation frequency of the libration angle, and is defined by equation (4), theequation of motion takes the form d ϕdτ + [1 + η k δ ([ τ ] − ∆ τ )] sin ϕ = 0 , (11)where δ ( x ) is the Dirac delta function. In this formulation, ∆ τ is the time interval betweenstochastic kicks, the dimensionless time variable is measured mod-∆ τ , and the amplitude of thekicks is given by η k = −
13 Ω ω (cid:18) ∆ JJ (cid:19) k = − √ µe (cid:18) ∆ JJ (cid:19) k . (12)The quantity η k is thus a stochastic variable that has a distribution of values inherited from thedistribution of values of [(∆ J ) /J ] k , which in turn is determined by the properties of the turbulence.The kicks in angular momentum can be either positive or negative, with no expected asymmetry,so that h η k i = 0. We can thus characterize the amplitude by the RMS (root-mean-square) of thevariable. Using the amplitudes expected from numerical MHD simulations ( § η rms = h η k i / ≈ . − . , (13)where we have used typical values of the mass ratio µ = 10 − and the eccentricity e = 0.1 to eval-uate the result. The typical expected value of the RMS fluctuation amplitude is thus η rms ∼ . t ≈ π/ Ω, which sets the correspondingdimensionless time interval to be ∆ τ ≈ πω / Ω. In general, both the amplitudes and the timeintervals between forcing kicks will vary from cycle to cycle according to their (respective) distri-butions. However, we can scale out the time variations (Appendix A; Adams & Bloch 2008) and 8 –characterize the turbulence with a well-defined (single-valued) forcing time. As a result, we presentthe remainder of this analysis using a single forcing time interval and focus on the effects of thedistribution of forcing amplitudes η k .To analyze the behavior of the stochastic pendulum, as defined above, it is useful to work interms of phase space variables. We thus rewrite the equation of motion into two parts, i.e., dϕdτ = V , and dVdτ = − [1 + η ] sin ϕ . (14)We then define the probability distribution function for the phase space variables, P ( ϕ, V ; t ), whichobeys the Fokker-Planck equation (MM04; Binney & Tremaine 1987) ∂P∂τ + V ∂P∂ϕ − sin ϕ ∂P∂V = D ϕ ∂ P∂V , (15)where the phase space diffusion constant D is set by the amplitude of the fluctuation spectrum.Specifically, we can write the diffusion constant in terms of previously defined variables, D = h η k i ∆ τ ≈ h η k i Ω8 πω . (16)Next we argue that the libration angle ϕ varies rapidly compared to the velocity V . This claim canbe verified by noting that dV /dτ ∼ η k so that V ∼ τ / (since the random variable η k implies arandom walk growth). However, the libration angle obeys the equation dϕ/dτ = V , which in turnimplies that ϕ ∼ τ / (see MM04 for further discussion). In the long term, the libration angle thusvaries faster than V and we can average the Fokker-Planck equation over the angle ϕ to obtain thetime evolution equation for the averaged probability distribution function p ( V ; τ ), i.e., ∂p∂τ = h sin ϕ i D ∂ p∂V . (17)This equation has the well-known solution p ( V ; τ ) = 1( πDτ ) / exp (cid:20) − V Dτ (cid:21) , (18)where we have defined an effective diffusion constant D = 2 Dh sin ϕ i . In the long time limit, ϕ ∼ τ / , and the angle fully samples the range [0 , π ], so that h sin ϕ i = 1 / D = D .Notice also that the velocity V grows as V ∼ τ / in the long time limit. Since the energy of thependulum grows large at long times, the kinetic term dominates the potential energy term, andthe energy E ≈ V /
2. With this substitution, the probability distribution function for the energytakes the simple form P ( E ; τ ) = (cid:18) πDτ (cid:19) / E − / exp (cid:20) − EDτ (cid:21) . (19)This result implies that the expectation value of the energy grows linearly with time in the longtime limit, i.e., h E i = D τ . (20) 9 – Fig. 1.— Mean time evolution of the energy of a stochastic pendulum with noise terms of various“colors”. For each curve/line, 10 realizations of the time evolution of the oscillator have beenaveraged together. The curves show the results for white noise (solid), pink noise (dotted), andbrown noise (dashed). The first two types of noise lead to energy evolution that is linear in time, h E i ∼ t , whereas brown noise leads to energy that increases h E i ∼ √ t . 10 –This linear time dependence can be understood in terms of a random walk of velocity, as shown inAppendix B.This time behavior is seen in numerical simulations of the stochastic pendulum equation.Figure 1 shows the time development of the mean energy, averaged over 10 realizations, for threedifferent types of fluctuations. The solid curve shows white noise, where the fluctuation amplitudeand the time interval ∆ t between kicks are uniformly distributed about well-defined mean values.The dotted curve shows the case of “pink” noise, in which the distribution of time intervals is takento have a decaying exponential form; this distribution leads to more time intervals near the lowerend of the spectrum. As long as the mean time interval and mean fluctuation amplitude remainthe same, and as long as both are uncorrelated with previous samples, the linear time dependencepredicted by equation (20) holds. For completeness, Figure 1 also shows the case of “brown” noise,in which the time intervals undergo a random walk (akin to brownian motion) so that a correlationof time intervals is present; in this case, the time development of the energy is slower with h E i ∝ t / (see also MM04). Thus, the long term evolution of an ensemble of stochastic pendulums obeys theanalytic expectations derived above, provided that (1) the fluctuations are independent, (2) thelong-time limit has been reached, and (3) the system can freely random walk back into a resonant(bound) state (from a higher energy, unbound state). Given the time evolution of the distribution function, we now consider the fraction of lowenergy states as a function of time. Here, states of the system with sufficiently low energy arebound and the pendulum oscillates — instead of circulates — so that the system is in resonance.We are thus finding an estimate of the fraction of systems that remain in resonance as a function oftime, where systems leave resonance due to exposure to turbulent forcing. Note that our treatmentthus far provides only an upper limit to the fraction of bound states because we allow the systemto freely random walk both in and out of resonance. The fraction of bound states derived in thissection thus represents the largest possible value. In the following section, we find a more accurateestimate by accounting for the reduced probability of systems returning to resonance after leaving.As formulated here, the energy E is the kinetic energy, so that the criterion for the resonanceto be undone (for the libration angle to circulate) is given (approximately) by E >
1. Note that thesystem must have
E > < E < h E i grows ever larger. This probability is adecreasing function of time and is given by P bound = (cid:18) πDτ (cid:19) / Z dE √ E exp (cid:20) − EDτ (cid:21) = 2 √ π Z z e − z dz = Erf( z ) , (21) 11 – Fig. 2.— Time evolution of the probability for the planetary system to stay in a mean motionresonance as a function of time. The solid curve shows the result of 10 realizations of the stochasticpendulum equation; in this set of numerical experiments, the system is allowed to come in and outof resonance freely. The dashed curve shows the expected asymptotic behavior of the system, i.e., asurvival probability of the form P bound ∝ t − / . The dot-dashed curve shows the survival probabilityfor when a one-way barrier is imposed so that systems that leave resonance (the pendulum becomesunbound) are not allowed to re-enter the bound configuration. The dotted curve shows the analyticapproximation to this probability evolution derived in the text. 12 –where z = (2 /Dτ ) / . Here Erf( z ) is the error function (e.g., Abramowitz & Stegun 1970), whichcan be expanded in the limit of small z as follows:Erf( z ) = 2 √ π (cid:18) z − z z − z
42 + . . . (cid:19) , (22)where small z values correspond to late times. To leading order, the probability that the planetarysystem remains in resonance is thus given by the relation P bound ≈ (cid:18) πDτ (cid:19) / = (cid:18) π h E i (cid:19) / = 8 η rms [Ω t ] / = 4 p /πη rms √ N orb , (23)where this expression is valid only for sufficiently late times when τ > /D . For planetary systems,the natural clock is set by the orbit time, so in the final equality we measure time using the numberof orbits N orb ≡ Ω t/ (2 π ). If we then insert the expected fluctuation amplitude η rms ≈ .
005 (eq.[13]), the above expression simplifies to the form P bound ≈ /N orb1 / , where N orb includes onlythe orbits (time) for which turbulence is active. For typical systems, where the number of orbits N orb = 10 − , this expression implies that the probability of remaining in resonance lies in therange P bound ≈ . − .
60. Keep in mind, however, that this derivation assumes that planetarysystems can freely random walk both in and out of resonance, so that equation (23) represents anupper limit to the fraction of surviving resonant systems. But even this upper limit demonstratesthat turbulence significantly reduces the expected number of resonant systems. In the more realisticcase, where systems can random walk out resonance more easily than they can random walk backinto resonance ( § § systems. The smooth dashed and dotted curves show the analytic predictions for the timedependence derived above.This estimate of the probability for remaining bound is approximate in several respects. First,we have not included the potential energy term ( ω cos ϕ ) in the calculation of the probabilitydistribution P ( E ; τ ). This approximation leads to an uncertainty of ∼ √ ϕ ), so this behavior is natural. In actualplanetary systems, however, other physical variables are relevant and must have the proper valuesto allow the planets to be in resonance. This difficulty implies that once a planetary systems leavesresonance, it will be unlikely to random walk back into resonance. This behavior, in turn, implies 13 –that the fraction of planetary systems remaining in resonance will fall below the fiducial law derivedhere. We consider this complication in the following sections. The discussion presented thus far assumes that the stochastic oscillator can diffuse in and outof a bound state, i.e., into resonance as well as out of resonance. For the simple model case, whichhas only one variable, this behavior is sensible. For real planetary systems, however, the conditionsof resonance are more complicated than the simple pendulum model. In particular, for highlyinteractive systems (e.g., the observed 2:1 resonance in GJ876; see § P bound ∝ τ − / , represents an upper limit on the fraction of systems that remain in resonance.In this section, we explore a more complicated model for the probability evolution that includes theone-way nature of this process. We also consider the intermediate case of a partially transmittingbarrier, where systems can return to resonance after leaving, but with reduced probability. We now consider a simple generalization of the time evolution of the probability. At any giventime, the fraction of systems that would remain bound in the absence of the one-way barrier isgiven by P bound from equation (23). We then assume that some fraction of these systems are“lost” at a well-defined rate, so that the time evolution of the fraction of bound states (in theabsence of the diffusion term) would be an exponential decay. This ansatz thus assumes that theremaining systems are distributed across the possible bound energy values, and that each systemhas some probability of leaving its bound state per unit time. With this assumption, the Fokker-Planck equation for the averaged probability distribution function p ( V ; τ ) now takes the modified(heuristic) form ∂p∂τ = D ∂ p∂V − γ (cid:18) πDτ (cid:19) / p , (24)where we assume that the bound fraction has the asymptotic form derived above and where theparameter γ encapsulates the details of this “decay process”. Note that a similar form arises(Caughey 1963) when dissipation is included in the Fokker-Planck equation (these connections 14 –should be explored further). Equation (24) can be solved to obtain the result: p ( V ; τ ) = 1( πDτ ) / exp (cid:20) − V Dτ (cid:21) exp " − γ (cid:18) τπD (cid:19) / . (25)With this complication, the fraction of systems that remain bound as a function of time now takesthe form P bound ≈ (cid:18) πDτ (cid:19) / exp " − γ (cid:18) τπD (cid:19) / . (26)Although heuristic, this model captures the basic time behavior for stochastic pendulums with aone-way barrier, as shown in Figure 2. For the simplest case of a constant diffusion constant, the solution for a one-way barrier can befound using separation of variables. In this case, the solution to the diffusion equation is assumedto have the general form p ( t, V ) = T ( t ) g ( V ); the diffusion equation then can be written as1 T dTdt D = 1 g d gdV = − Λ , (27)where Λ is the separation constant.The solutions g ( V ) have the form g n ( V ) = A n cos Λ n V , (28)where the subscript denotes one of a series of solutions that satisfy the boundary condition. Wewe invoke the boundary condition that g ( V = √
2) = 0, which implies that the distribution func-tion vanishes when the kinetic energy reaches E = V / n , i.e., Λ n = π √ n + 1 / , (29)where n is an integer. The full solution thus takes the form p ( t, V ) = ∞ X n =0 A n cos(Λ n V ) exp (cid:20) − D Λ n t (cid:21) . (30)For the initial condition p ( t = 0 , V ) = δ ( V ), where δ ( x ) is the Dirac delta function, the constantsare specified so that A n = 1 / √ . (31)The fraction of bound states at any time is then given by P bound = 2 Z √ p ( t, V ) dV = ∞ X n =0 − n π ( n + 1 /
2) exp (cid:20) − D π ( n + 1 / t (cid:21) . (32) 15 –This solution, though approximate, indicates that the long term evolution of the ensemble ofresonant states follows an exponential decay. For the simple solution of equation (30), each “mode”decays exponentially, albeit with a differing decay rate. In the long run, however, only the lowestorder mode survives, and the time evolution becomes purely exponential. The evolution of the actual stochastic pendulum differs from that described by a simple dif-fusion equation in several respects. Most importantly, the diffusion constant D is not really aconstant. In the formulation developed here, D ∝ sin ϕ , and we have taken the average value h sin ϕ i = 1/2. However, this limit only applies in the long time limit when most of the oscillatorsin the ensemble have large energy, so that their angles uniformly sample the possible values. Forbound states, by definition, the energies are low, and hence this limiting form is not strictly applica-ble. In particular, the typical amplitudes of the bound oscillators will be smaller than average, andthe true effective diffusion constant will be less than that of the long time limit. In addition to itslower value, the diffusion constant will also have some type of time dependence: In the beginning,when all of the oscillators in the ensemble have small amplitudes, D will be correspondingly small.As time increases, even the bound oscillators will have larger amplitudes and hence the evolutionis properly described by larger values of D (but still smaller than that of the long time limit). As shown above, the time evolution of the number of bound states decays exponentially forsystems with a strict one-way barrier, i.e., for the case in which systems that leave resonance arenot allowed to return. We now consider the case in which a given system has some probability ρ ret of re-entering a bound (resonant) state after leaving, where 0 ≤ ρ ret ≤
1. In this case, all ofthe systems are allowed to stay in play, but only a fraction of the systems can make the transitionfrom an unbound state to a bound state. As shown below, even when the fraction ρ ret ≪
1, thismodification makes a large qualitative change in the long-term behavior of the system.We first note that at late times, essentially all of the bound states in the original problem(without a barrier) are those that have returned to resonance after leaving. This claim followsdirectly from the above results: The fraction of bound states of all types decreases as ∼ t − / ,while the fraction of bound states that have never left resonance decreases as ∼ exp[ − γt ], wherethe decay rate γ depends on the details of the diffusion process. At late times, however, the power-law behavior wins, and hence most bound states are due to systems that have returned from higherenergy states. If the higher energy states are allowed to re-enter resonance with probability ρ ret ,then the fraction of systems in a resonant state will again follow the t − / law, but with a smallernormalization. Specifically, the normalization is reduced by the factor ρ ret and the long-term 16 –solution for the fraction of bound states becomes P bound ≈ ρ ret (cid:18) πDτ (cid:19) / . (33)This behavior is illustrated in Figure 3, which shows the time evolution for an ensemble ofstochastic pendulums. The upper two curves show the evolution for the case in which the systemscan freely re-enter a resonant (bound) state after leaving. The bottom two curves show the fractionof bound states for the case in which the systems are allowed to have only a 20% probability ofreturning to resonance after leaving. The solid and dashed curves show the results of the actualnumerical integrations, whereas the dotted curves show the analytic approximations derived above.Notice the good agreement between the predicted asymptotic forms and the numerical results. Thefluctuations of the numerically determined fractions about the analytic predictions have amplitudesthat are consistent with root- N noise; note that the numerical ensemble contains only 10 systemsand the number of bound states is much smaller by the end of the time interval shown in Figure 3.On a related note, one can consider the problem in which systems leaving a bound state havesome probability p X of never returning to resonance. Perhaps surprisingly, the long term behaviorfor this case is wildly different than for the previous case. Here, for any nonzero value of p X , thetime dependence of the number of bound states is always an exponential decay. As outlined above, the 2:1 resonance is generally the strongest, and planets in such a config-uration should survive longest in the presence of turbulence. For other resonances, the oscillationfrequencies are generally smaller so that the periods of libration are longer. For the case of the 3:1resonance, for example, ones finds (MD99) that ω ≈ − . µ Ω e . (34)The negative sign indicates that the stable equilibrium point of the oscillator occurs at ϕ = π rather than ϕ = 0. After defining φ ≡ ϕ − π , the equation of motion for φ becomes the same asbefore. The above frequency is thus lower than that of the (simplified) 2:1 resonance by a factor of ∼ (3 /e ) / . The potential energy of the oscillator is thus less deep by a factor (3 /e ) / , about 5–6for typical cases, so that the system is more easily bounced out of a resonant state by turbulence.In order to assess the probability of remaining in a bound (resonant) configuration, one canevaluate the equations derived in § /e ) / ∼
5. In other words, atany given time, the probability of a planetary system remaining in a mean motion resonance, inthe face of turbulence, is inversely proportional to the relevant value of ω . For typical values of 17 – Fig. 3.— Time evolution of the fraction of planetary systems that stay in a mean motion resonanceas a function of time. The solid and dashed curves show the result of 10 realizations of thestochastic pendulum equation. For the top (solid) curve, the system is allowed to come in and outof resonance freely; for the bottom (dashed) curve, the system has only a 20% chance of returningto a bound state after leaving. The two dotted curves shows the expected asymptotic behavior ofthe system, i.e., a survival probability of the form P bound ∝ t − / . 18 –the orbital elements of extrasolar planets, where e ∼ .
1, the probability of staying in a 3:1 meanmotion resonance is about 5 times smaller than for a 2:1 resonance.
The basic analytical model developed here contains three physical variables that describe theeffects of turbulent perturbations on mean motion resonance. The first variable is the natural oscil-lation frequency ω of the resonance ( § ∼
100 times longer than the orbital time scale of the planets, anddepends on the system properties in a well-known manner (see § t ≈ π/ Ω (or its dimensionless counterpart ∆ τ ≈ πω / Ω) over which theturbulent perturbations are re-established to provide an independent realization of the torques.The final variable is the amplitude of the perturbations. Since the torques and their correspondingchanges in angular momentum can be either positive or negative, and hence the mean vanishes, theperturbation amplitude can be characterized by the root mean square j rms = h (∆ J ) /J i / ; thisquantity also has a dimensionless counterpart η rms given by equations (12) and (13). Numericalsimulations of MHD turbulence indicate that the dimensionless amplitude is expected to be of order η rms ∼ .
005 (including the reduction factor due to gap clearing).For the simplest version of this model, where systems can freely random walk back into aresonant state, the three variables ( ω , ∆ t, η rms ) determine the long term evolution of the ensemble.In particular, the expectation value of pendulum energy increases linearly with time accordingto equation (20). The probability of the system staying (bound) in resonance decreases as thesquare root of time according to equation (23). Note that both of these results depend only on thedimensionless time τ = ω t and a dimensionless diffusion constant D = η rms2 / ∆ τ . In other words,the action of the turbulent fluctuations can be described by a diffusion constant that determineshow the variables of the system – considered here as a nonlinear oscillator – random walk throughphase space.In more realistic versions of the model, where systems cannot easily return to a resonantstate, the long term behavior results in a lower probability of remaining bound. As a result, thesimple result described above provides an upper limit on the expected number of resonant statesas a function of time. In the opposite limit in which systems that leave resonance can neverreturn, the number of bound states decays exponentially, with varying decay rates, as indicatedby equations (26) and (30). For the case in which systems can re-enter resonance after leaving,but with reduced probability, the number of bound states decreases as t − / as before, but withreduced normalization, as shown by equation (33). 19 –
3. NUMERICAL SIMULATIONS
In this section, we consider an ensemble of numerical integrations motivated by an observedextrasolar planetary system. For this exploration, we adopt the outer planets c ( P ∼
60 d) and b( P ∼
30 d) of the Gliese 876 planetary system (Marcy et al. 2001). These planets are observed to liedeep in the 2:1 mean motion resonance, with the angles θ and θ librating with small amplitudes(Laughlin & Chambers 2001, Rivera et al. 2005). Indeed, Gliese 876 b and c represent, by far, themost convincing case of an extrasolar planetary system in mean motion resonance. Among the 25multiple-planet systems that have been discovered to date, 2:1 mean motion resonances have alsobeen tentatively identified for HD 82943 b and c (Gozdziewski & Maciejewski 2001, Mayor et al.2004, Lee et al. 2006), HD 128311 b and c (Vogt et al. 2005), and HD 73526 b and c (Tinney etal. 2005). In each of these three cases, however, the libration widths of the resonant arguments arevery uncertain, and for HD 128311 and HD 73526, the best dynamical fit to the radial velocity datashows only a single argument in libration. Following the discovery of 55 Cancri c and d (Marcyet al. 2002), it was suggested by several authors (e.g., Ji et al. 2003) that 55 Cancri b and c arepossibly participating in a 3:1 mean motion resonance. Self-consistent 6-body fits to the 55 Cancridata set published by Fischer et al. (2007), however, show no evidence that 55 Cancri b and c arein 3:1 resonance. In this set of numerical experiments, we start an ensemble of planetary systems with the obser-vationally determined orbital elements of GJ876, so that the system begins in a 2:1 mean motionresonance, and then integrate the system forward in time. The integrations include additionalstochastic fluctuations, which would be present if the circumstellar disk that formed the planetswere still present. We then monitor how easily the system is knocked out of its resonant configu-ration. The ensemble of numerical integrations reported here are thus analogous to those done inthe previous section with the stochastic pendulum. In this case, however, we integrate the full 18phase space variables (6 orbital elements for each planet and 6 more for the star) instead of onlyone variable ( ϕ ) for the pendulum model. For this system, integrating the motion of the star isimportant because its mass is relatively small, only 0.32 M ⊙ , whereas the masses of the planetsare 0.79 and 2.53 m J , and hence the system is highly interactive (Laughlin & Chambers 2001).As a result, the GJ876 system is as far from the idealized (one variable) model of the stochasticpendulum as any planetary system observed to date; as a result, any agreement found between thesimple model and the full integrations can be considered robust.The starting orbital elements are taken to be those observed (Rivera et al. 2005) at the presentday; specifically, the epoch for the initial conditions is JD 2449679.6316, the date for which the firstpublished Lick Observatory radial velocity was made. The systems are then integrated into thefuture for 2 × years. Since the planets in this particular system have relatively short periods(60.830 days and 30.455 days for the outer and inner planets, respectively), the number of orbits is For an up-to-date χ -ranked list of self-consistent fits to the 55 Cnc data sets see http://oklo.org?p=257 and thelinks therein.
20 –large – about 2 . × orbits of the inner planet. The integrations are performed with a Bulirsch-Stoer (B-S) integration scheme (e.g., Press et al. 1992) to obtain high accuracy with reasonablecomputational speed. In the absence of turbulent forcing, the system can be integrated forward intime for millions of years and is found to stay in resonance, where we use the following definitionsfor the libration angles: θ = λ − λ + 2( ̟ − ̟ ) and θ = λ − λ + ( ̟ − ̟ ) , (35)where the λ j are the mean longitudes and the ̟ j are the longitudes of periapse of the two planets.Note these angles are defined to be a linear combination of the standard ones (e.g., Lee & Peale2002), and are chosen because they provide a cleaner separation of the apsidal libration and themean motions (J. Chambers, private communication). We measure the degree of resonance bymonitoring the maximum amplitude of these libration angles over a given interval of time. Forthis purpose, we use a fixed monitoring time of 100 years, which is much longer than the period ofthe oscillations of the resonance, about 3.5 years, so that many oscillation periods are included inthe determination of the maximum. This maximum libration angle is then tracked as a functionof time. For the sake of definiteness, we consider the systems to be bound (in resonance) whenthe maximum angle remains less than θ ∗ = 90 degrees; when the maximum angle exceeds θ ∗ =90 degrees, we consider the resonance to be compromised. Note that the exact number of bound(resonant) states as a function of time depends on the choice of the angular threshold value θ ∗ ;however, the general trends (shown below) are qualitatively similar for any choice of threshold anglein the range 90 ≤ θ ∗ ≤ §
2. In terms of the simplependulum model of the system, variations of the stochastic forcing times could be scaled out of theproblem (Appendix A) and would effectively add width to the distribution of forcing strengths. Inthese experiments, we take the size of the velocity perturbations to have the form ∆ v = v A ξ , wherethe amplitude v A = 0.0032 km/s, and where ξ is a uniformly distributed random variable on theinterval − ≤ ξ ≤
1. (This amplitude value was chosen to be a round number, namely 2 × − ,in code units, where G = 1, mass is measured in M ⊙ , and time is measured in years.) Both the ˆ x and ˆ y components of the velocity are perturbed (independently) in this manner, but the verticalvelocity component is left unchanged. Including both components of the velocity perturbation, wefind that the fractional change in speed, and hence angular momentum, has amplitude j rms ∼ − ,consistent with the benchmark value of equation (9).The results from one ensemble of simulations are shown in Figures 4 and 5. Figure 4 showsthe time evolution of the maximum libration angle (as defined above) for five different trials, i.e.,five different realizations of the turbulent fluctuations. Note that the systems can re-enter boundstates – defined here to be maximum libration angles less then θ ∗ = 90 degrees – after leaving. 21 – Fig. 4.— Time evolution of the resonance angles for a collection of planetary systems based on theobserved system GJ876. The systems are subjected to turbulent forcing as described in the text.The five curves show the first resonance angle θ as a function of time (given here in years) for thefive different realizations of the turbulent fluctuations. 22 – Fig. 5.— Time evolution of the fraction of bound states, as measured by the maximum excursionsof the resonance angles, for a collection of planetary systems based on the observed system GJ876.The systems are subjected to turbulent forcing as described in the text. The solid and dot-dashedcurves show the fraction of systems that remain in resonance as a function of time (in years) forthe first and second resonance angles, respectively. The dotted curve shows the fraction of boundstates for the first resonance angle that would result if systems that leave resonance could neverreturn to a bound state. For reference, the dashed curve shows the power-law behavior P b ∝ t − / expected for a simple stochastic pendulum. Note that the survival probability P b ≪ .
01 if weextrapolate this result out to time scales of ∼ N fluctuations are ∼ ∼ P bound ∝ t − / for reference. Note that the time evolution of the fraction ofbound states falls well below the prediction of a pure power-law with P bound ∼ t − / . The longterm behavior of the number of bound states seems to follow a steeper power-law, but does reach afull exponential decay during the time window studied (which would be expected if systems neverreturn to resonance).As an order of magnitude check on the time scales, note that the fractional change in velocityper stochastic kick has amplitude ∆ v ∼ − v orb , where v orb is the orbital speed of the inner planet.As outlined in §
2, however, the energy associated with the potential well of the resonance is muchsmaller than that of orbit. If we denote the change in velocity required to leave resonance as v res ,then v res ∼ . v orb , which in turn implies ∆ v ∼ . v res . If ∆ v accumulates as a random walk,the system thus requires ∼ stochastic forcing kicks to account for enough energy to compromiseresonance, and this number translates into a time scale of ∼
4. DISCUSSION AND CONCLUSION
The main finding of this paper is that mean motion resonance in planetary systems is readilycompromised through the action of turbulent fluctuations in circumstellar disks. If MRI operatesand the accompanying turbulent torques are common during the epoch of planetary migration,then planetary systems in mean motion resonances are predicted to be rare. Specifically, even forthe most favorable case in which the system can freely random walk in and out of resonance, only asmall percentage of the solar systems that produce pairs of planets in a resonant configuration willremain in resonance over typical disk lifetimes of ∼ t − / (see Fig. 2). This law holds both forthe case in which systems can freely re-enter resonance, and for the case in which they re-enter res-onance with a given (fixed) probability, although the normalization is lower for the latter case (eq.[33]). If systems that leave resonance can never return to a bound state, then the time evolutionof the fraction of bound states is a decaying exponential ( § t − / power-law behavior of the simplest pendulum system andthe decaying exponential form that applies for systems that can never re-enter resonance afterbecoming unbound. The numerical simulations show that systems can indeed random walk backinto a resonant state after leaving, consistent with the finding that the number of bound statesdoes not decay exponentially. On the other hand, the return to resonance is relatively rare, inparticular more unlikely than for the (one variable) stochastic pendulum, a result that is consistentwith the fraction of bound states falling below the t − / law. This latter power-law behavior —the solution to the the long term behavior of a stochastic pendulum – thus provides an upper limiton the probability of resonant states as a function of time.The most important astronomical implication of this work is a quantitative estimate of therarity of solar systems surviving in a resonant configuration. Our results can be summarized bythe following expression, which represents our working estimate for the fraction of surviving bound(resonant) states: P bound ≈ ρ ret η rms (cid:18) πN orb (cid:19) / ≡ CN orb1 / , (36)where ρ ret is the probability of returning to resonance, η rms is set by the amplitude of the fluctua-tions, and N orb is the number of orbits over which the turbulence is active. In the second equality wehave defined the dimensionless quantity C = (32 /π ) / ρ ret /η rms . For typical turbulent amplitudes, η rms ≈ . C ≈ ρ ret . Given that returning to resonance requires somewhat specialconditions, we expect ρ ret to be small, i.e., ρ ret ≪
1. This expectation is borne out in our numericalsimulations ( § P bound decreases rapidly with time. As a result, we expect that 25 –the constant C will be of order ∼
10 (perhaps as large as ∼
50) so that P bound ≈ /N orb1 / .If turbulence is present for a substantial fraction of the disk lifetime, then N orb = 10 − and P bound = 0.001 – 0.01, i.e., mean motion resonances will be rare.In addition to the prediction that resonant states are rare, this work has other astronomicalimplications: For the particular case of the Gliese 876 system, if no eccentricity damping is assumedfor the epoch of planet formation and migration, then the current eccentricities and libration widthsof the system suggest that the system migrated only a small distance (8% of the initial semi-majoraxes) while in resonance (Lee & Peale 2002). This result is consistent with any turbulence in thedisk of the GJ876 system having had little time or ability to provide substantial perturbations onthe resonant arguments. On the other hand, the observed HD 128311, HD 82943 and HD73526planetary systems are all consistent with having received significant exposure to turbulent per-turbations, but not enough to have destroyed their librations completely. Once a large sample ofmultiple-planet systems has been assembled, we should gain further insight from the distribution ofobserved resonant widths, in addition to the observed fraction of planets in mean motion resonance.In particular, if planets can random walk back into the resonance after leaving, as indicated by bothour numerical integrations and model equations, the observational sample should show a definitepreference for systems with large libration widths. This trend is exactly what we observe for thetenuous resonances in the three systems other than Gliese 876.This analysis also outlines the basic physical mechanism by which turbulence affects meanmotion resonance. We have modeled the process through a pendulum equation (which representsthe resonance) with the addition of stochastic forcing (which represents the turbulence). Thisformulation applies specifically to the class of systems for which the turbulent fluctuations provideforcing kicks that are fully independent, i.e., where both the amplitude of the fluctuations andthe time intervals are not correlated. In this case, the results are largely independent of thedistributions, and depend primarily on the expectation values of the amplitudes and time intervalsof the turbulent forcing. These latter quantities set the value of the effective diffusion constant(for an analogous problem, see Fatuzzo & Adams 2002, which considers ambipolar diffusion withturbulent fluctuations in molecular cloud cores).As an intuitive check on the results presented herein, consider the following heuristic argument.In a stochastic pendulum, the system accumulates angular momentum as a random walk, so thechange in angular momentum after N step steps is given by N step (∆ J ) k , where (∆ J ) k is the typicalamplitude of the angular momentum fluctuations. In order of magnitude, the amount of angularmomentum in the resonant system (when bound) is given by ( J orb ω / Ω) . The combination of thesetwo expressions indicates that the criterion for compromising resonance can be written in the form N step ≥ [(∆ J ) k /J orb ] − ( ω / Ω) . In order of magnitude, (∆ J ) k /J orb ∼ − and ω / Ω ∼ − , andhence the required number of steps is of order N step ∼ . Indeed, when the number of stochasticsteps (measured here by the number of orbits N orb ) exceeds this limit, the probability of remainingin a bound state decreases (see eqs. [23] and [36]). 26 –Notice also that the effective diffusion constant D ∝ (∆ J ) k / (∆ t ) k . In order to set the propervalue of the constant D , one must choose the time interval to be long enough so that the fluctu-ations are independent and short enough that the changes in angular momentum (∆ J ) increaselinearly with time (where we now suppress the index k ). Suppose, for example, one were tochoose a longer time interval to sample the fluctuations, say (∆ t ) ≫ (∆ t ) , where (∆ t ) is thetime interval for optimal independent sampling. The change in angular momentum over onesampling interval would then increase as a random walk (even within the time interval), i.e.,(∆ J ) = (∆ J ) [(∆ t ) / (∆ t ) ] / . The effective value of the diffusion constant then scales accord-ing to D ∝ (∆ J ) / (∆ t ) = (∆ J ) [(∆ t ) / (∆ t ) ] / (∆ t ) = (∆ J ) / (∆ t ) . We thus recover the correctvalue, provided that the changes in angular momentum are computed properly. In this paper, wehave taken the time interval ∆ t to be four times the orbital period of the inner planet (twice theperiod of the outer planet), consistent with the calculations of LSA (see also Nelson 2005). Theabove considerations imply that the diffusion constant D is inversely proportional to ∆ t , and thefraction of bound states scales as D − / , so that our quoted results are only weakly dependent onany uncertainty in the specification of ∆ t .Although our general conclusion is robust – namely, that turbulence easily compromises meanmotion resonances – a good deal of additional work should be done. The parameter space availablefor these planetary systems is enormous and further exploration should be carried out for differentsolar system architectures. In particular, the degree to which the full numerical integration agreewith the simple model equations should be considered, including a determination of the solar systemproperties required for consistency. This treatment assumes that the planets have already formed,and have already migrated to their (nearly) final radial locations, and then considers the effectsof turbulence on mean motion resonances. In actuality, planets can be subjected to turbulentfluctuations as they form, and as they migrate inward and become locked into resonance. Theseprocesses of planet formation, migration, resonance locking, and turbulent forcing should thusbe studied concurrently. In addition, the differences between the results of our full numericalintegrations (which include all 18 phase space variables of the system) and those of the idealizedstochastic pendulum (which has only one degree of freedom) poses a mathematically interestingissue for further study.This work was initiated during a sabbatical visit of FCA to U. C. Santa Cruz; we wouldlike to thank both the Physics Department at the University of Michigan and the AstronomyDepartment at the University of California, Santa Cruz, for helping to make this collaborationpossible. We thank REU student Jeff Druce for performing supporting calculations and thankan anonymous referee for many insightful comments that improved and clarified the paper. Thiswork was supported in part by the Michigan Center for Theoretical Physics. FCA is supported byNASA through the Origins of Solar Systems Program via grant NNX07AP17G. AMB is supportedby the NSF through grants CMS-0408542 and DMS-604307. GL is supported by the NSF throughCAREER Grant AST-0449986, and by the NASA Planetary Geology and Geophysics Programthrough grant NNG04GK19G. 27 – A. RESCALING FOR VARIABLE TIME INTERVALS
In this Appendix, we show that if the stochastic forcing impulses occur at irregular timeintervals, then the time variable can be rescaled to make the stochastic forcing kicks occur atregular (periodic) intervals. As a compensating result, however, the distribution of the forcingstrengths becomes wider and the natural oscillation frequency of the system takes on a distributionof values.We begin by writing the stochastic pendulum equation in the form d ϕdt + (cid:2) ω + q k Q ( λ k t ) (cid:3) sin ϕ = 0 , (A1)where Q is a periodic function and q k determines the forcing strength for the kth time interval; thevariable λ k has unit mean, but allows for different time intervals for the stochastic forcing terms.Notice that we generally use periodic Dirac-delta functions to specify the forcing time, but thisargument can be generalized to include any periodic function Q . If we re-scale the time variablefor each cycle according to t → λ k t , (A2)then the stochastic pendulum equation takes the form λ k d ϕdt + (cid:2) ω + q k Q ( t ) (cid:3) sin ϕ = 0 , (A3)where the time variable appearing in this equation is the rescaled one. If we then re-scale theremaining parameters ( ω , q k ) according to the relations ω → ω /λ k ≡ ω k and q k → q k /λ k ≡ ˜ q k , (A4)the stochastic pendulum equation has the same form as that for the case of perfectly periodic forcingterms. This result is thus the analog of Theorem 1 of Adams & Bloch (2008) for the stochasticHill’s equation. Notice that in this case the forcing strength q k acquires a new (generally wider)distribution to accommodate the variation in forcing periods, and, the natural oscillation frequency ω becomes a random variable. B. ALTERNATE DERIVATION OF THE LINEAR TIME DEPENDENCE FORTHE MEAN ENERGY
In this Appendix, we present an alternate derivation of the result that the energy of a stochas-tically driven pendulum tends to increase linearly with time (in the limit of long times). This resultis thus consistent with equation (20) in the text.We begin by considering the pendulum equation with a periodic delta-function forcing term,as described in the text ( § ϕ itself must be continuous at the moment of forcing, but 28 –the derivative of ϕ obeys a jump condition that can be written in the form dϕdτ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + = dϕdτ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − − q k sin ϕ k , (B1)where we have denoted the value of the angle at the kth cycle boundary as ϕ k . Since the angle ϕ must be a continuous function across the boundary, the potential energy is also a continuousfunction, and hence the total energy experiences discrete jumps of the form∆ E = − q k sin ϕ k dϕdτ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − + 12 q k sin ϕ k . (B2)Over long times, the first term in the above equation averages to zero, whereas the second termaverages to h q k i /
4. The mean energy thus grows as E ≈ h q k i τ τ , (B3)where ∆ t is the time interval of the stochastic forcing. Not only does this argument reproduce thelinear growth of the mean energy, as found using the Fokker-Planck analysis in the text, but it alsoresults in the same coefficient if we identify the effective diffusion constant D = h q k i / (∆ τ ). REFERENCES
Abramowitz, M., & Stegun, I. A. 1970, Handbook of Mathematical Functions (New York: Dover)Adams, F. C., & Bloch, A. M. 2008, SIAM J. Appl. Math., 68, 947Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214Beaug´e, C., Michtchenko, T. A., Ferraz-Mello, S. 2006, MNRAS, 365, 1160Binney, J., & Tremaine, S. 1987, Galactic Dynamics (Princeton: Princeton Univ. Press)Butler, P. R., et al. 2006, ApJ, 646, 502Caughey, T. K. 1963, J. Acoustic Soc. Amer., 35, 1683Fatuzzo, M., & Adams, F. C. 2002, ApJ, 570, 210Fatuzzo, M., Adams, F. C., & Melia, F. 2006, ApJ, 653, L49Fischer, D. A., et al. 2007, ArXiv e-prints, 712, arXiv:0712.3917Gammie, C. F. 1996, ApJ, 457, 355Goldreich, P., & Sari, R. 2003, ApJ, 585, 1024 29 –Go´zdziewski, K., & Maciejewski, A. J. 2001, ApJ, 563, L81Ji, J., Kinoshita, H., Liu, L., & Li, G. 2003, ApJ, 585, L139Holman, M. J., & Murray, N. W. 1996, AJ, 112, 1278Johnson, E. T., Goodman, J., & Menou, K. 2006, ApJ, 647, 1413Laughlin, G., & Chambers, J. E. 2001, ApJ, 551, 109Laughlin, G., Steinacker, A., & Adams, F. C. 2004, ApJ, 608, 489 (LSA)Lee, M. H. 2004, ApJ, 611, 517Lee, M. H., & Peale, S. J. 2002, ApJ, 567, 596Lee, M. H., Butler, R. P., Fischer, D. A., Marcy, G. W., & Vogt, S. S. 2006, ApJ, 641, 1178Mallick, K., & Marcq, P. 2004, J. Phys. A, 37, 4769 (MM04)Marcy, G. W., Butler, R. P., Fischer, D., Vogt, S. S., Lissauer, J. J., & Rivera, E. J. 2001, ApJ,556, 296Marcy, G. W., Butler, R. P., Fischer, D., Laughlin, G., Vogt, S. S., Henry, G. W., & Pourbaix, D.2002, ApJ, 581, 1375Mayor, M. et al. 2001, ESO Press Release 07/01Mayor, M. et al. 2004, A & A, 415, 391Moorhead, A. V., & Adams, F. C. 2005, Icarus, 178, 517Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge: Cambridge Univ.Press), (MD99)Nelson, R. P. 2005, A&A, 443, 1067Nelson, R. P., & Papaloizou, J.C.B. 2003, MNRAS, 339, 993 (NP03)Papaloizou, J.C.B., & Terquem, C. 2006, Rep. Prog. Phys., 69, 119Peale, S. J. 1976, ARA&A, 14, 215Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical Recipes inFORTRAN: The art of scientific computing (Cambridge: Cambridge Univ. Press)Quillen, A. C. 2006, MNRAS, 365, 1367Rivera, E. J., et al. 2005, ApJ, 634, 625 30 –Tinney, C. G., Butler, R. P., Marcy, G. W., Jones, H. R. A., Laughlin, G., Carter, B. D., Bailey,J. A., & O’Toole, S. 2006, ApJ, 647, 594Vogt, S. S., Butler, R. P., Marcy, G. W., Fischer, D. A., Henry, G. W., Laughlin, G., Wright, J. T.,& Johnson, J. A. 2005, ApJ, 632, 638