Evolution of Planetary Systems with Time Dependent Stellar Mass Loss

Observations indicate that intermediate mass stars, binary stars, and stellar remnants often host planets; a complete explanation of these systems requires an understanding of how planetary orbits evolve as their central stars lose mass. Motivated by these dynamical systems, this paper generalizes in two directions previous studies of orbital evolution in planetary systems with stellar mass loss: [1] Many previous treatments 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 time dependent. The mass loss rate can be increasing or decreasing, but the stellar mass always decreases monotonically. For this class of models, we develop analytic approximations to specify the final orbital elements for planets that remain bound after the epoch 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 the asymptotic limit where the stellar mass vanishes. [2] We consider the chaotic evolution for two planet systems with stellar mass loss. Here we focus on a model consisting of analogs of Jupiter, Saturn, and the Sun. By monitoring the divergence of initially similar 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 in the range 5 - 10 Myr; we find that the Lyapunov time decreases with increasing stellar mass loss 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 dynamical evolution that occurs in solar systems with stellar mass loss.


Introduction
Solar systems orbiting other stars display a diverse set of architectures and motivate further studies concerning the dynamics of planetary systems. Part of the richness of this dynamical problem arises from the intrinsic complexity of N-body systems, even in the absence of additional forces (Murray & Dermott 1999). The ledger of physical behavior experienced by such systems is enormous, and includes mean motion resonances, secular interactions, and sensitive dependence on the initial conditions (chaos). Additional complications arise from additional forces that are often present: During early stages of evolution, circumstellar disks provide torques that influence orbital elements, and turbulent fluctuations act on young planets. Over longer time scales, solar systems are affected by tidal forces from both stars and planets, and by general relativistic corrections that lead to orbital precession. Another classic problem in solar system dynamics concerns planetary orbits around central stars that are losing mass (Gyldén 1884, Jeans 1924; see also Hadjidemetriou 1963Hadjidemetriou , 1966. Although this issue has received some recent attention (see below), this paper expands upon existing work in two main directions. Recent work often focuses on the particular case of constant 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 numerical simulations to consider more realistic, time-dependent stellar mass loss). In this paper, for single planet systems, we extend existing calculations to account for time dependence of the mass loss rates and obtain a number of analytic results. For systems with two or more planets, we also show that the Lyapunov exponents, which determine the time scales for chaotic behavior, depend on the time scales for mass loss. As outlined below, these two results can account for a great deal of the possible behavior in solar systems where the central star loses mass.
A number of previous studies have considered planetary dynamics for host stars that are losing mass. For our own Solar System, long term integrations have been carried out to study the fate of the planets in light of mass loss from the dying Sun (Duncan & Lissauer 1998). Recent related work estimates an effective outer boundary r B to the Solar System (due to stellar mass loss) in the range r B = 10 3 − 10 4 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 mass loss alone (Villaver & Livio 2009). In the future of our own system, Earth is likely to be engulfed by the Sun (Schröder & 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 experience significant mass loss (Bear & Soker 2011, Spiegel & Madhusudhan 2012. For planets orbiting stars that are losing mass, a more general treatment of the dynamics has been carried out for both single planet systems (Veras et al. 2011) and multiple planet systems (Veras & Tout 2012;Voyatzis et al. 2013); these studies provide a comprehensive analysis for the particular case of constant mass loss rates. In addition to causing planets to become unbound, stellar mass loss can drive orbital evolution that leads to unstable planetary systems surrounding the remnant white dwarfs remaining at the end of stellar evolution (Debes & Sigurdsson 2002). Indeed, observations indicate that white dwarfs can anchor both circumstellar disks ) and planetary systems ); many white dwarf atmospheres contain an excess of heavy elements Jura 2003), which is assumed to be a signature of accretion of a secondary body (a planet or asteroid). Finally, mass loss in binary star systems can lead to orbital instability, allowing planets to change their host star ; for additional related work, see also , Moeckel & Veras 2012). This paper builds upon the results outlined above. Most previous studies have focused on stellar mass 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 obtain general analytic results that apply to a wide class of mass loss functions. The first goal of this paper is to study single planet systems where the stellar mass loss rate varies with time. As the system loses mass, the semimajor axis of the orbit grows, and the planet becomes unbound for critical values of the stellar mass fraction m f = M f /M 0 * . For systems that become unbound, we find the critical mass fraction m f as a function of the mass loss rate and the form of the mass loss function. In other systems with mass loss, the orbit grows but does not become unbound. In these cases, we find the orbital elements at the end of the mass loss epoch, again as a function of the mass loss rate and the form of the mass loss function. For initially circular orbits and slow mass loss (time scales much longer than the initial orbital period), the critical mass fraction and/or the final orbital elements are simple functions of parameters that describe the mass loss rate. For initial orbits with nonzero eccentricity, however, the outcomes depend on the initial orbital phase. In this latter case, the allowed values of the critical mass fraction m f (or the final orbital elements) take a range of values, which we estimate herein.
Next we consider the effects of additional planets on the results described above. If the planets are widely spaced, they evolve much like individual single planet systems. However, if the planets are sufficiently close together so that planet-planet interactions are important, the systems are generally chaotic. As a secondary goal, this paper estimates the Lyapunov exponents for twoplanet systems with stellar mass loss. For the sake of definiteness, we focus on planetary systems containing 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 Lyapunov exponent) is proportional to the mass loss time scale. As a result, by the time the star has lost enough mass for the planets to become unbound, the planets have begun to erase their initial conditions through chaos. For systems that evolve far enough in time, one can use the semi-analytic results derived for single planet systems with initially circular orbits (see above) as a rough estimate of 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 sensitive dependence on initial conditions, one cannot unambiguously predict the value of ξ f required for an planet to become unbound.
For completeness, we note that the problem of planetary orbits with stellar mass loss is analogous to the problem of planetary orbits with time variations in the gravitational constant G. For single planet systems (the pure two-body problem), the gravitational force depends only on the (single) product GM * , so that the two problems are equivalent (e.g., Vinti 1974). However, for the case with time varying gravitational constant, the product GM * could increase with time. Current experimental limits show that possible variations occur on time scales much longer than the current age of the universe (see the review of Uzan 2003), so that these effects only (possibly) become important in the future universe (Adams & Laughlin 1997). We also note that when considering time variations of the constants, one should use only dimensionless quantities, in this case the gravitational fine structure constant α G = Gm 2 P /c (e.g., Duff et al. 2002).
This paper is organized as follows. We first present a general formulation of the problem of planetary orbits with stellar mass loss in Section 2 and then specialize to a class of models where the mass loss has a specific form (that given by equations [13] and [14]). These models include a wide range of behavior for the time dependence of the mass loss rates, including constant mass loss rates, exponential mass loss, and mass loss rates that decrease quickly with time; these results are described in Section 3. In the following Section 4, we consider two planet systems and calculate the Lyapunov time scales for a representative sample of mass loss functions. Next we apply these results to representative astronomical systems in Section 5. The paper concludes, in Section 6, with a summary of our results and a discussion of their implications.

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 system is preserved and the total angular momentum is a constant of motion. This constraint is explicitly satisfied in the analytical solutions that follow. For the numerical solutions, this property is used as a consistency check on the numerical scheme. On the other hand, the total energy of the system is not conserved because the total mass decreases with time (equivalently, the system no longer exhibits time reversal symmetry).
General forms for the equations of motion with variable stellar mass are presented in many previous papers (from Jeans 1924to Veras et al. 2011. Some of the subtleties of the various approaches are outlined in Hadjidemetriou (1963). In this section and the next we specialize to systems with a single planet and focus on the case where the planet mass is much smaller than the stellar mass, M p ≪ M * . The specific angular momentum J can be written in the form where a 0 is the starting semimajor axis and M 0 * is the initial mass of the star. Equation (1) can be taken as the definition of the angular momentum parameter η. For a starting circular orbit η = 1, whereas eccentric orbits have η = 1 − e 2 < 1, where e is the initial orbital eccentricity. The radial equation of motion can be written in the dimensionless form where η is constant and where we have defined a dimensionless (total) mass The dimensionless radial coordinate ξ = r/a 0 and the dimensionless time variable is given in units of Ω −1 = (a 3 0 /GM 0 * ) 1/2 .
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 an analogous one-body problem, the equation of motion describes the orbit of the reduced mass. In the version of the problem with mass loss represented by equation (2), the motion is also that of a reduced mass (e.g., Jeans 1924, MacMillan 1925, Hadjidemetriou 1963. In this treatment, the mass loss 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 * → 0. For finite planetary masses, the scaling between the two-body problem and the equivalent single body problem changes quantities by O(M P /M * ). In applications of interest, however, our choice of mass loss functions (and their uncertainties) provides the greatest degree of approximation -much larger 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 semimajor axis is unity (by definition) and the initial conditions require that the starting radial coordinate ξ 0 lies in the range 1 − e ≤ ξ 0 ≤ 1 + e, where the starting eccentricity e is given by η = 1 − e 2 . Note that choosing the value of ξ 0 is equivalent to specifying the starting phase of the orbit (up to a sign). The initial energy E 0 = −1/2 by definition and the initial (radial) velocityξ 0 is given bẏ The initial velocity can be positive or negative, where the choice of sign completes the specification of the starting phase of the orbit.

Change of Variables
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 through the expression 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 forṁ 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 dimensionless variables, this finding implies that the function should vary over a limited range. We thus change the dependent variable from ξ to f and write the equation of motion in the forṁ 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 and dependent variables.
The leading coefficient in equation (8) represents an important quantity in the problem: Note that the time scale for mass loss is given by u/u and the orbital time scale is given by u 2 f 3/2 (this latter time scale is the inverse of the orbital frequency, and is shorter than the orbital period by a factor of 2π). The ratio λ of these two fundamental time scales is given by The leading coefficient in equation (8) is thus λ 2 , the square of the ratio of the orbital time scale to the time scale for mass loss. For small values of λ 2 , the mass loss time is long compared to the orbit time, and the orbits are expected to be nearly Keplerian; for larger λ 2 , the star loses a significant amount of mass per orbit and a Keplerian description is no longer valid. For the former case, where mass loss is slow compared to the orbit time, we can use the parameter λ to order the terms in our analytic estimates.
In addition to the coefficient λ 2 , given by the ratio of time scales, another important feature of equation (8) is the index β appearing within the square brackets, where The index β encapsulates the time dependence of the mass loss. This paper will focus on model equations 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 system takes the form The energy has a starting value E = −1/2, by definition, and increases as mass loss proceeds. If and 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 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 defined according to a ∝ |E| −1 , is also monotonically increasing.

Mass Loss Functions
Next we want to specialize to the class of mass loss functions where β = constant. The defining equation (10) for the mass loss index can be integrated to obtain the forṁ 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 rate has the formṁ = −γm (2−β) .
In addition to simplifying the equation of motion, this form for the mass loss function is motivated by stellar behavior, as discussed below. The dimensionless parameter γ is defined to be the ratio of the initial orbital time scale to the initial mass loss time scale. Specifically, if we define τ = (M * /Ṁ * ) 0 , then γ is given by where M 0 * is the stellar mass and a 0 is the semimajor axis at t = 0. For typical orbits, a 0 = 1 -100 AU, so that we expect the parameter γ to be small, often in the range 10 −7 < ∼ γ < ∼ 10 −4 .
The mass loss rate of stars is often characterized by the physically motivated forṁ whereṀ 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 same power-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 or negative, depending on how the stellar luminosity and radius vary with mass during the different phases of mass loss (see Hurley et al. 2000 for a detailed discussion). For example, if we consider main-sequence stars, the stellar cores adjust quickly enough that the luminsoity obeys the standard 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/2 ≤ q ≤ 1). For main-sequence stars we thus obtain the scaling lawṁ ∼ −m αm , where the index α m = p + q − 1 is predicted to lie in the range 2.5 ≤ α m ≤ 3; the corresponding mass loss index lies in the range −1 ≤ β ≤ −1/2. Next we consider stars on the first giant branch or the asymptotic giant branch. In this phase of stellar evolution, mass loss occurs from an extended stellar envelope, but the luminosity is produced deep within the stellar core. As the star loses mass, the core and hence the luminosity remains relatively constant, whereas the radius scales approximately as R * ∼ M −1/3 * (Hurley et al. 2000). For this case, one obtains the scaling lawṁ ∼ −m −4/3 , with a mass loss index β = 10/3. In general, foṙ m ∼ −m αm , the mass loss index β = 2 − α m . As these examples show, the mass loss index can take on a wide range of values −1 ≤ β ≤ 4.
To fix ideas, we consider the time dependence for mass loss functions that are often used. For a constant mass loss rate, the most common assumption in the literature, the index β = 2, and the mass evolution function has the form The value β = 2 marks the boundary between models where the mass loss rate accelerates with time (β > 2) and those that decelerate (β < 2). For the case of exponential time dependence of the stellar mass, the index β = 1, and the mass loss function has the form The value β = 1 marks the boundary between models where the system reaches zero stellar mass in a finite time (β > 1), and those for which the mass m → 0 only in the limit u → ∞. For the case with index β = 0, which represents an important test case, the mass evolution function becomes m(t) = 1 1 + γt and u(t) = 1 + γt .
For β = 0, analytic solutions are available (see Section 3.2), which inform approximate treatments for more general values of the index β. Finally, the case where β = −1 plays a defining role (see Section 3.1) and corresponds to the forms The value β = −1 marks the boundary between models where the planet becomes unbound at finite stellar mass (β > −1) and those for which the planet becomes unbound only in the limit m → 0 or u → ∞ (β < −1). In general, for constant β = 1, the time dependence of the mass takes the form The particular case β = 1 results in the decaying exponential law of equation (18).

Equation of Motion with Constant Index β
For constant values of the mass loss index β, the equation of motion reduces to the form where the ratio of time scales λ is given by By writing the equation of motion in the form (22), we immediately see several key features of the solutions: When the parameter λ ≪ 1, the left-hand side of equation (22) is negligible, and the equation of 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 traces through its nearly elliptical path (see below). Nonetheless, this behavior is often seen in numerical studies of planetary systems with stellar mass loss (e.g., see Debes & Sigurdsson 2002). Orbital evolution with λ ≪ 1 is often called the "adiabatic regime". We note that this terminology is misleading, however, because "adiabatic" refers to evolution of a thermodynamic system at constant energy (heat), whereas the systems in question steadily gain energy through stellar mass loss (the gravitational potential becomes less negative).
When the parameter λ ≫ 1, the left-hand side of equation (22) dominates, and the solutions for f (u) take the form of power-laws with negative indices. In this regime, the equation of motion approaches the form so that the function f (u) has power-law solutions with indices p given by the quadratic equation The general form for the solution f (u) in this regime is thus After the solutions enter this power-law regime, the energy can quickly grow and the planet can become unbound. To illustrate this behavior, consider the differential equation (12) for the orbital energy. We first consider the regime where λ ≪ 1 and the function f is nearly constant. For the benchmark case f = 1, the equation can be integrated to obtain As long as f ≈ 1, the energy remains negative and the planet remains bound, except in the limit u → ∞. Now let λ > 1 so that the solutions enter into the power-law regime. If we let the solution have the form f = A/u, for u > u c , the differential equation for the energy can be integrated to obtain where the subscript c denotes the reference point where the solutions enters into the power-law regime. Since Au c ∼ u 2 c and |E c | ∼ u 2 c /2, the energy quickly becomes positive once the power-law regime is reached. This argument indicates that the planet becomes unbound when the time scale ratio λ is of order unity (as seen in previous studies). The subsequent subsections provide further verification of this finding.
In order for the solution to make the transition from f ≈ constant to the power-law solutions that cause the orbits to become unbound, the ratio of time scales λ must grow with time. However, growth requires that β > −1 (see equation [23]). We can understand this requirement as follows: The orbital time scale P varies with u (and hence time) according to P ∼ u 2 f 3/2 . Since f is nearly constant, this relation simplifies to the form P ∼ u 2 . The time scale for mass loss τ is given by τ = u/u, which has the form τ ∼ u 1−β from equation (13). As a result, when β = −1, the orbit time has the same dependence on stellar mass as the mass loss time scale, so that the ratio λ is nearly constant as the star loses mass. For β < −1, the ratio λ of time scales decreases with time, and the system grows "more stable".

Results for Single Planet Systems with Stellar Mass Loss
This section presents the main results of this paper for single planet systems with a central star that loses mass. First, we consider mass loss index β = −1, which marks the critical value such that systems with β > −1 become unbound at finite values of the stellar mass, whereas systems with β < −1 only become unbound in the limit m → 0. Next we consider mass loss index β = 0; in this case, the solutions can be found analytically, and these results guide an approximate analytic treatment of the general case, which is addressed next. We also consider the limiting case where stellar mass loss takes place rapidly.

The Transition Case
Here we consider systems where the mass loss index β = −1, which corresponds to the transition value between cases where the ratio λ of time scales grows with time (β > −1) and those where the ratio decreases with time (β < −1). In this regime, the equation of motion reduces to the form The equation of motion can be simplified further by making the change of (independent) variable so that the equation of motion becomes This version of the equation of motion (first considered by Jeans 1924) contains no explicit dependence on the independent variable w, so that the equation can be integrated to obtain the expression where E is a constant that plays the role of energy. Note that we have chosen the sign such that E > 0 and that E = 1 for initially circular orbits. In order for the function f (w) to have oscillatory solutions, the fourth order polynomial must be positive between two positive values of f . In order for this requirement to be met, the parameters must satisfy the inequalities 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 solutions occur.

Systems with Vanishing Mass Loss Index
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 where E is a constant of integration. The parameter E plays the role of energy for the orbit problem where 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 we have 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ξ = 0 and the energy constant has the value where the eccentricity e = (1 − η) 1/2 . In general, the initial value f 0 = ξ 0 can lie anywhere in the range 1 − e ≤ f 0 ≤ 1 + e, and the energy constant has the general form where the choice of sign is determined by whether the planet is initially moving outward (+) or inward (−). With the energy constant E specified, the turning points for the function f are found to be If we consider the function f to play the role of the radial coordinate, then equation (38) defines analogs of the semimajor axis a * and eccentricity e * , which are given by For a given starting value f 0 , the effective eccentricity is given by the expression Note that the effective eccentricity e * of the function f is larger than the initial eccentricity e of the 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 and written in the form f df If we integrate this equation from one turning point to the other, the change in mass of the system the same for every cycle, i.e., where E is given by equation (37). Following standard procedures, we can find the solution for the orbit shape, which can be written in the form The orbit equation thus takes the usual form, except that the original eccentricity e is replaced with the effective eccentricity e * and the effective "radius" variable (f ) scales with the mass/time variable u = 1/m.
For this mass loss function (with β = 0), we can find a simple relationship between the value of the 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 of the 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 is thus given by Here, f is evaluated when the planet becomes unbound. In this case, however, the value of f is constrained to lie in the range f 1 ≤ f ≤ f 2 , where the turning points are given by equation (38).
For small γ, the orbit oscillates back and forth between the turning points many times before the planet becomes unbound. The final value of f is thus an extremely sensitive function of the starting orbital phase. This extreme sensitivity is not due to chaos, and can be calculated if one knows the exact orbital phase at the start. In practice, however, the final value of f can be anywhere in the These results were obtained through numerical integration of equation (22). Note that the range of allowed values for the time scale ratio λ f is much larger than for the case with β = 0, whereas the 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 points of the orbit appraoch f 1,2 = 1. Using f = 1 and η = 1 in equation (44), we find λ = 1.

Limit of Rapid Mass Loss
If we now consider the case where the mass loss is rapid, so that the equation of motion has solution (26) throughout the evolution, we can fix the constants A and B by applying the initial conditions. Since f = ξ/u and u = 1 at the start of the epoch, f (1) = ξ 1 , where ξ 1 is the starting value of the orbital radius. By definition, the semimajor axis is unity, and the starting orbital eccentricity is given by e 2 = 1 − η. The starting radius thus lies in the range which is equivalent to 1 − e ≤ ξ 1 ≤ 1 + e. The derivative f ′ = df /du is given by In the regime of interest where γ ≫ 1, the second term is small compared to the first. In this limit, f ′ (1) = −ξ 1 , where ξ 1 lies in the range indicated by equation (45). The constants A and B are thus determined to be B = 0 and A = ξ 1 , so that the solution has the simple form The energy of the orbit, by definition, starts at E(u = 1) = E 1 = −1/2, and the energy obeys the differential 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), We can then read off the value of u f , and hence the mass m f , where the energy becomes positive and the planet becomes unbound, i.e., 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 ξ 1 rather than the starting semimajor axis, which is unity (notice also that this condition is equivalent to that given by equations [46][47][48] in Veras et al. 2011). For circular orbits, we must have ξ 1 = 1, so that planets become unbound when the stellar mass decreases by one half (as expected).
For cases where the mass loss is rapid, but the planet remains bound, we can find the orbital properties 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., 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, where ξ 2,1 = 1 ± e and η = 1 − e 2 . After mass loss has taken place, the new (dimensionless) stellar mass is m ∞ , and the first integral of the equation of motion can be written in the forṁ where the final constant term takes into account the change in (dimensionless) energy at the moment of mass loss. The radial position at t = 0 is ξ 0 ; since the planet is initially in a bound elliptical orbit, the radial coordinate must lie in the range 1 − e ≤ ξ 0 ≤ 1 + e.
The energy E f of the new orbit is thus given by The energy E f is negative, and the orbit is bound, provided that the remaining mass m ∞ > 1−ξ 0 /2. This condition is thus consistent with equation (49), which defines the mass scale at which orbits become unbound in the limit of rapid mass loss. In terms of the energy E f , the turning points of the new orbit take the form We can then read off the orbital elements for the new (post-mass-loss) orbit, i.e., where the new orbital energy E f is given by equation (53).

