EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 533 (September 2011) A&A, 533 (2011) A2 Full HTML
Free Access
Issue
A&A
Volume 533, September 2011
Article Number A2
Number of page(s) 9
Section Celestial mechanics and astrometry
DOI https://doi.org/10.1051/0004-6361/201016199
Published online 11 August 2011

© ESO, 2011

1. Introduction

A classical problem within the realm of orbital stability studies for theoretical and observed planets in stellar binary systems is the circular restricted 3-body problem (CR3BP) (e.g., Szebehely 1967; Roy 2005). The CR3BP describes the motion of a body of negligible1 mass moving in the gravitational field of the two massive primaries considered here to be two stars. The stars move in circular orbits about the center of mass and their motion is not influenced by the third body, the planet. Furthermore, the initial velocity of the planet is assumed in the same direction as the orbital velocity of its host star, which is usually the more massive of the two stars.

The study of planetary orbital stability is timely for various astronomical reasons. First, although most planets are found in wide binaries, various cases of planets in binaries with separation distances of less than 30 AU have also been identified (e.g., Patience et al. 2002; Eggenberger et al. 2004; Eggenberger & Udry 2010, and references therein). Second, many more cases of planets in stellar binary systems are expected to be discovered noting that binary (and possibly multiple) stellar systems occur in high frequency in the local Galactic neighbourhood (Duquennoy & Mayor 1991; Lada 2006; Raghavan et al. 2006). Moreover, orbital stability studies of planets around stars, including binary systems, are highly significant in consideration of ongoing and future planet search missions (e.g., Catanzarite et al. 2006; Cockell et al. 2009).

The CR3BP has been studied in detail by Dvorak (1984), Dvorak (1986), Rabl & Dvorak (1988), Smith & Szebehely (1993), Pilat-Lohinger & Dvorak (2002), and Mardling (2007). It has also been the focus of the previous papers in this series. In Paper I, Eberle et al. (2008) obtained the planetary stability limits using a criterion based on Jacobi’s integral and Jacobi’s constant. The method, also related to the concept of Hill stability (Roy 2005), showed that orbital stability can be guaranteed only if the initial position of the planet lies within a well-defined limit determined by the mass ratio of the stellar components. Additionally, the stability criterion was found to be also related to the borders of the “zero velocity contour” (ZVC) and its topology assessed by using a synodic coordinate system regarding the two stellar components.

In Paper II, Eberle & Cuntz (2010) followed another theoretical concept based on the hodograph eccentricity criterion. This method relies on an approach given by differential geometry that analyzes the curvature of the planetary orbit in the synodic coordinate system. The centerpiece of this method is the evaluation of the effective time-dependent eccentricity of the orbit based on the hodograph in rotating coordinates as well as the calculation of the mean and median values of the eccentricity distribution. This approach has been successful in mapping quasi-periodicity and multi-periodicity for planets in binary systems. It has also been tested by comparing its theoretical predictions with work by Holman & Wiegert (1999) and Musielak et al. (2005) in regard to the extent of the region of orbital stability in binary systems of different mass ratios.

Previously, the work by Eberle & Cuntz (2010) dealt with the assessment of short-time orbital stability, encompassing time scales of 103 yr or less. One of the findings was the identification of a quasi-periodic region in the stellar mass and distance (μ, ρ0) parameter space (see Sect. 2.1). Due to the relatively short time scales considered in the previous study, there is a strong motivation to revisit the onset of orbital instability using longer time scales and different types of methods. The focus of this paper is the analysis of orbital stability by Lyapunov exponents, which are among the most commonly used numerical tools for investigating chaotic behaviour of different dynamical systems (e.g., Hilborn 1994). The exponents have already repeatedly been used in orbital mechanics studies of the Solar System (e.g., Lissauer 1999; Murray & Holman 2001). For example, Lissauer (1999) discussed the long-term stability of the eight Solar System planets, while also considering previous studies. He concluded that the Solar System is most likely astronomically stable, in the view of the limited life time of the Sun; however, the orbits of Pluto and of many asteroids may become unstable soon after the Sun becomes a white dwarf.

