## Abstract

The orbits of two individual planets in two known binary star systems, γ Cephei and HD 196885 are numerically integrated using various numerical techniques to assess the chaotic or quasi-periodic nature of the dynamical system considered. The Hill stability (HS) function which measures the orbital perturbation of a planet around the primary star due to the secondary star is calculated for each system. The maximum Lyapunov exponent (MLE) time series are generated to measure the divergence/convergence rate of stable manifolds, which are used to differentiate between chaotic and non-chaotic orbits. Then, we calculate dynamical mean exponential growth factor of nearby orbits (MEGNO) maps from solving the variational equations along with the equations of motion. These maps allow us to accurately differentiate between stable and unstable dynamical systems. The results obtained from the analysis of HS, MLE and MEGNO maps are analysed for their dynamical variations and resemblance. The HS test for the planets shows stability and quasi-periodicity for at least 10 Myr. The MLE and the MEGNO maps have also indicated the local quasi-periodicity and global stability in a relatively short integration period. The orbital stability of the systems is tested using each indicator for various values of planet inclinations (ipl ≤ 25°) and binary eccentricities. The reliability of HS criterion is also discussed based on its stability results compared with the MLE and MEGNO maps.

## INTRODUCTION

The discovery of extra solar planets has been growing substantially since the first planet, 51 Pegasi b, was detected almost two decades ago (Mayor & Queloz 1995). Since then 872 extra solar planets have been confirmed as of 2013 April 25.1 Near half of solar-type stars (Duquennoy & Mayor 1991; Raghavan et al. 2006) and a third of all stars in the Galaxy (Raghavan et al. 2010) are in a binary or multistar system with 40 planets confirmed in such systems (Desidera & Barbieri 2007). The confirmation of the existence of planets in binaries has raised a new astrophysical challenge which includes the study of long-term orbital stability of such planets. The ultimate destiny of exoplanet research, including observations from the Kepler mission,2 is to detect a planet on a stable orbit within the habitable zone of the host star (Borucki et al. 1997; Koch et al. 2007; Borucki et al. 2008, 2010). The degree of stability is largely governed by the planet's semimajor axis, eccentricity and orbital inclination. Orbital long-term stability is believed to be a necessary condition for life to develop. The study of orbital stability furthers one self-evident aim of mankind which is to find an answer to the century old question, ‘Are we alone in the Milky Way Galaxy?’

By using different stability criteria the question of stability has been addressed by many others in the past. While studying the Trojan-type orbits around Neptune, Zhou, Dvorak & Sun (2009) showed that the inclination of orbits can be as high as 60° while maintaining orbital stability. Several authors (Szebehely 1980; Szebehely & McKenzie 1981; Szenkovits & Makó 2008) calculated orbital stability of planets by using several techniques that include the integrals of motion, zero velocity surfaces (ZVS), and a Hill stable region that is mapped by a parameter space of orbital radius and mass ratio, μ, for a coplanar circular restricted three body problem (CRTBP). Quarles et al. (2011) has also used the maximum Lyapunov exponent (MLE) to determine the orbital stability or instability for the CRTBP case. The stability limits were defined based on the values of MLE that are dependent on the mass ratio μ of the binaries and the initial distance ratio ρ0 of the planet. Other chaos indicator techniques such as mean exponential growth factor of nearby orbits (MEGNO) maps have also been used to study the dynamical stability of irregular satellites (Hinse et al. 2010) and extrasolar planet dynamics (Goździewski & Maciejewski 2001; Goździewski et al. 2001). The MEGNO criterion is known to be efficient in distinguishing between chaotic and quasi-periodic initial conditions within a dynamical phase space.

Knowing the orbital stability of planets is a crucial step for further studies of planetary systems. In order to probe a planet for its habitability, there exists a primary requirement that the system be orbitally stable. In this paper, we have used three stability indicators, HS, MLE and MEGNO maps, in the study of the orbital dynamics of planets in the selected stellar binaries. We have also compared all the results from three different methods in order to determine if HS can be considered a reliable, efficient tool in the stability analysis of exoplanets.

This paper is outlined as follows. In Section 2, we discuss the basic theory of the analysis tools. In Section 3, we present our numerical methods and the results of the stability analysis from different chaos indicators followed by discussion. Finally, we will conclude in Section 4 with a brief overview of our results.

## THEORY

### Basic definitions and equations

For the motion of a planet of mass m1 around a star of mass m2 in an orbital plane, the initial velocity was calculated using the time derivative of the position matrix given by Murray & Dermott (1999). While calculating the initial conditions we used the values of the parameters (a, e, i Ω, ω and f) whenever they have been observationally determined.

Our particular interest is on the stellar binaries that are less than or equal to 25 au apart. For the stars with greater than 25 au separations, the effects of the secondary star on the planet would not be significant, especially while considering the intent of our present study.

The list of planets in the binaries and their orbital parameters are given in Tables 1 and 2. The initial setup in our simulations is in a barycentric coordinate system with the appropriate placement of Star B (left) and Star A (right) relative to the barycentre. The positive x-axis is taken to be the reference. The true anomaly of both the binary and planet are assumed to be equal to zero.

Table 1.

Orbital parameters of γ Cephei (Neuhauser et al. 2007).