Systems with General Mass Loss Indices
In order to address the general case, we first change variables according to the ansatz After substitution, the equation of motion becomes where the subscripts denote derivatives with respect to the new variable x. The ratio of time scales is now given by We can integrate the differential equation (57) to obtain the implicit form where E > 0 has the same meaning as before. To move forward, we define the integral quantity so that Note that J = O(λ 2 ), which means that J will be negligible for most of the mass loss epoch (see Appendix A). The energy of the system can be written in the form At the start of the evolution (where u = 1 and x = 1), the energy E = −1/2 by definition. Using this specification, we can find the value of the integration constant E, which takes the form where f 0 (= ξ 0 ) 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 Since J = O(λ 2 ), the condition |J| ≪ E holds for most times. As a result, working to leading order, we can set J = 0 in equation (64) and recover analogs to the orbital solutions found earlier in Section 3.2 (for a related result, see Radzievskii & Gel'Fgat 1957; for an alternate approach, see Rahoma et al. 2009). The only difference is that the dependent (time-like) variable u is replaced with 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 horizontal lines) 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 and grows with time. The simple oscillatory behavior for f (u) ceases near the point where the time scale ratio λ becomes of order unity. Note that the function f falls outside the boundary marked by the turning points near u = 12 in the Figure. Next we consider the time evolution of the energy of the orbit. After some rearrangement, the energy from equation (62) can be rewritten in the form Since J = O(λ 2 ), it is often convenient to work in the limit J → 0 where the energy becomes This form for the energy shows why the product am of the effective semimajor axis and the mass is slowly varying: In dimensionless units, the energy E = −m/2a so that The product am is thus nearly constant as long as the time scale ratio λ is small, and the departure is of order λ. When λ is small, the orbit cycles through many turning points before the mass changes substantially, so that the average of the above equation becomes so that the average of am is constant to second order.

Number of Cycles:
If we ignore J for now and integrate equation (59) over one cycle, we obtain The integral on the LHS gives us π/E. If we integrate over N cycles we obtain the expression where we assume that m = 1 at the start. The total number of possible cycles occurs when m N → 0, so that Since E ∼ 1 and π(β + 1) ∼ 10, whereas γ ≪ 1, we expect the number of cycles N T ∼ 1/(10γ) to often be large.
The Last Cycle: The above analysis (if we continue to work in the regime J ≪ 1) suggests that the last cycle occurs when the right hand side of equation (69) is no longer large enough to balance the left hand side, which is the same for each cycle. As a result, a minimum mass must be left in the star in order for the orbit to complete a cycle (in the function f ). This condition can be written in the form m c = πγα which holds for α = 1 + β > 0. At the point when the mass falls below this threshold, the time scale ratio is given by λ = (Ef ) 3/2 /(πα). Since Ef ∼ 1 and πα ∼ 3 − 10, the system crosses the threshold 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 scale ratio λ f (evaluated at the moment that planet becomes unbound), we find the condition which can then be written in the form 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 result is written in terms of the variable x instead of u.
The orbital eccentricity, calculated the usual way, oscillates with time with an increasing amplitude of oscillation (e.g., see Figure 3.4). As shown here, however, the function f executes nearly Keplerian behavior, with nearly constant turning points, where this statement is exact in the limit J → 0. The oscillation of eccentricity, although technically correct, is misleading. The turning points of the orbit (in the original radial variable ξ) are strictly increasing functions of time. The oscillation in e arises because the orbits are not ellipses, and, in part, because the period of the orbits in ξ are not the same as the period of the orbits in f . As a result, the oscillations in the calculated eccentricity do not imply that the near-elliptical shape of the orbit is varying between states of greater and lesser elongation. Instead, these oscillations imply that if the mass loss stops and the orbit once again becomes an ellipse, the value of the final eccentricity of that ellipse 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 the planet remains bound after the epoch of stellar mass loss. In this case, we want to estimate the orbital elements of the planet. We are mostly interested in specifying the orbital elements at the end of the mass loss epoch, but we can also evaluate them at any time while the star continues to lose mass. Suppose that the orbit passes through N turning points of the function f during the mass loss epoch. The orbit will then complete a partial cycle so that the final value of f lies between the turning points, f 1 ≤ f ≤ f 2 . In the ideal case, where we have complete information describing both the starting orbital elements and the final value of the stellar mass, and where the mass loss function 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 ≫ 1 (see equation [71]), and it is unlikely that stellar mass loss can be exactly described by a model with constant β 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. Additional planets, or other perturbations, increase this uncertainty (see Section 4). In this case of incomplete information, we can write the mean value of the energy (averaged over the cycles) in the form where we have replaced f with the value 1/E of its effective semimajor axis, and where the remaining term averages to zero. This estimate for the final energy has an uncertainty -a range of possible values -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 where e * is defined by equations (39) and (40). With the energy specified, the final value a f of the semimajor axis is given by The expected value of the energy E f is given by equation (75), but it can take on any value in the range defined by equation (76). Similarly, the final value e f of the orbital eccentricity is given by 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 corresponding range of possible eccentricities.
The expressions derived above for the final orbital elements are expected to be valid provided that the time scale ratio λ is small compared to unity (and hence |J| ≪ 1). The time evolution of the orbital elements is illustrated in Figure 3.4 for a representative system with mass loss index β = 2 and mass loss parameter γ = 10 −4 . Numerical integration of the full equation of motion (solid curves) show that the semimajor axis and eccentricity both oscillate and (on average) grow with time. (Note that the Figure plots log[a].) The values of the elements (a, e) calculated from the average energy (from equation [75]) provide a good approximation to the mean evolution of the orbital elements (see the central dotted curves in Figure 3.4). Furthermore, using the range of allowed energy calculated from equation (76), we can calculate upper and lower limits to the expected behavior of the semimajor axis and eccentricity (shown as the upper and lower dotted curves). 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, the limits on the energy allow for the planet to become unbound, and the upper limit for the semimajor axis approaches infinity. The approximation scheme thus breaks down at this point.

Numerical Results
The equations of motion can be numerically integrated to find the value ξ f of the radial coordinate when the system becomes unbound (when the energy becomes positive). For the case of exponential mass loss, β = 1, the result is shown in Figure 3.5 as a function of the mass loss parameter γ (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.9 curve is shown only for γ > 0.003; for smaller values of γ, the curve oscillates more quickly as a function of γ, and the curve would appear as a solid black band in the figure. The bottom panel shows the value of λ, the ratio of the orbital period to the mass loss time scale, evaluated when the system becomes unbound. As expected, the parameter λ is of order unity when the system energy becomes positive and the planet becomes unbound. For the case of circular orbits at the initial epoch (η = 1), the value of λ ∼ 1.3 for small γ. For starting orbits with nonzero eccentricity, the value 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 β. . The nearly monotonic curve shows the result for the case of circular starting orbits; the oscillating curve shows the result for angular momentum parameter η = 0.9 (which corresponds to starting eccentricity e = √ 0.1 ≈ 0.316 . . .). Bottom panel shows the value of the parameter λ when the planet becomes unbound, for both circular starting orbits (smooth curve) and eccentric orbits (η = 0.9; oscillating curve). All orbits are started at periastron. circular orbits) as a function of γ: where c 0 is a constant. For c 0 = 0.74212, the error is less than 0.4% for mass loss indices in the range 0 ≤ β ≤ 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 γ < ∼ 0.1. For larger values of γ > ∼ 0.1, the final value u f → 2. Equation (79) describes the final mass values (for initially circular orbits) for realistic mass 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). Note that mass loss on the AGB takes place through a series of pulses, but this complication primarily affects 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 parameter γ. 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 becomes unbound is expected: In physical terms, this result means that the mass loss time scale has become shorter than the orbital period, so that the potential well provided by the star is changing fast enough 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 mass loss index β. Here the time scale ratios are evaluated at the moment when the planet becomes unbound. Results are shown for the limiting case of small γ (from Figure 3.5 we see that the time scale ratio λ approaches a constant value as γ → 0). All of the values are of order unity; for the particular case where β = 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 fitting function of the form log λ f = c 1 β + c 2 β 2 , where c 1 and c 2 are constants, we obtain a good fit with the values c 1 = 0.21658 and c 2 = 0.04102. The fitting function is shown as the dashed curve in Figure 3.5, and is almost indistinguishable from the numerically integrated solid curve.