The numerical procedure of computing the Lyapunov exponents has been developed by Wolf et al. (1985) based on previous work by Benettin et al. (1980). Hence, the main objective of the present paper is to establish the Lyapunov exponent criterion and verify it by performing the power spectra analysis (e.g., Hilborn 1994) as well as by comparing the obtained results to those presented in Papers I and II. Our newly established criterion is then used to investigate the stability of planetary orbits in stellar binary systems (approximated here as the CR3BP) of different mass and distance ratios. The methodology of our study of orbital stability based on the Lyapunov exponent criterion is presented and discussed in Sect. 2. Detailed model simulations are described in Sect. 3, and our conclusions are given in Sect. 4.

2. Theoretical approach

2.1. Basic equations

In the following, we consider the so-called coplanar CR3BP (e.g., Szebehely 1967; Roy 2005), which we will define 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 × 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.

Following the conventions described by Eberle et al. (2008), we write the equations of motion in terms of the parameters μ and ρ0 with μ and α = 1 − μ being related to the ratio of the two stellar masses m1 and m2 (see below). Moreover, R0 denotes the initial distance of the planet from its host star, the more massive of the two stars with mass m1, whereas D denotes the distance between the two stars, allowing us to define ρ0. In addition, we use a rotating reference frame, which also gives rise to Coriolis and centrifugal forces. The equations of motion are given as where The variables in the above equations describe the position of the planet, which in essence constitutes a test particle. Its position is defined in Cartesian coordinates {x,y,z}. We denote the time derivative or velocity of a coordinate using the dot notation . We also represent the set of second order differential equations, the equations of motion, by a set of first order differential equations. The velocity is defined by the coordinates {u,v,w} whose time derivatives are the accelerations. By defining the mass ratio and using normalized coordinates, we can define the distances {r1,r2} with reference to the location of the stars in the rotating coordinate system.

We enumerate the co-linear Lagrange points in the synodic frame by the order of which the ZVC opens. The point between the stars opens first; therefore, we denote it as L1. The point to the left of the star that does not host the planet opens second; thus, we denote it L2. The point to the right of the star hosting the planet opens third and it is denoted L3. The two Trojan Lagrange points which are of lesser importance to our study are L4 and L5 (see Fig. 1).

2.2. Lyapunov exponents

A fundamental difference between stable and unstable planetary orbits is that two nearby trajectories in phase space will diverge following a power law (usually linear) for the former and exponentially for the latter. The parameter that is used to measure this rate of divergence is called Lyapunov exponent as it was originally introduced by Lyapunov (1907); see also Baker & Gollub (1990). A dynamical system with n degrees of freedom is represented in 2n phase space; thus, to fully determine the stability of the system all 2n Lyapunov exponents must be calculated. The Lyapunov exponents are the most commonly used tools to determine the onset of chaos and chaotic behaviour of both dissipative (e.g., Musielak & Musielak 2009, and references therein) and non-dissipative systems of orbital mechanics (e.g., Lissauer 1999,and references therein). The positive Lyapunov exponents measure the rate of divergence of neighbouring orbits, whereas negative exponents measure the convergence rates between stable manifolds. For dissipative dynamical systems the sum of all Lyapunov exponents is less than 0 (e.g., Musielak & Musielak 2009); however, for Hamiltonian (non-dissipative) systems the sum is equal to 0 (e.g., Hilborn 1994).

Specific applications of the Lyapunov exponents to the circular restricted 3-body problem were discussed by many authors, including Gonczi & Froeschlé (1981), Jefferys & Yi (1983), Lecar et al. (1992), Milani & Nobili (1992), Smith & Szebehely (1993), Murray & Holman (2001), and others. Some of these authors also considered the so-called Lyapunov time, which measures the e-folding time for the divergence of nearby trajectories. It should be noted that the Lyapunov time is not well-defined for cases when close encounters between a planet and one of the stars occur or when the planet is ejected from the system (Lissauer 1999). We shall treat such cases with special caution in this paper.

thumbnail Fig. 1

Locations of the Lagrange points L1, L2, L3, L4, and L5 as used in the CR3BP. Label P denotes the primary star and label S denotes the secondary star.

thumbnail Fig. 2

Case study for the initial planetary distance ratio ρ0 = 0.355 with the planetary orbit in the synodic coordinate system, Lyapunov spectrum, and power spectrum. The stellar mass ratio is μ = 0.3.

thumbnail Fig. 3

Case study for the initial planetary distance ratio ρ0 = 0.474 with the planetary orbit in the synodic coordinate system, Lyapunov spectrum, and power spectrum. The stellar mass ratio is μ = 0.3.

thumbnail Fig. 4

