Stability of planetary, single M dwarf, and binary star companions to Kepler detached eclipsing binaries and a possible five body system

The ‘Kepler Eclipsing Binary Star Catalog’ contains more than 2,000 eclipsing binary stars that have been observed during the Kepler mission (Prša et al., 2011; Slawson et al., 2011). The high precision observations from Kepler enable eclipse time studies to be performed where variations in the eclipse times of binary stars can be used to detect third bodies (e.g. Getley et al., 2017; Borkovits et al., 2016). Binary stars that have orbits aligned with the Earth will eclipse each other and detached and isolated binary stars should have eclipses that occur at predictable intervals. Plots of observed eclipse times (O) minus the calculated eclipse times (C), or O-C plots, may show variations from these predicted intervals. If these variations are also periodic it may be the result of a third body orbiting the binary stars (Beuermann et al., 2010).

When performing an eclipse timing study on the eclipsing binaries contained in the Kepler catalogue, several O-C diagrams were found where the values begin to decrease, or increase, and then suddenly and rapidly reverse direction and change sign, i.e., eclipses that occur earlier than expected change to later than expected, or vice versa. The O-C curves for these systems then rapidly reverse sign again, or "flip-flop" (see Fig.  1 for a visual example). The secondary eclipse O-C curve is out of phase with the primary eclipse O-C curve by a half orbital period. Examples of these "flip-flop" systems can be seen in Borkovits et al. (2016). Most of these systems also appear to have eclipse depth variations with differing magnitudes for each system. These systems all have similar reported eccentricities of the eclipsing binary as well as the highly eccentric orbit of the reported third body/outer companion. For the purposes of this paper, eclipsing binary is defined as the primary and secondary stars that eclipse each other i.e. producing the eclipses seen in the system O-C diagrams while third body/outer companion refers to one (or more) additional bodies orbiting the eclipsing binary.

The "flip-flop" features of the O-C diagrams and the high eccentricities raise the question of the dynamical stability of the systems and whether the systems with the reported configurations are stable. The dynamical stability of systems are important as outer companions in unstable orbits may result in the outer companion being ejected from the system within a short time period. However, stable orbits suggest the outer companion will remain within the system and are, therefore, more likely to exist as described and be observed (Horner et al., 2012a; Horner et al., 2012b). If an outer companion is stable for a range of configurations then the outer companion is more likely to exist as any detection errors won’t have a dramatic effect on the determination of the stability of the system.

The aims of this study are to: perform a dynamical stability analysis on the systems found with highly eccentric binary star orbits and extremely high eccentric outer companion orbits; report on the source of the "flip-flop" effect and the stability of the systems KIC5255552, KIC5653126, KIC5731312, KIC7670617, KIC7821010, KIC8023317, KIC10268809, KIC10296163, KIC11519226, KIC11558882 and KIC12356914 with the proposed third bodies; comment on the likelihood of these proposed third bodies existing; and comment on the likelihood of more of these "flip-flop" systems existing that continue to go undetected.

The Kepler data were used to produce O - C diagrams for detached eclipsing binaries to study eclipse timing variations. We created a C++ program, called bet or binary eclipse timings, to determine the mid-eclipse times of as many primary and secondary eclipses in the Kepler detached binary systems as possible (see Getley et al., 2017). bet is based on the software transit analysis package (Gazak et al., 2012) which uses the analytic formulae from Mandel & Agol (2002). The analytic formulae describe a system of two objects, using parameters including orbital period, radius ratio of the two objects, mid-eclipse time, orbital inclination and eccentricity, during various points throughout an orbit. The O-C diagrams of the Kepler systems shown in this paper were created using bet and found to contain rapid variations with the primary and secondary eclipse O-C curves out of phase.

rebound is an N-body integrator with Python and C implementations (Rein & Spiegel, 2015). Systems of bodies are able to be set up and integrated over time to estimate the orbital characteristics, such as semi-major axis and eccentricity, at various intervals. By simulating the positions and the evolution of the estimated orbital characteristics of a system over a long time period we can determine if the proposed system is in a stable orbit (allowing it to have been observed) or if it is in an unstable orbit and likely to eject one or more of the bodies. Eclipse times were obtained from the simulation and an O-C diagram produced to make sure that the distinctive characteristics of the actual O-C diagrams were present. The systems with these orbital characteristics were also integrated for $10^{6}$ years. These same systems were then integrated again 40 times for $10^{4}$ years with random values for the mean longitude, argument of pericentre and longitude of ascending node of the orbit of the outer companion and eclipsing binary. The purpose of the random values was to see if the third bodies were stable in this very specific configuration or if third bodies were stable for a range of configurations.

