A dynamical stability study of Kepler Circumbinary Planetary systems with one planet
Abstract
To date, 17 circumbinary planets have been discovered. In this paper, we focus our attention on the stability of the Kepler circumbinary planetary systems with only one planet, i.e. Kepler16, Kepler34, Kepler35, Kepler38, Kepler64 and Kepler413. In addition to their intrinsic interest, the study of such systems is an opportunity to test our understanding of planetary system formation and evolution around binaries. The investigation is done by means of numerical simulations. We perform numerical integrations of the full equations of motion of each system with the aim of checking the stability of the planetary orbit. The investigation of the stability of the above systems consists of three numerical experiments. In the first one we perform a long term (1Gyr) numerical integration of the nominal solution of the six Kepler systems under investigation. In the second experiment, we look for the critical semimajor axis of the six planetary orbits, and finally, in the third experiment, we construct two dimensional stability maps on the eccentricitypericentre distance plane. Additionally, using numerical integrations of the nominal solutions we checked if this solutions were close to the exact resonance.
keywords:
(Stars:) binaries: general, Celestial mechanics, Planets and satellites: dynamical evolution and stability, Methods: numerical1 Introduction
Before 1992, the only planets we knew were the ones of our own Solar System. To date, more than one thousand and five hundred planets are known to revolve around stars other than the Sun. Many of the new planetary systems exhibit widely characteristics from our Solar System (e.g. medium or high eccentricities and the presence of Jupiter sized bodies close to the star) and consequently this has led to a major revision the planetary formation and evolution theories. Modern observations give values of up to 70 for the frequency of multiple stellar systems in our Galaxy (KiselevaEggleton and Eggleton 2001 and references therein), which makes the study of the process of planet formation and evolution in such systems, an important field of study.
Nowadays, over seventy binary systems have been discovered harbouring planets (e.g. see Roell et al. 2012) and even more are expected to be found. The interest of the scientific community, regarding, besides planet formation, various other issues related to planets in stellar binaries, such as habitability (e.g. Eggl et al. 2012, Eggl et al. 2013, Kane & Hinkel 2013), has greatly increased. Although many questions still remain, significant effort has been channeled into formation studies in stellar binaries, which is also a test field for theories of planetary formation around single stars. Planet formation in binaries is currently a hotly debated topic that puzzles the scientific community, especially as in some of the discovered systems the circumstellar planet is very close to its parent star (e.g. the Cephei system, Hatzes et. al 2003, for which the planet semimajor axis is 2 AU and the companion star has semimajor axis 20 AU and eccentricity 0.36). Paardekooper et al. (2012), Meschiari (2012), Pelupessy & Portegies Zwart (2013), Thebault & Haghighipour (2014) may provide the reader with relevant information regarding current issues of planet formation in stellar binaries.
Most of the planets that have been discovered in stellar binaries are in circumstellar orbits. The first planet orbiting both stars of a stellar binary was found in the pulsar system PSR B162026 (Backer et al. 1993) with a minimum mass between 1.5 and 3.5 M. More recently, a system with two circumbinary planets was found to orbit HW Virginis (Lee, J. W. et al. 2009) and another pair of planets was detected around the eclipsing binary NN Serpentis (Qian et al. 2009, Beuermann et al. 2010). A planet was found in another eclipsing binary, DP Leonis (Qian et al. 2010). In the same year, a planet orbiting DT Virginis, also known as Ross 458 AB, was detected by direct imaging (Goldman et al. 2010, Burgasser et al. 2010). At around the same time (Chavez et al. 2012), the light curve of the cataclysmic variable FS Aurigae was found to show the presence of a third body. In that study, it was found that the most plausible explanation to the observed variations in the light curve of the system was the presence of a third body on a wide orbit with a mass of 48 M, but a mass as small as 2 M was not ruled out.
The first circumbinary planet in the Kepler catalog was Kepler16b (Doyle et al. 2011); this circumbinary planet was the first one to be confirmed by two different detection methods, i.e. by the transit method and by RV measurememnts. As a consequence, its parameters are better constrained compared to the rest of the circumbinary planets and the same can be said for the rest of the circumbinary planets that have now been confirmed by Kepler. A new planet was found to orbit NY Virginis in 2012 (Qian et al. 2012) and shortly after, another planet (possibly accompanied by an additional one) was reported to orbit the RR Cae system, an eclipsing binary containing white dwarfs (Qian et al. 2012b). More circumbinary planets were found by Kepler: Kepler34b and Kepler35b (Welsh et al. 2012), Kepler38b (Orosz et al. 2012a), Kepler47b and c (Orosz et al. 2012b, Kostov et al. 2013), Kepler64b (Kostov et al. 2013, Schwamb et al. 2013) and finally, the most recent of such Kepler discoveries has been Kepler413b (Kostov et al. 2014).
Previously, there has been some work regarding the stability of those systems. Doyle et al. (2011) integrated numerically the bestfitting solution for Kepler16b for two million years and they noticed no significant excursions in orbital distance that might lead to instability. The stability of the Kepler16 system was also investigated by Jaime et al. (2012). In the context of determining areas of stable orbits around stellar binaries based on the concept of ’invariant loops’ (e.g. Pichardo et al. 2005), they found that the Kepler16b planet was lying in the stable region around the stellar binary, but not far away from the unstable area. In addition, Popova and Shevchenko (2013) investigated the stability of the Kepler16, systems by constructing stability maps based on the Lyapunov exponent concept and on an ’escapecollision’ criterion, acording to which an orbit was classified as stable if the distance between the planet and one of the stars did not become less than 10 AU or greater than 10AU. They found that the planet Kepler16b turns out to be just outside the chaotic area, among some ’teeth’ of chaotic motion corresponding to certain mean motion resonances. Additionally, they found similar results when applied their method to Kepler34 and Kepler35.
Welsh et al. (2012), along with the announcement of the discovery of Kepler34b and Kepler35b, performed a series of numerical integrations in order to assess the stability of the two planetary systems. Using the nominal orbital solutions, they performed direct Nbody integrations for 10 Myrs and found no indication of instability. They also carried out 1000 simulations of systems with different masses and orbital parameters for 1 Myrs, without finding, again, any sign of instability. Moreover, they integrated an ensemble of a few thousand threebody systems, each consistent with the observed masses and orbital parameters, except that they varied the semimajor axis of the planet with the aim of finding the minimum semimajor axis for which the system was not flagged as unstable (i.e. change of the initial value of the planetary semimajor axis by more than 50%) during the 10000 yr simulations. Besides Kepler34b and Kepler35b, that was also done in the same article for Kepler16b. Regarding Kepler38b, Orosz et al. (2012a) tested the stability of the system by performing direct Nbody simulations and found the critical period to be 75 days. Kostov et al. (2013) and (2014) used the MEGNO method to investigate the stability of the Kepler64 and the Kepler413 systems respectively. Carrying out numerical simulations with an integration time of 2700 yrs (which corresponded to binary periods) for the former and 200000 days (which corresponded to about binary periods) for the latter, and after they constructed dynamical maps, they concluded that the bestfit orbital parameters located the planets in a quasiperiodic region in (a, e) space, rendering the orbital solution to be plausible from a dynamical point of view. Regarding Kepler64, that conclusion seems to be confirmed by Schwamb et al. 2013, who carried out numerical integrations of the system using the best fit physical and orbital parameters and found the system to be stable on gigayear timescales, without giving any further information though.
In this paper, we focus our attention on the long term stability and dynamical evolution of the Kepler planetary systems with only one circumbinary planet. The investigation is done by means of numerical simulations. The paper is structured as follows: in section 2, we describe the method we use and the setup of our numerical experiments. The results of the various simulations conducted are presented in section 3. In section 4 we discuss the implications of our results and present concluding remarks.
2 Numerical modelling and Method
In this paper, we focus on the long term evolution and stability of the six singleplanet circumbinary Kepler planetary systems. In the context of the three body problem, we perform numerical integrations of the full equations of motion. Our investigation is split into three numerical experiments.
In the first experiment, we used the nominal orbital solution of the six systems, as given in the relevant literature, and we integrated their motion for 1 Gyr. This a reasonable amount of time to choose for our purposes as it is a considerable fraction of the time that the stars will remain on the main sequence and it is also much longer compared to almost all previous studies mentioned above, as most of the latter investigated the stability of the systems for a rather limited time interval ( Myrs). All orbits were assumed coplanar, as the systems under investigation have a mutual inclination of the order of just a few degrees and using a three dimensional model would have an insignificant effect on the outcome (e.g. see Doolin and Blundell 2011, Georgakarakos 2013).
The initial conditions for our simulations were derived from the orbital elements of the systems that were published in the relevant papers. More specifically, for Kepler16, Kepler34 and Kepler35, we used the data that appear in ’Supplementary Table 1’ in ’Supplementary Information’ of Welsh et al. (2012). For Kepler38, we obtained the values of all parameters from Tables 5, 6 and 7 that appear in Orosz et al. (2012a), using the column ’Best Fit’ in all cases (for the planetary mass we used the maximum proposed value). For the case of Kepler64 we used two set of values since there are two published articles for this system. The first set appears in Tables 2 and 3 in Kostov et al. (2013) (again, for the planetary mass, we used the maximum proposed value), while the second set appears in Tables 7 and 8 in Schwamb et al. (2013). Finally, for Kepler413 we used the ’Best Fit’ column in Tables 3 and 4 in Kostov et al. (2014). The mean longitudes for Kepler16, Kepler34 and Kepler35, are given in the corresponding articles, while for the case of Kepler38, Kepler64 (Kostov solution) and Kepler413 the time of pericentre passage of the binary T is given and it can be used to calculate the mean longitude by using Eqs. (2.39) and (2.53) from Murray and Dermott (1999). For the case of the Schwamb solution of Kepler64, we used the same time of pericentre passage of the binary that appears in Kostov et al. (2013), since that parameter does not appear in the corresponding article. In order to compute the missing mean longitudes, we used the fact that the times on which the planet transits the primary are given for all systems, therefore we took the time of the first observed primary transit and we evolved the inner binary from the its passage of the pericentre epoch to that given time. Since the planet has to transit the primary, it means that the Earth, the planet and the primary should be in the same line. The transit times for Kepler38, Kepler64 (Kostov), Kepler64 (Schwamb), Kepler413 were taken from Table 1 of Orosz et al. (2012), Table 3 of Kostov et al. (2013), Table 3 of Schwamb et al. (2013) and Table 5 from Kostov et al. (2014) respectively. Tables 1 and 2 give the orbital elements and masses of the systems under investigation.
System  (AU)  (deg)  (deg)  (deg)  T (JD)  

