On the stability of planetary orbits in binary star systems I. The S-type orbits

Many exoplanets are discovered in binary star systems in internal or in circumbinary orbits. Whether the planet can be habitable or not depends on the possibility to maintain liquid water on its surface, and therefore on the luminosity of its host stars and on the dynamical properties of the planetary orbit. The trajectory of a planet in a double star system can be determined, approximating stars and planets with point masses, by solving numerically the equations of motion of the classical three-body system. In this study, we analyze a large data set of planetary orbits, made up with high precision long integration at varying: the mass of the planet, its distance from the primary star, the mass ratio for the two stars in the binary system, and the eccentricity of the star motion. To simulate the gravitational dynamics, we use a 15th order integration scheme (IAS15, available within the REBOUND framework), that provides an optimal solution for long-term integration. In our data analysis, we evaluate if an orbit is stable or not and also provide the statistics of different types of instability: collisions with the primary or secondary star and planets ejected away from the binary star system. Concerning the stability, we find a significant number of orbits that are only marginally stable, according to the classification introduced by Musielak et al. 2005. For planets of negligible mass, we estimate the critical semi-major axis $a_c$ as a function of the mass ratio and the eccentricity of the binary, in agreement with the results of Holman and Wiegert 1999. However, we find that for very massive planets (Super-Jupiters) the critical semi-major axis decrease in some cases by a few percent, compared to cases in which the mass of the planet is negligible.