Case study for the initial planetary distance ratio ρ0 = 0.595 with the planetary orbit in the synodic coordinate system, Lyapunov spectrum, and power spectrum. The stellar mass ratio is μ = 0.3.

thumbnail Fig. 5

Case study for the initial planetary distance ratio ρ0 = 0.290 with the planetary orbit in the synodic coordinate system, Lyapunov spectrum, and power spectrum. The stellar mass ratio is μ = 0.5.

thumbnail Fig. 6

Case study for the initial planetary distance ratio ρ0 = 0.370 with the planetary orbit in the synodic coordinate system, Lyapunov spectrum, and power spectrum. The stellar mass ratio is μ = 0.5.

thumbnail Fig. 7

Case study for the initial planetary distance ratio ρ0 = 0.400 with the planetary orbit in the synodic coordinate system, Lyapunov spectrum, and power spectrum. The stellar mass ratio is μ = 0.5.

The previously obtained results clearly show that the Lyapunov exponents can be calculated for the case of the CR3BP considered in this paper, for which we have 2n = 6. This requires a state vector for the system containing 6 elements. Details about the nature of the state vector are given in Appendix A. From the equations of motion, the Jacobian J can be determined, which will be the foundation of how the Lyapunov exponents will be determined. In addition, the nature of this Jacobian will elude to certain properties of the Lyapunov exponents.

For Hamiltonian systems the trace of the Jacobian should equal zero, TrJ = 0. This requires that the CR3BP should have either all zero diagonal elements or an even amount of 3 positive/negative diagonal elements. As a result this forces the sum of the Lyapunov exponents to be also zero. Then we should expect some symmetry in the spectrum of the Lyapunov exponents. Since this is indeed the case, Lyapunov exponent spectra shown in this paper will only present the positive Lyapunov exponents and omit the negative Lyapunov exponents as they do not give any additional information about the system. It is the convention that positive Lyapunov exponents indicate that both dissipative (e.g., Hilborn 1994) and non-dissipative (e.g., Ozorio de Almeida 1990) systems are chaotic. Here the end value of the Lyapunov exponents will be used to make the distinction between chaotic and non-chaotic orbits.

It is well known that the largest Lyapunov exponent is sufficient to determine the magnitude of the chaos in a system (e.g., Wolf et al. 1985). Therefore, this study will consider only the effect of the maximum Lyapunov exponent because it will show the greatest effect of chaos on the system. This motivates us to also examine how large the maximum Lyapunov exponent can be before introducing enough chaos to affect the stability of the CR3BP.

3. Results and discussion

We performed simulations for stellar mass ratios from μ = 0.0 to 0.5 in increments of 0.01. A Yoshida sixth order symplectic integration scheme was used (e.g., Yoshida 1990). As a measure of the precision of the integration scheme, we note that a time step of ϵ = 10-4 yr is used for the individual cases. However, for producing plots of the entire parameter space we adopt a time step of ϵ = 10-3 yr as the smaller time step did not noticeably enhance the quality of the plots. The order of error in energy for each time step was 10-14 and 10-10, respectively. We performed simulations for different time limits, which range from 10 to 105 yr in increments of powers of 10. We also performed case studies with time limits of 106 and 107 yr, although we do not expect significant changes to occur at those longer time scales. By using different time scales we can estimate when certain phenomena occur and ascertain how they will affect other regions over longer periods of time.

We display runs for selected initial conditions (i.e., starting distances ρ0 of the planet from the stellar primary) corresponding to μ = 0.3 in Figs. 2 to 4 and to μ = 0.5 in Figs. 5 to 7. In Figs. 2 to 7, we first present the planetary orbits in the synodic coordinate system (X*, Y*). Secondly in each figure, the Lyapunov spectra of the simulations are shown using a logarithmic scale for the Y*-axis. Lastly in each figure, the time series power spectra of the simulations are shown with their normalized amplitudes. The power spectra were obtained through the usage of a FFT subroutine in Matlab that uses the X*-component separation distance as a function of time and converts the output frequencies to periods. The selected starting distance ratios ρ0 of the planet are 0.355, 0.474, 0.595 for μ = 0.3 and 0.290, 0.370, 0.400 for μ = 0.5. Note that ρ0 indicates the relative initial distance of the planet given as R0/D, with D = 1 AU as the distance between the two stars and R0 as the initial distance of the planet from its host star, the primary star of the binary system.