Kepler16 
0.68970.20255  0.22431  0.15944  90.3401  96.536  92.3520  – 
Kepler34 ^{a}  1.04791.0208  0.22882  0.52087  89.8584  71.4359  300.1970  – 
Kepler35 ^{a}  0.88760.8094  0.17617  0.1421  90.4238  86.5127  89.1784  – 
Kepler38 
0.9490.249  0.1469  0.1032  89.265  268.680  62.96  2454971.66790 
Kepler64 (K) 
1.470.37  0.1769  0.204  87.59  214.3  159.6  2454973.862 
Kepler64 (S) 
1.3840.386  0.1744  0.2117  87.360  217.6  162.6  2454973.862 
Kepler413 
0.8200.5423  0.10148  0.0365  87.332  279.74  356.5  2454973.230 
System  ()  (AU)  (deg)  (deg)  (deg)  (deg)  

Kepler16 
0.333  0.7048  0.00685  90.0322  41.2971  0.003  106.51 
Kepler34 ^{a}  0.22  1.0896  0.182  90.355  82.0928  1.74  106.5 
Kepler35 ^{a}  0.127  0.60345  0.042  90.76  63.4349  1.24  136.4 
Kepler38 
0.384  0.4644  0.032  89.446  32.8285  0.012  274.7 
Kepler64 (K) 
5  0.642  0.1  90.0  105.0  0.0  264.43 
Kepler64 (S) 
0.532  0.634  0.0539  90.022  348.0  0.89  275.4 
Kepler413 
0.21  0.3553  0.1181  90.022  94.60  3.139  283.1 
In the second numerical experiment we searched for the critical planetary semimajor axis, i.e. the smallest planetary semimajor axis for which the planet would be on a stable orbit in the vicinity of the stellar binary, just before the planetary orbit becomes unstable for the first time as we approach the binary. Starting with the orbital elements of each system the same as in Tables 1 and 2, we varied only the planetary semimajor axis with the aim of finding the value for which the planetary orbit would become unstable for the first time as we moved closer to the stellar binary. These simulations were carried out for years (longer than any other similar study). Making a reasonable choice for the initial semimajor axis, it was then reduced by a step of 0.01AU until the planetary orbit became unstable. For each value of the semimajor axis, the planet was started at eight different initial positions, i.e. the initial planetary true anomaly was given values from to with a spacing of . If the planetary orbit became unstable at one or more of the initial positions, then, for that value of the semimajor axis, the planetary orbit was classified as unstable.
Finally, in the third experiment, we searched for unstable areas on the plane ( being the initial planetary eccentricity and being the initial planetary pericentre distance) for each planetary system. We constructed a grid of 10191 models, corresponding to the same number of points, with ranging from 0 to 0.9 and making a reasonable choice for our values by having a look at the stellar apocentre and planetary pericentre distances. As in the search for the critical semimajor axis of the nominal solution, the planet was started at eight different values of the true anomaly. When one of those initial positions led to an unstable planetary orbit, then the specific pair was classified as unstable and we moved on to the next pair of values. Additionally, we used the results of the numerical integrations performed for the nominal solutions in order to check whether or not the nominal solutions were close to the exact resonance, as defined in Murray & Dermott (1999), Eq. (8.22).
In order to perform our numerical simulations, we used the symplectic integrator with time transformation developed by Mikkola (1997), specially designed to integrate hierarchical triple systems. The code uses standard Jacobi coordinates, i.e. it calculates the relative position and velocity vectors of the inner and outer orbit at every time step. Then, by using standard two body formulae, we computed the orbital elements of the two binaries.
All bodies were treated as point masses and only gravitational interactions at the Newtonian level were considered. Effects such as stellar finite size and general relativity were ignored at this point. The approximation of treating the bodies as point masses, in particular for the inner binary, should be taken with some caution. Tidal deformation of the stars as well as general relativity effects might affect the dynamics of the system (e.g. see Soderhelm 1984; Kiseleva, Eggleton & Mikkola 1998, Eggleton, Kiseleva & Hut 1998, Borkovits, ForgácsDajka & Regály 2004). These effects will be discussed in a the next section.
3 Results
3.1 Long term results
First, we present the results from our long term numerical simulations. As we stated previously, the nominal orbital solutions of six Kepler circumbinary systems were integrated forward in time over 1 Gyr. During our simulations, none of the systems demonstrated any sign of instability for the initial conditions used. Figure 1 shows the variation of the planetary semimajor axis and eccentricity for the Kepler64 system (Kostov et al. 2013 solution). For this system, the minimun and maximum semimajor axis values were 0.638 AU and 0.657 AU respectively, while for the eccentricity, the minimum and maximum values were 0.06 and 0.18. Therefore no significant excursions in semimajor axis and eccentricity that might lead to instability are noticed in the two plots.
3.2 Critical semimajor axis
Besides checking the long term stability of the six circumbinary systems, we also found the critical planetary semimajor axis for each system with the method described in the previous section. In addition, we computed the critical semimajor axis for the six Kepler systems under investigation by using the results of some earlier work: the empirical formula of Holman & Wiegert (1999) and the semianalytical stability criterion of Mardling & Aarseth (2001).
Holman and Wiegert (1999) investigated the stability of Ptype and Stype orbits in stellar binary systems. They performed numerical simulations of particles on initially circular and prograde circumstellar or circumbinary orbits, in the binary plane of motion and with different initial orbital longitudes. The binary mass ratio was taken in the range , the binary eccentricity was in the range and the integrations lasted for binary periods. If a particle survived the whole integration time at all initial longitudes, then the system was classified as stable. For circumbinary orbits, and by using a least squares fit to their data, they obtained an expression for the critical semimajor axis for the planet, , given by:
(1)  
where is the binary semimajor axis, is the binary eccentricity and .
Mardling and Aarseth (1999, 2001) approached the problem of stability of a hierarchical triple system by noticing that stability against escape in the three body problem was analogous to stability against chaotic energy exchange in the binarytides problem. The way energy and angular momentum are exchanged between the two orbits of a stable (unstable) hierarchical triple system is similar to the way they are exchanged in a binary undergoing normal (chaotic) tideorbit interaction. Having that in mind, they derived a semianalytical formula for the critical value of the outer pericentre distance :
(2) 
where , is the outer binary eccentricity and is the inner semimajor axis. If , then the system is considered to be stable. The above formula is valid for prograde and coplanar systems and it applies to escape of the outer body. The constant was determined empirically and it was found to be 2.8. A heuristic correction (with being the mutual inclination in degrees) is applied for noncoplanar orbits, to account for the increased stability. Also, the criterion ignores a weak dependence on the inner eccentricity and inner mass ratio. Finally, as stated in Mardling and Aarseth (2001), the formula holds for , beyond which all unstable systems for the planets can suffer exchange; however the formula has not been tested for systems with planetary masses (Aarseth 2004), probably because the authors were mainly interested in using the formula in star cluster simulations.
System  Nominal  Numerical  Holman & Wiegert  Mardling & Aarseth  Published 