For the rebound models used, the value for the longitude of the ascending node for the eclipsing binary (i.e. $\Omega_{binary}$) was fixed to 90°. We found when it was fixed to 0° that although the "flip-flop" features of the O-C diagrams still occurred, the primary and secondary eclipse O-C curves were in phase rather than out of phase as seen within the real O-C diagrams. The value for the longitude of the ascending node for the outer companion (i.e. $\Omega_{companion}$) was set such that $\Omega_{companion}=\Omega_{binary}+\Delta\Omega$.

The individual masses for the primary and secondary star were calculated using the sum of the masses in Borkovits et al. (2016) and the temperature of the systems. Making the assumption that the primary star significantly dominates the temperature of the system we can search for the corresponding mass of a star at that temperature from Pecaut & Mamajek (2013)¹¹1 With additional details from http://www.pas.rochester.edu/~emamajek/EEM_dwarf_UBVIJHK_colors_Teff.txt . This becomes the estimate for the mass of the primary star. Using either a calculated mass ratio or a sum of masses of the primary and secondary star with the estimated primary star mass, calculating the estimated mass of the secondary star becomes trivial. Finally, we compare the J-H colour/magnitude difference of the system with the estimate for the primary star in order to perform a check on the assumption that the primary star significantly dominates the system. We tested this process for estimating masses against Kepler systems with known masses for the primary and secondary stars, Kepler-16 (Doyle et al., 2011), Kepler-34 and 35 (Welsh et al., 2012), Kepler-38 (Orosz et al., 2012b), Kepler-47 (Orosz et al., 2012a), all with at least one confirmed planet. Our estimates for the primary and secondary masses agree with the reported masses within $\sim$10% or less. We also tested against systems with no confirmed outer companions, KIC9851142 (Çakırlı, 2015) and KIC1571511 (Ofir et al., 2012), and found our mass estimates agreed with the reported masses within $\sim$20% or less.

The first systems to be selected for the dynamical stability study were KIC5255552, KIC5731312, KIC7670617, KIC10268809, KIC12356914 as these systems were identified as part of our own eclipse time study of the Kepler eclipsing binaries that had matching entries in Borkovits et al. (2016). These systems all contained a unique "flip-flop" feature or sudden period change in their O-C diagrams as seen in Fig.  1. The inferred properties of these systems were compared to see what all the systems had in common. The systems were found to have binary eccentricities ranging between $\sim 0.25$ and $\sim 0.42$ and third bodies with eccentricities of at least 0.385. The rest of the systems in Borkovits et al. (2016) were checked to see if there were any other systems that matched these criteria. Finally, the O-C diagrams of the systems were visually compared to find other possible candidates. The complete list of systems and their orbital properties can be found in Table  1 and Table  2. Two systems, KIC4055092 and KIC9715925, were found to match the selection criteria, however these systems were not a part of the dynamical stability study as both of these systems have mass estimates for the primary star that exceed the mass estimates for the total system. Another two systems, KIC6794131 and KIC7177553, were also possible candidates for the dynamical stability study, however accurate values for $m_{a+b}$ were not obtainable from Borkovits et al. (2016). As such, reliable models in rebound were unable to be made for these four systems and they were not included in the study.

The systems in Table  2 are listed separately due to the long period nature of the outer companions. These third bodies all have periods longer than the window of Kepler’s observations, and so, while models and fits can give us an indication of the properties and type of third bodies located within the systems, the margin of error in the values are likely too great to make firm conclusions. We can expect any estimate of the orbital period to be a lower limit due to the uncertainty involved in observing a system for less than one complete orbital period. The outer companion mass is likely to be an upper limit as lower masses are more detectable at longer periods (Watson & Marsh, 2010). For those systems with orbital periods less than the period of Kepler’s observations the values listed in Table  1 are likely to be accurate with a smaller margin of error.

The primary and secondary masses for the systems listed in Table 1 and Table 2 were calculated as described and are listed in Table 3. A number of the systems in Table  1 and Table  2 have outer companion masses that are almost as large, or even larger, than one or both of the stars in the binary system. If these third bodies significantly contribute to the flux of the system then the individual mass for the primary star would be larger than estimated and the mass for the secondary star would be lower (although the total mass of the binary system would be unaffected).

With the systems set up in rebound and integrated, plots are produced showing eccentricity vs time and semi-major axis vs time. By considering these plots we are able to view the evolution of the system over the defined period and determine whether any object is likely to be ejected from the system. For example, by considering the change in semi-major axis we can tell if an outer companion stays within the system or is moving further away from the binary stars and being ejected out of the system.

The light curves for the systems with inclinations of close to 90° were also visually inspected to look for any additional eclipsing events. Additional eclipsing events are a direct way of confirming the existence of additional bodies and may provide additional information about the characteristics and orbital properties of any additional bodies.