Using the method of Lyapunov exponents, we are able to verify and extend the methods described by Eberle et al. (2008) and Eberle & Cuntz (2010). Absolute orbital stability can be more rigorously shown through the Lyapunov exponent method. It occurs when , where represents the point at which the initial condition results in a ZVC that opens at L1. For larger values of ρ0, stability is not guaranteed due the behaviour discussed by Eberle et al. (2008); this result is also consistent with our Lyapunov exponent criterion for orbital stability.

In this paper, the main criterion for orbital stability is the Lyapunov exponent criterion, which is based on the maximum Lyapunov exponent. From a theoretical point of view, an orbit is stable when all Lyapunov exponents are exactly zero. Obviously, this “perfect” criterion for orbital stability will be very difficult to reach numerically because it would require an infinite simulation time. Based on our finite simulation times, we obtain 3 positive and 3 negative Lyapunov exponents, and their sum becomes close to zero within the limits of our numerical simulations. Hence, in order to determine orbital stability numerically, we must look at the values of the three positive Lyapunov exponents at the beginning and at the end of our simulations. By comparing these values, we determine whether the exponents decrease in time, and if so what is the rate of their decrease, or whether they stay approximately constant in time.

Using this information, we are able to identify the maximum Lyapunov exponent and adopt it as our primary indicator of orbital stability. We classify an orbit as stable if the initial value of the maximum Lyapunov exponent is below a certain threshold and if it decreases at a “reasonable rate” (see below); otherwise, the orbit is classified as unstable. In our plots of Lyapunov spectra (see the second panels in Figs. 2 to 7), we present all three positive Lyapunov exponents to show their values and study how they change in time.

Table 1

Errors in the Jacobi constant (JCE) for models of μ = 0.3 at different time intervals.

Table 2

Maximum Lyapunov exponent study for the models of μ = 0.3 at different time intervals.

Table 3

Maximum Lyapunov exponent study for the models of μ = 0.5 at different time intervals.

In addition to the Lyapunov exponent criterion for orbital stability, we also use the so-called orbital energy criterion, which requires that the kinetic energy should not exceed twice the value of the potential energy. This is evaluated in our numerical simulations during each time step. A failure of this criterion would imply a break in conservation of the Jacobi constant as detailed in Szebehely (1967) and shown numerically in Table 1. It should also be noted that the cases presented in Table 1 are also presented in Figs. 2 to 4. There are two important points associated with this criterion. First, our systems are Hamiltonian, which means that their energy must be conserved. On the other hand, the fact that the energy conservation of the planet may be broken because we are neglecting the effects of the third mass on the larger masses has already been discussed by Szebehely (1967). Second, even if the energy is only approximately conserved in our numerical simulations, we still request that the Jacobi constant remains constant, which is a required stability condition for CR3BP. To be consistent with this criterion (see also Paper I), we stop our numerical simulation once changes (even small) in the Jacobi constant occur; such cases are depicted in Tables 2 and 3.

Alike in the previous paper by Eberle et al. (2008), we are able to identify the same regions of orbital stability, instability, and domains of quasi-periodic orbital stability. A key difference, however, is that the orbit diagrams depicted in Figs. 2 to 7 have been simulated for 105 yr, whereas the corresponding figures in the previous paper have been simulated for only 103 yr. The power spectra are shown to indicate how we determined the region of quasi-periodic orbital stability including the correct magnitude.

We begin by classifying Figs. 2, 3, 5, and 7 as stable candidate configurations. Three of the cases show similar behaviour in the maximum Lyapunov exponent as shown in Figs. 3b, 5b, and 7b. These Lyapunov spectra have a common trend by starting at a maximum value on the order of 10-1 and converging, albeit slowly, to smaller orders of ten. Figure 3 shows a somewhat different behaviour compared to the other cases in its class of stable candidates. This case establishes a quasi-periodic orbit, which illustrates a resonance in the power spectrum and reveals the same trend of stability in the Lyapunov exponent spectrum. Figure 2 shows an additional variance in behaviour due to the elevated nature of the maximum Lyapunov exponent. This shows that a limited amount of chaos exists in the system while remaining stable for the full simulation time. This case also demonstrates a resonance in the power spectrum.

In contrast, Fig. 4 demonstrates a case of instability. The Lyapunov spectrum shows a different nature than that of the preceding cases. In Fig. 4b it begins at a maximum value greater than 1 and converges to a value between 10-1 and 1. By establishing the preceding cases as stable cases, we can contrast the final values of the corresponding maximum Lyapunov exponents. In the unstable case of Fig. 4b, it is two orders of magnitude greater than in the stable cases of Figs. 3b, 5b, and 7b. However, in comparison with Fig. 2b we find only a difference of one order of magnitude.