Introduction
The first possible observational evidence of a Jupiterlike planet in the binary star system gamma Cephei was reported by Campbell et al. (1988), based on the measurements of the radial velocity variations in a sample of stars. However, due to possibility that the low signal was due to chromospheric activities of the star, the existence of a planet in this system was subsequently questioned (Walker et al. 1992). It was thanks to more accurate measurements that the true nature of the gamma Cephei planet was finally confirmed (Hatzes et al. 2003), fifteen years after the first observation. It is worth noticing that the small distance (about 20 AU) between the two stars in gamma Cephei has as a consequence a complex dynamics of the planetary orbit, whose stability is not guaranteed. Assuming the system as a classic Newtonian, isolated, three-body problem, no general analytical solution exists, also in the special case where the mass of the planet is negligible, i.e. for the so called restricted three body problem. A recent review of the three-body problem in the context of both historical and modern developments is presented by Musielak and Quarles (2014). To date, the catalogue 1 maintained by Schwarz et al. (2016) reports 103 binary and 26 multiple confirmed star systems hosting planets, and 28 candidates. The vast majority of multiple systems hosting planets are made up by three stars.
In some case, the binary (or multiple) star system hosts more than one planet.
Orbits of planets in binary systems are traditionally classified into three categories (Dvorak 1986): i) the Planet-type (P-type) external orbits around both stars in the binary, ii) the Satellite type (S-type) internal orbits, around one of the two stars, and iii) the Librationtype (L-type) orbits, corresponding to librations around the Lagrangian equilibrium points L 4 and L 5 , which are stable when the stellar mass ratio m 1 /(m 1 + m 2 ) is less than ∼ 0.04 (assuming m 1 < m 2 ). While L-type orbits are not normally of interest for exoplanets in binary systems, P-type and S-type orbits are of relevant astronomical impact and deserves to be deeply studied in their characteristics. For a recent review of planet formation and dynamical evolution in binary systems see Marzari and Thebault (2019).
The stability problem of planetary orbits in binary star systems have been investigated by many authors, with different methods and in different schemes. After pioneering work by Henon (1968Henon ( , 1969, Hénon and Guyot (1970), and Szebehely (1980), numerical approaches have been done by Dvorak (1986) (P-orbits in restricted assumptions), Rabl and Dvorak (1988) (Sorbits in restricted assumptions), and Pilat-Lohinger and Dvorak (2002) (S-orbits with less restrictions). Dvorak (1986) and Rabl and Dvorak (1988) numerically integrated a set of orbits by mean of a scheme based on Lie-series, as developed by Delva (1985) whose drawback is that of being limited to infinitesimal mass for the third-body, in what the primaries motion is assumed as the analytical solution of the two-body problem. The work done by Pilat-Lohinger and Dvorak (2002) using Bulirsh-Stoer algorithm, does not encounter in principle this limitation, but it was applied by the authors, again, to zero mass planets. Later on Holman and Wiegert (1999), hereafter referred to as H99, using a modified symplectic mapping technique initially developed by Wisdom and Holman (1991) to study the secular behavior of planetary motion around a dominating mass, investigated the planet orbital stability extending significantly with respect to previous works the time basis of integration; the authors assume that the mass of the planet is negligible. A straightforward Runge-Kutta-Fehlberg numerical technique was adopted by Musielak et al. (2005) who dealt with the stability of initially circular planetary orbits in both S-and P-configurations in the field of two stars revolving each other in circular motion. The Frequency Map Analysis (FMA) method (Laskar et al. 1992) was for the first time applied by Turrini et al. (2004Turrini et al. ( , 2005 to investigate the dynamical stability of the giant planet in the γ Cephei binary system. The FMA method was also applied by Marzari and Gallina (2016), finding a significant number of stable planetary orbits in binaries beyond the critical distance from the primary star as evaluated by H99. Fatuzzo et al. (2006), by performing a huge set of numerical integrations of Earth-like planetary orbits in binary systems, report the statistical distribution of the survival times of planets in these systems. The authors conclude that in principle some unstable planets could have a sufficiently long survival time in the habitability zone. In systems with multiple planets orbiting around a single star, the minimum distance beyond which there are no close encounters is determined by the Hill stability limit (Marchal and Bozis 1982;Gladman 1993). All these studies and results stimulated more specific works aimed on observational implications, including considerations on the habitability regions (see, for instance Cuntz 2014Cuntz , 2015. Just a few planetary systems are detected in star cluster, and so far only one (PSR B1620-26 b) in a globular cluster. The survivability and characteristics of planetary systems in star clusters was recently investigated by Cai et al. (2018Cai et al. ( , 2019; van Elteren et al. (2019); Stock et al. (2020). Among the observations of planetary systems in binaries, the case of ν Octantis is enigmatic since the planet is located outside the zone of orbital stability, as estimated by H99 for prograde orbits. To solve this puzzle, Eberle and Cuntz (2010) suggested that the planet moves in a stable retrograde orbit. In a system (K2-290 ) composed by three stars, retrograde motions with respect to the primary star's spin have been recently observed (Hjorth et al. 2021).
The secular motion has been considered in different works. In general, the classical analysis of secular motion is not accurate for objects in highly eccentric orbits; Michtchenko and Malhotra (2004) introduce a semi-numerical approach to study the secular motion of two massive planets in co-planar orbits on high eccentricity (0.1 -0.6) orbits. Heppenheimer (1978) introduced an analytical approximation for the three-body problem, providing a description of the secular orbital evolution for a planet in a binary star system. To assess the range in parameter space in which secular motion models provide accurate results, Andrade-Ines et al. (2016) compared predictions for secular motion based on first and second order analytical models with Nbody simulations. A semi-analytical correction to the Heppenheimer (1978) solution was recently presented by Andrade-Ines and Eggl (2017); this correction provides a quite accurate description of the S-type planetary secular motion in a given range of parameters. Quarles et al. (2020) recently performed a numerical study of the stability of S-type orbits following an approach similar to that of Holman and Wiegert (1999) enlarging their initial parameter space. In particular, Quarles et al. (2020) probe the role of inclined planetary orbits by testing 4 inclinations for eccentric binaries, and 8 inclinations for circular binaries, integrating the motion of Earth-mass particles with a symplectic scheme. With respect to Quarles et al. (2020), we probe a wider planet mass parameter space by simulating the motion of test particles and planets with 1 and 30 Jupiter-mass. Moreover, we provide a deeper characterization of the trajectories by recording, for unstable orbits, ejections or collisions with one of the two stars, and, for stable orbits by evaluating the average distance of the planet from the unperturbed orbit, i.e. the orbit that the planet would have without the secondary star.
The paper is organized as follows. In Section 2 we define the dataset made up with a huge number planetary orbits that we have numerically integrated at different initial conditions, and we describe the data analysis. In Section 3 we report on the results. We end up with the discussion in Section 4. Fig. 1 Distribution of the masses for the exoplanets discovered in binary star systems (from exoplanets.org). Although the majority of the planets have a small mass, a significant fraction of them has a mass equal to the mass of Jupiter or more.

Code and dataset
We approach the problem of the planet orbital stability in a numerical way, by integrating the orbits under different initial conditions. We consider a system made up by two star A (called primary) and B (called secondary) in motion around their center of mass, with the planet in a co-planar orbit with zero initial eccentricity. With unperturbed orbit we will call the orbit around the primary star that the planet would have without the presence of the secondary star; in our configuration the unperturbed orbit is circular. Two examples of planetary orbits are shown in Fig. 2 (top panel). We notice that the trajectory of the planet fills a ring nearby the unperturbed orbit, whose "thickness" is a consequence of the gravitational field of the secondary star.
Many ordinary differential equations (ODEs) integration methods can be applied to calculate the orbits of a gravitational N-body system (Aarseth 2003). As test particles do not interact with each other, H99 simulate large numbers of single planet systems in one single integration. For finite mass planets this approach is not possible and each simulation must be managed individually. H99 use a mixed approach for the numerical integration of the orbits. When the planet is near the primary star, the symplectic integration scheme Wisdom and Holman (1991) is applied. This technique is accurate for motion around a single star, but also in the case of a binary star system it provides a reasonable approximation, the gravity of the primary star being dominant. Far from this condition, a Bulirsh-Stoer (BS ) integration method is applied; this algorithm combines a fairly accuracy at a relatively little computational cost (Press et al. 2007). To achieve a very accurate integration of the orbits, in this study we will use the highprecision integration scheme IAS15 (Rein and Spiegel 2015) throughout all the duration of the motion. With IAS15, a 15th-order integrator, the systematics errors are kept below machine precision for long term integration over at least 10 9 orbital time scale. IAS15 is an integration options available in REBOUND Rein and Liu (2012), a modern open source code for gravitational dynamics. REBOUND is written in standard C99; an easy to use and convenient python wrapper is also provided.
The planetary system is made up by two stars: primary star A, and the secondary star B, in elliptical motion in the center of mass frame of the two stars. The distance of the planet from the primary star plays a crucial role on the stability. Generally speaking, we expect that the farther the planet is from the primary star, the stronger is the perturbation due to the secondary star, so that the orbit can eventually become unstable. The unperturbed circular orbit is characterized by the initial (transverse) speed of the planet v P given by: where r is the initial distance of the planet from the primary star, m P is the mass of the planet and m A is the mass of the primary star. Following H99, we consider as free parameters the mass ratio and the eccentricity of the binary star system, and the initial semi-major axis of the planet orbit. The mass ratio is defined as: For equal mass stars, µ = 0.5. The smaller the µ the lower the gravitational perturbation due to the secondary star. As in the study of H99, in our dataset the mass ratio of the binary µ goes from 0.1 to 0.9 (with ∆µ = 0.1) and the eccentricity e from 0 to 0.8 (with ∆e = 0.1). The semi-major axis of the binary is equal to 1 AU, and the binary initial phase is at periapse or at apoapse. Unlike the study of H99, where only planets with a negligible mass are considered, we also study the motion of planets with finite non negligible mass. We made numerical integration for a planet of infinitesimal mass and of mass equal to 1 and 30 Jupiter masses, which is a reasonable choice, taking into account the mass distribution of observed for planets in binary systems (Fig. 1). By putting all these degrees of freedom together, our dataset is made by more than ten thousand orbits. Concerning the masses of the planets, Stevens and Gaudi (2013)  The integration time extension is an important parameter of the numerical study. Actually, long-term integration provides a greater confidence about the stability of a planetary orbit, but this obviously implies a cost in term of computational time. The duration of the integration should be a trade-off between the duration of the numerical integration and the accuracy of the results. Following the studies of Holman and Wiegert (1999), an integration time equal to 10 000 star orbits it is sufficient to obtain reliable results.

Analysis of the orbits
We classify the simulated orbits on the basis of their stability according to the following definitions. We consider an orbit as unstable if (i) the planet will collide in a finite time with one of the two stars or (ii) it is ejected from the binary system. We catch a collision if the star-planet distance reduces to a solar radius, an ejection if the planet moves away from the center of mass of the binary system by more than 100 AU. In these cases we record the orbit lifetime (aka survival time), that is the time elapsed from the beginning of the orbit to the collision or the escape. In this classification, we do not explicitly refer to chaos, i.e. that sensitivity to initial conditions that makes impossible an accurate long term prediction of the position and velocity of the system particles. Orbits that we classify as stable, can have chaotic behavior. In any case, all the stable orbits we find in our sample are well confined in a ring surrounding the unperturbed orbit (Fig.  2). In general, stable orbits are characterized by how much they deviate from the unperturbed orbit; the orbits close to the primary star are almost circular, and little affected by the gravity of the secondary star. A simple way to characterize the stable orbits is through the difference ∆r of the distance r of the planet from the primary star with respect to the initial value mediated along the orbit: where the average is made over the entire simulated orbit of the planet, so T corresponds to 10 000 orbital periods the binary. Hereafter this mean is referred simply as ∆r. Of course, ∆r = 0 only if there is no interaction with the secondary star (m B = 0). Musielak et al. (2005) classified the planetary orbits in binary on the basis of the value of the above parameter (Eq. 3): (1) stable if ∆r ≤ 5%, (2) marginally stable if 5% < ∆r ≤ 35%, (3) unstable if ∆r > 35%. The 5% threshold for stability is motivated by some studies (e.g. Kasting et al. (1993), Underwood et al. (2003)) which show that this limit is required for the Earth to remain in the habitable zone and therefore allow the evolution of life.

Results
We report in this section how the stability depend on the mass of the planet, the initial semi-major axis of the planetary orbit, and on the mass ratio and eccentricity of the binary star system.

Collisions and ejections
We introduce with some examples 2 the statistics of possible collisions of the planet with one of the two stars of the binary and of ejections from the binary system. As Fig. 2 Top left and right panels: a stable (∆r = 3.9%), and a marginally stable (∆r = 7.9%) orbit (see Section 2.2). Bottom left, right panels: a collision of the planet with the secondary star, and an ejection of the planet away from the binary. The red circle is the orbit of the secondary star in the frame centered in the primary, the blue circle is the orbit of the unperturbed planetary motion.

Fig. 3
A peculiar unstable orbit of a planet of 30 Jupiter masses. At first, the planet moves away from the main star to be captured by the gravitational field of the secondary star and ejected away from the system. The orbit of the secondary star is significantly modified into an elliptic orbit. Fig. 4 Synopsis of stability of the studied orbits for three different values of the planet mass (0, 1, 30 Jupiter mass) in a binary with µ = 0.5, e = 0.5. ST indicate stable and marginally stable orbits, CA and CB indicate orbits leading to a collision with the primary or secondary star. Orbits leading the planet to be ejected away from the binary star system are labelled with EJ. Error bars represent 1 σ confidence intervals. a demo, we built a real-time animation based on a 4th order Runke-Kutta integrator 3 . The two orbits shown on the top left and right panel of Fig. 2 are classified as stable and marginally stable, on the basis of their estimated ∆r (see Section 2.2). As the full planetary orbit is too long to be displayed effectively, the plot is obtained integrating the final part of motion within a time interval T = 30 000 days. We notice that despite ∆r of the two orbits differ significantly, the planetary orbit shows similar characteristics.
The bottom left panel of Fig 2 shows a collision with the secondary star. The orbit is the result of the numerical integration from the initial time to the star collision. The bottom right panel of Fig. 2 shows a orbit where the planet is ejected away from the binary. It is interesting to notice that, before escaping, the planet moves in an outer orbit. This leads to the hypothesis that, under certain conditions, a transition from an internal (S-type) orbit to an external (P-type) orbit could occur. An interesting orbit is shown in Fig. 3. In this case, before being ejected away from the binary system, the massive planet scatters with the secondary star, producing a significant perturbation of the stellar orbit. As the consequence of the interaction planetstar, the star moves from a substantially circular orbit to a wider elliptical orbit. Figure 4 reports the statistics of planet-star collision and ejection events for a binary with e = 0.5 and µ = 0.5. In this particular case, the number of collisions and ejections do not depend on the planetary mass.

Stability of the orbits
For a given binary, the planetary motion depends mainly on the initial semi-major axis a of the planetary orbit. In general, a planet with a small a, i.e. revolving close to the primary star, has a stable orbit. We apply here the idea, introduced by H99, of critical semimajor axis a c as threshold from stability to instability, to evaluate how orbital stability depends on the mass ratio µ and the eccentricity e of the binary star system. We notice that a c is well defined only if the boundary between the region of stability (where a ≤ a c ) and instability where (a > a c ) is sharp. In Tables 1 and 2 we report the results for a binary star system with e = 0.5 and µ = 0.5; in Table 1 the binary is initially at periapse, in Table 2 the binary is initially at apoapse. Our choice of considering eight values (0, π/4, π/2, (3/4)π, π, (5/4)π, (3/2)π, (7/4)π, 2π) for the initial longitude Planet of negligible mass, µ = 0.5, e = 0.  Table 1 Characterization of the planetary orbits for a given value of mass ratio µ = 0.5, and eccentricity of the binary e = 0.5. a is the initial major semi-axis or the planetary orbit. The calculation is made for 8 equispaced longitudes (P1-P8) with the binary initially at periapse. Each cell in this table is associated with a complete simulation of a planetary orbit. The critical semi-major axis is ac = 0.12. For stable orbits (a ≤ 0.12) the value of ∆r is reported. For unstable orbits, we use a code to indicate whether the instability is due to a collision with the primary or secondary star (CA, CB) or if the planet is ejected away (EJ). The survival time of the orbit in units of 10 binary periods is quoted in brackets.
Planet of negligible mass, µ = 0.5, e = 0.5 Binary Initially at Apoapse a (AU)  Table 2 Characterization of the planetary orbits for a given value of mass ratio µ = 0.5, and eccentricity of the binary e = 0.5 (see Table 1). The calculation is made for 8 equispaced longitudes (A1-A8) with the binary initially at apoapse.  Table 3 Critical semi-major axis ac for a planet mass (m) of 1, 30 Jupiter mass (MJ), and for a planet of negligible mass (test particle), as a function of the mass ratio µ and the eccentricity e of the binary. The text is in boldface when the ac estimation differs from the value obtained for a test particle.  of the planet allows a robust estimate of a c = 0.12, that is the same value obtained by H99. We characterize the stable orbits by the value of ∆r (see end of Sect. 2). For stable orbits, in some cases we estimate a value of ∆r greater than 5%, which define a marginally stable orbits, according to the classification introduced for the first time by Musielak et al. (2005). For unstable orbits, we evaluate whether i) there is a collision with the primary star (CA) or ii) with the secondary star (CB), or iii) if the planet is ejected away from the binary star system (EJ). In the particular case of a binary with µ = 0.5 and e = 0.5 (Tables 1 and 2), most of the instabilities correspond to collisions with the primary star, but there are also ejections and a few collisions with the secondary star. We obtain different survival times than H99, likely due to the different, and much more accurate, integration scheme we adopted. We find some stable orbits beyond the critical semi-major axis (a c = 0.12): one is quoted in Table 1 with a = 0.13, four are quoted in Table 2 with a = 0.13, 0.14. We do not find mean motion resonances in these orbits. The region of stability (a < a c ), in which all the orbits are stable at any longitude, is well defined for all the binaries systems that we considered. The dependence of the critical semi-major axis a c on the mass ratio µ (from 0.1 to 0.9) and eccentricity e (from 0 to 0.8) of the binary star system for a planet of negligible mass and equal of 1 and 30 Jupiter mass is reported in Table 3. The results for planets of negligible mass are in agreement with H99. In most cases, we find that the value of a c does not depend on the mass of the planet, with the exception of some entries, that show different values of a c , especially for 30M J . However this result is not reported in Table 3 with a high statistical significance. The large error is due to the fact that a c values are using a ∆a c step equal to 0.01 (see Table  1 and 2). In order to obtain a greater resolution, we recalculated a c in four cases using a smaller ∆a c step ( Table 4). The new outcomes confirm a significant dependence for planets of mass equal to 30M j . We find that a mass of the planet equal to M j does not produce significant effects in our simulations, with the exception of the (e, µ) = (0.2, 0.50) configuration, where however the significance is quite low. In all the 81 binary configurations we find that a fraction of the orbits are stable beyond the critical semi-major axis. Considering a narrow neighborhood of a c , the number of stable orbits strongly depends on the (µ, e) value. In the case reported in Table 1 and 2 (µ = 0.5, e = 0.5, a c = 0.12) there are 4 stable orbits with a = 0.13; the fraction of stable orbits in this case is 4 over 16 (25%). The mean number of stable orbits beyond a c over all binary configurations does not depend significantly on the mass of the planet. Averaging over all binary configurations, the fraction of stable orbits in a narrow neighborhood of a c is 55%.

Discussion
In our numerical study, we simulated more than ten thousand S-type orbits of planets in a binary star system, using a high precision, 15th order, integration scheme. We estimated the critical semi-major axis, which defines the size of the stable regions, as a function of the mass ratio and orbital eccentricity of the binary star system hosting the planet, finding results in agreement with H99 for planets with negligible mass.
We also studied the motion for planets with a mass of 1 and 30 Jupiter. This is a reasonable choice, because it covers the observational data. Indeed, looking at the distribution of the masses measured in binary systems ( Fig. 1), we note that the majority of the planets have a mass equal to about two hundredths of the Jupiter mass, with a distribution of masses extending up to 16.1 Jupiter mass (Liu et al. 2008). Considering all the planetary systems (i.e. including those around single stars), the largest mass reported in the data archive (exoplanets.org) is equal to 22.6 Jupiter mass (Sato et al. 2010); in this extreme case the small mass tertiary around the binary falls in the brown-dwarf regime. Marzari and Gallina (2016), using the FMA method to select stable orbits, found a significant numbers stable planetary orbits in binaries beyond this critical distance from the primary star. In general, this result does not disagree with our simulations, since we find in many cases stable orbits beyond the critical distance, as shown for example in Tables 1 and 2 for two stars of equal mass and eccentricity equal 0.5. For very massive objects (30M j) we find differences in the a c estimation of the order of 3% compared with test (negligible mass) particles (see Table 4). It is worth noting that we find these dependence on the planetary mass for low binary mass ratio µ. This is expected, since in these cases the mass of the secondary star is smaller than that of the primary star, so where the mass of the planet plays a more important role.
For all the stable orbits we estimated the spread (∆r) of the distance of the planet from the primary star with respect to its initial value, showing that a fraction of these orbits is only marginally stable, according with the classification of Musielak et al. (2005). Our data analysis also produces statistics of planet-star collisions and planetary ejections away from the binary, in the cases of planet mass equal to 0, 1, 30 Jupiter masses. In the particular case of stars of equal mass and moving on an e = 0.5 orbit, for all planet masses studied here, the largest fraction of unstable orbits result into a collision of the planet onto the primary star, followed by ejection and collisions onto the secondary. A planetstar collision might result in the pollution of stars with planetary material, which, in principle, would be detectable by accurate spectroscopic measurements. The measurements of iron abundance for a sample of 23 wide binaries reported by Desidera et al. (2004) seem to rule out a contamination with a mass of iron greater than 1 Earth Mass during the Main Sequence lifetime of the stars. The instability of planetary orbits in double star systems with a final fate as ejection of the planet away from the binary system has an relevance in increasing the number of freely-floating planets as detected by gravitational microlensing (Sumi et al. 2011).
In terms of the parameter space, we made some assumptions. We considered binary systems hosting a single planet moving on an unperturbed circular co-planar orbit around the primary star. For the data analysis, we implemented a code based on standard Python tools (numpy, pandas). A joint application of the high-order integration code and our data analysis tools on a High Performance Computing (HPC) platform will allow a significant extension of the space of parameters to investigate. An interesting approach to analyse the large amount of output data that will be produced is based on Machine Learning techniques (see e.g. Tamayo et al. (2016)).
An interesting perspective is to investigate our results in the scenario of the subresonances theory. Indeed, the planet experiences resonant perturbations from the companion star. Starting from the observation that planet may be dislodged from its host star if it is simultaneously affected by two or more resonances, Mudryk and Wu (2006) found that overlap between subresonances lying within mean-motion resonances can explain the boundary of orbital stability within binary systems.