Lyapunov Exponents for Two Planet Systems with Mass Loss
The discussion thus far has focused on single planet systems, whereas many solar systems contain multiple planets. In order to see how multiple planets affect orbital evolution during mass loss, we generalize the treatment to study systems consisting of two planets and a central star with 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) are considered. As a starting point, we fix the planetary masses and initial orbital elements (a, e) to  those of Jupiter and Saturn, and set the initial stellar mass to M 0 * = 1.0 M ⊙ . We also restrict the orbits 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 is generalized later.
As long as the system suffers no close encounters, the orbit of an individual planet is similar to that described by the variable mass two-body problem. An example of the evolution of the osculating orbital elements (a, e) for our benchmark system (see above) is shown in Figure 4 for a mass loss time scale of 10 5 yr. As each planet orbits in an outward spiral, the semimajor axis increases approximately exponentially in time (in inverse proportion to the stellar mass). The eccentricity oscillates rapidly on orbital time scales, and more slowly on secular time scales (as the planets exchange angular momentum), but remains close to its starting value until stellar mass loss 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 amount of mass is lost, the orbital elements a → ∞ and e → 1, and planets can become unbound. Notice that at this point, the stellar mass is only a few percent of the initial value; as a result, this scenario is rather artificial for stars like our Sun, which are only expected to lose about half of their initial masses. However, larger stars lose a greater fraction of their original masses. For example, a star with initial mass M 0 * ≈ 8M ⊙ is expected to end its life as a white dwarf with roughly ∼ 15% of its 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 small enough separation so that orbital crossings can occur. In this regime, chaos dominates and the orbital elements evolve in a less predictable manner. An example is shown in Figure 4 for the same parameters used in Figure 4, but with the initial eccentricity of the inner planet (Jupiter) increased to e = 0.3 (a typical value for the current exoplanet sample). Although the system initially exhibits nearly periodic behavior, by the time t = τ this stable evolution has been compromised.
Studies of dynamical systems generally use the maximum Lyapunov exponent as an indication of the level of chaos present (e.g., Lichtenberg & Lieberman 1983;Strogatz 1994). If a system is chaotic, two nearby trajectories in phase space initially differing by a small amount δ 0 should diverge according to The Lyapunov time is thus τ ly = 1/Λ. The long-term dynamical stability of the solar system has been explored in the absence of stellar evolution (Batygin & Laughlin 2008;Laskar 2008), and current estimates of the Lyapunov time for the solar system (while the mass of the Sun remains constant) are τ ≈ 5 Myr (Sussman & Wisdom 1992), but this value decreases when stellar mass loss is introduced, as demonstrated below.
Here we determine the Lyapunov times as a function of the mass loss time scale via numerical integrations. We define a "real" system along with a "shadow" system where the initial conditions of the shadow system differ by a small amount δ 0 . By integrating both systems simultaneously and  Fig. 9.-Osculating semimajor axis and orbital eccentricity for a pair of planets orbiting an initially solar-mass star with mass loss time scale τ = 10 5 years. Planets have masses, initial semimajor axis and eccentricities of Jupiter and Saturn. The orbital elements evolve in a roughly predictable manner, with the semimajor axes increasing smoothly and the eccentricities oscillating on secular time scales, but remaining relatively constant until the star has lost the majority of its initial mass. monitoring the quantity δ(t), we can calculate the divergence of the neighboring trajectories and then estimate the maximum Lyapunov exponent. Since a three-body system restricted to a plane consists of 12 phase space variables, there is some choice in defining the quantity δ. For the sake of definiteness, we define δ according to 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 complicated behavior, the functions δ(t) will vary for effectively equivalent cases. As a result, for each system of interest, we run 1000 cases with both a "real" and a "shadow" system. To extract the Lyapunov exponent, we can either average together the 1000 runs to construct a single function δ(t) and use the result to find the exponent, or, we can find the exponent from each of the 1000 individual cases and then find the average exponent. Both schemes produce the same values; here we present results for the former case. Figure 4 shows an example of the time evolution of δ(t) for different values of the mass loss time scale τ . After an initial period of transient growth (roughly delimited by t < ∼ 0.3τ ), the divergence metric δ(t) increases exponentially with time and the Lyapunov exponent can be obtained by finding the slope of the line defined by ln δ = ln δ 0 + Λt.
For each value of the mass loss time scale τ , the maximum Lyapunov exponent was calculated as outlined above. Since the maximum possible separation between the reference and shadow systems is 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 initial interval of transient behavior but before saturation occurs. In most cases, Λ was calculated from the time-series data for times τ /3 ≤ t ≤ τ ; this time interval is delimited by the vertical dashed lines 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" Lyapunov time of ∼ 10 Myr (the value obtained without mass loss), and the curves of divergence saturate before t = 0.3τ . In this case, Λ was calculated only for time series data with t < 10 Myr. Notice that 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. To estimate this uncertainty, we have calculated the Lyapunov time values for a wide range of choices for 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 loss time τ . We performed the analysis described above for each of the two planets in the system separately and averaged the results. The horizontal dotted line -included here for referencecorresponds to our numerically determined Lyapunov time for the model solar system in the absence of stellar mass loss. This value is in relatively good agreement with previous calculations for the complete Solar System (Sussman & Wisdom 1992), but differs slightly because our model considers only two of the four giant planets. Note that as τ → ∞, the Lyapunov time approaches the dotted line, i.e., the value expected with no mass loss. The solid curve shows the calculated values of the Lyapunov time, whereas the dashed line indicates the least-squares fit (where the fit was taken over Black curve is τ = 10 4 ; blue is τ = 10 5 ; purple is τ = 10 6 ; and red is τ = 10 7 (yr). After a period of initial growth, the trajectories diverge exponentially, indicated by the linear shape of the latter portions of the graphs. The regions between the dashed lines were used to calculate the Lyapunov exponent. Note that the time variable has been scaled by the mass loss time scale τ .
range of mass loss time scales τ ≤ 10 7 ). Notice that the slope of this line is close to unity. More specifically, we obtained τ ly ∼ τ p where p = 0.99.
Note that this fitted line cannot be meaningfully extrapolated below τ = 10 2 − 10 3 . In this regime, the mass loss time scale τ becomes comparable to the orbital periods of the planets, and the dynamics 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 the difference between nearby trajectories. An especially compelling option is to use the semimajor axis a, because unlike the physical distance between the reference and shadow trajectories, there is no upper limit on this quantity. The previous calculation was thus repeated using Our results are similar, however, which indicates that the Lyapunov times do not depend sensitively on the choice of δ.
Next we would like to ensure that the (nearly) linear relation between the Lyapunov time and mass loss time is not an artifact of the exponential (β = 1) mass loss law that was chosen. Toward this end, we have explored two additional functional forms for the mass loss law: The first used vanishing mass loss index β = 0 (see equation [19]), whereas the second used a constant mass loss rate with β = 2 (see equation [17]). For both of these mass loss functions, we obtained nearly the same 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 in Figure 4, are thus nearly identical, independent of the index β of the mass loss function. This finding 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, we generally expect such solar systems to be only moderately influenced by chaos. In order for chaos to fully erase initial conditions for a dynamical system, nearby trajectories in phase space must diverge for several Lyapunov times. For example, in order for a starting uncertainty of 1 • in position angle to grow into 360 • , one needs ∼ 6 Lyapunov times, or, about 6 mass loss time scales. For exponential mass loss, e.g., this time interval would result in the star losing 99.75% of its initial mass. Since most stars do not lose such a large percentage of their mass, the effects of chaos are not expected to completely erase the initial conditions of these systems. Nonetheless, chaos will partially erase the initial conditions; this trend will affect our ability to predict the phase of the orbit at the end of mass loss and will thus introduce uncertainty into predictions of the final (post-mass-loss) orbital elements (see Section 3.4).