thumbnail Fig. 8

Depiction of the maximum Lyapunov exponent (colour coded) for various mass ratios μ and initial conditions ρ0 using a linear scale (left) or a logarithmic scale (right) where  − n corresponds to 10 − n. Crosses denote cases where the simulation was terminated due to the planet being captured by one of the stars or being ejected from the system. Green, red, and black curves represent the initial conditions that cause the ZVC to open at L1, L2, and L3, respectively.

Figure 6 has been examined in a similar manner. Having already classified Figs. 5 and 7 as stable configurations, we emphasize that they reveal similar trends in the maximum Lyapunov exponent as given by Fig. 3. However, Fig. 6 shows a different orbit diagram and is described by a maximum Lyapunov exponent that gives the outcome of instability. This case conveys a much noisier power spectrum along with a maximum Lyapunov exponent on the order of 1 or greater. Finer detail is illustrated in Tables 1 and 2. Particularly Table 2 shows a boundary in the maximum Lyapunov exponent where values near 0.1 and smaller for the first 100 yr tend toward stability.

Figure 8 conveys the bigger picture for the overall parameter space. It represents contour plots of the maximum Lyapunov exponent with respect to the (μ, ρ0) parameter space in linear and logarithmic scale. The crosses depict initial conditions for runs that prematurely ended due to the orbital energy criterion. This demonstrates where the regions of instability occur as reflected by the respective colour coded scale of the maximum Lyapunov exponent. Comparing Fig. 8 (left) to the previous result by Eberle et al. (2008), it can be seen that the main regions of stability remain the same. Other noteworthy features include that the instability islands present near ρ0 = 0.4 have grown as the simulation time has been increased as well as the existence of a plateau near ρ0 = 0.48. The colour scale in Fig. 8 (left) has a maximum colour of dark red at a maximum Lyapunov exponent of 0.15. Some regions appear to be coloured black; however, this does not correspond to the adopted colour scale as it is caused by the finite contour line width due to the close proximity of the lines. Therefore, the plateau has a peak value of a maximum Lyapunov exponent near 0.15. Considering the diagrams in Fig. 3, we conclude that a region of quasi-periodicity exists on this plateau.

Figure 8 (right), the contour plot in logarithmic scale, is shown to display some of the finer details pertaining to structure inside the stability regions. The blue-green coloured regions demonstrate areas of possible stable chaos that are hidden in Fig. 8 (left). Some other smaller contours are also produced in the stability region, but the colour coding in general shows little change, which is partially due to the chosen spacing between contours. The average difference in height between these levels is  − 0.15; on the logarithmic scale, it indicates a change by a factor of 1.41 in the maximum Lyapunov exponent.

Inspecting Fig. 8 also allows a comparison with previous results, particularly results of Paper I. There it was shown that absolute orbital stability occurs when , where represents the point at which due to the initial condition the ZVC opens at L1. This condition is represented by the line of C1 in Fig. 8. It is evident that the Lyapunov exponent criterion is almost perfectly consistent with this condition. Additionally, as discussed in Paper I and references therein, no stability can be expected if because now the ZVC is open at L1, L2, and L3. This finding is also almost perfectly reflected through the behaviour of the maximum Lyapunov exponent as a orbital stability criterion. Both panels of Fig. 8 show no stability regions on the right side of C3, which is reflective of the opening of the ZVC at Lagrange points L1, L2, and L3.

4. Conclusions

We present a detailed case study of the CR3BP with analyses through the usage of the Lyapunov exponents. We are able to characterize stability limits based on the value of the maximum Lyapunov exponent at the end of each simulation. Cases where the maximum Lyapunov exponent exceeds a value of 0.15 indicate that the planet will experience an event that causes the orbital velocity to decrease. These events include intersecting a Lagrange point or the ZVC. After such an event, the planet will experience a series of near misses (or collision) with one of the stars in the binary system leading to overall instability. Chaos theory and the concept of Lyapunov time prevent us from predicting exactly when the planet will be ejected. Using the Lyapunov time as a measure of the length of predictive time, we can show that this relationship is proportional to the inverse of the maximum Lyapunov exponent. Using our critical value of 0.15, the unstable systems will lose predictability within  yr or less. After that time, it is unknown when the planet will be ejected as this could take many multiples of the Lyapunov time to occur.