The Kepler "flip-flop" systems appear visually unique upon the first consideration of their O-C diagram (Fig.  1). The primary and secondary eclipses O-C variations are out of phase with each other, and there are sharp and rapid "flip-flops" indicating eclipses rapidly transitioning from earlier than expected to later than expected (or vice versa). An example of a simulated model’s O-C diagram can be seen in Fig.  2. The simulated O-C diagram shows the same out of phase and rapid variations that can be seen in the actual O-C diagrams from observed data.

The models from rebound allowed us to produce visual representations of the bodies and their orbits within the systems found in Table  1 and Table  2. By producing visual representations of the binary star orbits (Fig.  3), animating the binary star and outer companion orbits and the inclination evolution of the systems (Fig.  4) we were able to determine that all systems with the "flip-flop" O-C variations exhibit similar behaviour/orbital configurations as described in Section 2. The binary stars are locally bound together and both orbit and exhibit apsidal precession around the centre of mass of the entire system. The period of the eclipsing binary apsidal precession around the centre of mass appears to be the same as the orbital period of the outer companion, likely due to the dynamical interactions between the outer body, primary and secondary stars. The third bodies orbit the centre of mass opposite the binary stars. The orbits of the binary stars and the system as a whole are provided as animations available as additional supplementary material online. The models provide clarity on the orbits of the bodies within the system and explain the features seen in the O-C diagrams.

The results of integrating the systems for $10^{6}$ years can be seen in Figures  5(a) to  5(f). All of these systems were found to be stable over $10^{6}$ years. The eccentricities of the eclipsing binary combined with the high eccentricities of the outer companion do not appear to compromise the long-term stability of the systems. While the eccentricities of the objects in the systems varied over differing time-scales and by differing amounts, the semi-major axis remained relatively constant and, therefore, the outer companions remained within each system. As illustrated by Fig.  5(d) and Fig.  5(f), while the eccentricity of the binary stars can vary significantly, this did not necessarily translate to a major change in eccentricity of the outer companion or the semi-major axis of the system. The systems were also found to be stable for $10^{4}$ years when random values were used for the mean longitude, argument of pericentre and longitude of pericentre of the outer companion and the eclipsing binary. This increases the likelihood of the outer companions existing as slight changes or deviations from the proposed orbital properties still produced stable orbits.

The out of phase variations in the O-C diagrams for primary and secondary eclipses are likely the result of apsidal motion (Zasche et al., 2015). The apsidal motion and rapid eclipse time transitions are features that appear in all of the O-C diagrams of the models when an outer companion as described in Table  1 or Table  2 is present. The light curves of some of these systems also show significant eclipse depth variations. The eclipse depth variations are likely due to the dynamics of the system at play due to apsidal and nodal precession (Kane et al., 2012), and the evolution of the inclination in the system over time. Apsidal motion and nodal precession are illustrated in the simulated orbits in Fig.  3 while inclination evolution over time for a system can be seen in Fig.  4. Inclination evolution does not necessarily only change the depth of the eclipses seen but also whether we see the eclipses at all. For example, secondary eclipses for KIC11558882 are not initially seen in the light curve but begin to appear around 800 days (BJD - 2,454,833) and remain for the rest of the observing window (Fig.  6).

The Kepler mission viewed these systems for approximately 1400 days (Conroy et al., 2014), and it is fortunate that the observation period of Kepler coincided with the point in the outer companion’s orbit that results in the sudden "flip-flop" nature of the period changes. For third bodies that have orbital periods greater than 1400 days, part of the orbit will be unobserved and the "flip-flop" effect potentially missed. The greater the orbital period of the outer companion, the greater the chance of missing this dynamical effect in the observations. The sudden period changes are so rapid, some occurring over approximately 100 days, that even an orbital period of $\sim$1700 days could result in this system characteristic going undetected in the Kepler data.

The set of orbital properties within a system jointly influence the potential for transits or eclipses to be seen in the light curve. The probability of a transit occurring decreases as the orbital period increases  (Kane & von Braun, 2008) so while KIC10268809, for example, has inclinations that may indicate the possibility of transits (84° and 94° for the binary stars and outer companion respectively), the very long orbital period of the outer companion results in transits being unlikely to occur. Extra events can be seen in the light curve of KIC5255552, indicating that transits occur, and there are also additional eclipses that a third body may not account for, thus indicating the possibility of a quadruple system (Zhang, J. et al., 2018). None of the other systems considered in this study have definite or clear additional events occurring within the light curve, however it is possible KIC11519226 contains an additional eclipse (described in section 4.4). The equation for the probability of a third body transit/eclipse being seen from Earth is

where a is the semi-major axis, e is the eccentricity and $\omega$ is the longitude of periastron of the third body and and the orientation of the orbit of the third body is assumed to be random (Charbonneau et al., 2006). Using equation 1 and the mean radius of stars from Pecaut & Mamajek (2013) with the masses and other orbital characteristics in Table  1 we can calculate the probabilities of seeing transits from the systems with third bodies. We find that the probability of extra events occurring in KIC5255552, KIC8023317 and KIC11519226 to be less than 1% and that the extra events seen in KIC5255552 must be due to an extremely fortuitous occurrence.