(AU)  (AU)  (AU)  (AU)  (AU)  
Kepler16  0.7048  0.67  0.65  0.64  0.59 (W) 
Kepler34  1.0896  1.00  0.84  0.87  0.88 (W) 
Kepler35  0.60345  0.52  0.50  0.53  0.49 (W) 
Kepler38  0.4644  0.43  0.39  0.43  0.37 (O) 
Kepler64 (K13)  0.642  0.65  0.53  0.58  (K13) 
Kepler64 (S)  0.634  0.58  0.52  0.53  0.57 (S) 
Kepler413  0.3553  0.31  0.26  0.35  (K14) 
The results regarding the critical semimajor axis can be found in Table 3, which, besides the values given by the two formulae just mentioned, also presents some relevant results by Welsh et al. (2012), Orosz et al. (2012a) and Schwamb et al. (2013).
Generally, our simulations gave us a larger critical semimajor axis than all previous results in all but one case (for Kepler413 the Mardling & Aarseth criterion leads to a critical semimajor axis greater than our value). However, this does not necessarily mean that our results are inconsistent with those of the other authors. The main reason for that is the definition of the critical semimajor axis: we found the semimajor axis for which the system was stable and immediately after that the system became unstable (although there may also be smaller values of the semimajor axis for which the system becomes stable again), while for example Holman and Wiegert (1999) looked for the smallest semimajor axis for which the system was stable. The definition of the critical axis in Welsh et al. (2012) is also similar to that. It is worth pointing out that for Kepler16 Jaime et al. (2012), using the invariant loops criterion of Pichardo et al. (2005, 2008), found a value of 0.63 AU. Again, the value we found for this system (0.67 AU) is greater than that.
Another factor that may explain part of the discrepancy between our value for the critical semimajor axis and the rest of the results is the definition of stability on which those results were based. It is a factor that may have a significant effect on any such study (e.g. see Georgakarakos 2008). For instance, the stability criterion in Welsh et al. (2012) would be expected to produce more conservative estimates. However, the fact that their integration time was only the one tenth of our simulation time could possibly cancel out that effect. Similarly, Holman’s and Wiegert’s integration time is only a fraction of ours, but as we said, in that case, it is the definition of the critical axis that explains the discrepancy between the results.
Finally, our results for Kepler16, Kepler34 and Kepler35 seem to be in agreeement with the results obtained from Popova & Shevchenko (2013).
3.3 Stability maps
In this subsection, we present the results from the third numerical experiment we conducted, which was the construction of stability maps for all planetary systems under investigation. In the way described in section 2, we found stable and unstable areas in the plane for each system and the results can be seen in Figures 28. The areas denoted by black dots show the unstable planetary orbits, while the white area represents the stable ones. The location of the nominal solution for each system is denoted by an empty black circle and also, indicated by straight grey lines, we mark the location of various mean motion resonances between the binary system and the planet. The resonances are of the type k:1, with k ranging from 4 to 14 for Kepler16, from 5 to 14 for Kepler34, from 4 to 14 for Kepler35, from 3 to 14 for Kepler38, from 4 to 14 for Kepler64 (both solutions) and from 3 to 14 for Kepler413. As we stated earlier, in this experiment, each point on the grid was integrated for years. Since we deal with systems that consist of two stars and a planet, a stability simulation may require a very long integration timescale, as sometimes instability signs may not become evident fast, but they become noticeable after a significant number of secular periods (e.g. see Georgakarakos 2013). However, for the six Kepler systems, the planetary secular periods of motion is of the order of a few decades and therefore our year integration time should be sufficient for investigating the stability of the systems.
Fig. 2 shows the stability map of the Kepler16 system. The planet seems to reside in a stable area, close to a region of unstable orbits associated with the 5:1 and 6:1 mean motion resonances. Our map also is consistent with the corresponding one of Popova & Shevchenko (2013). There are some areas in the map denoted by light grey dots. These correspond to unstable planetary orbits which were integrated for 1 Myr, ten times longer than the rest of the simulations. This was done in order to get an idea of how a longer integration time could possibly affect our stability results and maps. Apparently, only a limited area was added to the map, mainly for higher eccentricity values.
Figure 3 shows the stability map for Kepler34b. The planet seems to be far enough from the unstable area, lying between the 10:1 and 11:1 mean motion resonances. The unstable area in our map is larger than the corresponding plot of Popova and Shevchenko (2013), probably due to the longer integration time we used.
Regarding Kepler35, as seen from Fig. 4, the planet is located in the stable area on the plane, between the 6:1 and 7:1 resonances. As before, our map has a more extensive unstable area than the corresponding map of Popova and Shevchenko (2013), especially for intermediate and high eccentricities.
The Kepler38 stability map is shown in Fig. 5. The planet is clearly sitting in a stable area between two prongs of unstable planetary orbits which correspond to the 5:1 and 6:1 mean motion resonances.
For the Kepler64 system we have two figures, Figure 6, which is based on the solution of Kostov et al. (2013) and Figure 7,which is based on the solution of Schwamb et al. (2013). For both models, the nominal solution for the planetary orbit is very close to the unstable points associated with the 7:1 mean motion resonance. However, the stability maps have been based on an upper limit for the planetary mass. It is possible that a smaller planetary mass would result in a different stability map (e.g. Georgakarakos 2013). We would like to point out here that, as it is most likely, the planet did not form at its present location, which means that it possibly had to go through areas of instability during its migration phase. However, the crossing of the planet through the unstable areas during its inward migration phase would still be possible due to the facts that the migration timescale could be shorter than the instability timescale and that eccentricity damping forces can lower the planetary eccentricity during that phase. Of course, the previous remark is valid not only for Kepler64, but also for the rest of the systems we investigate here.
Finally, the last planetary system that we model here is Kepler413. Fig. 8 shows the stability map for the system. The orbit of the planet is lying in a stable area between the 6:1 and 7:1 mean motion resonances.
3.4 Checking if the nominal solutions are in resonance
As we saw earlier, all planets investigated here are close to a mean motion commensurability with the stellar binary and therefore we decided to check whether the systems were close to the exact resonance, using the definition as it appears in Murray and Dermott (1999) subsection 8.9.1. For that reason, we explored all possible resonant angles. For example, Kepler16b sits between the 5:1 and 6:1 MMR and hence, in this case, the potentially resonant angles are given by for the 5:1 MMR or for the 6:1 MMR, where is the mean longitude of the stellar binary, is the mean longitude of the planet, is the longitude of the pericentre of the stellar binary and is the longitude of the pericentre of the planet. The integer numbers and must satisfy the relationship for the 5:1 MMR and for the 6:1 MMR (Murray & Dermott 1999). We found that all the possible resonant angles of the systems were circulating at a high frequency and therefore not close to the exact resonance. The two most interesting exceptions were Kepler34 and Kepler35. Regarding Kepler34, although the planet is closer to the 10:1 MMR than to the 11:1, we checked both possibilities. We found that the nominal solution is close to the exact resonance associated with the 10:1 MMR. Fig. 9 shows the time evolution of the resonant angle . Both circulation and libration features can be noticed in that figure. Similarly, for Kepler35, which is between the 6:1 MMR and the 7:1 MMR, we found that the important angle was . Its time evolution can be seen in Fig. 10. Again, the plot exhibits both circulation and libration features.
4 Discussion and Summary
4.1 Effect of tidal evolution and general relativity
As it was stated earlier, the previous simulations were done under the assumption that the bodies were point masses. However, when the binary separation is of the order of a few stellar radii, tidal friction may become important, particularly if a system undergoes evolutionary phases of high eccentricity. Tidal dissipation can affect the orbital evolution of the stellar binary, which in turn, may have a significant effect on the planetary orbit. In addition, the orbit may also be affected by the nonspherical shape that the stars acquire due to the tidal effects and due to their rotation around their axes. Finally, general relativistic effects may also play a role in the stellar binary orbital evolution, as the extra frequency due to the precession of the longitude of the pericentre may have an effect on the orbital evolution.
From the systems we considered, Kepler34 experiences the closest pericentre distance between the stellar components in terms of their stellar diameter and for that reason it is chosen to see whether the above mentioned effects may become important in the evolution of the system (Kepler34 is the largest masscombined system with 2.0687 M and has the largest eccentricity, both important factors for GR). More specifically, and , which yield a pericentre distance of 0.10963 AU, which is approximately ten times the diameter of the largest star of the system. For comparison, the Kepler413 has and that give us a pericentre distance of 0.0978 AU, which is approximately 13.5 times the diameter of the largest star on the system. Additionally, the combined mass of the stars in Kepler413 is 1.3623 M
We integrated numerically the equations of motion (Eggleton et al. 1998, Prodan & Murray 2012) by using a BulirschStoer integrator (Press et al. 1996). The initial conditions used were: semimajor , eccentricity , longitude of pericentre , , , radius of the first star , radius of the second star , spin period for both stars 16 days, tidal Love number and tidal dissipation factor (the last two are adopted values for sunlike stars). The presence of the planet was neglected and the integration time was set to 9 Gyrs.
As expected, the results showed that the binary orbit tended to shrink and to get circularised, i.e. there was a clear decrease in the values of the semimajor axis and eccentricity. The eccenticity reduces by about after 8 Gyrs, while the semi major axis exhibits a drop of around in the same period of time. Finally, the argument of pericentre circulates with a period of around yrs. A graphical representation of some of the results can be found in Fig. 11. Since the above timescales are much longer than the secular ones, the evolution of the planetary orbit should not be significantly affected.
Hence, based on our results of Kepler34, we conclude that tidal friction, general relativistic effects and the deformation of the stars due to rotation are not very important in modelling the systems we are studying here.
4.2 Summary
We have investigated the stability of the six Kepler singleplanet circumbinary systems, i.e. Kepler16, Kepler34, Kepler35, Kepler38, Kepler64 and Kepler413 by integrating numerically the full equations of motion of each system.The orbits were assumed to be coplanar and our investigation was split into three experiments.
In the first experiment, we integrated the nominal solution of each system, as given in the corresponding papers, for 1Gyr, a timescale that was much longer than all previous work by other authors ( 10Myrs), except for the Kepler64 system (Schwamb et al. 2013), for which the authors used a similar timescale to ours. Our simulations showed that all systems were stable for the time that were integrated.
In the second numerical experiment we searched for the critical planetary semimajor axis (as defined earlier). The integration was done for years which was longer than any previous timescale for a similar experiment. Besides comparison of our results with the corresponding results of Welsh et al. (2012) and Orosz et al. (2012), we also compared our values for the critical semimajor axis with results based on the stability criteria of Holman & Wiegert (1999) and Mardling & Aarseth (2001). Generally, our simulations gave larger values compared to the rest of the results by other authors.
In the third experiment, we constructed stability maps on the plane. The simulations were done for years, which was at least 10 times longer that some similar past results. Moreover, this analysis was done for the first time for Kepler38. The simulations showed that all planets reside in stable areas between instability prongs along the location of mean motion resonances. However, the position of Kepler64b is very close to unstable orbits assosiated with the 7:1 mean motion resonance.
Additionally, we decided to use the data from the numerical integrations we performed for the nominal solutions in order to check whether the systems were close to the exact resonance and we found that Kepler34 and Kepler35 were the most interesting cases in that respect. This result may be related to the dynamical evolution of those two systems during their formation stage. In the context of planet formation around single stars, numerical simulations of the process of gas driven migration suggest that, as the gas in a protoplanetary disk is dispersed, giant planets can be found locked in a multiresonant configuration (Masset & Snellgrove 2001, Lee & Peale 2002, Morbidelli & Crida 2007, Morbidelli et al. 2007). Whether this resonancelocking occurs also in the gaseous protoplanetary disks around binary stars, has been studied by Nelson 2006 and Pierens & Nelson (2008, 2013). Our results suggest that, for some reason, Kepler34 and possibly Kepler35 demonstrate some features of that resonance trapping. The fact that we do not see similar features in the orbital motion of the rest of the systems we investigated, may be related to the further evolution of the systems due to the interaction of the planets with a remnant planetesimal disk, which can drive the planets out of the resonant condition, as is believed to have happened in the solar system (Malhotra 1995, for a review see Levison et al. 2007). This again may suggest a different evolutionary path for Kepler34 (and possibly Kepler35) compared to the rest of the systems in our sample, an interesting issue to be studied in future contributions.
The current reasearch therefore builds up on previous evidence that these systems are stable. Whether this scenario may change due to new discoveries such as more planetary companions or revised orbital solutions, only more observations will tell. Our results may be important for issues such as habitability, as the determination of habitable zones or looking for Earthlike planets in the Kepler systems under investigation can not ignore the presence of the already known gas giant planets. In addition, our results, and especially our stability diagrams, could also prove very helpful when investigating planet formation scenarios and especially late stage formation, when Lunarsized embryos develop into fully formed planets. As the stability limit for a planetary or an even smaller mass body in a circumbinary orbit appears to have a weak dependence on the mass of the planetary body (e.g. Georgakarakos 2013), our stability diagrams can provide useful information regarding where a planetary body may form under the gravitational perturbations of the six Kepler stellar binaries of this study.
5 Acknowledgements
We thank the referee for the useful comments and corrections. CC would lile to thank to FIME, UANL and Mexican SEP program PROMEP for their financial support on the development of this research. MRR acknowledges support from UNAMDGAPAPAPIIT project 115413 . HA thanks CONACyT Research Project 179662 and UNAMDGAPA project 108914.
Footnotes
 pagerange: A dynamical stability study of Kepler Circumbinary Planetary systems with one planet–References
 pubyear: 2014
 We use the ’Supplementary Table 1’ in Supplementary Information’ of Welsh et al. (2012)
 We use Tables 5, 6 and 7 of Orosz et al. (2012a) using ’Best Fit’
 We use Tables 2 and 3 of Kostov et al. (2013)
 We use Tables 7 and 8 of Scwamb et al. (2013)
 We use the ’Best Fit’ collumn in Tables 3 and 4 of Kostov et al. (2014)
 We use the ’Supplementary Table 1’ in Supplementary Information’ of Welsh et al. (2012)
 We use Tables 5 and 6 of Orosz et al. (2012a) using ’Best Fit’
 We use Tables 2 and 3 of Kostov et al. (2013)
 We use Tables 7 and 8 of Scwamb et al. (2013)
 We use the ’Best Fit’ collumn in Tables 3 and 4 of Kostov et al. (2014)
References

Prodan S., Murray N., 2012, ApJ, 747, id. 4, 15 pp.
Qian, S.B., Dai, Z.B., Liao, W.P. et al. 2009, ApJL, 706, L96
Qian S.B., Liao W.P., Zhu L.Y., Dai Z.B. 2010, ApJL, 708, L66
Qian S.B., DZhu L.Y., Dai Z.B., FernándezLajús E., Xiang F.Y., He J.J., 2012, ApJL, 745, L23
Qian S.B., Liu L., Zhu L.Y., FernándezLajús E., Baume G. L., 2012, MNRAS Let., 422, 24
Roell T., Neuhäuser R., Seifahrt A., Mugrauer M., 2012, A&A, 542, id. A92, 10 pp.
Schwamb M. E. et al., 2013, ApJ, 762, id. 127, 21 pp.
Soderhjelm S. 1984, A&A, 141, 232
Thébault P., Haghighipour N. 2014, preprint (arXiv:1406.1357)
Welsh W. F., Orosz J. A., Carter J. A. et al. 2012, Nature, 481, 475