The method also shows evidence of a region of quasi-periodicity. This is a region where the maximum Lyapunov exponent remains near the critical value without exceeding that value. This is shown for cases simulated for 105 yr. Based on this time scale, we can conclude that our quasi-periodic plateau represents a region of stable or quasi-stable chaos. This is due to the Lyapunov times being much less than the simulation run time. The chaos is shown through the value of the maximum Lyapunov exponent, whereas the stability is shown through the motion of the planet over longer time scales. The lack of near miss events and encounters with regions of decreasing velocity prohibits the planet from becoming orbitally unstable.

Comparisons to previously established criteria for stability show that our results are consistent with previously obtained stability limits (e.g., Holman & Wiegert 1999; Musielak et al. 2005; Cuntz et al. 2007; Eberle et al. 2008; Eberle & Cuntz 2010). In Eberle et al. (2008), Paper I, the onset of instability was related to the topology of the ZVC, whereas in Eberle & Cuntz (2010), Paper II, it was shown that the onset of orbital instability occurs when the median of the effective eccentricity distribution exceeds unity. Both results are consistent with the findings of Paper III, i.e., the inspection of the maximum Lyapunov exponent. However, the latter offers the advantage to link the study of orbital instability for the CR3BP (and any subsequent generalization, if available) to chaos theory, including the evaluation of different types of chaos.

Although our results have been obtained for the special case of the CR3BP, we expect that it may also be possible to augment our findings to planets in generalized stellar binary systems. Desired generalizations should include studies of the elliptical restricted 3-body problem (ER3BP) (e.g., Pilat-Lohinger & Dvorak 2002; Szenkovits & Makó 2008) as well as of planets on inclined orbits (Dvorak et al. 2007), noting that especially the former have important applications to real existing systems in consideration of the identified stellar and planetary orbital parameters (Eggenberger & Udry 2010).

1

Negligible mass means that although the body’s motion is influenced by the gravity of the two massive primaries, its mass is too low to affect the motions of the primaries.

Acknowledgments

This work has been supported by the U.S. Department of Education under GAANN Grant No. P200A090284 (B.Q. and J.E.), the Alexander von Humboldt Foundation (Z.E.M.) and the SETI institute (M.C.).