γ Cephei Ab
Mass 1.4 M 0.362 M 1.6 MJ
Semimajor Axis (a19.02 au  1.94 au
Eccentricity (e0.4085  0.115
Argument of Periapsis (ωp0° 180° 94°
γ Cephei Ab
Mass 1.4 M 0.362 M 1.6 MJ
Semimajor Axis (a19.02 au  1.94 au
Eccentricity (e0.4085  0.115
Argument of Periapsis (ωp0° 180° 94°
Table 2.

Orbital parameters of HD 196885 (Chauvin et al. 2007).

HD 196885  Ab
Mass 1.33 M 0.45 M 2.98 MJ
Semimajor Axis (a21 au  2.6 au
Eccentricity (e0.42  0.48
Argument of Periapsis (ωp0° 180° 93$${^{\circ}_{.}}$$
HD 196885  Ab
Mass 1.33 M 0.45 M 2.98 MJ
Semimajor Axis (a21 au  2.6 au
Eccentricity (e0.42  0.48
Argument of Periapsis (ωp0° 180° 93$${^{\circ}_{.}}$$

For two primaries in elliptic orbits moving about their barycentre the dynamical system of a third smaller mass can be written as the first-order differential equations. The equations of motion are given below (Szebehely 1967; Szenkovits & Makó 2008)

(1)
\begin{eqnarray} \begin{array}{lllllll} \displaystyle x^{\prime } &=& \displaystyle u & \quad u^{\prime } &=& \displaystyle 2v + {1\over (1+e\cos {f})}\\ &&&&& \times\;\displaystyle \left[x-{\alpha (x+\mu ) \over r_1^3} - {\mu (x-1+\mu ) \over r_2^3}\right], \nonumber \\ \displaystyle y^{\prime } &=& \displaystyle v & \quad v^{\prime } &=& \displaystyle -2u +{y\over (1+e\cos {f})}\left[1-{\alpha \over r_1^3} - {\mu \over r_2^3}\right], \nonumber \\ \displaystyle z^{\prime } &=& \displaystyle w & \quad w^{\prime } &=& \displaystyle -z + {z\over (1+e\cos {f})}\left[1-{\alpha \over r_1^3} - {\mu \over r_2^3}\right],\end{array} \end{eqnarray}
where
(2)
\begin{eqnarray} \mu &=& \frac{m_2} {(m_1+m_2)},\nonumber \\ \alpha &=& 1-\mu ,\nonumber \\ r_1^2 &= &{(x-\mu )}^2 + y^2 + z^2,\nonumber \\ r_2^2 &=& {(x+\alpha )}^2 + y^2 + z^2. \end{eqnarray}

The Jacobi constant (C0) for an initial state (x0, y0, z0, f0) is given by

(3)
\begin{eqnarray} C_0 = {x_0^2 + y_0^2 + {2{(1-\mu )}\over {r_1}} + {2{\mu }\over {r_2}} \over 1+e\cos (f_0)} -\dot{x}_0^2 - \dot{y}_0^2 - \dot{z}_0^2. \end{eqnarray}

In equation (1), the variables represent the velocity of a test particle (planet) in Cartesian coordinates (x, y, z). The distances r1 and r2 are defined in terms of mass ratio, normalized coordinates and the position of the stars within a pulsating–rotating coordinate system.3

### Lyapunov exponents

For stable planetary orbits, the two nearby trajectories in phase space will converge and for unstable orbits, the trajectories diverge exponentially. The rate of divergence is measured by using the method of Lyapunov exponents (Lyapunov 1907). Wolf et al. (1985) developed a numerical method of computing the Lyapunov exponents in fortran following the earlier works by Benettin et al. (1980).

Lyapunov exponents are commonly used because they give the measure of an attractor of a dynamical system as it converges or diverges in phase space. The positive Lyapunov exponents measure the rate of divergence of neighbouring orbits, whereas negative exponents measure the convergence rates between stable manifolds (Tsonis 1992; Ott 1993). The sum of all Lyapunov exponents is less than zero for dissipative systems (Musielak & Musielak 2009) and zero for non-dissipative (Hamiltonian) systems (Hilborn & Sprott 1994). Lyapunov exponents for the CRTBP have been calculated previously (Gonczi & Froeschle 1981; Murray & Holman 2001) and similar methods have been used in the elliptic restricted three body problem (ERTBP). In this work, we have calculated the Lyapunov exponents for the N-body problem where N = 3. In order to calculate the Lyapunov exponents a dynamical system with n degrees of freedom is represented in a 2n phase space. Then the state vectors (2n) containing six elements are used to calculate the Lyapunov exponents. The details on the calculation of Jacobian J from the equations of motion can be found in Quarles et al. (2011).

For a Hamiltonian system (see above) to be stable, the sum of all the Lyapunov exponents should be zero. In order to numerically meet such a criterion, a simulation of the system would require an impractically long period of time. Within the limits of simulation, the sum of all six exponents must remain numerically close to zero (∼10−10). We have used the largest positive Lyapunov exponent to determine the magnitude of the chaos while ignoring the lesser positive and the negative exponents as they do not provide significant additional information about the evolving system. The positive MLE is known to indicate a chaotic behaviour in both dissipative (Hilborn & Sprott 1994) and non-dissipative (Ozorio de Almeida 1990) systems. Chaos can be proven up to the integration time if a given chaos indicator has converged to an unstable manifold.

### Hill stability

Hill (1878a,b) developed the equations of motion for a particle around the primary mass in presence of a nearby secondary mass. The purpose of the Hill equations was to calculate the orbital perturbation of the particle due to the secondary mass. Later the idea was further developed and used in the study of orbital stability of planets (Szebehely 1967; Walker & Roy 1981; Marchal & Bozis 1982).

The significant radial gravitational influence of the secondary mass reaches as far as the Lagrange points, L1 and L2, forming the Hill sphere (Hill 1878a,b). The contour lines within the sphere are the zero velocity curves. After measuring a particle's position and velocity, a constant of motion relation can be implemented (Szebehely 1967; Murray & Dermott 1999) 2Uv2 = CJ, where v is the velocity, U is the generalized potential, and CJ is the constant of integration called the Jacobi constant. When the velocity of the particle is zero, 2U = CJ, a contour represents a ZVS and the motion of a particle within such a surface is considered Hill stable.

The measure of Hill stability (HS) for the ERTBP, S(f), is given by a parameter-dependent potential, $$\Omega ({\boldsymbol x},f)$$, where f is the true anomaly and Ccr is the value of the Jacobi constant at the Hill radius or the Lagrange point L1 (Szenkovits & Makó 2008)

(4)
\begin{eqnarray} && \!\!\!\!\Omega (x,y,z,f) =\displaystyle{1 \over 2}\left[x^2 + y^2 - ez^2\cos (f)\right] + {1 - \mu \over r_1^2} + {\mu \over r_2^2} \nonumber \\ &&+\displaystyle {1 \over 2} \mu (1-\mu ), \nonumber \\ && \!\!\!\!S(f) =\displaystyle {1 \over C_{cr}}\left[2\Omega ({\boldsymbol x},f)-\left[1+e\cos (f)\right]v^2\right]-1. \end{eqnarray}

Using the orbital parameters obtained from the numerical integration the potential, $$\Omega ({\boldsymbol x},f)$$, is calculated to obtain the HS function S(f). Although the HS function depends on the true anomaly, it can also be represented as a time series. We have implemented this representation in our results concerning the HS function. When the measure of S(f) of a planet is positive then we have the indication of quasi-periodic orbits and its motion is Hill stable. But when the measure of S(f) is negative then the planet enters the instability region.

### The MEGNO chaos indicator

The MEGNO criterion was first introduced by Cincotta & Simó (1999, 2000) and Cincotta, Giordano & Simó (2003) and found widespread applications in dynamical astronomy (Goździewski et al. 2001; Goździewski, Breiter & Borczyk 2008; Hinse et al. 2008, 2010; Goździewski & Migaszewski 2009; Frouard, Vienne & Fouchard 2011; Compère, Lemaître & Delsate 2012; Kostov et al. 2013). The MEGNO (usually denoted as 〈Y〉) formalism has the following mathematical properties. In general, MEGNO has the parametrization 〈Y〉 = α × t + β (see references above). For a quasi-periodic initial condition, we have α ≃ 0.0 and β ≃ 2.0 (or 〈Y〉 → 2.0) for t → ∞ asymptotically. If the orbit is chaotic, then 〈Y〉 → λt/2 for t → ∞. Here λ is the MLE of the orbit. In practice, when generating our MEGNO maps, we terminate a given numerical integration of a chaotic orbit when 〈Y〉 > 12.0. Quasi-periodic orbits have |〈Y〉 − 2.0| ≤ 0.001.

We used the mechanic4 software (Słonina, Goździewski & Migaszewski 2012) optimized to N-body code to calculate the orbits of the given masses and the MEGNO maps on a multi-CPU computing environment. Typically we allocated 60 CPUs for the calculation of one map considering a typical grid of (500 × 300) initial conditions in (a, e) space. The numerical integration of the equations of motion and the associated variational equations (Mikkola & Innanen 1999) are based on the odex integration software (Hairer, Norsett & Wanner 1993) which implements a Gragg–Bulirsch–Stoer algorithm. The MEGNO indicator is calculated from solving two additional differential equations as outlined in Goździewski et al. (2001).

MEGNO and the MLE have a close relation where both indicators provide the magnitude of the exponential divergence of orbits. Froeschlé, Lega & Gonczi (1997) introduced the fast Lyapunov indicator (FLI), which exhibits the least dependence on initial conditions. Mestre, Cincotta & Giordano (2011) showed that MEGNO and FLI are related to each other. FLI is used to detect weak chaos and is considered a faster means to determine the same characteristics as MLE. Recently, Maffione et al. (2011) compared various chaos indicators including FLI and MEGNO. The MEGNO technique and FLI are considered to be in the same class of chaos detection tools (Morbidelli 2002), and we have chosen the MEGNO technique to compare against the HS criterion. More on mathematical properties of MEGNO and its relationship with the Lyapunov exponents can be found in Hinse et al. (2010).

## RESULTS AND DISCUSSION

### Numerical simulation

To establish the HS criterion and calculate the MLE, we numerically simulated each of the planets in the stellar binaries using a Yoshida sixth-order symplectic and a Gragg–Bulirsch–Stoer integration scheme (Yoshida 1990; Hairer et al. 1993; Grazier et al. 1996). A stepping of ϵ = 10−3 yr/step was used in each case to have a better measure of the precision of the integration scheme. The error in energy was calculated at each step which falls in the range of 10−14–10−10 during the total integration period. Numerical simulations were completed for a million years to calculate the MLE and 10 Myr to calculate the HS. MEGNO maps are calculated using 100 000 yr per initial condition. Chaos, quasi-perodicity, or regular motion can be shown up to the integration time; however, the long-term evolution of these systems can only be proven for chaos where the other types of motion are inferred due to the moderate presence or lack of chaos.

The purpose of simulating for 10 Myr was to display the evolution of eccentricity and semimajor axis without applying any indicator tools to provide a consistent comparison in order to analyse the orbital behaviour and to establish the full effectiveness of the indicators. The time series plots for the systems, γ Cephei and HD 196885, are relatively constant with minor oscillations for 10 Myr (Fig. 1). In these cases, we found that the eccentricity of the giant planets is oscillating with a constant amplitude. For example, Figs 1(a) and (c) demonstrate oscillations from 0 to 0.1 and 0.4 to 0.5 in values of eccentricity for γ Cephei and HD 196885, respectively. The amplitude of the oscillation changed with a different choice of initial conditions. As a result, specific choices can minimize the oscillation amplitude and can render the simulation for γ Cephei to be in closer agreement with previous studies by Haghighipour (2006). One such initial condition involved the choice of eccentricity. If e = 0 for the planetary orbit initially, we observed that the amplitude of oscillation is minimum and consistent with Haghighipour (2006) while the amplitude increases with larger initial e values.

Figure 1.

Variation of orbital elements for the giant planets in γ Cephei and HD 196885 (ipl = 0.0) simulated for 1 × 107 yr.

Figure 1.

Variation of orbital elements for the giant planets in γ Cephei and HD 196885 (ipl = 0.0) simulated for 1 × 107 yr.

### MLE: Indicator analysis

The MLE time series for the simulated planets in the stellar binaries are given in Figs 2 and 3. The MLE is plotted using a logarithmic scale along the y-axis and a linear scale along the time-axis. We obtained six Lyapunov exponents from our simulation among which three are negative and three are positive. We inspected the first three positive LEs and found the magnitude of the largest value which is used for our purpose of establishing the stability of a system. In Fig. 2, we have taken the MLE as the primary indicator of the orbital stability. For a given initial condition, the MLE must quickly drop below a cutoff value similar to Quarles et al. (2011) and decrease at a constant rate where we can determine it as stable or unstable.

Figure 2.

Lyapunov exponent time series for the giant planet in γ Cephei with respect to variations in ebin and ipl. The coplanar cases show (a),(b) and (c) considered initial values of ebin = 0.4085, 0.6 and 0.65, respectively. (d) illustrates the time series with the same ebin as in (a) but with different initial values for ipl = 7°, 15° and 25° as indicated by the legend. Note: the time axis of (c) has been truncated to the point of instability.

Figure 2.

Lyapunov exponent time series for the giant planet in γ Cephei with respect to variations in ebin and ipl. The coplanar cases show (a),(b) and (c) considered initial values of ebin = 0.4085, 0.6 and 0.65, respectively. (d) illustrates the time series with the same ebin as in (a) but with different initial values for ipl = 7°, 15° and 25° as indicated by the legend. Note: the time axis of (c) has been truncated to the point of instability.

Figure 3.

Lyapunov exponent time series for the giant planet in HD 196885 with respect to variations in ebin and ipl. The coplanar cases show (a), (b) and (c) considered initial values of ebin = 0.42, 0.46 and 0.6, respectively. (d) illustrates the time series with the same ebin as in (a) but with different initial values for ipl = 7°, 15° and 25° as indicated by the legend. Note: the time axis of (c) has been truncated to illustrate the nature of the instability.

Figure 3.

Lyapunov exponent time series for the giant planet in HD 196885 with respect to variations in ebin and ipl. The coplanar cases show (a), (b) and (c) considered initial values of ebin = 0.42, 0.46 and 0.6, respectively. (d) illustrates the time series with the same ebin as in (a) but with different initial values for ipl = 7°, 15° and 25° as indicated by the legend. Note: the time axis of (c) has been truncated to illustrate the nature of the instability.

The MLE indicates stability for both of the considered systems in the coplanar case, as expected from the observations. This is demonstrated by the values of the MLE which start on the order of 1 and slowly converge down by orders of 10 on a logarithmic scale. The MLE for γ Cephei is calculated using three values of initial binary eccentricity: first at the nominal observational value (Fig. 2a), secondly within the semichaotic region (Fig. 2b) and third within a region of chaos (Fig. 2c). The ebin values assumed for the semichaotic and chaotic regions were obtained from MEGNO maps (Section 3.4). While varying the value of ebin, the orbital inclination was kept constant at 0°. For the first case (ebin = 0.4085), Fig. 2a shows the MLE slowly decreasing to −13 in a million years. The MLE also shows a similar trend for the second case (ebin = 0.6). Hence, considering this decreasing trend and the nature of Lyapunov exponents (Section 2.2), the results reflect the outcome of orbital stability for the planet for our two choices of ebin. For the third choice of ebin = 0.65, the MLE time series began displaying instabilities within 2.25 × 105 yr.

The MLE for γ Cephei Ab (Fig. 2d) is also calculated using three additional values of initial orbital inclination (ipl = 7°, 15°, 25°). While varying the initial planetary inclination, ebin was set to the nominal value. Irrespective of our considered initial values in inclination the MLEs of the planet remained unaffected where all cases demonstrate a decreasing trend to −13 in a million years, reflecting the criterion for orbital stability.

Similar to the case of γ Cephei, we calculated the MLE for the planet in HD 196885 for three values of binary eccentricity. For the nominal observational value of ebin (Fig. 3a) and ebin chosen to be at the semichaotic region (Fig. 3b), the MLEs clearly display stable orbits. But, as the ebin was raised to 0.6 the MLE time series displayed an orbital instability within 1 × 103 yr (Fig. 3c), as also shown by the MEGNO result (see Fig. 9c).

The MLE time series for HD 196885 Ab, Fig. 3 is distinct for all considered initial values of the planet's inclination, ipl, and the orbit of the planet exists in a stable configuration. However, when using our four choices of ipl values, the MLE time series demonstrates that the planetary orbits are more perturbed when it is positioned at 7° and 15°. The MLE values drop quickly to −7 and levels off indicating chaotic behaviour but does not lead to instability. This may be due to the effects of a near resonance behaviour. However, the stability criterion is not violated, thus this system could exhibit long-term chaotic motions akin to Pluto in the Solar system. For ipl = 0° and 25° cases, the MLE values appear to be converging slower up to −10 for the first few thousand years than the other considered cases and indicate stable orbits.

### Establishing the Hill stability criterion

The HS time series for the planets in γ Cephei and HD 196885 are shown in Figs 4– 7. Starting with the coplanar case and the nominal value of binary eccentricity, we have studied the cases when ipl = 7°, 15° and 25°. We have also calculated the HS functions for two other choices of ebin. The choices of our ebin and ipl values are similar to the ones that we made to calculate the MLE time series and MEGNO maps for both of the binary systems.

Figure 4.

The HS function for the giant planet in γ Cephei for different cases of binary eccentricity (for ipl = 0°) simulated for 1 × 107 yr (see Figs 2a, b and c). Note: the time axis of (c) has been truncated to illustrate the nature of the instability (shown by the vertical arrow) and the point of ejection.

Figure 4.

The HS function for the giant planet in γ Cephei for different cases of binary eccentricity (for ipl = 0°) simulated for 1 × 107 yr (see Figs 2a, b and c). Note: the time axis of (c) has been truncated to illustrate the nature of the instability (shown by the vertical arrow) and the point of ejection.

Figure 5.

The HS function for the giant planet in γ Cephei for different cases of initial planetary inclination (for ebin = 0.4085) simulated for 1 × 107 yr (see LE time series Fig. 2d).

Figure 5.

The HS function for the giant planet in γ Cephei for different cases of initial planetary inclination (for ebin = 0.4085) simulated for 1 × 107 yr (see LE time series Fig. 2d).

Figure 6.

The HS function for the giant planet in HD 196885 for different cases of binary eccentricity (for ipl = 0°) simulated for 1 × 107 yr (see Figs 2a, b and c). Note: the time axis of (c) has been truncated to illustrate the nature of the instability and the point of ejection.

Figure 6.

The HS function for the giant planet in HD 196885 for different cases of binary eccentricity (for ipl = 0°) simulated for 1 × 107 yr (see Figs 2a, b and c). Note: the time axis of (c) has been truncated to illustrate the nature of the instability and the point of ejection.

Figure 7.

The HS function for the giant planet in HD 196885 for different cases of initial planetary inclination (for ebin = 0.42) simulated for 1 × 107 yr (see LE times series Fig. 3d).

Figure 7.

The HS function for the giant planet in HD 196885 for different cases of initial planetary inclination (for ebin = 0.42) simulated for 1 × 107 yr (see LE times series Fig. 3d).

The measure of HS for γ Cephei Ab stays positive for ebin = 0.4085 and ebin = 0.6 throughout the integration period (10 Myr) (Figs 4a,b). The oscillations on average are small and positive which reflect the condition for HS criterion established by Szenkovits & Makó (2008). Their calculation of HS is positive and constantly increasing in time for a million years, which is consistent with our result. However, the plots are limited to a million years in the calculations of Szenkovits & Makó (2008) which makes the periodicity for long-term oscillations for the planets to remain unclear. For the nominal case, seven spikes in the HS are noted, which indicates that the planet shows a quasi-periodic motion every 1.4 Myr. When ebin is increased to 0.6 the period of HS function decreases and occurs every 2.5 Myr. This periodicity is caused by the gravitational perturbation on the planet due to the secondary star near the pericentre of the binary. As the ebin is increased, the effects of secondary star on the planet is minimized because the time of interaction has been significantly reduced due to Kepler's second law, thus also reducing the period of HS function.

The HS criterion has shown that the orbital stability of the planet decreases as we increase the ebin of γ Cephei. When the ebin was set to 0.65, the planet starts to display orbital instability at ∼2.25 × 105 yr, a point where the HS function starts a negative trend in the time series (a cutoff point shown by a vertical arrow, see Fig. 4c). The planet was ejected from the system at ∼6 × 105 yr. The results supplement our previous instability prediction made from MLE analysis.

We also found that the measure of HS for γ Cephei Ab stays positive for all of our choices of ipl. The periodicity however decreases with increase in ipl values (Figs 5a, b and c). The increment in orbital inclination of the planet minimizes the effect due to the secondary star's close approach to the planet.

Similarly, the HS time series stays positive for the planet in HD 196885 binary system for the cases when ebin is set at 0.42 (observational value) and 0.46 (Figs 6a,b). The planet displays a quasi periodic motion every 1.1 Myr for our first choice of ebin. But we found insignificant change in the periodicity for the second choice of ebin. The periodicity was not affected because of the minute difference in ebin (0.04).

When the ebin of HD 196885 was set to 0.6, within the presumed unstable region, the HS time series demonstrates that the stability was lost within 5 × 103 yr (Fig. 6c). For this system, the ejection of the planet also occurs at ∼5 × 103 yr of simulation time. The ejection occurs in such a short period that the HS function was not efficient enough to predict the instability.

The stability of HD 196885 Ab was not affected by our choice of orbital inclination. The system demonstrates stable orbits for ipl values as high as 25°. The HS time series remains positive reflecting the established stability criterion from Szenkovits & Makó. Just like the planet in γ Cephei, the period of HS function for HD 196885 Ab decreases with the increase in ipl. It would be interesting to exploit the nature of HS at higher ipl values when the system hits the Kozai resonances, but we have limited our choice of ipl at 25° in this paper.

### Analysis of MEGNO maps

The dynamical maps of MEGNO are generated using a resolution of (500 × 300) producing 150 000 initial conditions in eccentricities and semimajor axes for the respective planets within the selected binaries. In Figs 8 and 9 MEGNO maps for various binary eccentricities were simulated for 1 × 105 yr. The cross-hair in each subplot represents the osculating orbit of planet Ab for the respective binary (see Tables 1 and 2). The colour bar on top of each map indicates the strength in the value of MEGNO(〈Y 〉). The blue colour denotes regions of quasi-periodicity and the yellow indicates regions of chaos.

Figure 8.

MEGNO maps for the planet in γ Cephei at ipl = 0° and various binary eccentricities simulated for 1 × 105 yr.

Figure 8.

MEGNO maps for the planet in γ Cephei at ipl = 0° and various binary eccentricities simulated for 1 × 105 yr.

Figure 9.

MEGNO maps for the planet in HD 196885 (ipl = 0°) and various binary eccentricities simulated for 1 × 105 yr.

Figure 9.

MEGNO maps for the planet in HD 196885 (ipl = 0°) and various binary eccentricities simulated for 1 × 105 yr.

For different eccentricity values of the binary in γ Cephei, Fig. 8, the MEGNO indicator shows a clear distinction between quasi-periodic and chaotic regions. Within the observational value, ebin ∼ 0.4085 (Fig. 8a), the planet unmistakably demonstrates stable orbits. When the binary eccentricity, ebin = 0.2 (Fig. 8b), is decreased the cross-hair is completely inside the quasi-periodic region, hence increasing the orbital stability. Conversely, as the eccentricity of the binary orbit is increased, ebin=0.6, the location of the chaotic mean-motion resonances (yellow spikes at constant semimajor axis) are shifted to lower semimajor axes of the planet (Fig. 8c), hence decreasing the orbital stability. Studies done by Haghighipour (2006) for ebin = 0.20–0.45 at interval of 0.05 and ipl = 0° also show that the planet in γ Cephei demonstrates stable orbits.

We also tested the stability of the planet in γ Cephei for the given values of orbital elements by generating a global map of ebin versus ipl, where ebin ranges from 0.0 to 1.0 and ipl ranges from 0° to 25° (Fig. 10). The dynamics of the planet was found unchanged for the ebin as high as 0.6 and any choice of ipl ≤ 25°. The system seems to show chaotic behaviour at ebin = 0.6 and low ipl values but our earlier investigation (Fig. 8c) indicates the cross-hair at the semichaotic region. We investigated a finer resolution window concerning this region and determined that the initial condition existed on the border of a chaotic region. Nonetheless, this confirms that for our choice of orbital inclination of the planet with the binary orbit, the planet is stable but may display chaos within the long term (>10 Myr). Moreover, our results are consistent with those obtained by Haghighipour (2006) where he showed the stable configuration of the system for all values of the planet's orbital inclination less than 40°.

Figure 10.

MEGNO maps for the planet in (a) γ Cephei showing variation in the (ebin, ipl) plane for the nominal (apl, epl) values and (b) HD 196885 showing variation in the (ebin, ipl) plane for its respective (apl, epl) values simulated for 1 × 105 yr.

Figure 10.

MEGNO maps for the planet in (a) γ Cephei showing variation in the (ebin, ipl) plane for the nominal (apl, epl) values and (b) HD 196885 showing variation in the (ebin, ipl) plane for its respective (apl, epl) values simulated for 1 × 105 yr.

Fig. 9(a) shows a locally stable region for the planet in HD 196885 system. The cross-hair in the map is located right at the edge of the stable region. A small change in semimajor axis of the planet can divert the planet towards the chaotic region loosing global stability. Similarly the location of cross-hair indicates that the quasi-periodicity increases with a decrease in the ebin, Fig. 9(b), while the system becomes chaotic with increment in ebin, Fig. 9(c).

Fig. 10(b), displays the case with the previously observed values of ebin (Fig. 9) in a global map of ebin versus ipl for the given values of orbital parameters of the system. It is found that the dynamics of HD 196885 Ab is unaffected by the choices of ebin ≤ 0.46 and ipl ≤ 25°. Giuppone et al. (2012) showed that the system is mostly stable when the planetary orbit is nearly coplanar or highly inclined orbits near the Lidov–Kozai equilibrium point, (ipl = 47°). We did not test the system for ipl > 25°, however, our result is consistent to the previous work with respect to MLE and HS time series of the respective inclination range.

### Evolution of osculating parameters

Considering the orbital parameters of γ Cephei and HD 196885 from the observations, the evolution of argument of periapsis (ω) and longitude of ascending node (Ω) are studied. In Figs 11(a) and 12(a), we have shown the difference between longitude of periapsis between the planet and the binary (where ϖ = ω + Ω). The longitude of periapsis oscillates and is bound within 0° to 360° which exhibits an apocentric libration phenomenon (Murray & Dermott 1999).

Figure 11.

(a) Difference between the osculating longitude of pericentres (ϖ) for γ Cephei for 5 × 104 yr. (b), (c) and (d) Eevolution of osculating values (e cos ϖ, e sin ϖ) for different orbital inclinations. (e) Evolution of osculating values (i cos Ω, i sin Ω) for different orbital inclination.

Figure 11.

(a) Difference between the osculating longitude of pericentres (ϖ) for γ Cephei for 5 × 104 yr. (b), (c) and (d) Eevolution of osculating values (e cos ϖ, e sin ϖ) for different orbital inclinations. (e) Evolution of osculating values (i cos Ω, i sin Ω) for different orbital inclination.

Figure 12.

(a) Difference between the osculating longitude of pericentres (ϖ) for HD 196885 for 5 × 104 yr. (b), (c) and (d) Evolution of osculating values (e cos ϖ, e sin ϖ) for different orbital inclinations. (e) Evolution of osculating values (i cos Ω, i sin Ω) for different orbital inclination.

Figure 12.

(a) Difference between the osculating longitude of pericentres (ϖ) for HD 196885 for 5 × 104 yr. (b), (c) and (d) Evolution of osculating values (e cos ϖ, e sin ϖ) for different orbital inclinations. (e) Evolution of osculating values (i cos Ω, i sin Ω) for different orbital inclination.

Plots of the osculating values, e cos ϖ and e sin ϖ, for the planet's nominal choice of parameters indicated no change in the osculating parameters of both systems. When the orbital inclination was increased (up to 25°) the precession rates were noticeably increased (see Figs 11b, c, d and 12b, c, d), an effect similar to period-doubling bifurcation (Ott 1993) leading to chaos. This change is not noticed directly in LE time series for γ Cephei (Fig. 2d) but we do observe variation (different converging rate) in the LE time series for HD 196885 (Fig. 3d) and variation (decreasing periodicity) in HS time series for both of the systems (Figs 5 and 7).

Figs 11(e) and 12(e) display plots of the osculating parameters, i cos Ω and i sin Ω, for different cases of orbital inclination. For the coplanar case, we would have a point–circle at the origin and the circle would grow larger for larger values of ipl, as seen in these diagrams. The osculating parameters with respect to inclination indicates regular motion in the case of γ Cephei. But in the case of HD 196885 significant variations occur with increases in inclination (Fig. 12e).

### Reliability of hill stability function

Following the work done by Szenkovits & Makó (2008), we used the HS criterion to study the orbital stability of planets in stellar binaries (γ Cephei and HD 196885). We also used the previously well-known chaos indicators, MLE and the MEGNO maps to study the same systems. The MLE (Figs 2a and 3a) and the cross-hair in the quasi periodic region of MEGNO map (Figs 8a and 9a) demonstrate a stable configuration for given parameters in each of the considered systems. The positive value and non-decreasing global trend of HS time series (Figs 4a and 6a) also provide the necessary evidence of stable system. The HS time series indicates that the systems are quasi-periodic and the periodicity decreases as the planet's inclination and binary eccentricity is increased. Since the calculation for MLE is limited to 1 Myr, we are unable to notice any periodicity in the time series as seen in the HS time series but it does also indicate trends towards stability. Similarly, MEGNO maps do not demonstrate any changes in stability when the respective planet's inclination increased up to 25°, nonetheless, it shows that the planets in both of the system exist within a stable configuration.

The stability of the planets is tested for highly eccentric binary orbits using MLE and MEGNO maps (Figs 2c and 3c and 10). These figures clearly demonstrate unstable orbits for chosen ebin values. Similar tests done by using HS time series (Figs 4c and 6c) produce similar results regarding stability of the system.

The HS time series supplements two of our results from MLE and MEGNO for the cases when orbits are periodic and chaotic, thus establishing itself as a reliable stability criterion. For the planet in γ Cephei (see Figs 2c and 4c), the HS function indicated an instability within the planetary orbit on an equal simulation time-scale compared to MLE (see Section 3.2). However, it took a slightly longer simulation time to predict instability as compared to MLE for the case of HD 196885. For the system where the planet ejects in relatively short time (∼5 × 103 yr in this case) MLE is more efficient to predict regular or chaotic dynamical systems and the HS function seems to be an inefficient indicator (Figs 3c and 6c).

There are limitations to the definition of HS as well. Most notably HS is defined for S-type orbits only. Considering this limitation, it may be best suited in future work to the study of hypothetical or observed moons around gas giants as observed S-type planet populations remain modest (Kipping et al. 2012). However, the populations of multiplanet single star systems are rapidly increasing with the use of Transit Timing Variations (see Lissauer et al. 2011, 2013 for recent results) and the HS criterion would provide adequate stability determinations for similar systems.

## CONCLUSIONS

We have applied various chaos indicator techniques in order to study the dynamics of S-type extrasolar planets in the binaries that are less than 25 au apart. With the time series obtained from the MLE, the HS function, and the maps from the MEGNO indicator, we have shown that both the systems exist within a stable configuration for given parameters. Using these chaos indicators, we have also tested the orbital stability of the system for various choices of binary eccentricity and planet's orbital inclination.

The calculated MLE and HS time series for the planets in both systems for different values of planet's orbital inclination, ipl = 0°, 7°, 15°, 25°, and given ebin value display the orbital stability of the system. For the given values of orbital parameters it is found that the planet in γ Cephei/HD 196885 can maintain stability for ebin as high as 0.6/0.46. Similar studies made by using the global MEGNO maps of ebin versus ipl for given orbital parameters concur with our results from MLE and HS and demonstrate that the planet in γ Cephei/HD 196885 can maintain stability for ipl ≤ 25° and ebin ≤ 0.6/0.46.

The MEGNO chaos indicator has been effective in determining the quasi-periodic regions. The location of eccentricity–semimajor axis cross-hairs in the MEGNO maps for the planets in γ Cephei and HD 196885 systems (Figs 8a and 9a) are located well inside the quasi periodic region (blue) and in the teeth between the chaotic and quasi periodic regions, respectively. This resembles the global and local stability of the planets. These planets do not survive if the binary orbits are highly eccentric (see Figs 10a,b).

The HS time series for a planet is successfully measured using the numerical integration of orbital parameters and conversion into a potential related to the ERTBP. The measure of the HS of the planets in γ Cephei and HD 196885 has shown changes in periodicity with the variation of ipl and ebin values. These periodicities are believed to have originated due to the differences in the time of interaction on the planet due to the secondary star near time of periastron passage for the binary.

A concise study of evolution of osculating parameters shows that they have insignificant variation in both of the systems for the nominal case. But, as we increase the ipl values up to 25° the osculating values (e cos ϖ and e sin ϖ) evolve causing the precession rate in γ Cephei and HD 196885 to increase. In the set of osculating parameters (i cos Ω and i sin Ω) the inclination varies with a very small amplitude while the longitude of ascending node rotates in the full range [0°, 360°]. In case of HD 196885 these parameters display significant evolution which provides possible evidence towards chaotic behaviour in this regime.

Aside from the dynamical analysis of planets in stellar binaries, we are able to successfully test the reliability of HS against the results obtained from both MLE and MEGNO. Direct comparison of stability shows that the HS test can be set as one of the three stringent criteria in the study of stable/unstable nature of a planetary orbit. Our results show that the HS indicator is comparable, in the context of determining the orbital stability, to other well-known indicators. Like MLE and MEGNO, it has consistently predicted the stable/unstable nature of planets in binaries. Calculations show that the HS function predicts instability of orbits in comparable time with LE (case of γ Cephei). Also, unlike MEGNO maps, HS time series cannot produce definite quasi-periodic, chaotic or decreasing stability regions without considering a map of a larger parameter space.

SS and BQ would like to thank department of Physics UT Arlington, Zdzislaw Musielak and Manfred Cuntz for their continuous support and guidance. TCH gratefully acknowledges financial support from the Korea Research Council for Fundamental Science and Technology (KRCF) through the Young Research Scientist Fellowship Programme and financial support from Korea Astronomy and Space Science Institute (KASI) grant number 2013-9-400-00. Numerical computations were partly carried out using the SFI/HEA Irish Centre for High-End Computing (ICHEC) and the PLUTO computing cluster at the Korea Astronomy and Space Science Institute. Astronomical research at the Armagh Observatory is funded by the Northern Ireland Department of Culture, Arts and Leisure (DCAL). We would like to thank Cezary Migaszewski for his useful comments and suggestions which have brought a lot of improvements to this work.

1
www.exoplanet.eu (Schneider et al. 2011)
3
Although the masses were integrated using an N-Body representation, we have included the ERTBP equations of motion to illustrate the necessary transformations to arrive at the pulsating–rotating coordinate system (see Szebehely 1967 for full details) used in our analysis.

## REFERENCES

Benettin
G.
Galgani
L.
Giorgilli
A.
Strelcyn
J.-M.
Meccanica
,
1980
, vol.
15
pg.
9

Borucki
W.
, et al.  .
Sun
Y.-S.
Ferraz-Mello
S.
Zhou
J.-L.
Proc. IAU Symp. 249, Finding Earth-size Planets in the Habitable Zone: the Kepler Mission.
,
2008
Cambridge
Cambridge Univ. Press
pg.
17

Borucki
W. J.
Koch
D.
Basri
G.
Batalha
N.
Brown
T.
Caldwell
D.
Sci
,
2010
, vol.
327
pg.
977

Borucki
W. J.
Koch
D. G.
Dunham
E. W.
Jenkins
J. M.
Soderblom
D.
ASP Conf. Ser. Vol. 119, The Kepler Mission: A Mission To Detennine The Frequency Of Inner Planets Near The Habitable Zone For A Wide Range Of Stars
,
1997
San Francisco
Astron. Soc. Pac.
pg.
153

Chauvin
G.
Lagrange
A.-M.
Udry
S.
Mayor
M.
A&A
,
2007
, vol.
475
pg.
723

Cincotta
P.
Simó
C.
Celest. Mech. Dyn. Astron.
,
1999
, vol.
73
pg.
195

Cincotta
P. M.
Giordano
C. M.
Simó
C.
Physica D Nonlinear Phenom.
,
2003
, vol.
182
pg.
151

Cincotta
P. M.
Simó
C.
A&AS
,
2000
, vol.
147
pg.
205

Compère
A.
Lemaître
A.
Delsate
N.
Celest. Mech. Dyn. Astron.
,
2012
, vol.
112
pg.
75

Desidera
S.
Barbieri
M.
A&A
,
2007
, vol.
462
pg.
345

Duquennoy
A.
Mayor
M.
A&A
,
1991
, vol.
248
pg.
485

Froeschlé
C.
Lega
E.
Gonczi
R.
Celest. Mech. Dyn. Astron.
,
1997
, vol.
67
pg.
41

Frouard
J.
Vienne
A.
Fouchard
M.
A&A
,
2011
, vol.
532
pg.
A44

Giuppone
C. A.
Morais
M. H. M.
Boué
G.
Correia
A. C. M.
A&A
,
2012
, vol.
541
pg.
A151

Gonczi
R.
Froeschle
C.
Celest. Mech.
,
1981
, vol.
25
pg.
271

Goździewski
K.
Maciejewski
A. J.
ApJ
,
2001
, vol.
563
pg.
L81

Goździewski
K.
Migaszewski
C.
MNRAS
,
2009
, vol.
397
pg.
L16

Goździewski
K.
Bois
E.
Maciejewski
A. J.
Kiseleva-Eggleton
L.
A&A
,
2001
, vol.
378
pg.
569

Goździewski
K.
Breiter
S.
Borczyk
W.
MNRAS
,
2008
, vol.
383
pg.
989

Grazier
K. R.
Newman
W. I.
F.
Goldstein
D. J.
Kaula
W. M.
BAAS
,
1996
, vol.
28
pg.
1181

Haghighipour
N.
ApJ
,
2006
, vol.
644
pg.
543

Hairer
E.
Norsett
S. P.
Wanner
G.
Springer Series in Computational Mathematics: Solving Ordinary Differential Equations I, Nonstiff Problems
,
1993
2nd edn
Berlin
Springer-Verlag
Hilborn
R. C.
Sprott
J. C.
Am. J. Phys.
,
1994
, vol.
62
pg.
861

Hill
G. W.
MNRAS
,
1878a
, vol.
38
pg.
192

Hill
G. W.
Astron. Nachr.
,
1878b
, vol.
91
pg.
251

Hinse
T. C.
Michelsen
R.
Jørgensen
U. G.
Goździewski
K.
Mikkola
S.
A&A
,
2008
, vol.
488
pg.
1133

Hinse
T. C.
Christou
A. A.
Alvarellos
J. L. A.
Goździewski
K.
MNRAS
,
2010
, vol.
404
pg.
837

Kipping
D. M.
Bakos
G. Á.
Buchhave
L.
Nesvorný
D.
Schmitt
A.
ApJ
,
2012
, vol.
750
pg.
115

Koch
D.
, et al.  .
Hartkopf
W. I.
Guinan
E. F.
Harmanec
P.
Proc. IAU Symp. 240, The Kepler Mission and Eclipsing Binaries.
,
2007
Cambridge
Cambridge Univ. Press
pg.
236

Kostov
V. B.
McCullough
P. R.
Hinse
T. C.
Tsvetanov
Z. I.
Hébrard
G.
Díaz
R. F.
Deleuil
M.
Valenti
J. A.
ApJ
,
2013
, vol.
770
pg.
52

Lissauer
J. J.
, et al.  .
Nat
,
2011
, vol.
470
pg.
53

Lissauer
J. J.
, et al.  .
ApJ
,
2013

preprint (arXiv:e-prints)
Lyapunov
A.
Ann. Fac. Sci. Toulouse
,
1907
, vol.
9
pg.
203

Maffione
N. P.
Darriba
L. A.
Cincotta
P. M.
Giordano
C. M.
Celest. Mech. Dyn. Astron.
,
2011
, vol.
111
pg.
285

Marchal
C.
Bozis
G.
Celest. Mech.
,
1982
, vol.
26
pg.
311

Mayor
M.
Queloz
D.
Nat
,
1995
, vol.
378
pg.
355

Mestre
M. F.
Cincotta
P. M.
Giordano
C. M.
MNRAS
,
2011
, vol.
414
pg.
L100

Mikkola
S.
Innanen
K.
Celest. Mech. Dyn. Astron.
,
1999
, vol.
74
pg.
59

Morbidelli
A.
Modern Celestial Mechanics : Aspects of Solar system dynamics. Taylor & Francis, London
,
2002
Murray
C. D.
Dermott
S. F.
Solar System Dynamics.
,
1999
Cambridge
Cambridge Univ. Press
Murray
N.
Holman
M.
Nat
,
2001
, vol.
410
pg.
773

Musielak
Z. E.
Musielak
D. E.
Int. J. Bifurcation Chaos
,
2009
, vol.
19
pg.
2823

Ott
E.
Chaos in Dynamical Systems.
,
1993
Cambridge
Cambridge Univ. Press
Ozorio de Almeida
A. M.
Hamiltonian Systems.
,
1990
Cambridge
Cambridge Univ. Press
Quarles
B.
Eberle
J.
Musielak
Z. E.
Cuntz
M.
A&A
,
2011
, vol.
533
pg.
A2

Raghavan
D.
Henry
T. J.
Mason
B. D.
Subasavage
J. P.
Jao
W.-C.
Beaulieu
T. D.
Hambly
N. C.
ApJ
,
2006
, vol.
646
pg.
523

Raghavan
D.
, et al.  .
ApJS
,
2010
, vol.
190
pg.
1

Schneider
J.
Dedieu
C.
Le Sidaner
P.
Savalle
R.
Zolotukhin
I.
A&A
,
2011
, vol.
532
pg.
A79

Słonina
M.
Goździewski
K.
Migaszewski
C.
Arenou
F.
Hestroffer
D.
Proc. Orbital Couples: Pas de Deux in the Solar system and the Milky Way, Mechanic: A New Numerical MPI Framework for the Dynamical Astronomy.
,
2012
Paris
Obervatoire de Paris
pg.
125

Szebehely
V.
Theory of Orbits. The Restricted Problem of Three Bodies.
,
1967
New York
Szebehely
V.
Celest. Mech.
,
1980
, vol.
22
pg.
7

Szebehely
V.
McKenzie
R.
Celest. Mech.
,
1981
, vol.
23
pg.
3

Szenkovits
F.
Makó
Z.
Celest. Mech. Dyn. Astron.
,
2008
, vol.
101
pg.
273

Tsonis
A.
Chaos: From Theory to Applications
,
1992
Plenum Press, New York
Walker
I. W.
Roy
A. E.
Celest. Mech.
,
1981
, vol.
24
pg.
195

Wolf
A.
Swift
J. B.
Swinney
H. L.
Vastano
J. A.
Physica D Nonlinear Phenom.
,
1985
, vol.
16
pg.
285

Yoshida
H.
Phys. Lett. A
,
1990
, vol.
150
pg.
262

Zhou
L.-Y.
Dvorak
R.
Sun
Y.-S.
MNRAS
,
2009
, vol.
398
pg.
1217