Applications
To illustrate the efficacy of the results found in the previous sections, we consider two representative astronomical problems. For solar type stars, we find an effective outer edge of the solar system, i.e., the boundary between planetary bodies that remain bound after the epoch of mass loss and those that escape (Section 5.1). Next we consider planets that remain bound to white dwarfs, and find their final orbital elements (Section 5.2).

Outer Boundary of the Solar System
Suppose that a solar-type star, with initial mass M 0 * = 1 M ⊙ , loses some portion of its mass over a time interval ∆t = 1 Myr, so that the fraction m f remains afterward. Since solar type stars lose 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 and arbitrary index β > 1, the mass loss parameter γ is given by where a 0 is the initial semimajor axis of the orbit. For initially circular orbits, we can use the results of Section 3.5 to find the conditions for which planets become unbound. Using equation (79) to define 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 become unbound during the mass loss epoch. First we define the constant The critical semimajor axis a c is then given by We can thus find the critical semimajor axis a c for any given value of the index β and the remaining mass fraction m f . Note that the resulting value of a c is a sensitive function of the index β. To leading order, a c ∼ (34,000 AU) m 2(1+β)/3 f . If we use m f = 1/2, the value expected for the Sun, then the critical value of the semimajor axis a c ≈ 8500, 5350, and 3370 AU for indices β = 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. For initially 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 the mean values found here.
Note that in the limit β → ∞, we formally obtain a c → 0 for any value of the remaining mass fraction m f < 1. However, this formulation of the problem breaks down before that limit is reached: For large values of β, even though the time interval is fixed, the mass loss rate accelerates rapidly 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.3 apply, and circular orbits become unbound (remain bound) for mass fraction m f < 1/2 (m f > 1/2).

