1 - Department of Physics, Science Hall, University of Texas at Arlington, Arlington, TX 76019-0059, USA
2 - Institut für Theoretische Astrophysik, Universität Heidelberg, Albert Überle Str. 2, 69120 Heidelberg, Germany
3 - Kiepenheuer-Institut für Sonnenphysik, Schöneckstr. 6, 79104 Freiburg, Germany

Abstract
Aims. We study the onset of orbital instability for a small object, identified as a planet, that is part of a stellar binary system with properties equivalent to the restricted three body problem.
Methods. Our study is based on both analytical and numerical means and makes use of a rotating (synodic) coordinate system keeping both binary stars at rest. This allows us to define a constant of motion (Jacobi's constant), which is used to describe the permissible region of motion for the planet. We illustrate the transition to instability by depicting sets of time-dependent simulations with star-planet systems of different mass and distance ratios.
Results. Our method utilizes the existence of an absolute stability limit. As the system parameters are varied, the permissible region of motion passes through the three collinear equilibrium points, which significantly changes the type of planetary orbit. Our simulations feature various illustrative examples of instability transitions.
Conclusions. Our study allows us to identify systems of absolute stability, where the stability limit does not depend on the specifics or duration of time-dependent simulations. We also find evidence of a quasi-stability region, superimposed on the region of instability, where the planetary orbits show quasi-periodic behavior. The analytically deduced onset of instability is found to be consistent with the behavior of the depicted time-dependent models, although the manifestation of long-term orbital stability will require more detailed studies.

Key words: stars: binaries: general - celestial mechanics - stars: planetary systems

More than thirty years have passed since Heppenheimer (1974) presented a preliminary theory of planet formation in binary systems, leading to the conclusion that planet formation is virtually impossible in binary systems with separation distances of less than 30 AU. Today, observations have been reported indicating that planets occur in more than twenty binary systems, as well as in several triple star systems (Eggenberger et al. 2004; Patience et al. 2002; Eggenberger & Udry 2007). Although most planets are found in wide binaries, four cases of planets in binaries with separation distances between 20 and 25 AU have also been identified, which are: GJ 86 (Queloz et al. 2000; Lagrange et al. 2006), $\gamma$ Cep (Hatzes et al. 2003; Neuhäuser et al. 2007), HD 41004A (Zucker et al. 2004), and HD 196885A (Correia et al. 2008). The spectral type of the main stellar component ranges between F8V and K2V. These findings are consistent with previous theoretical results showing that planets can successfully form in binary (and possibly multiple) stellar systems (e.g., Kley 2001; Quintana et al. 2002), known to occur in high frequency in the local Galactic neighborhood (Raghavan et al. 2006; Duquennoy & Mayor 1991; Lada 2006).

Bonavita & Desidera (2007) performed a statistical analysis of binaries and multiple systems concerning the frequency of hosting planets, leading to the conclusion that there is no significant statistical difference between binary systems and single stars. That planets in binary systems are now considered to be relatively common is also implied by the recent detection of debris disks in various main-sequence stellar binary systems by the Spitzer Space Telescope (e.g., Trilling et al. 2007, and references therein). Trilling et al. observed 69 main-sequence binary star systems and found emission in excess of predicted photospheric flux levels for approximately 9% and 40% of these systems at 24 and 70 $\mu$ m, respectively, interpreted as being caused by protoplanetary dust.

In the past few decades, significant progress has been made in the study of the stability of planetary orbits in stellar binary systems. Most of these studies focused on S-type systems, where the planet is orbiting one of the stars, and the second star is considered a perturbator. Conversely, P-type orbits lie well outside the binary system, where the planet essentially orbits the center of mass of both stars. Early results were obtained by Hénon & Guyot (1970), Szebehely & McKenzie (1981), Dvorak (1986,1984), Dvorak et al. (1989), Benest (1988,1993,1989,1996), Kubala et al. (1993), and Holman & Wiegert (1999), who explored the stability of planets in binary systems with different mass ratios and eccentricities. More results have been given by David et al. (2003) and Musielak et al. (2005). More recently, Mudryk & Wu (2006) and Barnes & Greenberg (2007) focused on the role of resonances in the ejection of planets, whereas Fatuzzo et al. (2006) presented a detailed statistical analysis of ejection times in response to a large number of configurations, concerning both the planet and the companion star.

In previous work (Stuit 1995; Musielak et al. 2005), we studied the stability of S-type orbits in stellar binary systems, and deduced the orbital stability limits for planets. These limits were found to depend on the stellar mass ratios, a result that will be explored further in this paper. We considered initially circular planetary orbits and classified them by using the following criteria: $0\% \leq$ truly stable $\leq$ 5%, $5\% <$ stable $\leq 15\%$ , $15\% <$ marginally stable $<35\%$ , and unstable $\geq 35\%$ , where the percentage refers to the orbital variability with respect to the initial distance between the primary star and the planet specified at the beginning of the calculations. The stability limit of $5\%$ was motivated by this limit being required for Earth to remain within the conservative habitable zone around the Sun (e.g., Kasting et al. 1993; Underwood et al. 2003). The limit of $35\%$ was based on our studies that showed that test planets outside that limit were not confined to their binary system.

Our paper is structured as follows. In Sect. 2, we present the theoretical approach adopted in our study. Case studies for planets in selected stellar binary systems are given in Sect. 3. Section 4 presents our conclusions.

2 Theoretical approach

2.1 Definition

In the following, we consider the so-called coplanar circular restricted three body problem (e.g., Szebehely 1967; Roy 2005), which is defined as follows. Two stars are in circular motion about their common center of mass and their masses are much larger than the mass of the planet. In our case, the planetary mass is assumed to be $1 \times 10^{-6}$ of the mass of the star it orbits; also note that the planetary motion is constrained to the orbital plane of the two stars. Moreover, it is assumed that the initial velocity of the planet is in the same direction as the orbital velocity of its host star, which is the more massive of the two stars, and that the starting position of the planet is to the right of the primary star (3 o'clock position), along the line joining the binary components (see Fig. 1).

2.2 Basic equations

The origin of the coordinate system is chosen at the center of mass of the two binary stars. Therefore, we find

Next we describe the equations of motion in fixed (sidereal) coordinates loosely following Szebehely (1967), Danby (1988), and Roy (2005). This consists in computating the orbital positions of stars 1 and 2 as a function of time t, assuming that the stars move in circular motion. These positions are denoted as X₁(t), Y₁(t), X₂(t) and Y₂(t), respectively. Furthermore, we consider a small object, a planet, with mass m₃ subjected to the gravity of both stars. The distance of the planet to star i = 1, 2 is given as

$\displaystyle \ddot{X}_3 = -{GM_1\over {r_1}^3}{(X_3-X_1)} -{GM_2\over{r_2}^3}{(X_3-X_2)}$
$\displaystyle \ddot{Y}_3 = -{GM_1\over {r_1}^3}{(Y_3-Y_1)} -{GM_2\over{r_2}^3}{(Y_3-Y_2)}.$			(4)

$\displaystyle \ddot{x}_3 = -{\left({{\alpha\over{r_1}^3}{(x_3-\mu\cos{\tau})} + {\mu\over{r_2}^3}}{(x_3+\alpha\cos{\tau})}\right)}$
$\displaystyle \ddot{y}_3 = -{\left({{\alpha\over{r_1}^3}{(y_3-\mu\sin{\tau})} + {\mu\over{r_2}^3}}{(y_3+\alpha\sin{\tau})}\right)}$			(5)

$\displaystyle r_1(t) = {[(x_3(t)-\mu\cos{\tau})^2+(y_3(t)-\mu\sin{\tau})^2]}^{1\over 2}$
$\displaystyle r_2(t) = {[(x_3(t)+\alpha\cos{\tau})^2+(y_3(t)+\alpha\cos{\tau})^2]}^{1\over 2}.$			(6)

2.3 Jacobi's integral and constant

The next step consists in transforming the equations of motion into a rotating (synodic) coordinate system that rotates along with the binary stars, which simplifies the restricted three body problem. The binary stars remain fixed in such a coordinate system. In addition, fixed equilibrium points (Lagrange points) can be found at which the small mass rotates in sync with the binary stars.

In the synodic coordinate system, with all distances labeled by an asterisk (*), we find

$\displaystyle {\partial\overline{\phi}\over {\partial {x^*_3}}} = \ddot{x}^*_3-... ...{{x^*_3-{\alpha\over{r_1}^3}{(x^*_3-\mu)} - {\mu\over{r_2}^3}}{(x^*_3+\alpha)}}$

$\displaystyle {\partial\overline{\phi}\over {\partial {y^*_3}}} \ = \ \ddot{y}^... ...dot{x}^*_3 = {\left({1-{\alpha\over{r_1}^3} - {\mu\over{r_2}^3}}\right){y^*_3}}$

(7)

$\displaystyle r_1(t) = {[(x^_3(t)-\mu)^2+{y^_3(t)}^2]}^{1\over 2}$
$\displaystyle r_2(t) = {[(x^_3(t)+\alpha)^2+{y^_3(t)}^2]}^{1\over 2}.$			(8)

$\begin{displaymath} \overline{\phi} = {1\over 2} \left[\left({x^*_3}^2+{y^*_3}^2\right)+2\left({\alpha \over r_1}+{\mu\over r_2}\right)\right] , \end{displaymath}$

(9)

$\begin{displaymath} {1\over2}\left(\dot{x}^{*2}_3+\dot{y}^{*2}_3\right) = \overline{\phi}-{\overline{C}\over2} \end{displaymath}$

(10)

For the small mass, the dimensionless distance from the center of mass and the synodic velocity are

$\begin{displaymath} {v^*}^2 = 2\overline{\phi} - \overline{C} = {r}^2+2\left({\alpha \over r_1}+{\mu\over r_2}\right)-\overline{C}. \end{displaymath}$

(12)

The constant $\overline{C}$ can be set by the system parameters $\mu$ and $\alpha$ , noting that $\alpha = 1-\mu$ , and the initial conditions. Jacobi's integral constitutes a relationship between the velocity and position of the planet that remains unchanged during the time-dependent development of the system. It can also be shown that Jacobi's integral at the L₄ and L₅ equilibrium positions can be utilized to specify the Jacobi constant C, given as $C = \overline{C} + \mu(1-\mu)$ . This results in an absolute minimum at L₄ and L₅ for any $\mu$ , which is useful when comparing the characteristic behavior of C for various initial conditions.

The system parameters $\mu$ and $\alpha$ and the initial conditions, allowing v^*₀ to be set, can also be used to formulate an expression for the Jacobi constant C, which is

$\begin{displaymath} C = 2\left[\alpha\left({{{\rho_0}^3+2}\over {2\rho_0}} \righ... ...ft({(1+\rho_0)^3+2}\over {2(1+\rho_0)}\right)\right]-{v^*_0}^2 \end{displaymath}$

(13)

The initial velocity is set by adding the initial velocity of the primary star $R_1\omega$ to the velocity that would result in a circular orbit about a stationary star. This allows the Jacobi constant to be expressed as

$\begin{figure} \par\includegraphics[width=6.5cm,height=6.5cm]{09758F2.eps}\end{figure}$	Figure 2: Locations of the Lagrange points L₁, L₂, L₃, L₄, and L₅ in synodic coordinates $X^{\ast }$ and $Y^{\ast }$ .
Open with DEXTER

2.4 The Jacobi constants at L $_{\sf 1}$ , L $_ {\sf 2}$ , and L $_{\sf 3}$

If we consider the situation when the net force acting on the planet is zero, including fictitious centrifugal forces due to the rotation of the non-inertial synodic reference frame, then the planet will have zero acceleration in that frame, resulting in an equilibrium position. There is a total of five solutions. The five equilibrium points are commonly known as Lagrange points L₁, L₂, L₃, L₄, and L₅ (see Fig. 2). Following Danby (1988), L₁ is located between the primary and secondary stars, and L₂ and L₃ are located to the left of the secondary star and to the right of the primary star, respectively; note that in some textbooks these three Lagrange points are defined differently (e.g., Roy 2005).

The relevance of the points L₁, L₂, and L₃ to the problem considered here is that the values of the Jacobi's constants, C₁, C₂, and C₃, corresponding to these points are useful for establishing the stability of planetary orbits in the system. According to Danby (1988) and Roy (2005), the constants obey the following relation: C₃ < C₂ < C₁. Specifically, if C < C₁, then a planet initially orbiting the primary star may be captured by the secondary star. If C < C₂, then the planet is free to escape from the system through point L₂. Finally, if C < C₃, then the planet may also escape from the system through point L₃.

In the following, we calculate C₁, C₂, and C₃, and use them to determine $\rho_0^{(1)}$ , $\rho_0^{(2)}$ , and $\rho_0^{(3)}$ (see Table 1), which represent critical values of the planet's relative initial distance $\rho _0$ in correspondence to the different types of planetary orbital stability discussed above. These values of $\rho_0^{(i)}$ will only depend on the stellar mass ratio $\mu$ (see Table 1); i.e., the $\rho_0^{(i)}$ increase when $\mu$ decreases. The planet's absolute stability is guaranteed only if $\rho_0 < \rho_0^{(1)}$ ; in this case, the planet moves about the primary star in a periodic orbit. The criteria for planetary orbits to be unstable are either $\rho_0 > \rho_0^{(2)}$ or $\rho_0 > \rho_0^{(3)}$ , or both. An interesting case occurs if $\rho_0 > \rho_0^{(1)}$ and $\rho_0 < \rho_0^{(2)}$ . The existence of periodic orbits in this case depends on the value of C (e.g., Szebehely 1967; Roy 2005; Danby 1988, and references therein).

In this paper, we describe case studies for binary systems of different mass ratios $\mu$ . Our focus is to obtain insight into planetary orbital stability, especially the transition to unstable orbits.

2.5 Tests of computer code

The accuracy of the computer code utilized for pursuing our simulations is checked by integrating the equations of motion for a small body placed near the triangular L₄ point. A mass ratio of $\mu = 0.001$ is assumed considering that the triangular Lagrange points are only stable for $\mu \leq {\mu}_0 = 0.0385 = ({1-{\sqrt{69} / 9}}) / 2$ (Szebehely 1967). If the small mass is placed near the stable point, it will oscillate about that point. Even if in the simulation it could be placed exactly at the stable point, an impossible task owing to the truncation errors inherit in any digital variable, the Runge-Kutta integrator would still introduce an error proportional to the size of the time step $\tau$ of $O({\tau}^5)$ for each step, but $O({\tau}^4)$ over multiple steps. This means that shorter time steps will keep the body more closely near the stable point, but will extend the running time of the simulation.

$\begin{figure} \par\includegraphics[width=12cm,clip]{aa09758F3.eps}\end{figure}$	Figure 3: For a mass ratio of $\mu = 0.5$ , runs for four different distance ratios are shown, which are ${\rho _0} = 0.200$ , 0.300, 0.425, and 0.450 (solid lines). We also show the zero velocity contours (dash-dotted lines) indicating the permitted regions of motion and the Lagrange points (crosses).
Open with DEXTER

$\begin{figure} \par\includegraphics[width=12cm,clip]{aa09758F4.eps}\end{figure}$	Figure 4: Same as Fig. 3, but now for a mass ratio of $\mu = 0.3$ and distance ratios of ${\rho _0} = 0.277$ , 0.400, 0.474, and 0.610.
Open with DEXTER

$\begin{figure} \par\includegraphics[width=12cm,clip]{aa09758F5.eps}\end{figure}$	Figure 5: Same as Fig. 3, but now for a mass ratio of $\mu = 0.1$ and distance ratios of ${\rho _0} = 0.395$ , 0.461, 0.634, and 0.670.
Open with DEXTER

In Table 2, it is shown how different values of $\tau$ affect the precision and running time of the simulation. The approximate scale of the motion about the stable point is given by the quantity $\epsilon$ . Here $T_{\rm sim}$ is the physical time of the simulation, and $T_{\rm run}$ the computational time of the simulation to run on a Pentium 4 CPU Dual 3 GHz with 1 GB RAM. For shorter time steps, the deviation from the equilibrium is less; however, the time it takes to run the simulation increases by the same magnitude.

As depicted in Table 2, $T_{\rm sim}$ is 50 yr corresponding to five orbits of the binary system. This is the same simulation time as used for the bulk of simulations undertaken as part of a broad qualitative scan of the system for different mass ratios and initial conditions. We used a time step of $5 \times 10^{-7}$ yr step^-1 because it allowed a quick survey of the parameter space with an error accumulation of the motion of the small planet given as $\epsilon/D \simeq 1\%$ . However, this may be slightly misleading because the amount of error accumulation is highly dependent on the stability of the respective point. The error accumulation in the simulations will most likely be greater for cases in which instability occurs.

Since this test of accuracy was done only for a known stable point, we also investigated a few simulations on quasi-stable orbits with the same time steps as used in the 50 year simulations, as well as with time steps an order of magnitude smaller for simulation times of up to 1000 years corresponding to 100 binary orbits with no discernable difference. Moreover, the validity of our method is also confirmed by the fact that one of the case studies (i.e., $\mu = 0.1$ , ${\rho_0} = 0.461$ ) has previously been discussed by Szebehely (1967) who found very similar results to those given in this paper.

3 Case studies of planetary motion

In Figs. 3, 4, and 5, we illustrate the transition from stability to instability by progressively increasing the value of $\rho _0$ , i.e., the relative distance of the planet from the primary star, for binary systems with a fixed mass ratio $\mu$ , given as $\mu = 0.5$ , 0.3, and 0.1, respectively; see Cuntz et al. (2007) and Eberle et al. (2008) for results concerning $\mu = 0.2$ . While the initial distance of the planet from its host star is increased, the initial velocity of the planet is set assuming initially circular orbits. Moreover, the changes in the potential and kinetic energies of the planet due to increased values of $\rho _0$ widen the domain for the possible planetary motion, as reflected by the increase in space between the zero velocity curves. If $\rho _0$ is increased beyond $\rho_0^{(1)}$ , $\rho_0^{(2)}$ , and $\rho_0^{(3)}$ , the zero velocity contour opens at L₁, L₂, and L₃, respectively, which dramatically affects the orbital stability of the planet.

Figure 3 focuses on models with $\mu = 0.5$ , where both stars have equal mass. Here we study the cases of ${\rho _0} = 0.200$ , 0.300, 0.425, and 0.450. For ${\rho _0} = 0.200$ , the planet's absolute stability is guaranteed since $\rho_0 < \rho_0^{(1)}$ with $\rho_0^{(1)}$ given as 0.251 (see Table 1). For the case of ${\rho_0} = 0.300$ , an extreme case of orbital instability is encountered, as the planet follows a highly intricate quasi-chaotic orbit; obviously, the planet is orbiting each star in a highly unpredictable manner. A very interesting case occurs for ${\rho_0} = 0.425$ . Even though planetary orbital stability is not guaranteed, the planet nevertheless follows a quasi-periodic orbit, at least for the very short time of our orbital simulations (50 yr). For that particular value of $\rho _0$ , the zero velocity contour is close to opening at both L₂ and L₃, and the planet is now on the brink of being able to escape from the system entirely. For ${\rho_0} = 0.450$ , instability is observed again. However, even though it would be possible for the planet to leave the system due to its energetics, it happens that the planet is captured by the secondary star after 45.9 yr. We also did simulations for $\mu = 0.4$ (not shown here), and found a behavior due to changes in $\rho _0$ that is very similar to $\mu = 0.5$ .

Model simulations for $\mu = 0.3$ are depicted in Fig. 4. Here we study the cases of $\rho _0$ = 0.277, 0.400, 0.474, and 0.610. In accord to our previous discussion, orbital stability is found for ${\rho _0} = 0.277$ (see Table 1). For the other three cases, orbital stability is no longer secured, which means that unstable orbits are expected. However, for ${\rho_0} = 0.400$ and 0.474, the planetary orbit follows quasi-periodic behavior, as also found for $\mu = 0.5$ and ${\rho_0} = 0.425$ (see Fig. 6). For ${\rho_0} = 0.610$ , a complicated planetary orbit emerges. Although the planet is in principle able to leave the binary system entirely, as implied by the topology of the zero velocity contour, it is found that it loops around L₃ between an elapsed time of 35 and 39 yr, and that it is captured by the secondary star after 45.8 yr.

$\begin{figure} \par\includegraphics[width=8cm,height=8cm]{09758F6.eps}\end{figure}$

Figure 6: Limits of stability for planetary orbits for different mass ratios $\mu$ . We show the results based on Jacobi's constant C₁ (solid line) (criterion of stability) and based on Jacobi's constant C₂ (dashed line) and C₃ (dash-dotted line). In addition, we indicate the model parameters $\mu$ and $\rho _0$ of the simulations depicted in Figs. 3 to 5. Here the filled and open boxes refer to case studies of orbital stability and instability, respectively. The gray area indicates the domain of quasi-stability, tentatively determined by our model simulations.

Open with DEXTER

As examples for low mass ratios, we study $\mu = 0.1$ (see Fig. 5) for the cases of ${\rho _0} = 0.395$ , 0.461, 0.634, and 0.670. Note that the difference between $\rho_0^{(1)}$ and $\rho_0^{(2)}$ is very small and, consequently, the initial conditions that result in zero velocity curves passing through L₁ and L₂ are close to each other. As the mass ratio approaches zero, the values of C₁ and C₂ approach each other much more quickly than C₃. This implies that even though the range of initial conditions for which $\rho_0 < \rho_0^{(1)}$ is wide, once $\rho_0^{(1)}$ is passed, the transition to significant instability occurs well prior to $\rho _0$ approaching $\rho_0^{(3)}$ . Our studies find orbital stability for ${\rho _0} = 0.395$ and quasi-periodicity for ${\rho_0} = 0.461$ . For ${\rho_0} = 0.461$ , this result is consistent with previous analytical work by Szebehely (1967), which constitutes an indirect validation of our adopted numerical and analytical methods. For ${\rho_0} = 0.634$ , the planet approaches the zero velocity contour after about 28 yr, and is caputured by the secondary star after 35.1 yr. For ${\rho_0} = 0.670$ , it is found that the planet passes by the secondary star after about 18 yr, and thereafter continues to escape from the system.

Our simulations indicate that, when the mass ratio is small, i.e., $\mu < 0.4$ , the planet appears to remain in a stable, yet peculiar, orbit about the more massive star, well beyond the point at which the zero velocity curve has first opened up. Instability seems to occur when the orbit of the planet comes very close to the zero velocity curve. In this case, it abruptly changes direction and falls toward the more massive star, whereupon it is catapulted in a new direction and the probability of capture increases. However, for a given parameter combination of ( $\mu$ , $\rho _0$ ), it is nonetheless difficult to gauge when this is going to occur because there is no obvious way of determining how close the planet has to come to the zero velocity curve for its trajectory to be deflected.

Based on our studies, highly intricate orbits are encountered for ${\rho_0} = 0.300$ and 0.450 in the case of $\mu = 0.5$ , for ${\rho_0} = 0.610$ in case of $\mu = 0.3$ and, furthermore, for ${\rho_0} = 0.634$ and 0.670 in the case of $\mu$ = 0.1. In these model simulations, the initial distance of the planet from its host star (i.e., the primary star of the binary system) is sufficiently large to be able to initiate orbital instability. This is indicated by the fact that the respective $\rho _0$ value of the model is above the $\rho _0$ value associated with the C₁, C₂, or C₃ limit, or two or more of these limits (see Fig. 6). In some other cases, quasi-periodic orbits are found, even though the planet is expected to follow an unstable orbit, based on the condition previously discussed. This phenomenon will be studied in more detail in Paper II of this series.

4 Conclusion

For the coplanar circular restricted three body problem, we developed stringent mathematical criteria that permit insight into the orbital stability of planets in a stellar binary system. This is accomplished by comparing the planet's relative initial distance $\rho _0$ to the critical values $\rho_0^{(1)}$ , $\rho_0^{(2)}$ , and $\rho_0^{(3)}$ , defined for a fixed stellar mass ratio $\mu$ . An adequate way of demonstrating the different types of behavior is to assess the topology of the zero velocity contour, defined by the $\mu$ and $\rho _0$ values of the system.

The orbital stability of a planet is ensured if its orbit is entirely encapsulated by the zero velocity contour, a result previous pointed out by Roy (2005) and others. The rationale of the zero velocity contour is that it limits the allowable region of the planet as dictated by its limited available kinetic energy. As part of this study, we illustrate the onset of planetary instability by sets of model simulations with different initial planetary starting positions for systems with a given mass ratio. The property of the planet to remain within the zero velocity contour is tantamount to orbital stability, although it does not necessarily imply stability in the sense of quasi-periodicity. Based on the very short timespan of our simulations, we encountered evidence that quasi-periodic orbits may even occur for planets outside the previously established stability region. The relationship between these two stability assessments will be investigated in more detail in Paper II.

Previous numerical case studies (e.g., David et al. 2003; Musielak et al. 2005; Holman & Wiegert 1999) show a behavior consistent with the theoretical prediction of a stringent stability limit, although some minor, but noticeable, differences exist (see discussion by Cuntz et al. 2007), which are typically attributable to minor shortcomings in the analytical fitting procedure owing to the time limitation of the numerical simulations. Note however that for generalized binary systems, other methods are required to establish long-term stability of planetary orbits. Important applications of our work include evaluating numerically deduced stability limits for cases where analyticalresults exist and obtaining insight into the onset of orbital instability, including routes to chaos.

The instability transition for the restricted 3-body problem

I. Theoretical approach

1 Introduction