References

  1. Baker, G. L., & Gollub, J. P. 1990, Chaotic Dynamics: An Introduction (Cambridge: Cambridge University Press) [Google Scholar]
  2. Benettin, G., Galgani, L., Giorgilli, A., & Strelcyn, J.-M. 1980, Meccanica, 15, 9 [NASA ADS] [CrossRef] [Google Scholar]
  3. Catanzarite, J., Shao, M., Tanner, A., Unwin, S., & Yu, J. 2006, PASP, 118, 1319 [NASA ADS] [CrossRef] [Google Scholar]
  4. Cockell, C. S., Léger, A., Fridlund, M., et al. 2009, Astrobiology, 9, 1 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  5. Cuntz, M., Eberle, J., & Musielak, Z. E. 2007, ApJ, 669, L105 [NASA ADS] [CrossRef] [Google Scholar]
  6. Duquennoy, A., & Mayor, M. 1991, A&A, 248, 485 [NASA ADS] [Google Scholar]
  7. Dvorak, R. 1984, Celest. Mech. Dyn. Astron., 34, 369 [Google Scholar]
  8. Dvorak, R. 1986, A&A, 167, 379 [NASA ADS] [Google Scholar]
  9. Dvorak, R., Schwarz, R., Sűli, Á., & Kotoulas, T. 2007, MNRAS, 382, 1324 [NASA ADS] [CrossRef] [Google Scholar]
  10. Eberle, J., & Cuntz, M. 2010, A&A, 514, A19 (Paper II) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Eberle, J., Cuntz, M., & Musielak, Z. E. 2008, A&A, 489, 1329 (Paper I) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Eggenberger, A., & Udry, S. 2010, in Planets in Binary Star Systems, ed. N. Haghighipour (New York: Springer), 19 [Google Scholar]
  13. Eggenberger, A., Udry, S., & Mayor, M. 2004, A&A, 417, 353 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  14. Fouchard, M., Lega, E., Froeschlé, Ch., & Froeschlé, C. 2002, Cel. Mech. Dyn. Astron., 83, 205 [Google Scholar]
  15. Froeschlé, C., Lega, E., & Gonczi, R. 1997, Cel. Mech. Dyn. Astron., 67, 41 [CrossRef] [Google Scholar]
  16. Gonczi, R., & Froeschlé, C. 1981, Cel. Mech. Dyn. Astron., 25, 271 [Google Scholar]
  17. Hilborn, R. C. 1994, Chaos and Nonlinear Dynamics (Oxford: Oxford University Press) [Google Scholar]
  18. Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621 [NASA ADS] [CrossRef] [Google Scholar]
  19. Jefferys, W. H., & Yi, Z.-H. 1983, Cel. Mech. Dyn. Astron., 30, 85 [Google Scholar]
  20. Lada, C. J. 2006, ApJ, 640, L63 [NASA ADS] [CrossRef] [Google Scholar]
  21. Lecar, M., Franklin, F., & Murison, M. 1992, AJ, 104, 1230 [NASA ADS] [CrossRef] [Google Scholar]
  22. Lissauer, J. J. 1999, Rev. Mod. Phys., 71, 835 [NASA ADS] [CrossRef] [Google Scholar]
  23. Lyapunov, M. A. 1907, Ann. Fac. Sci., University of Toulouse, 9, 203 [Google Scholar]
  24. Malhotra, R. 1998, in Solar System Formation and Evolution, ed. D. Lazzaro, R. Vieira Martins, S. Ferraz-Mello, & J. Fernandez (San Francisco: ASP), 149, 37 [Google Scholar]
  25. Mardling, R. A. 2007, in Dynamical Evolution of Dense Stellar Systems, ed. E. Vesperini, M. Giersz, & A. Sills, Proc. IAU Symp., 246, 199 [Google Scholar]
  26. Milani, A., & Nobili, A. M. 1992, Nature, 357, 569 [NASA ADS] [CrossRef] [Google Scholar]
  27. Murray, N., & Holman, M. 2001, Nature, 410, 773 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  28. Musielak, Z. E., & Musielak, D. E. 2009, Int. J. Bif. Chaos, 19, 2823 [CrossRef] [Google Scholar]
  29. Musielak, Z. E., Cuntz, M., Marshall, E. A., & Stuit, T. D. 2005, A&A, 434, 355; Erratum 2008, 480, 573 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Ozorio de Almeida, A. M. 1990, Hamiltonian Systems: Chaos and Quantization (Cambridge: Cambridge University Press) [Google Scholar]
  31. Patience, J., White, R. J., Ghez, A. M., et al. 2002, ApJ, 581, 654 [NASA ADS] [CrossRef] [Google Scholar]
  32. Pilat-Lohinger, E., & Dvorak, R. 2002, Celest. Mech. Dyn. Astron., 82, 143 [Google Scholar]
  33. Rabl, G., & Dvorak, R. 1988, A&A, 191, 385 [NASA ADS] [Google Scholar]
  34. Raghavan, D., Henry, T. J., Mason, B. D., et al. 2006, ApJ, 646, 523 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
  35. Roy, A. E. 2005, Orbital Motion (Bristol and Philadelphia: Institute of Physics Publ.) [Google Scholar]
  36. Smith, R. H., & Szebehely, V. 1993, Celest. Mech. Dyn. Astron., 56, 409 [NASA ADS] [CrossRef] [Google Scholar]
  37. Szebehely, V. 1967, Theory of Orbits (New York and London: Academic Press) [Google Scholar]
  38. Szenkovits, F., & Makó, Z. 2008, Celest. Mech. Dyn. Astron., 101, 273 [NASA ADS] [CrossRef] [Google Scholar]
  39. Tsonis, A. A. 1992, Chaos. From Theory to Applications (New York: Plenum Press) [Google Scholar]
  40. Wolf, A., Swift, J. B., Swinney, H. L., & Vastano, J. A. 1985, Physica D, 16, 285 [Google Scholar]
  41. Yoshida, H. 1990, Phys. Lett. A, 150, 262 [Google Scholar]

Appendix

A.1. Basic concepts and definitions

From stability analysis we can use a Jacobian matrix to show how a state vector x evolves in time (see Tsonis 1992, for details). The governing equations are

This differential equation has a general solution, which is . Assuming J has distinct eigenvalues, we can find a matrix U to diagonalize J to form a diagonal matrix D such as

By rewriting the above equation we obtain the more useful form given as