Orbital Elements for White Dwarf Planets
Now consider a progenitor star with initial mass M 0 * = 5 M ⊙ , which evolves into a white dwarf with mass M wd = 1 M ⊙ (Hurley et al. 2000); the final mass fraction m f = 1/5 (u f = 5). For purposes of illustration, we assume that the mass is lost over a single epoch that can be described by a single value of the mass loss index β, with a time scale ∆t = 1 Myr. For orbits with starting semimajor axes a and eccentricities e, we would like to know the final orbital elements, after the epoch of mass loss. For orbits with starting semimajor axis a in the range 1 − 100 AU (closer planets are often accreted by the star), the value of γ falls in the range γ = 10 −7 − 10 −4 . The mass loss parameter γ is thus small and nearly independent of the index β (see equation [85]). The time scale ratio λ = γu β+1 f 3/2 . For the systems of interest, the largest γ value is thus ∼ 10 −4 , the largest value of u = 5, and the largest value of β = 3; the largest value of the time scale ratio is thus λ ∼ 0.063, so that λ 2 ≤ 0.004. In the approximation scheme developed in Section 3.5, we have exact results when the integral J is small, where J = O(λ 2 ) ≤ 0.004 (see also Appendix A).
The initial conditions for a planetary orbit include not only the semimajor axis and eccentricity (a, e), but also the phase of the orbit at t = 0. 1 This latter quantity is specified by the initial value of the radial coordinate ξ 0 = f 0 , which lies in the range 1 − e ≤ ξ 0 ≤ 1 + e. With ξ 0 (f 0 ) 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 and forth between its turning points (given by equation [38]) while the system loses mass (analogous to 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 Note that we can ignore the second γ 2 term in the above expression. The expected final value of the (dimensionless) semimajor axis (from equation [77]) is given by and the corresponding expected value of the eccentricity (from equation [78]) is given by where we work to the same order of approximation as for a f (recall that e is the starting, pre-massloss value of the eccentricity). Since γ ∼ 10 −7 (a 0 /1AU) 3/2 , the correction term is small: 390,625 γ 2 ∼ 4× 10 −9 (a 0 /1AU) 3 ≤ 0.004 since a 0 ≤ 100 AU. The mean value of final semimajor axis is thus about 5 times the starting value, as expected, and the leading order correction has been quantified. The square of the eccentricity increases by a similar amount. Because of the range of possible orbital phases at the end of mass loss, the orbital elements can differ from these mean values according to the relations where we have used equation (76) and where e 2 * = e 2 + γ 2 (1 − e) 2 is the effective eccentricity of the function f (u) during the epoch of mass loss. Thus, the total (relative) width of the distribution of possible final orbital elements is thus ∼ 2500γe. Note that this range is often larger than the correction to the mean values. Consider a planet with starting semimajor axis a 0 = 100 AU and eccentricity e = 0.30. The mean value of the final semimajor axis is a f ≈ 502 AU, only 2 AU larger than the value suggested by the simple scaling law aM * ≈ constant that is often used. However, the range of possible values about this mean is about ∆a f = ±19 AU. Similarly, the final eccentricity has mean value e f ≈ 0.306, and the width of the range is about ±0.006.

Summary of Results
This paper has reexamined the classic problem of the evolution of planetary orbits in the presence of stellar mass loss. Although this issue has been addressed in previous studies (see Section 1), we generalize existing work to include a new analytic formulation for time-dependent mass loss and to determine Lyapunov time scales for multiple planet systems. In particular, we consider 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 number of 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 time scale for mass loss is significantly longer than the orbital period. By writing the equation of motion 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 a strictly increasing function of time, so that the semimajor axis (defined via a ∼ 1/|E|) increases monotonically (whereas the orbital radius ξ oscillates in and out).
[2] Previous literature often claims that the orbital eccentricity remains constant during the early phases of stellar mass loss. In contrast, this work shows that the eccentricity oscillates back and forth between well defined limits during the phase of mass loss, and that the amplitude of these oscillations grow with time (one example is shown in Figure 3.4). Moreover, the upper and lower limits of the eccentricity range can be calculated analytically using equations (75 -78). Note that these oscillations in the eccentricity, while technically correct, result from assigning orbital elements (which describe ellipses) to orbital paths that are not elliptical. The actual orbit expands with time, and, for example, the outer turning point of the orbit increases monotonically (it does not oscillate).
[3] In the limit of rapid mass loss, λ → ∞, we obtain analytic solutions that describe orbits for the entire class of mass loss functions (see Section 3.3). The condition for the planet becoming unbound is given by equation (49). For planets that remain bound, the new orbital elements are given by equation (55).
[4] Not all mass loss functions lead to planets becoming unbound (except, of course, in the extreme case where the stellar mass vanishes m → 0). The critical value of the mass loss index is β = −1, where systems with mass loss characterized by β < −1 only lose planets in the m → 0 limit. Note that systems with the transition value of the mass loss index β = −1 were first considered by Jeans (1924).
[5] For the particular, intermediate value of the mass loss index β = 0, we can find analytic expressions for the function f (u) and for the final values of the time scale ratio λ f when the planet becomes unbound (see Section 3.2). For arbitrary values of the mass loss index β = 0, this approach can be generalized to find analytic expressions for f (u) and the orbital elements of the planet (see Section 3.4). The resulting expressions are approximate, correct to order O(λ 2 ), and are thus accurate over most of the mass loss epoch.
[6] One way to characterize the dynamics of these systems is through the parameter λ f . We define λ to be the ratio of dimensionless mass loss rate to the orbital frequency, and λ f is the value at 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 starting with nonzero eccentricity, however, the parameter λ f can depart substantially from unity and varies significantly with β (compare Figures 3.2 and 3.2).
[7] For multiple planet systems, we find that the Lyapunov times decrease in the presence of stellar mass loss, so that chaos should play a larger role in planetary dynamics as stars leave the 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 loss time scale of τ ∼ 10 6 yr, the Lyapunov time is ∼ 2 − 3 × 10 5 yr. Three different forms for the stellar mass loss function have been considered, all yielding similar results, which suggests that this trend is robust.

Discussion
In spite of its apparent simplicity, the classic problem of planetary orbits with stellar mass loss is 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 this scheme, the planet is assigned the Keplerian orbital elements that it would have if it were moving within a purely Keplerian potential. In systems with stellar mass loss, however, this approach might be more misleading than illuminating (see also Radzievskii & Gel'Fgat 1957, Hadjidemetriou 1963). 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 oscillating eccentricity would seem to imply motion along an elliptical path, where the shape of the ellipse cycles back and forth between being more elongated and more round. However, this "oscillation" of the eccentricity is an artifact of assigning a Keplerian orbital element (e) to orbital motion that is not Keplerian. During the mass loss epoch, the orbit has inner and outer turning points, analogous to those of an ellipse. But the turning points of the actual orbit do not oscillate -the orbit smoothly spirals outward (e.g., see Figure 2 of Veras et al. 2011).
Another interesting complication arises: These orbits display a type of sensitive dependence on initial conditions even in the absence of chaos. For planets that remain bound after the mass loss epoch, the final orbital elements (a f , e f ) depend on the phase of the orbit, at both the start and the end of the mass loss epoch. However, the number of orbital cycles during mass loss is large, N c ∼ 1/γ ≫ 1, and a precision of one part in N c is necessary to specify the final orbital phase (given the initial phase). The orbit can be calculated to sufficient accuracy to specify the final orbital elements provided that the starting state, the duration of stellar mass loss, and the form of the mass loss function are all known well enough. In practice, however, real astronomical systems will not follow these particular forms to such high accuracy, so that the final orbital elements cannot be predicted 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).

Future Work
There are many opportunities for future work. Analytic studies can be taken into two directions. First, the formulation developed here can be applied to a wide range of astronomical systems, including predictions of orbital elements for planets that remain bound after stellar mass loss 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, multiple phases of mass loss, and higher order approximations for the conditions under which planets become unbound.
For the calculations of the Lyapunov time scales, this paper has focused on analogs of our own solar system, and has considered only the motion of Jupiter, Saturn, and the Sun, since these are the most gravitationally dominant bodies. However, future calculations should also include additional planets and a wider range of starting orbital elements. In addition, this paper has considered relatively short integrations, spanning at most only 10 − 100 Myr. Although stellar mass loss is not expected to continue for longer than 100 Myr, the increased semimajor axes of the planets will change the overall dynamics and the decreased stellar mass (relative to the planets) could lead to an increase in dynamical instabilities. As a result, longer integrations should be performed, using the lower (constant) stellar mass and the increased semimajor axes of the planets as input. The calculations of this paper are limited to coplanar planeteary systems; future work should explore the effects of different inclination angles. Finally, given the diversity exhibited in the observed sample of extrasolar planets, different planetary configurations should also be considered. These types of calculations will help us understand the long term fate of planetary systems in general and will help direct future observations.