In some systems, the sudden period flip in the O-C diagram may be the only indication of the presence of an outer companion. It is likely, given the large number of eclipsing binary stars observed with Kepler, that there are a number of systems that have been observed and classified as not containing an outer companion when in actuality the observations of Kepler haven’t been long enough to observe the effects of an outer companion. With only 11 systems displaying the "flip-flop" behaviour out of the more than 2000 Kepler eclipsing binary systems and almost half of the systems having an outer companion reported with greater than a $\sim$1400 day orbital period, it is likely that there are many more systems that have outer companions that remain undetected due to orbital configurations that did not result in notable O-C diagrams within the Kepler viewing window. The approximately 1400 day viewing window of Kepler will necessarily bias the detection results to systems that have outer companions with orbital periods of less than 1400 days. As Tables  1 and  2 contain a similar number of systems it is possible, if not likely, that the "flip-flop" characteristic seen in the O-C diagrams will exist in a wide range of systems that have already been observed but not during this "flip-flop" window.

All of the systems in Table  1 and Table  2 were integrated 40 times each for $10^{4}$ years with random initial values for the mean longitude, argument of pericentre and longitude of pericentre of the third body orbit and the eclipsing binary. While the random values can produce systems with O-C diagrams that vary significantly from the previously calculated values, the systems are still found to be stable. This exercise shows that even for a wide range of (though not necessarily all) orbital configurations systems with these mass and eccentricity values are likely to be stable.

In this study we used custom software bet based on transit analysis package or tap to perform an eclipse timing study on Kepler eclipsing binary stars. During the eclipse timing study we found systems that had O-C diagrams that displayed "flip-flop" or out of phase variations between the primary and secondary eclipse O-C curves and rapid period change variations. rebound was used to simulate these systems. The systems in Table  1 and Table  2 were chosen as they all exhibited a unique "flip-flop" effect within their O-C diagrams. Outer companions with the characteristics described all account for the features seen in the O-C diagrams such as the out of phase eclipse time variations and the "flip-flop" effect. With the systems simulated in rebound we then integrated the systems as described for $10^{6}$ years. We found that all systems were dynamically stable for at least $10^{6}$ years and, therefore, bodies in these orbital configurations are likely to be stable and observable. We also integrated these systems with random values for the mean longitude, argument of pericentre and longitude of pericentre of the third body and the binary star for $10^{4}$ years and found that the systems were stable for a wide range of orbital configurations. The evidence suggests the outer companions for the systems listed in Table  1 are an additional pair of stars in a binary configuration (KIC5255552, KIC5653126 and KIC11519226), a single M Dwarf star (KIC5731312 and KIC8023317) and a planet (KIC7821010).

We also suspect that a larger number of systems that have been observed would also show similar "flip-flop" characteristics if observed over longer or much longer time spans. However, due to the limits of the Kepler viewing window and large orbital periods estimated for the third bodies/outer companions the "flip-flop" effect continues to go undetected. As more and more systems are found with multiple bodies, the dynamical stability of the system as a whole is an important consideration when determining the likelihood of their existence. Of particular note is KIC7821010 which has a third body mass of $\sim$2.6 Jupiter masses. At $\sim$2.6 Jupiter masses it is well within planetary mass range and shows that even a relatively small mass can have large effects on the motion of its parent stars.

Other stand-out systems from this study include KIC5255552, where there are additional eclipses in the light curve (Zhang, J. et al., 2018) that may indicate the presence of a fourth star bound in a binary with the third star. A fifth body in the KIC5255552 system is a possibility and further observations of the system are crucial in determining the true nature of this system. While a triple star explanation cannot be ruled out for the systems KIC5653126 and KIC11519226, the photometric and dynamical analysis performed for this study suggests these systems are detached eclipsing binary stars with binary star companions.

Some of the systems presented, for example KIC11558882, cannot be reliably studied with ground based observations. The orbital period of the binary stars can be so great that observing eclipses to get meaningful data was only made possible with Kepler. Without space based observations these systems, and their O-C variations, may have continued to go undetected.

TESS (Ricker et al., 2014) is an all-sky survey of bright local stars with the ability of detecting planets with orbital periods of a few hours to a year or more. The launch of TESS provides more opportunities to locate comparable systems that are more local to the solar system and capable of follow up studies. With the launch of TESS and future projects we expect the number of systems that have similar characteristics to increase significantly.

This research has been supported by an Australian Government Research Training Program Scholarship. We would like to thank the referee and editor for their helpful feedback.

The data underlying this article will be shared on reasonable request to the corresponding author.