Evolution of Planetary Systems with Time Dependent Stellar Mass Loss
aa r X i v : . [ a s t r o - ph . S R ] M a r EVOLUTION OF PLANETARY SYSTEMSWITH TIME DEPENDENT STELLAR MASS LOSS
Fred C. Adams , , Kassandra R. Anderson , and Anthony M. Bloch Physics Department, University of Michigan, Ann Arbor, MI 48109 Astronomy Department, University of Michigan, Ann Arbor, MI 48109 Math Department, University of Michigan, Ann Arbor, MI 48109
ABSTRACT
Observations indicate that intermediate mass stars, binary stars, and stellar rem-nants often host planets; a complete explanation of these systems requires an under-standing of how planetary orbits evolve as their central stars lose mass. Motivated bythese dynamical systems, this paper generalizes in two directions previous studies oforbital evolution in planetary systems with stellar mass loss: [1] Many previous treat-ments focus on constant mass loss rates and much of this work is carried out numerically.Here we study a class of single planet systems where the stellar mass loss rate is timedependent. The mass loss rate can be increasing or decreasing, but the stellar massalways decreases monotonically. For this class of models, we develop analytic approx-imations to specify the final orbital elements for planets that remain bound after theepoch of mass loss, and find the conditions required for the planets to become unbound.We also show that for some mass loss functions, planets become unbound only in theasymptotic limit where the stellar mass vanishes. [2] We consider the chaotic evolutionfor two planet systems with stellar mass loss. Here we focus on a model consistingof analogs of Jupiter, Saturn, and the Sun. By monitoring the divergence of initiallysimilar trajectories through time, we calculate the Lyapunov exponents of the system.This analog solar system is chaotic in the absence of mass loss with Lyapunov time τ ly ≈ −
10 Myr; we find that the Lyapunov time decreases with increasing stellar massloss rate, with a nearly linear relationship between the two time scales. Taken together,the results of this paper help provide an explanation for a wide range of dynamicalevolution that occurs in solar systems with stellar mass loss.
Subject headings: planets and satellites: dynamical evolution and stability — planet-star interactions — stars: evolution — stars: mass loss — white dwarfs
1. Introduction
Solar systems orbiting other stars display a diverse set of architectures and motivate furtherstudies concerning the dynamics of planetary systems. Part of the richness of this dynamical 2 –problem arises from the intrinsic complexity of N-body systems, even in the absence of additionalforces (Murray & Dermott 1999). The ledger of physical behavior experienced by such systems isenormous, and includes mean motion resonances, secular interactions, and sensitive dependence onthe initial conditions (chaos). Additional complications arise from additional forces that are oftenpresent: During early stages of evolution, circumstellar disks provide torques that influence orbitalelements, and turbulent fluctuations act on young planets. Over longer time scales, solar systemsare affected by tidal forces from both stars and planets, and by general relativistic corrections thatlead to orbital precession. Another classic problem in solar system dynamics concerns planetaryorbits around central stars that are losing mass (Gyld´en 1884, Jeans 1924; see also Hadjidemetriou1963, 1966). Although this issue has received some recent attention (see below), this paper expandsupon existing work in two main directions. Recent work often focuses on the particular case ofconstant mass loss rates, although stellar mass loss rates typically vary with time; in addition,this recent work is primarily carried out numerically (note that Veras et al. 2013 use numericalsimulations to consider more realistic, time-dependent stellar mass loss). In this paper, for singleplanet systems, we extend existing calculations to account for time dependence of the mass lossrates and obtain a number of analytic results. For systems with two or more planets, we also showthat the Lyapunov exponents, which determine the time scales for chaotic behavior, depend on thetime scales for mass loss. As outlined below, these two results can account for a great deal of thepossible behavior in solar systems where the central star loses mass.A number of previous studies have considered planetary dynamics for host stars that are losingmass. For our own Solar System, long term integrations have been carried out to study the fateof the planets in light of mass loss from the dying Sun (Duncan & Lissauer 1998). Recent relatedwork estimates an effective outer boundary r B to the Solar System (due to stellar mass loss) inthe range r B = 10 − AU, where orbiting bodies inside this scale remain safely bound (Veras& Wyatt 2012). Planets orbiting more massive stars, which lose a larger percentage of their mass,have their survival threatened by possible engulfment during the planetary nebula phase (Villaver& Livio 2007, Mustill & Villaver 2012), and are more likely to become unbound due to stellar massloss alone (Villaver & Livio 2009). In the future of our own system, Earth is likely to be engulfedby the Sun (Schr¨oder & Connon Smith 2008), but planets in wider orbits are expected to survive.However, gaseous planets that escape engulfment are still subject to evaporation and can experiencesignificant mass loss (Bear & Soker 2011, Spiegel & Madhusudhan 2012). For planets orbiting starsthat are losing mass, a more general treatment of the dynamics has been carried out for both singleplanet systems (Veras et al. 2011) and multiple planet systems (Veras & Tout 2012; Voyatzis etal. 2013); these studies provide a comprehensive analysis for the particular case of constant massloss rates. In addition to causing planets to become unbound, stellar mass loss can drive orbitalevolution that leads to unstable planetary systems surrounding the remnant white dwarfs remainingat the end of stellar evolution (Debes & Sigurdsson 2002). Indeed, observations indicate that whitedwarfs can anchor both circumstellar disks (Melis et al. 2010) and planetary systems (Zuckermanet al. 2010); many white dwarf atmospheres contain an excess of heavy elements (Melis et al.2010; Jura 2003), which is assumed to be a signature of accretion of a secondary body (a planet or 3 –asteroid). Finally, mass loss in binary star systems can lead to orbital instability, allowing planetsto change their host star (Kratter & Perets 2012; for additional related work, see also Perets &Kratter 2012, Moeckel & Veras 2012).This paper builds upon the results outlined above. Most previous studies have focused on stellarmass loss rates that are constant in time, and most recent work has been carried out numerically.However, stars generally have multiple epochs of mass loss, the corresponding rates are not constant,and these solar systems span an enormous range of parameter space. It is thus useful to obtaingeneral analytic results that apply to a wide class of mass loss functions. The first goal of thispaper is to study single planet systems where the stellar mass loss rate varies with time. As thesystem loses mass, the semimajor axis of the orbit grows, and the planet becomes unbound forcritical values of the stellar mass fraction m f = M f /M ∗ . For systems that become unbound, wefind the critical mass fraction m f as a function of the mass loss rate and the form of the mass lossfunction. In other systems with mass loss, the orbit grows but does not become unbound. In thesecases, we find the orbital elements at the end of the mass loss epoch, again as a function of themass loss rate and the form of the mass loss function. For initially circular orbits and slow massloss (time scales much longer than the initial orbital period), the critical mass fraction and/or thefinal orbital elements are simple functions of parameters that describe the mass loss rate. For initialorbits with nonzero eccentricity, however, the outcomes depend on the initial orbital phase. In thislatter case, the allowed values of the critical mass fraction m f (or the final orbital elements) takea range of values, which we estimate herein.Next we consider the effects of additional planets on the results described above. If the planetsare widely spaced, they evolve much like individual single planet systems. However, if the planetsare sufficiently close together so that planet-planet interactions are important, the systems aregenerally chaotic. As a secondary goal, this paper estimates the Lyapunov exponents for two-planet systems with stellar mass loss. For the sake of definiteness, we focus on planetary systemscontaining analogs of the Sun, Jupiter, and Saturn, i.e., bodies with the same masses and (usually)the same orbital elements. We find that the time scale for chaos (the inverse of the Lyapunovexponent) is proportional to the mass loss time scale. As a result, by the time the star haslost enough mass for the planets to become unbound, the planets have begun to erase their initialconditions through chaos. For systems that evolve far enough in time, one can use the semi-analyticresults derived for single planet systems with initially circular orbits (see above) as a rough estimateof the conditions (e.g., the final value ξ f of the radius) required for a planet to become unbound.Multiple planets and nonzero initial eccentricity act to create a distribution of possible values (e.g.,for ξ f ) centered on these results. Since the two-planet systems are chaotic, and display sensitivedependence on initial conditions, one cannot unambiguously predict the value of ξ f required for anplanet to become unbound.For completeness, we note that the problem of planetary orbits with stellar mass loss is anal-ogous to the problem of planetary orbits with time variations in the gravitational constant G . Forsingle planet systems (the pure two-body problem), the gravitational force depends only on the 4 –(single) product GM ∗ , so that the two problems are equivalent (e.g., Vinti 1974). However, for thecase with time varying gravitational constant, the product GM ∗ could increase with time. Currentexperimental limits show that possible variations occur on time scales much longer than the cur-rent age of the universe (see the review of Uzan 2003), so that these effects only (possibly) becomeimportant in the future universe (Adams & Laughlin 1997). We also note that when consideringtime variations of the constants, one should use only dimensionless quantities, in this case thegravitational fine structure constant α G = Gm P /c ~ (e.g., Duff et al. 2002).This paper is organized as follows. We first present a general formulation of the problem ofplanetary orbits with stellar mass loss in Section 2 and then specialize to a class of models wherethe mass loss has a specific form (that given by equations [13] and [14]). These models include awide range of behavior for the time dependence of the mass loss rates, including constant mass lossrates, exponential mass loss, and mass loss rates that decrease quickly with time; these results aredescribed in Section 3. In the following Section 4, we consider two planet systems and calculatethe Lyapunov time scales for a representative sample of mass loss functions. Next we apply theseresults to representative astronomical systems in Section 5. The paper concludes, in Section 6, witha summary of our results and a discussion of their implications.
2. Model Equations for Orbits with Stellar Mass Loss
In this section we develop model equations for solar systems where the central star loses mass.We assume that mass loss takes place isotropically, so that the rotational symmetry of the systemis preserved and the total angular momentum is a constant of motion. This constraint is explicitlysatisfied in the analytical solutions that follow. For the numerical solutions, this property is usedas a consistency check on the numerical scheme. On the other hand, the total energy of the systemis not conserved because the total mass decreases with time (equivalently, the system no longerexhibits time reversal symmetry).General forms for the equations of motion with variable stellar mass are presented in manyprevious papers (from Jeans 1924 to Veras et al. 2011). Some of the subtleties of the variousapproaches are outlined in Hadjidemetriou (1963). In this section and the next we specialize tosystems with a single planet and focus on the case where the planet mass is much smaller than thestellar mass, M p ≪ M ∗ . The specific angular momentum J can be written in the form J = GM ∗ a η , (1)where a is the starting semimajor axis and M ∗ is the initial mass of the star. Equation (1) canbe taken as the definition of the angular momentum parameter η . For a starting circular orbit η = 1, whereas eccentric orbits have η = 1 − e <
1, where e is the initial orbital eccentricity. Theradial equation of motion can be written in the dimensionless form d ξdt = ηξ − m ( t ) ξ , (2) 5 –where η is constant and where we have defined a dimensionless (total) mass m ( t ) ≡ M ( t ) M ∗ . (3)The dimensionless radial coordinate ξ = r/a and the dimensionless time variable is given in unitsof Ω − = ( a /GM ∗ ) / .The goal of this work is to find general solutions to the problem where the dimensionless mass m ( t ) monotonically decreases with time. In the reduction of the standard two-body problem to ananalogous one-body problem, the equation of motion describes the orbit of the reduced mass. Inthe version of the problem with mass loss represented by equation (2), the motion is also that of areduced mass (e.g., Jeans 1924, MacMillan 1925, Hadjidemetriou 1963). In this treatment, the massloss functions (defined below) refer to the dimensionless scaled mass m ( t ). The model equation (2)is exact in the limit where the planet mass is small compared to the stellar mass, i.e., M P /M ∗ → O ( M P /M ∗ ). In applications of interest, however, our choice ofmass loss functions (and their uncertainties) provides the greatest degree of approximation – muchlarger than than the O ( M P /M ∗ ) difference between the reduced problem and the full problem.This paper thus focuses on the dimensionless problem of equation (2). The starting semimajoraxis is unity (by definition) and the initial conditions require that the starting radial coordinate ξ lies in the range 1 − e ≤ ξ ≤ e , where the starting eccentricity e is given by η = 1 − e . Notethat choosing the value of ξ is equivalent to specifying the starting phase of the orbit (up to asign). The initial energy E = − / ξ is given by˙ ξ = 2 ξ − η − ξ ξ = [(1 + e ) − ξ ][ ξ − (1 − e )] ξ . (4)The initial velocity can be positive or negative, where the choice of sign completes the specificationof the starting phase of the orbit. The equation of motion (2) is complicated because it contains an arbitrary function, namely m ( t ), that describes the mass loss history. On the other hand, the independent variable (time)does not appear explicitly. As a result, we may define a new effective “time” variable u throughthe expression u ≡ m , (5)where m = m ( t ). The generalized time variable u starts at u = 1 and is monotonically increasing.In terms of the variable u , the basic equation of motion (2) takes the equivalent form˙ u d ξdu + ¨ u dξdu = ηξ − uξ . (6) 6 –Next we note that both standard lore and numerical solutions (beginning with Jeans 1924)show that, in physical units, the product aM ≈ constant . In terms of the current dimensionlessvariables, this finding implies that the function f ≡ ξu = ξm (7)should vary over a limited range. We thus change the dependent variable from ξ to f and writethe equation of motion in the form˙ u u f (cid:20)(cid:0) u f ′′ + 2 uf ′ (cid:1) + u ¨ u ˙ u (cid:0) uf ′ + f (cid:1)(cid:21) = η − f , (8)where primes denote derivatives with respect to the variable u . Keep in mind that equation (8)is equivalent to the original equation of motion (2), with a change in both the independent anddependent variables.The leading coefficient in equation (8) represents an important quantity in the problem: Notethat the time scale for mass loss is given by u/ ˙ u and the orbital time scale is given by u f / (thislatter time scale is the inverse of the orbital frequency, and is shorter than the orbital period by afactor of 2 π ). The ratio λ of these two fundamental time scales is given by λ ≡ ˙ u u f . (9)The leading coefficient in equation (8) is thus λ , the square of the ratio of the orbital time scaleto the time scale for mass loss. For small values of λ , the mass loss time is long compared tothe orbit time, and the orbits are expected to be nearly Keplerian; for larger λ , the star loses asignificant amount of mass per orbit and a Keplerian description is no longer valid. For the formercase, where mass loss is slow compared to the orbit time, we can use the parameter λ to order theterms in our analytic estimates.In addition to the coefficient λ , given by the ratio of time scales, another important featureof equation (8) is the index β appearing within the square brackets, where β ≡ u ¨ u ˙ u . (10)The index β encapsulates the time dependence of the mass loss. This paper will focus on modelequations with constant β (such models have a long history, from Jeans 1924 to Section 2.2). The Orbital Energy:
For the chosen set of dimensionless variables, the energy E of the systemtakes the form E = 12 ˙ u (cid:0) uf ′ + f (cid:1) + η u f − u f . (11)The energy has a starting value E = − /
2, by definition, and increases as mass loss proceeds. Ifand when the energy becomes positive, the planet is unbound. Although the energy expression(11) appears somewhat complicated, the time dependence of the energy reduces to the simple form d E du = 1 u f . (12) 7 –Note that the derivative of the energy is positive definite, so that the energy always increases.Since the energy is negative and strictly increasing, the semimajor axis of the orbit, when definedaccording to a ∝ |E| − , is also monotonically increasing. Next we want to specialize to the class of mass loss functions where β = constant . The definingequation (10) for the mass loss index can be integrated to obtain the form˙ u = γu β , (13)where γ is a constant that defines the mass loss rate at the beginning of the epoch (when t = 0, m = 1, and u = 1). For a given (constant) value of the index β , the dimensionless mass loss ratehas the form ˙ m = − γm (2 − β ) . (14)In addition to simplifying the equation of motion, this form for the mass loss function is motivated bystellar behavior, as discussed below. The dimensionless parameter γ is defined to be the ratio of theinitial orbital time scale to the initial mass loss time scale. Specifically, if we define τ = ( M ∗ / ˙ M ∗ ) ,then γ is given by γ = 1 τ (cid:18) a GM ∗ (cid:19) / (15) ≈ . × − (cid:18) τ (cid:19) − (cid:16) a (cid:17) / (cid:18) M ∗ M ⊙ (cid:19) − / , where M ∗ is the stellar mass and a is the semimajor axis at t = 0. For typical orbits, a = 1 –100 AU, so that we expect the parameter γ to be small, often in the range 10 − < ∼ γ < ∼ − .The mass loss rate of stars is often characterized by the physically motivated form˙ M = − ˙ M C (cid:18) L ∗ L ⊙ (cid:19) (cid:18) R ∗ R ⊙ (cid:19) (cid:18) M ∗ M ⊙ (cid:19) − , (16)where ˙ M C is constant and depends on the phase of stellar evolution under consideration (Kudritzki& Reimers 1978, Hurley et al. 2000). Since the radius and luminosity depend on stellar mass(for a given metallicity), the physically motivated expression of equation (16) can take the samepower-law form as equation (14), which corresponds to a constant mass loss index (see equations[10] and [13]).Using the scaling law (16), the power-law index appearing in equation (14) can be positive ornegative, depending on how the stellar luminosity and radius vary with mass during the differentphases of mass loss (see Hurley et al. 2000 for a detailed discussion). For example, if we considermain-sequence stars, the stellar cores adjust quickly enough that the luminsoity obeys the standard 8 –mass-luminosity relationship L ∗ ∼ M p ∗ (where the index p ≈
3) and the mass-radius relationship R ∗ ∼ M q ∗ (where the index q typically falls in the range 1 / ≤ q ≤ m ∼ − m α m , where the index α m = p + q − . ≤ α m ≤
3; the corresponding mass loss index lies in the range − ≤ β ≤ − /
2. Nextwe consider stars on the first giant branch or the asymptotic giant branch. In this phase of stellarevolution, mass loss occurs from an extended stellar envelope, but the luminosity is produced deepwithin the stellar core. As the star loses mass, the core and hence the luminosity remains relativelyconstant, whereas the radius scales approximately as R ∗ ∼ M − / ∗ (Hurley et al. 2000). For thiscase, one obtains the scaling law ˙ m ∼ − m − / , with a mass loss index β = 10/3. In general, for˙ m ∼ − m α m , the mass loss index β = 2 − α m . As these examples show, the mass loss index cantake on a wide range of values − ≤ β ≤ β = 2, and themass evolution function has the form m ( t ) = 1 − γt and u ( t ) = 11 − γt . (17)The value β = 2 marks the boundary between models where the mass loss rate accelerates withtime ( β >
2) and those that decelerate ( β < β = 1, and the mass loss function has the form m ( t ) = exp[ − γt ] and u ( t ) = exp[ γt ] . (18)The value β = 1 marks the boundary between models where the system reaches zero stellar massin a finite time ( β > m → u → ∞ . For the casewith index β = 0, which represents an important test case, the mass evolution function becomes m ( t ) = 11 + γt and u ( t ) = 1 + γt . (19)For β = 0, analytic solutions are available (see Section 3.2), which inform approximate treatmentsfor more general values of the index β . Finally, the case where β = − m ( t ) = (1 + 2 γt ) − / and u ( t ) = (1 + 2 γt ) / . (20)The value β = − β > −
1) and those for which the planet becomes unbound only in the limit m → u → ∞ ( β < − β = 1, the time dependence of the mass takes the form m ( t ) = 1 u ( t ) = [1 − ( β − γt ] / ( β − . (21)The particular case β = 1 results in the decaying exponential law of equation (18). 9 – β For constant values of the mass loss index β , the equation of motion reduces to the form λ (cid:2) u f ′′ + (2 + β ) uf ′ + βf (cid:3) = η − f , (22)where the ratio of time scales λ is given by λ = γ u β +2 f . (23)By writing the equation of motion in the form (22), we immediately see several key features of thesolutions:When the parameter λ ≪
1, the left-hand side of equation (22) is negligible, and the equationof motion reduces to the approximate form f ≈ η = constant . This equality is only approximate,because the function f also executes small oscillations about its mean value as the orbit tracesthrough its nearly elliptical path (see below). Nonetheless, this behavior is often seen in numericalstudies of planetary systems with stellar mass loss (e.g., see Debes & Sigurdsson 2002). Orbitalevolution with λ ≪ λ ≫
1, the left-hand side of equation (22) dominates, and the solutionsfor f ( u ) take the form of power-laws with negative indices. In this regime, the equation of motionapproaches the form u f ′′ + (2 + β ) uf ′ + βf = 0 , (24)so that the function f ( u ) has power-law solutions with indices p given by the quadratic equation( p + 1)( p + β ) = 0 . (25)The general form for the solution f ( u ) in this regime is thus f ( u ) = Au + Bu β . (26)After the solutions enter this power-law regime, the energy can quickly grow and the planet canbecome unbound. To illustrate this behavior, consider the differential equation (12) for the orbitalenergy. We first consider the regime where λ ≪ f is nearly constant. For thebenchmark case f = 1, the equation can be integrated to obtain E = − u . (27)As long as f ≈
1, the energy remains negative and the planet remains bound, except in the limit u → ∞ . Now let λ > f = A/u , for u > u c , the differential equation for the energy can be integrated toobtain E = E c + 1 Au c (cid:16) − u c u (cid:17) , (28)where the subscript c denotes the reference point where the solutions enters into the power-lawregime. Since Au c ∼ u c and |E c | ∼ u c /
2, the energy quickly becomes positive once the power-lawregime is reached. This argument indicates that the planet becomes unbound when the time scaleratio λ is of order unity (as seen in previous studies). The subsequent subsections provide furtherverification of this finding.In order for the solution to make the transition from f ≈ constant to the power-law solutionsthat cause the orbits to become unbound, the ratio of time scales λ must grow with time. However,growth requires that β > − P varies with u (and hence time) according to P ∼ u f / . Since f is nearlyconstant, this relation simplifies to the form P ∼ u . The time scale for mass loss τ is given by τ = u/ ˙ u , which has the form τ ∼ u − β from equation (13). As a result, when β = −
1, the orbittime has the same dependence on stellar mass as the mass loss time scale, so that the ratio λ isnearly constant as the star loses mass. For β < −
1, the ratio λ of time scales decreases with time,and the system grows “more stable”.
3. Results for Single Planet Systems with Stellar Mass Loss
This section presents the main results of this paper for single planet systems with a centralstar that loses mass. First, we consider mass loss index β = −
1, which marks the critical value suchthat systems with β > − β < − m →
0. Next we consider mass loss index β = 0; inthis case, the solutions can be found analytically, and these results guide an approximate analytictreatment of the general case, which is addressed next. We also consider the limiting case wherestellar mass loss takes place rapidly. Here we consider systems where the mass loss index β = −
1, which corresponds to the tran-sition value between cases where the ratio λ of time scales grows with time ( β > −
1) and thosewhere the ratio decreases with time ( β < − γ (cid:2) u f ′′ + uf ′ − f (cid:3) = ηf − f . (29)The equation of motion can be simplified further by making the change of (independent) variable w ≡ log u , (30) 11 –so that the equation of motion becomes γ (cid:20) d fdw − f (cid:21) = ηf − f . (31)This version of the equation of motion (first considered by Jeans 1924) contains no explicit de-pendence on the independent variable w , so that the equation can be integrated to obtain theexpression γ (cid:18) dfdw (cid:19) = γ f + 2 f − ηf − E , (32)where E is a constant that plays the role of energy. Note that we have chosen the sign such that E > E = 1 for initially circular orbits. In order for the function f ( w ) to have oscillatorysolutions, the fourth order polynomial p ( f ) = γ f − Ef + 2 f − η (33)must be positive between two positive values of f . In order for this requirement to be met, theparameters must satisfy the inequalities ηE ≤
98 and γ ≤ ( E/ / . (34)The first inequality is always satisfied for the cases of interest. The second inequality in equation(34) determines the maximum value of the mass loss parameter γ for which oscillatory solutionsoccur. In the particular case where β = 0, the equation of motion can be simplified. In particular,the first integral can be taken analytically to obtain the form γ (cid:0) u f ′ (cid:1) = − ηf + 2 f − E , (35)where E is a constant of integration. The parameter E plays the role of energy for the orbit problemwhere the function f ( u ) plays the role of the radial coordinate. Although E is constant, the energy E of the physical orbit (where ξ is the radial coordinate) increases with time. Notice also that wehave adopted a sign convention so that E >
0. The value of E depends on the initial configuration.For the particular case where the orbit starts at periastron, for example, the initial speed ˙ ξ = 0and the energy constant has the value E = 1 − γ (1 − e ) , (36)where the eccentricity e = (1 − η ) / . In general, the initial value f = ξ can lie anywhere in therange 1 − e ≤ f ≤ e , and the energy constant has the general form E = 1 − γ f ± γ (cid:2) f − η − f (cid:3) / , (37) 12 –where the choice of sign is determined by whether the planet is initially moving outward (+) orinward ( − ). With the energy constant E specified, the turning points for the function f are foundto be f , = 1 ± [1 − ηE ] / E . (38)If we consider the function f to play the role of the radial coordinate, then equation (38) definesanalogs of the semimajor axis a ∗ and eccentricity e ∗ , which are given by a ∗ = 1 E and e ∗ = [1 − ηE ] / . (39)For a given starting value f , the effective eccentricity is given by the expression e ∗ = e + γ (1 − e ) f ∓ γ (1 − e ) (cid:0) f − η − f (cid:1) / . (40)Note that the effective eccentricity e ∗ of the function f is larger than the initial eccentricity e ofthe original orbit (before the epoch of mass loss). In particular, for a starting circular orbit e = 0,the effective eccentricity e ∗ = γ = 0. The integrated equation of motion (35) can be separated andwritten in the form f df [( f − f )( f − f )] / = E / γ duu . (41)If we integrate this equation from one turning point to the other, the change in mass of the systemthe same for every cycle, i.e., ∆ m = γπE / , (42)where E is given by equation (37). Following standard procedures, we can find the solution for theorbit shape, which can be written in the form ηf = η uξ = 1 + e ∗ cos θ . (43)The orbit equation thus takes the usual form, except that the original eccentricity e is replacedwith the effective eccentricity e ∗ and the effective “radius” variable ( f ) scales with the mass/timevariable u = 1 /m .For this mass loss function (with β = 0), we can find a simple relationship between the value ofthe time scale ratio λ and the value of f when the planet becomes unbound. To obtain this result,we insert the first integral from equation (35) into the general expression (11) for the energy E ofthe orbit and set E = 0. After eliminating the derivative f ′ , we can solve for the time scale ratio λ f as a function of the final value of f . When the planet becomes unbound, the time scale ratio isthus given by λ f = (2 f − η ) / ± (cid:0) f − η − f (cid:1) / f / . (44)Here, f is evaluated when the planet becomes unbound. In this case, however, the value of f isconstrained to lie in the range f ≤ f ≤ f , where the turning points are given by equation (38). 13 –For small γ , the orbit oscillates back and forth between the turning points many times before theplanet becomes unbound. The final value of f is thus an extremely sensitive function of the startingorbital phase. This extreme sensitivity is not due to chaos, and can be calculated if one knows theexact orbital phase at the start. In practice, however, the final value of f can be anywhere in therange f ≤ f ≤ f .Figure 3.2 shows the final values of the time scale ratio λ as a function of the final value of f = ξ/u = ξm . Curves are shown for nine values of the starting eccentricity, where e = 0.1, 0.2, . . . f corresponds to two possible values of the time scale ratio λ , one fororbits that are increasing in f and one for orbits that are decreasing in f at the time when theplanet becomes unbound. These two values of λ (for a given f ) correspond to the two branches ofthe solution given by equation (44).For comparison, Figure 3.2 shows the same plane of parameters for the final values of thetime scale ratio λ f and f f for planetary systems with constant mass loss rate (where β = 2).These results were obtained through numerical integration of equation (22). Note that the rangeof allowed values for the time scale ratio λ f is much larger than for the case with β = 0, whereasthe range of final values f f is somewhat smaller.One can also show that for circular orbits ( η = 1 and e = 0), the value of the time scale ratio λ = 1 when the planet becomes unbound. For circular orbits in the limit γ →
0, the turning pointsof the orbit appraoch f , = 1. Using f = 1 and η = 1 in equation (44), we find λ = 1. If we now consider the case where the mass loss is rapid, so that the equation of motion hassolution (26) throughout the evolution, we can fix the constants A and B by applying the initialconditions. Since f = ξ/u and u = 1 at the start of the epoch, f (1) = ξ , where ξ is the startingvalue of the orbital radius. By definition, the semimajor axis is unity, and the starting orbitaleccentricity is given by e = 1 − η . The starting radius thus lies in the range1 − p − η ≤ ξ ≤ p − η , (45)which is equivalent to 1 − e ≤ ξ ≤ e . The derivative f ′ = df /du is given by f ′ = − ξu + ξ ′ u = − ξu + 1 γu β +1 dξdt . (46)In the regime of interest where γ ≫
1, the second term is small compared to the first. In this limit, f ′ (1) = − ξ , where ξ lies in the range indicated by equation (45). The constants A and B arethus determined to be B = 0 and A = ξ , so that the solution has the simple form f ( u ) = ξ u . (47) 14 – Fig. 1.— Time scale ratio λ f as a function of the final value of f = f f for planetary systems with β = 0. For a given angular momentum η , specified by the starting eccentricity, the allowed values of λ f form closed curves in the plane as shown. Curves are shown for a range of starting eccentricity,from e = 0.1 (inner curve) to e = 0.9 (outer curve). 15 – Fig. 2.— Time scale ratio λ f as a function of the final value of f = f f for planetary systemswith β = 2 (constant mass loss rates). For a given angular momentum η , specified by the startingeccentricity, the allowed values of λ f form closed curves in the plane as shown. Curves are shownfor a range of starting eccentricity, from e = 0.1 (inner curve) to e = 0.9 (outer curve). Comparewith Figure 3.2 (and note the change of scale). 16 –The energy of the orbit, by definition, starts at E ( u = 1) = E = − /
2, and the energy obeys thedifferential equation (12). Combining the solution of equation (47) with the differential equation(12) for energy, we can integrate to find the energy as a function of u (equivalent to time or mass), E = −
12 + 1 ξ (cid:18) − u (cid:19) . (48)We can then read off the value of u f , and hence the mass m f , where the energy becomes positiveand the planet becomes unbound, i.e., m f = 1 u f = 1 − ξ . (49)Note that this critical value of the mass depends on the orbital phase of the planet within its orbit,i.e., the result depends on ξ rather than the starting semimajor axis, which is unity (notice alsothat this condition is equivalent to that given by equations [46–48] in Veras et al. 2011). For circularorbits, we must have ξ = 1, so that planets become unbound when the stellar mass decreases byone half (as expected).For cases where the mass loss is rapid, but the planet remains bound, we can find the orbitalproperties for the post-mass-loss system. Consider the limiting case where the star has initial mass m = 1, and loses a fraction of its mass instantly so that it has a final mass m ∞ , i.e., m ( t ) = m ∞ + (1 − m ∞ ) H ( − t ) , (50)where H is the Heaviside step function. The mass loss thus occurs instantaneously at t = 0. For t <
0, the solutions to the orbit equation (2) have the usual form,˙ ξ = 2 ξ − ηξ − ξ − ξ )( ξ − ξ ) ξ , (51)where ξ , = 1 ± e and η = 1 − e . After mass loss has taken place, the new (dimensionless) stellarmass is m ∞ , and the first integral of the equation of motion can be written in the form˙ ξ = 2 m ∞ ξ − ηξ − − m ∞ ) ξ , (52)where the final constant term takes into account the change in (dimensionless) energy at the momentof mass loss. The radial position at t = 0 is ξ ; since the planet is initially in a bound ellipticalorbit, the radial coordinate must lie in the range 1 − e ≤ ξ ≤ e . The energy E f of the neworbit is thus given by 2 E f = − − m ∞ ) ξ . (53)The energy E f is negative, and the orbit is bound, provided that the remaining mass m ∞ > − ξ / E f , the turning points ofthe new orbit take the form ξ ± = m ∞ ± (cid:2) m ∞ − |E f | η (cid:3) / |E f | . (54)We can then read off the orbital elements for the new (post-mass-loss) orbit, i.e., a f = m ∞ |E f | and e f = q − |E f | η/m ∞ , (55)where the new orbital energy E f is given by equation (53). In order to address the general case, we first change variables according to the ansatz x = u α where α = β + 1 . (56)After substitution, the equation of motion becomes γ α x (cid:2) x f xx + 2 xf x (cid:3) + γ βx f = ηf − f , (57)where the subscripts denote derivatives with respect to the new variable x . The ratio of time scalesis now given by λ = γ x f . (58)We can integrate the differential equation (57) to obtain the implicit form γ α (cid:2) x f x (cid:3) + 2 γ β Z x x f f x dx = − ηf + 2 f − E , (59)where
E > J ≡ γ β Z x x f f x dx , (60)so that γ α (cid:2) x f x (cid:3) + J = 2 f − η − Ef f . (61)Note that J = O ( λ ), which means that J will be negligible for most of the mass loss epoch (seeAppendix A). The energy of the system can be written in the form u E = 12 γ x ( αxf x + f ) + η f − f . (62)At the start of the evolution (where u = 1 and x = 1), the energy E = − / E , which takes the form E = 1 − γ f ± γ (cid:0) f − η − f (cid:1) / , (63) 18 –where f (= ξ ) is the starting value of the function (radial variable).At an arbitrary time during the epoch of mass loss, we can write the derivative f x = df /dx inthe form γα (cid:2) x f x (cid:3) = ± f (cid:0) f − η − Ef − J f (cid:1) / . (64)Since J = O ( λ ), the condition | J | ≪ E holds for most times. As a result, working to leadingorder, we can set J = 0 in equation (64) and recover analogs to the orbital solutions found earlierin Section 3.2 (for a related result, see Radzievskii & Gel’Fgat 1957; for an alternate approach, seeRahoma et al. 2009). The only difference is that the dependent (time-like) variable u is replacedwith x = u β +1 . As a result, the turning points for the function f ( x ) will be given by equation (38)and the orbital elements for f ( x ) are given by equation (39).The basic behavior of the orbit is illustrated by Figure 3.4. The function f , plotted here versus u as the solid black curve, oscillates between the turning points (marked by the red horizontallines) given by equation (38). The radial coordinate (here log ξ is plotted as the dotted blue curve)oscillates also, but grows steadily. The eccentricity of the orbit (green dashed curve) also oscillates,but grows with time. Finally, the time scale ratio λ (magenta dot-dashed curve) also oscillates andgrows with time. The simple oscillatory behavior for f ( u ) ceases near the point where the timescale ratio λ becomes of order unity. Note that the function f falls outside the boundary markedby the turning points near u = 12 in the Figure.Next we consider the time evolution of the energy of the orbit. After some rearrangement, theenergy from equation (62) can be rewritten in the form2 u E = − E − J (65) ± γx (cid:0) f − η − Ef − J f (cid:1) / + γ x f . Since J = O ( λ ), it is often convenient to work in the limit J → u E = − E ± γx (cid:0) f − η − Ef (cid:1) / + γ x f + O ( J ) . (66)This form for the energy shows why the product am of the effective semimajor axis and the massis slowly varying: In dimensionless units, the energy E = − m/ a so that( am ) − = − E u = E + O ( λ ) . (67)The product am is thus nearly constant as long as the time scale ratio λ is small, and the departureis of order λ . When λ is small, the orbit cycles through many turning points before the masschanges substantially, so that the average of the above equation becomes h ( am ) − i = h− E u i = E + O (cid:0) λ (cid:1) , (68)so that the average of am is constant to second order. 19 – Number of Cycles:
If we ignore J for now and integrate equation (59) over one cycle, we obtain γαE / Z f df ( f − f ) / ( f − f ) / = Z dxx = 1 x − x . (69)The integral on the LHS gives us π/E . If we integrate over N cycles we obtain the expression γαN π = E / [1 − m αN ] , (70)where we assume that m = 1 at the start. The total number of possible cycles occurs when m N → N T = E / πγα = E / πγ ( β + 1) . (71)Since E ∼ π ( β + 1) ∼
10, whereas γ ≪
1, we expect the number of cycles N T ∼ / (10 γ ) tooften be large. The Last Cycle:
The above analysis (if we continue to work in the regime J ≪
1) suggeststhat the last cycle occurs when the right hand side of equation (69) is no longer large enough tobalance the left hand side, which is the same for each cycle. As a result, a minimum mass must beleft in the star in order for the orbit to complete a cycle (in the function f ). This condition can bewritten in the form m c = (cid:16) πγαE / (cid:17) /α , (72)which holds for α = 1 + β >
0. At the point when the mass falls below this threshold, the timescale ratio is given by λ = ( Ef ) / / ( πα ). Since Ef ∼ πα ∼ −
10, the system crosses thethreshold so that f cannot complete a cycle just before the time scale ratio λ reaches unity. Final States:
If we set the energy equal to zero and replace the variable x with the time scaleratio λ f (evaluated at the moment that planet becomes unbound), we find the condition(2 f − η ) / = λ f f / ± (cid:0) f − η − Ef − J f (cid:1) / , (73)which can then be written in the form λ f = (2 f − η ) / ± (cid:0) f − η − Ef − J f (cid:1) / f / . (74)This expression is thus a generalization of that obtained for the special case with β = 0 (see equation[44] and Figure 3.2). The differences are that we have included the extra term J and that the resultis written in terms of the variable x instead of u .The orbital eccentricity, calculated the usual way, oscillates with time with an increasingamplitude of oscillation (e.g., see Figure 3.4). As shown here, however, the function f executesnearly Keplerian behavior, with nearly constant turning points, where this statement is exact inthe limit J →
0. The oscillation of eccentricity, although technically correct, is misleading. Theturning points of the orbit (in the original radial variable ξ ) are strictly increasing functions of 20 – Fig. 3.— Evolution of the orbit during the epoch of stellar mass loss. In this example, the mass lossfunction has index β = 2, corresponding to a constant mass loss rate. The other parameters are γ =10 − , e = 0.3, and f = 1 (initially going inward). The black curve shows the function f ( u ) = ξ/u ;the red horizontal lines mark the analytically determined turning points of the function. The bluedotted curve shows the evolution of the radial coordinate ξ (plotted here as log [ ξ ]). The magentadot-dashed curve shows the evolution of the time scale ratio λ . Finally, the green dashed curveshows the eccentricity e of the orbit. 21 –time. The oscillation in e arises because the orbits are not ellipses, and, in part, because the periodof the orbits in ξ are not the same as the period of the orbits in f . As a result, the oscillationsin the calculated eccentricity do not imply that the near-elliptical shape of the orbit is varyingbetween states of greater and lesser elongation. Instead, these oscillations imply that if the massloss stops and the orbit once again becomes an ellipse, the value of the final eccentricity of thatellipse oscillates with the ending time of the mass loss epoch. Orbital Elements during and after the Mass Loss Epoch:
Next we consider the case where theplanet remains bound after the epoch of stellar mass loss. In this case, we want to estimate theorbital elements of the planet. We are mostly interested in specifying the orbital elements at theend of the mass loss epoch, but we can also evaluate them at any time while the star continues tolose mass. Suppose that the orbit passes through N turning points of the function f during themass loss epoch. The orbit will then complete a partial cycle so that the final value of f lies betweenthe turning points, f ≤ f ≤ f . In the ideal case, where we have complete information describingboth the starting orbital elements and the final value of the stellar mass, and where the mass lossfunction is is described exactly by a model with constant index β , we can calculate the final value f f . In practice, we will often have incomplete information: The number of cycles is generally large, N ≫ β for a precise number N of cycles (and then stops abruptly). As result,we are unlikely to know where the final value of f lies between the turning points. Additionalplanets, or other perturbations, increase this uncertainty (see Section 4). In this case of incompleteinformation, we can write the mean value of the energy (averaged over the cycles) in the form2 u E = − E + γ x E , (75)where we have replaced f with the value 1 /E of its effective semimajor axis, and where the remainingterm averages to zero. This estimate for the final energy has an uncertainty — a range of possiblevalues — due to the lack of knowledge of where the planet lies in its orbit during the final cycle.This range of energy is given by the form∆ EE = ± γx (2 f − η − Ef ) / u E (76)= ± γxe ∗ E / − γ x /E / , where e ∗ is defined by equations (39) and (40). With the energy specified, the final value a f of thesemimajor axis is given by a f = − u E f . (77)The expected value of the energy E f is given by equation (75), but it can take on any value in therange defined by equation (76). Similarly, the final value e f of the orbital eccentricity is given by e f = 1 + η (2 u E f ) ≈ − ηE + ηγ x /E , (78) 22 –where the second (approximate) equality holds for the mean value of energy given by equation (75).Since the energy E f can have a range of values, given by equation (76), the orbit has a correspondingrange of possible eccentricities.The expressions derived above for the final orbital elements are expected to be valid providedthat the time scale ratio λ is small compared to unity (and hence | J | ≪ β = 2 and mass loss parameter γ = 10 − . Numerical integration of the full equation of motion(solid curves) show that the semimajor axis and eccentricity both oscillate and (on average) growwith time. (Note that the Figure plots log[ a ].) The values of the elements ( a, e ) calculated fromthe average energy (from equation [75]) provide a good approximation to the mean evolution ofthe orbital elements (see the central dotted curves in Figure 3.4). Furthermore, using the rangeof allowed energy calculated from equation (76), we can calculate upper and lower limits to theexpected behavior of the semimajor axis and eccentricity (shown as the upper and lower dottedcurves). Note that the solutions for a and e oscillate back and forth between these limiting curves.In this example, the planet becomes unbound near u = 28. Prior to that epoch, near u ≈
20, thelimits on the energy allow for the planet to become unbound, and the upper limit for the semimajoraxis approaches infinity. The approximation scheme thus breaks down at this point.
The equations of motion can be numerically integrated to find the value ξ f of the radialcoordinate when the system becomes unbound (when the energy becomes positive). For the caseof exponential mass loss, β = 1, the result is shown in Figure 3.5 as a function of the mass lossparameter γ (top panel). The figure shows curves for initially circular orbits ( η = 1, smooth curve)and for nonzero starting eccentricity ( η = 0.9, rapidly oscillating curve). Note that the η = 0.9curve is shown only for γ > . γ , the curve oscillates more quickly as afunction of γ , and the curve would appear as a solid black band in the figure. The bottom panelshows the value of λ , the ratio of the orbital period to the mass loss time scale, evaluated when thesystem becomes unbound. As expected, the parameter λ is of order unity when the system energybecomes positive and the planet becomes unbound. For the case of circular orbits at the initialepoch ( η = 1), the value of λ ∼ γ . For starting orbits with nonzero eccentricity, thevalue of λ takes on a range of values, but remains of order unity. For the case shown ( η = 0.9,oscillating curve), the parameter λ varies between about 1/2 and 2.The above trend holds over a range of values for the mass loss index β . Figure 3.5 shows thevalue of the time/mass variable u = 1 /m when the planet becomes unbound as a function of themass loss parameter γ . Results are shown for mass loss indices in the range 0 ≤ β ≤
3. For allvalues of the index β , the curves become nearly straight lines in the log-log plot for small values of γ , which indicates nearly power-law behavior of the form u f ∼ γ − p , where the index p = 1 / ( β + 1).We can find a simple fitting function for the final value u f of the inverse mass variable (for initially 23 – Fig. 4.— Evolution of orbital elements during the epoch of stellar mass loss. In this example,the mass loss function has index β = 2. The other parameters are γ = 10 − , e = 0 .
3, and f = 1 (initially going inward). The solid curves show the semimajor axis (log a ; top) and orbitaleccentricity ( e ; bottom), calculated from numerical integration of the equation of motion. For eachorbital element, the three dotted curves show the average value, the upper limit, and the lowerlimit, as calculated from the analytic expressions given by equations (75 – 78). 24 – Fig. 5.— Radial coordinate of the planet at the moment when the system becomes unbound, shownhere as a function of the mass loss parameter γ for exponential mass loss (top panel). The nearlymonotonic curve shows the result for the case of circular starting orbits; the oscillating curve showsthe result for angular momentum parameter η = 0.9 (which corresponds to starting eccentricity e = √ . ≈ . . . . ). Bottom panel shows the value of the parameter λ when the planet becomesunbound, for both circular starting orbits (smooth curve) and eccentric orbits ( η = 0.9; oscillatingcurve). All orbits are started at periastron. 25 –circular orbits) as a function of γ : u f = (cid:20) c β β (cid:21) γ − / (1+ β ) , (79)where c is a constant. For c = 0.74212, the error is less than 0.4% for mass loss indices in the range0 ≤ β ≤
4. Power-law fits resulting from equation (79) are shown as the dashed lines in Figure 3.5.The power-law fits provide a good approximation for γ < ∼ .
1. For larger values of γ > ∼ .
1, the finalvalue u f →
2. Equation (79) describes the final mass values (for initially circular orbits) for realisticmass loss prescriptions: Most of the stellar mass is usually lost on the asymptotic giant branch(AGB), where the mass loss function has β ∼ constant (see Fig. 13 of Veras et al. 2011). Notethat mass loss on the AGB takes place through a series of pulses, but this complication primarilyaffects tides (Mustill & Villaver 2012), rather than the overall mass loss profile.A related result in shown in Figure 3.5, which plots the values of the ratio λ of time scales,evaluated at the moment when the planet becomes unbound, as a function of the mass loss param-eter γ . In the limit of small γ , the time scale ratio λ approaches a constant value (of order unity).The finding that λ has a value of order unity (in the limit of small γ ) when the planet becomesunbound is expected: In physical terms, this result means that the mass loss time scale has becomeshorter than the orbital period, so that the potential well provided by the star is changing fastenough that the orbital motion does not average it out.The limiting values of the time scale ratio λ are shown in Figure 3.5 as a function of the massloss index β . Here the time scale ratios are evaluated at the moment when the planet becomesunbound. Results are shown for the limiting case of small γ (from Figure 3.5 we see that the timescale ratio λ approaches a constant value as γ → β = 0, the final value of the time scale ratio λ = 1. Since this function λ f ( β )is useful for analysis of orbits in systems losing mass, we provide a simple fit. If we choose a fittingfunction of the form log λ f = c β + c β , (80)where c and c are constants, we obtain a good fit with the values c = 0.21658 and c = 0.04102.The fitting function is shown as the dashed curve in Figure 3.5, and is almost indistinguishablefrom the numerically integrated solid curve.
4. Lyapunov Exponents for Two Planet Systems with Mass Loss
The discussion thus far has focused on single planet systems, whereas many solar systemscontain multiple planets. In order to see how multiple planets affect orbital evolution during massloss, we generalize the treatment to study systems consisting of two planets and a central starwith decreasing mass. Such a 3-body configuration represents a crude model for our Solar System,where the motions of only the three most dominant objects (Jupiter, Saturn and the Sun) areconsidered. As a starting point, we fix the planetary masses and initial orbital elements ( a, e ) to 26 –
Fig. 6.— Value of the time/mass variable u = 1 /m at the moment when the system becomesunbound, shown here as a function of the mass loss parameter γ , for varying values of the index β .All of the curves correspond to circular starting orbits with η = 1 (which corresponds to startingeccentricity e = 0). The curves correspond to values of β = 0 (top curve) to β = 3 (bottom curve).Dashed (black) lines correspond to the fitting function of equation (79). 27 – Fig. 7.— Value of the time scale ratio λ evaluated at the moment when the system becomesunbound, shown here as a function of the mass loss parameter γ , for varying values of the index β . All of the curves correspond to circular starting orbits η = 1 (which corresponds to startingeccentricity e = 0). The curves correspond to values of β = − . -1 0 1 2 3 Fig. 8.— Value of the time scale ratio λ evaluated at the moment when the system becomesunbound, shown here as a function of the mass loss index β which defines the time dependence ofstellar mass loss. These values correspond to the limit of small mass loss parameter γ → e = 0. The dashed curve, which is nearlyidentical to the solid curve, shows a simple fit to the function λ f ( β ), as described in the text. Asshown in the text, λ f = 1 for the particular case β = 0. 29 –those of Jupiter and Saturn, and set the initial stellar mass to M ∗ = 1.0 M ⊙ . We also restrict theorbits to a plane, thereby reducing the number of phase space variables from 18 to 12. To start,the mass loss function is taken to be an exponential model with index β = 1, although this law isgeneralized later.As long as the system suffers no close encounters, the orbit of an individual planet is similarto that described by the variable mass two-body problem. An example of the evolution of theosculating orbital elements ( a, e ) for our benchmark system (see above) is shown in Figure 4 fora mass loss time scale of 10 yr. As each planet orbits in an outward spiral, the semimajor axisincreases approximately exponentially in time (in inverse proportion to the stellar mass). Theeccentricity oscillates rapidly on orbital time scales, and more slowly on secular time scales (as theplanets exchange angular momentum), but remains close to its starting value until stellar massloss has taken place for a few e-folding times. The product of the semimajor axis and stellar mass( aM ∗ ) is approximately constant until a few e-folding times have elapsed. After a critical amountof mass is lost, the orbital elements a → ∞ and e →
1, and planets can become unbound. Noticethat at this point, the stellar mass is only a few percent of the initial value; as a result, this scenariois rather artificial for stars like our Sun, which are only expected to lose about half of their initialmasses. However, larger stars lose a greater fraction of their original masses. For example, a starwith initial mass M ∗ ≈ M ⊙ is expected to end its life as a white dwarf with roughly ∼
15% ofits original mass (where the final mass fraction depends on the stellar metallicity).The evolution of the orbital elements can differ dramatically if the planets reach a smallenough separation so that orbital crossings can occur. In this regime, chaos dominates and theorbital elements evolve in a less predictable manner. An example is shown in Figure 4 for the sameparameters used in Figure 4, but with the initial eccentricity of the inner planet (Jupiter) increasedto e = 0 . t = τ this stable evolution has been compromised.Studies of dynamical systems generally use the maximum Lyapunov exponent as an indicationof the level of chaos present (e.g., Lichtenberg & Lieberman 1983; Strogatz 1994). If a systemis chaotic, two nearby trajectories in phase space initially differing by a small amount δ shoulddiverge according to δ ( t ) = δ e Λ t with Λ > . (81)The Lyapunov time is thus τ ly = 1 / Λ. The long-term dynamical stability of the solar systemhas been explored in the absence of stellar evolution (Batygin & Laughlin 2008; Laskar 2008), andcurrent estimates of the Lyapunov time for the solar system (while the mass of the Sun remainsconstant) are τ ≈ δ . By integrating both systems simultaneously and 30 – a ( AU ) t/τ e Fig. 9.— Osculating semimajor axis and orbital eccentricity for a pair of planets orbiting an initiallysolar-mass star with mass loss time scale τ = 10 years. Planets have masses, initial semimajoraxis and eccentricities of Jupiter and Saturn. The orbital elements evolve in a roughly predictablemanner, with the semimajor axes increasing smoothly and the eccentricities oscillating on seculartime scales, but remaining relatively constant until the star has lost the majority of its initial mass. 31 – a ( AU ) t/τ e Fig. 10.— Same as Figure 4, but with the initial eccentricity of Jupiter increased to e = 0 .
3. Thisincrease in eccentricity allows for orbital crossings and increases chaotic behavior. 32 –monitoring the quantity δ ( t ), we can calculate the divergence of the neighboring trajectories andthen estimate the maximum Lyapunov exponent. Since a three-body system restricted to a planeconsists of 12 phase space variables, there is some choice in defining the quantity δ . For the sakeof definiteness, we define δ according to δ = p ( x r − x s ) + ( y r − y s ) , (82)where ( x, y ) are Cartesian coordinates for a planet’s location, and where the subscripts r and s refer to the “real” and “shadow” trajectories respectively. Since chaotic systems display complicatedbehavior, the functions δ ( t ) will vary for effectively equivalent cases. As a result, for each systemof interest, we run 1000 cases with both a “real” and a “shadow” system. To extract the Lyapunovexponent, we can either average together the 1000 runs to construct a single function δ ( t ) and usethe result to find the exponent, or, we can find the exponent from each of the 1000 individual casesand then find the average exponent. Both schemes produce the same values; here we present resultsfor the former case. Figure 4 shows an example of the time evolution of δ ( t ) for different valuesof the mass loss time scale τ . After an initial period of transient growth (roughly delimited by t < ∼ . τ ), the divergence metric δ ( t ) increases exponentially with time and the Lyapunov exponentcan be obtained by finding the slope of the line defined by ln δ = ln δ + Λ t .For each value of the mass loss time scale τ , the maximum Lyapunov exponent was calculated asoutlined above. Since the maximum possible separation between the reference and shadow systemsis finite (using the definition of δ in equation [82]), the curves of growth eventually saturate. Thus,to extract the Lyapunov exponent Λ, we want to measure the curves of divergence after the initialinterval of transient behavior but before saturation occurs. In most cases, Λ was calculated fromthe time-series data for times τ / ≤ t ≤ τ ; this time interval is delimited by the vertical dashedlines in Figure 4. An exception was made for the case of extremely slow mass loss, however, where τ = 100 Myr. For this scenario, the time scale for mass loss is longer than the “natural” Lyapunovtime of ∼
10 Myr (the value obtained without mass loss), and the curves of divergence saturatebefore t = 0 . τ . In this case, Λ was calculated only for time series data with t <
10 Myr. Noticethat the curves shown in Figure 4 are not perfectly straight in the region between the dashed lines;this curvature introduces some uncertainty in the specification of the Lyapunov time scales. Toestimate this uncertainty, we have calculated the Lyapunov time values for a wide range of choicesfor the time intervals. This procedure implies an uncertainty of a few percent.Figure 4 shows our numerical values of the Lyapunov times τ ly as a function of the mass losstime τ . We performed the analysis described above for each of the two planets in the systemseparately and averaged the results. The horizontal dotted line – included here for reference –corresponds to our numerically determined Lyapunov time for the model solar system in the absenceof stellar mass loss. This value is in relatively good agreement with previous calculations for thecomplete Solar System (Sussman & Wisdom 1992), but differs slightly because our model considersonly two of the four giant planets. Note that as τ → ∞ , the Lyapunov time approaches the dottedline, i.e., the value expected with no mass loss. The solid curve shows the calculated values of theLyapunov time, whereas the dashed line indicates the least-squares fit (where the fit was taken over 33 – t/τ −5−4−3−2−10123 l n [ δ ] Fig. 11.— Curves showing the divergence of two trajectories separated by a small distance δ fordifferent values of the mass loss time scale τ . Black curve is τ = 10 ; blue is τ = 10 ; purpleis τ = 10 ; and red is τ = 10 (yr). After a period of initial growth, the trajectories divergeexponentially, indicated by the linear shape of the latter portions of the graphs. The regionsbetween the dashed lines were used to calculate the Lyapunov exponent. Note that the timevariable has been scaled by the mass loss time scale τ . 34 –range of mass loss time scales τ ≤ ). Notice that the slope of this line is close to unity. Morespecifically, we obtained τ ly ∼ τ p where p = 0 . . (83)Note that this fitted line cannot be meaningfully extrapolated below τ = 10 − . In this regime,the mass loss time scale τ becomes comparable to the orbital periods of the planets, and thedynamics of even single-planet systems becomes complicated (see the previous section).As a consistency check, we also explored other choices for the metric δ that measures thedifference between nearby trajectories. An especially compelling option is to use the semimajoraxis a , because unlike the physical distance between the reference and shadow trajectories, there isno upper limit on this quantity. The previous calculation was thus repeated using δ = | a r − a s | . (84)Our results are similar, however, which indicates that the Lyapunov times do not depend sensitivelyon the choice of δ .Next we would like to ensure that the (nearly) linear relation between the Lyapunov time andmass loss time is not an artifact of the exponential ( β = 1) mass loss law that was chosen. Towardthis end, we have explored two additional functional forms for the mass loss law: The first usedvanishing mass loss index β = 0 (see equation [19]), whereas the second used a constant mass lossrate with β = 2 (see equation [17]). For both of these mass loss functions, we obtained nearly thesame power-law relation for the Lyapnunov time scale versus the mass loss time scale, i.e., τ ly ∼ τ p ,where p = 0 .
98 and p = 1 .
01 using equations (19) and (17), respectively. The results, shown inFigure 4, are thus nearly identical, independent of the index β of the mass loss function. Thisfinding suggests that this power-law trend is robust.Since the Lyapunov time scale is found to be comparable to the time scale for mass loss, wegenerally expect such solar systems to be only moderately influenced by chaos. In order for chaos tofully erase initial conditions for a dynamical system, nearby trajectories in phase space must divergefor several Lyapunov times. For example, in order for a starting uncertainty of 1 ◦ in position angleto grow into 360 ◦ , one needs ∼ τ (yr) τ l y ( y r ) Fig. 12.— Calculated Lyapunov times τ ly versus mass loss time scales τ for a star losing massexponentially in time (mass loss index β = 1). The calculated quantities are shown (solid curve)along with a least-squares fit for τ ≤ yr (dashed line). As τ → ∞ , the Lyapunov timeapproaches that of the solar system with constant stellar mass ( τ ly ∼ τ (yr) τ l y ( y r ) Fig. 13.— Lyapunov time scale τ ly versus mass loss time scale τ for systems where the star losesmass through three different decay laws. From top to bottom, the blue curve shows the results for β = 0; the black curve shows the results for exponential mass loss ( β = 1); and the purple curveshows the results for constant mass loss rate ( β = 2). For all three cases, the square symbols showthe numerically calculated time scales and the dashed lines show a least-squares fit. These threeexamples show similar behavior, which indicates that the linear relationship between Lyapunovtime scales and mass loss time scales is largely independent of the particular mass loss formula. 37 –
5. Applications
To illustrate the efficacy of the results found in the previous sections, we consider two repre-sentative astronomical problems. For solar type stars, we find an effective outer edge of the solarsystem, i.e., the boundary between planetary bodies that remain bound after the epoch of massloss and those that escape (Section 5.1). Next we consider planets that remain bound to whitedwarfs, and find their final orbital elements (Section 5.2).
Suppose that a solar-type star, with initial mass M ∗ = 1 M ⊙ , loses some portion of its massover a time interval ∆ t = 1 Myr, so that the fraction m f remains afterward. Since solar type starslose most of their mass while they are on either the Red Giant Branch or Asymptotic Giant Branch(see Hurley et al. 2000), we expect β to be in the range 1–3. For a fixed time interval ∆ t andarbitrary index β >
1, the mass loss parameter γ is given by γ = 1 β − h − m β − f i t (cid:18) a GM ∗ (cid:19) / (85)= 1 . × − β − h − m β − f i (cid:16) a (cid:17) / , where a is the initial semimajor axis of the orbit. For initially circular orbits, we can use the resultsof Section 3.5 to find the conditions for which planets become unbound. Using equation (79) todefine the critical value of γ , and equating the resulting value to the expression from equation (85),we can solve for the critical semimajor axis a c , such that orbits with larger values of a becomeunbound during the mass loss epoch. First we define the constant A β ≡ ( β − / (cid:20) c β β (cid:21) β ) / . (86)The critical semimajor axis a c is then given by a c = (∆ t ) / ( GM ∗ ) / A β m βf h − m β − f i / (87) ≈ (34 , A β m βf h − m β − f i / . We can thus find the critical semimajor axis a c for any given value of the index β and the remainingmass fraction m f . Note that the resulting value of a c is a sensitive function of the index β . Toleading order, a c ∼ (34,000 AU) m β ) / f . If we use m f = 1/2, the value expected for the Sun, 38 –then the critical value of the semimajor axis a c ≈ β = 2,3, and 4 (these values for a c are in general agreement with the results of Veras & Wyatt 2012).For planets that escape, the corresponding velocities are small, in the range 0.3 – 0.5 km/s. Forinitially eccentric orbits, planets will become unbound for a wide range of starting semimajor axes,depending on the initial phase of the orbit (see Section 3); however, this range is centered on themean values found here.Note that in the limit β → ∞ , we formally obtain a c → m f <
1. However, this formulation of the problem breaks down before that limit isreached: For large values of β , even though the time interval is fixed, the mass loss rate acceleratesrapidly so that most of the mass is lost near the end of the time interval. Stellar mass is thus lost(effectively) through a step function in the limit of large β . In this limit, the results of Section 3.3apply, and circular orbits become unbound (remain bound) for mass fraction m f < / m f > / Now consider a progenitor star with initial mass M ∗ = 5 M ⊙ , which evolves into a whitedwarf with mass M wd = 1 M ⊙ (Hurley et al. 2000); the final mass fraction m f = 1/5 ( u f = 5). Forpurposes of illustration, we assume that the mass is lost over a single epoch that can be describedby a single value of the mass loss index β , with a time scale ∆ t = 1 Myr. For orbits with startingsemimajor axes a and eccentricities e , we would like to know the final orbital elements, after theepoch of mass loss. For orbits with starting semimajor axis a in the range 1 −
100 AU (closerplanets are often accreted by the star), the value of γ falls in the range γ = 10 − − − . Themass loss parameter γ is thus small and nearly independent of the index β (see equation [85]). Thetime scale ratio λ = γu β +1 f / . For the systems of interest, the largest γ value is thus ∼ − , thelargest value of u = 5, and the largest value of β = 3; the largest value of the time scale ratio isthus λ ∼ . λ ≤ . J is small, where J = O ( λ ) ≤ a, e ), but also the phase of the orbit at t = 0. This latter quantity is specified by the initial valueof the radial coordinate ξ = f , which lies in the range 1 − e ≤ ξ ≤ e . With ξ ( f ) determined,the integration constant E is given by equation (63).In the limit J →
0, the equation of motion shows that the function f oscillates back andforth between its turning points (given by equation [38]) while the system loses mass (analogousto the evolution depicted in Figure 3.4). If the duration of the mass loss epoch is specified exactly(equivalently, if m f = 1 /u f is known exactly), then we can determine the final value of the function A full specification would also include three additional angles, e.g., the longitude of periastron, the longitude ofthe ascending node, and the inclination angle, but we can orient the coordinate system to eliminate this complication.
39 – f f = f ( x f ), where x f = u αf . Let N c denote the number of complete cycles that that function f ( u )executes during the mass loss phase, where a cycle is defined as motion from one turning point tothe other (a half orbit). In addition, the orbit must turn through a partial cycle from its startingvalue f to the first turning point f j (where j = 1,2) and must turn through another partial cyclefrom the final turning point f k (where k = 1,2) to the final value f f . After integrating the equationof motion (61), we thus obtain E / παγ (cid:2) − m αf (cid:3) = ( δN ) + N c + ( δN ) f , (88)where the integration constant E is given by equation (37) and we have defined( δN ) ≡ ± Eπ Z f j f f df ( f − f ) / ( f − f ) / , (89)and ( δN ) f ≡ ± Eπ Z f f f k f df ( f − f ) / ( f − f ) / , (90)where the ± signs are chosen to keep the integrals positive. Note that equations (88 – 90) completelyspecify the final value of the function f f . After solving these equations for f f , we can use equation(66) to find the final value E f of the orbital energy, and then use equation (77) and (78) to find thesemimajor axis and eccentricity of the orbit after mass loss has ended.Although the procedure described above is exact in principle (subject to the approximation J → δN ) and ( δN ) f are less than unityby definition, whereas N c (and the left-hand side of the equation) is much larger, i.e., N c ∼ − for the orbits considered here. To find the value f f necessary to specify the phase of the final orbit,this approximate description for mass loss must be correct to better than 1 part in 10 (10 )for orbits with starting a = 100 AU (1 AU). It is unlikely that the mass loss function for a realastronomical system obeys this simple model to such a high degree of fidelity. As result, eventhough we have an exact solution for the model equation, we cannot predict with certainty thefinal phase of the orbit for realistic systems.Given the uncertainty outlined above, post-mass-loss orbits can be described in terms of theexpected values of the energy E f , semimajor axis a f , and eccentricity e f , as well as the possibleranges of values for the orbital elements (given the range of phases). For the sake of definiteness,we assume that the orbit starts at periastron (at the start of the mass loss epoch) and that themass loss index β = 3. The orbit has starting angular momentum η = 1 − e . The integrationconstant E is then given by E = 1 − γ (1 − e ) (from equation [36]) and the energy E f at the endof the mass loss epoch (from equation [75]) becomes E f = 150 (cid:8) − , γ + (1 − e ) γ + O ( γ ) (cid:9) . (91) 40 –Note that we can ignore the second γ term in the above expression. The expected final value ofthe (dimensionless) semimajor axis (from equation [77]) is given by a f = 5 (cid:2) − , γ (cid:3) − , (92)and the corresponding expected value of the eccentricity (from equation [78]) is given by e f = e + 390 , ηγ , (93)where we work to the same order of approximation as for a f (recall that e is the starting, pre-mass-loss value of the eccentricity). Since γ ∼ − ( a / / , the correction term is small: 390,625 γ ∼ × − ( a / ≤ .
004 since a ≤
100 AU. The mean value of final semimajor axisis thus about 5 times the starting value, as expected, and the leading order correction has beenquantified. The square of the eccentricity increases by a similar amount. Because of the range ofpossible orbital phases at the end of mass loss, the orbital elements can differ from these meanvalues according to the relations∆ a f a f = ∆( e f ) e f = ∆ E f E f ≈ ± γe ∗ E / , (94)where we have used equation (76) and where e ∗ = e + γ (1 − e ) is the effective eccentricity of thefunction f ( u ) during the epoch of mass loss. Thus, the total (relative) width of the distributionof possible final orbital elements is thus ∼ γe . Note that this range is often larger than thecorrection to the mean values. Consider a planet with starting semimajor axis a = 100 AU andeccentricity e = 0.30. The mean value of the final semimajor axis is a f ≈
502 AU, only 2 AU largerthan the value suggested by the simple scaling law aM ∗ ≈ constant that is often used. However, therange of possible values about this mean is about ∆ a f = ±
19 AU. Similarly, the final eccentricityhas mean value e f ≈ . ± .
6. Conclusion6.1. Summary of Results
This paper has reexamined the classic problem of the evolution of planetary orbits in thepresence of stellar mass loss. Although this issue has been addressed in previous studies (seeSection 1), we generalize existing work to include a new analytic formulation for time-dependentmass loss and to determine Lyapunov time scales for multiple planet systems. In particular, weconsider a class of model equations where the mass loss index β is constant (see equation [10]),which allows for a wide range of time dependence for the mass loss rates and allows for a numberof new results to be obtained analytically. Our main results can be summarized as follows:[1] Previous numerical studies show that planetary orbits often obey the approximate law am ≈ constant , where a is the semimajor axis, and where this approximation holds as long as the 41 –time scale for mass loss is significantly longer than the orbital period. By writing the equation ofmotion in the form given by (22) and (23), we show analytically why this law holds (see equations[65 – 68]). In addition, the differential equation (12) for the energy E shows that the energy is astrictly increasing function of time, so that the semimajor axis (defined via a ∼ / |E| ) increasesmonotonically (whereas the orbital radius ξ oscillates in and out).[2] Previous literature often claims that the orbital eccentricity remains constant during theearly phases of stellar mass loss. In contrast, this work shows that the eccentricity oscillates backand forth between well defined limits during the phase of mass loss, and that the amplitude ofthese oscillations grow with time (one example is shown in Figure 3.4). Moreover, the upper andlower limits of the eccentricity range can be calculated analytically using equations (75 – 78). Notethat these oscillations in the eccentricity, while technically correct, result from assigning orbitalelements (which describe ellipses) to orbital paths that are not elliptical. The actual orbit expandswith time, and, for example, the outer turning point of the orbit increases monotonically (it doesnot oscillate).[3] In the limit of rapid mass loss, λ → ∞ , we obtain analytic solutions that describe orbitsfor the entire class of mass loss functions (see Section 3.3). The condition for the planet becomingunbound is given by equation (49). For planets that remain bound, the new orbital elements aregiven by equation (55).[4] Not all mass loss functions lead to planets becoming unbound (except, of course, in theextreme case where the stellar mass vanishes m → β = −
1, where systems with mass loss characterized by β < − m → β = − β = 0, we can find analyticexpressions for the function f ( u ) and for the final values of the time scale ratio λ f when the planetbecomes unbound (see Section 3.2). For arbitrary values of the mass loss index β = 0, this approachcan be generalized to find analytic expressions for f ( u ) and the orbital elements of the planet (seeSection 3.4). The resulting expressions are approximate, correct to order O ( λ ), and are thusaccurate over most of the mass loss epoch.[6] One way to characterize the dynamics of these systems is through the parameter λ f . Wedefine λ to be the ratio of dimensionless mass loss rate to the orbital frequency, and λ f is the valueat the moment when the planet becomes unbound. For initially circular orbits, the parameter λ f is always of order unity, and approaches a constant value in the limit of small mass loss rates γ ;further, the value of this constant varies slowly with varying β (see Figure 3.5). For orbits startingwith nonzero eccentricity, however, the parameter λ f can depart substantially from unity and variessignificantly with β (compare Figures 3.2 and 3.2).[7] For multiple planet systems, we find that the Lyapunov times decrease in the presence ofstellar mass loss, so that chaos should play a larger role in planetary dynamics as stars leave the 42 –main sequence. In fact, the Lyapunov time scale is proportional to (but somewhat shorter than)the mass loss time scale over a range of conditions (see Figures 4 and 4). For a typical mass losstime scale of τ ∼ yr, the Lyapunov time is ∼ − × yr. Three different forms for the stellarmass loss function have been considered, all yielding similar results, which suggests that this trendis robust. In spite of its apparent simplicity, the classic problem of planetary orbits with stellar mass lossis dynamically rich. Here we discuss two issues that have been highlighted by this present work:The use of osculating orbital elements is a standard way to describe planetary motion. In thisscheme, the planet is assigned the Keplerian orbital elements that it would have if it were movingwithin a purely Keplerian potential. In systems with stellar mass loss, however, this approachmight be more misleading than illuminating (see also Radzievskii & Gel’Fgat 1957, Hadjidemetriou1963). Consider, for example, the osculating eccentricity of the orbit. As the star loses mass,the eccentricity oscillates (Figure 3.4) with an amplitude that grows with time. This oscillatingeccentricity would seem to imply motion along an elliptical path, where the shape of the ellipsecycles back and forth between being more elongated and more round. However, this “oscillation” ofthe eccentricity is an artifact of assigning a Keplerian orbital element ( e ) to orbital motion that isnot Keplerian. During the mass loss epoch, the orbit has inner and outer turning points, analogousto those of an ellipse. But the turning points of the actual orbit do not oscillate – the orbit smoothlyspirals outward (e.g., see Figure 2 of Veras et al. 2011).An alternate description of the dynamics can be constructed: Before the onset of mass loss, theorbit is a Keplerian ellipse and can be described by the usual orbital elements ( a, e ), which providethe initial conditions for the next stage (along with the phase of the orbit at t = 0). During themass loss epoch, the orbit is not an ellipse, but the scaled radial function f = ξm = ξ/u follows anearly Keplerian trajectory. For the particular mass loss function with index β = 0, this analogy isexact. In this case, the function f obeys the same equation of motion (35) as the radial coordinate ξ in the Kepler problem (albeit with a nonconventional time variable where dt ∼ du/u ). Thescaled function f has turning points (equation [38]), an effective semimajor axis a ∗ = 1 /E (seeequation [39]), and an effective eccentricity e ∗ (equation [40]). In the general case where β = 0,the scaled function f executes nearly Keplerian motion: Here the equation of motion for f has aKeplerian form up to a correction of order O ( J ) = O ( λ ), where λ ≪ a ∗ , e ∗ )characterizing the function f are constant, to O ( λ ), while the star loses mass. After the epochof mass loss, the planet (if it remains bound) once again enters into a Keplerian orbit, now withelements ( a f , e f ). The problem can thus be described in terms of three sets of orbital elements:before ( a, e ), during ( a ∗ , e ∗ ), and after ( a f , e f ) mass loss. But the orbital elements during mass loss( a ∗ , e ∗ ) correspond to orbits of the scaled function f (the orbit in physical space is not Keplerian). 43 –Another interesting complication arises: These orbits display a type of sensitive dependenceon initial conditions even in the absence of chaos. For planets that remain bound after the massloss epoch, the final orbital elements ( a f , e f ) depend on the phase of the orbit, at both the startand the end of the mass loss epoch. However, the number of orbital cycles during mass loss is large, N c ∼ /γ ≫
1, and a precision of one part in N c is necessary to specify the final orbital phase (giventhe initial phase). The orbit can be calculated to sufficient accuracy to specify the final orbitalelements provided that the starting state, the duration of stellar mass loss, and the form of themass loss function are all known well enough. In practice, however, real astronomical systems willnot follow these particular forms to such high accuracy, so that the final orbital elements cannot bepredicted with certainty. Instead, the expectation value of the orbital elements can be predicted,along with the range of possible variation about their mean values (see also Section 5.2). There are many opportunities for future work. Analytic studies can be taken into two di-rections. First, the formulation developed here can be applied to a wide range of astronomicalsystems, including predictions of orbital elements for planets that remain bound after stellar massloss and predictions of the conditions required for bodies to become unbound. On the other hand,the analytic treatment can be developed further to include more general mass loss functions, mul-tiple phases of mass loss, and higher order approximations for the conditions under which planetsbecome unbound.For the calculations of the Lyapunov time scales, this paper has focused on analogs of our ownsolar system, and has considered only the motion of Jupiter, Saturn, and the Sun, since these are themost gravitationally dominant bodies. However, future calculations should also include additionalplanets and a wider range of starting orbital elements. In addition, this paper has consideredrelatively short integrations, spanning at most only 10 −
100 Myr. Although stellar mass loss isnot expected to continue for longer than 100 Myr, the increased semimajor axes of the planets willchange the overall dynamics and the decreased stellar mass (relative to the planets) could lead toan increase in dynamical instabilities. As a result, longer integrations should be performed, usingthe lower (constant) stellar mass and the increased semimajor axes of the planets as input. Thecalculations of this paper are limited to coplanar planeteary systems; future work should explorethe effects of different inclination angles. Finally, given the diversity exhibited in the observedsample of extrasolar planets, different planetary configurations should also be considered. Thesetypes of calculations will help us understand the long term fate of planetary systems in general andwill help direct future observations. 44 –
Acknowledgments
We would like to thank Jake Ketchum, Kaitlin Kratter, Dimitri Veras, and Eva Villaver foruseful discussions. This work was supported by NSF grant DMS-0806756 from the Division ofApplied Mathematics, NASA grant NNX11AK87G (FCA), and NSF grants DMS-0907949 andDMS-1207693 (AMB).
A. Bounds on the J Integral
In the limit | J | ≪
1, we have a completely analytic description of the dynamics. It is thususeful to place bounds on the integral J , which can be done as follows: First write J = γ βI where I = 2 Z x x f f x dx . (A1)Note that since the function f is of order unity, the order of the integral J is given by J = O (cid:0) γ x (cid:1) = O (cid:0) λ (cid:1) . (A2)Next we integrate by parts to obtain I = x f − f − Z x xf dx . (A3)The second integral can be written2 Z x xf dx = 2 h f i Z x xdx = h f i ( x − , (A4)where we have invoked the mean value theorem. The function f varies between turning points sothat f ≤ f ≤ f . We thus have bounds I ≤ x ( f − f ) + f − f ≤ x ( f − f ) , (A5)and I ≥ x ( f − f ) + f − f ≥ x ( f − f ) . (A6)As a result, we have the bound | I | ≤ x ( f − f ) = x ( f + f )( f − f ) = x a ∗ e ∗ . (A7)We thus obtain the desired bound on J , i.e., | J | ≤ βγ x ( f + f )( f − f ) = γ x βa ∗ e ∗ . (A8)In order to evaluate this bound, we need expressions for the turning points f and f , or, equiva-lently, the effective semimajor axis a ∗ and eccentricity e ∗ . As derived in the text, we have approx-imations for these quantities, where these expressions are exact in the limit J →
0. 45 –
REFERENCES
Adams, F. C., & Laughlin, G. 1997, Rev. Mod. Phys., 69, 337Batygin, K., & Laughlin, G. 2008, ApJ, 683, 1207Bear, E., & Soker, N. 2011, MNRAS, 414, 1788Debes, J. H., & Sigurdsson, S. 2002, ApJ, 572, 556Diacu, F., & Selaru, D. 1998, J. Math. Phys., 39, 6537Duff, M. J., Okun, L. B., & Veneziano, G. 2002, JHEP, 0203, 023Duncan, M. J., & Lissauer, J. J. 1998, Icarus, 134, 303Gyld´en, H. 1884, AN, 109, 1Hadjidemetriou, J. D. 1963 Icarus, 2, 440Hadjidemetriou, J. D. 1966 Icarus, 5, 34Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543Jeans, J. H. 1924, MNRAS, 85, 2Jura, M., 2003, ApJ, 584, L91Kratter, K. M., & Perets, H. B. 2012, ApJ, 753, 91Kudritzki, R. P., & Reimers, D. 1978, A&A, 70, 227Laskar, J. 2008, Icarus, 196, 1Lichtenberg, A. J., & Lieberman, M. A. 1983, Regular and Chaotic Dynamics (Springer-Verlag,New York)MacMillan, W. D. 1925, MNRAS, 85, 904Melis, C., Jura, M., Albert, L., Klein, B., & Zuckerman, B. 2010, ApJ, 722, 1078Moeckel, N., & Veras, D. 2012, MNRAS, 422, 831Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Princeton: Princeton Univ. Press)Mustill, A. J., & Villaver, E. 2012, ApJ, 761, 121Perets, H. B., & Kratter, K. M. 2012, arXiv:1203.2914Prieto, C., & Docobo, J. A. 1997, A&A, 318, 657 46 –Radzievskii, V. V., & Gel’Fgat, B. E. 1957, Soviet Astronomy, 1, 568Rahoma, W. A., Abd El-Salam, F. A. & Ahmed, M. K. 2009, JApA, 30, 187Schr¨oder, K.-P., & Connon Smith, R. 2008, MNRAS, 386, 155Spiegel, D. S., & Madhusudhan, N. 2012, ApJ, 756, 132Strogatz, S. H. 1994, Nonlinear Dynamics and Chaos (Addison-Wesley, MA)Sussman, G. J., & Wisdom, J. 1992, Science, 257, 56Uzan, J.-P. 2003, Rev. Mod. Phys., 75, 403Veras, D., Wyatt, M. C., Mustill, A. J., Bonsor, A., & Eldridge, J. J. 2011, MNRAS, 417, 2104Veras, D., & Tout, C. A. 2012, MNRAS, 422, 1648Veras, D., & Wyatt, M. C. 2012, MNRAS, 421, 2969Veras, D., Mustill, A. J., Bonsor, A., & Wyatt, M. C. 2013, MNRAS, in pressVillaver, E., & Livio, M. 2007, ApJ, 661, 1192Villaver, E., & Livio, M. 2009, ApJ, 705, 81Vinti, J. P. 1974, MNRAS, 169, 417Voyatzis, G., Hadjidemetriou, J. D., Veras, D., & Varvoglis, H. 2013, MNRAS, in pressZuckerman, B., Melis, C., Klein, B., Koester, D., & Jura, M. 2010, ApJ, 722, 725