Using the multiplication theorem, we can develop relations to obtain the eigenvalues of J yielding

Furthermore, we find

This can be generalized from J to a function such as

From our general solution, we can assume . Therefore, we find

A volume of perturbations in phase space will be conserved if |dete(TrJ)t| = 1 or TrJ = 0 for Hamiltonian systems. A dissipative system will have |dete(TrJ)t| < 1 or TrJ < 0. The eigenvalues λ1,λ2,···   ,λn are the Lyapunov exponents of the flow or the logarithms of eigenvalues of J.

A.2. Numerical determination of Lyapunov exponents

From these equations we can define 6 dimensional tangent vectors xi and their derivatives where i = 1,...,6 as

Also we can define a Jacobian matrix, J, given as

A.3. Wolf method

Now we describe the scheme by which Wolf et al. (1985) determine the Lyapunov spectrum and its application to the CR3BP. The Wolf method follows a basic algorithm. The first step is to initialize a state vector of 6 elements. Then, the tangent vectors need to be initialized to some value. We choose to have all tangent vectors to be unit vectors for simplicity. This means that the elements {x1,y2,z3,u4,v5,w6} = 1 and all other elements will be equal to zero.

The next step consists of a loop that will use a standard integrator akin to the Runge-Kutta schemes to determine how the state and tangent vectors will change within a time step. This will continue for an adequate number of steps so that the tangent vectors become oriented along the flow. When this has been accomplished, it is necessary to perform a Gram-Schmidt renormalization (GSR) to orthogonalize the tangent space. Thereafter, we take the logarithm of the length of each tangent vector to obtain the Lyapunov exponents. We then continue the loop.

A.4. Procedure

We use the first tangent vector, x1, to define the basis for the GSR process. Thus, the first step consists in normalizing this vector. With the vectors of the new orthonormal set of tangent vectors, denoted by primes (′), we find From this new set of tangent vectors, the Lyapunov exponents will be calculated considering the lengths of each vector. Therefore, we find

.

All Tables

Table 1

Errors in the Jacobi constant (JCE) for models of μ = 0.3 at different time intervals.

In the text
Table 2

Maximum Lyapunov exponent study for the models of μ = 0.3 at different time intervals.

In the text
Table 3

Maximum Lyapunov exponent study for the models of μ = 0.5 at different time intervals.

In the text

All Figures

thumbnail Fig. 1

Locations of the Lagrange points L1, L2, L3, L4, and L5 as used in the CR3BP. Label P denotes the primary star and label S denotes the secondary star.
In the text
thumbnail Fig. 2

Case study for the initial planetary distance ratio ρ0 = 0.355 with the planetary orbit in the synodic coordinate system, Lyapunov spectrum, and power spectrum. The stellar mass ratio is μ = 0.3.
In the text
thumbnail Fig. 3

Case study for the initial planetary distance ratio ρ0 = 0.474 with the planetary orbit in the synodic coordinate system, Lyapunov spectrum, and power spectrum. The stellar mass ratio is μ = 0.3.
In the text
thumbnail Fig. 4

Case study for the initial planetary distance ratio ρ0 = 0.595 with the planetary orbit in the synodic coordinate system, Lyapunov spectrum, and power spectrum. The stellar mass ratio is μ = 0.3.
In the text
thumbnail Fig. 5

Case study for the initial planetary distance ratio ρ0 = 0.290 with the planetary orbit in the synodic coordinate system, Lyapunov spectrum, and power spectrum. The stellar mass ratio is μ = 0.5.
In the text
thumbnail Fig. 6

Case study for the initial planetary distance ratio ρ0 = 0.370 with the planetary orbit in the synodic coordinate system, Lyapunov spectrum, and power spectrum. The stellar mass ratio is μ = 0.5.
In the text
thumbnail Fig. 7

Case study for the initial planetary distance ratio ρ0 = 0.400 with the planetary orbit in the synodic coordinate system, Lyapunov spectrum, and power spectrum. The stellar mass ratio is μ = 0.5.
In the text
thumbnail Fig. 8

Depiction of the maximum Lyapunov exponent (colour coded) for various mass ratios μ and initial conditions ρ0 using a linear scale (left) or a logarithmic scale (right) where  − n corresponds to 10 − n. Crosses denote cases where the simulation was terminated due to the planet being captured by one of the stars or being ejected from the system. Green, red, and black curves represent the initial conditions that cause the ZVC to open at L1, L2, and L3, respectively.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.