Dynamical analysis of the circumprimary planet
in the eccentric binary system HD 59686
Abstract
We present a detailed orbital and stability analysis of the HD 59686 binarystar planet system. HD 59686 is a singlelined moderately close (AU) eccentric () binary, where the primary is an evolved K giant with mass and the secondary is a star with a minimum mass of . Additionally, on the basis of precise radial velocity (RV) data a Jovian planet with a minimum mass of , orbiting the primary on a nearly circular Stype orbit with and AU, has recently been announced. We investigate large sets of orbital fits consistent with HD 59686’s radial velocity data by applying bootstrap and systematic gridsearch techniques coupled with selfconsistent dynamical fitting. We perform longterm dynamical integrations of these fits to constrain the permitted orbital configurations. We find that if the binary and the planet in this system have prograde and aligned coplanar orbits, there are narrow regions of stable orbital solutions locked in a secular apsidal alignment with the angle between the periapses, , librating about . We also test a large number of mutually inclined dynamical models in an attempt to constrain the threedimensional orbital architecture. We find that for nearly coplanar and retrograde orbits with mutual inclination , the system is fully stable for a large range of orbital solutions.
=1 \fullcollaborationNameThe Friends of AASTeX Collaboration
1 Introduction
The first Doppler surveys looking for extrasolar planets were focused on finding Solar system analogs and usually avoided known binary stars with semimajor axes 200 AU (see Eggenberger & Udry, 2010; Thebault & Haghighipour, 2014). As a result, the number of known binary systems with planets orbiting around one of the components (in circumstellar or Stype orbits), or orbiting around both stars (in circumbinary or Ptype orbits, see Dvorak, 1986) is still relatively low when compared to planets orbiting single stars. To date^{1}^{1}1http://www.univie.ac.at/adg/schwarz/multiple.html we know of 50 Stype planets which are part of wide binaries separated by at least 50–1000 AU (Roell et al., 2012) and 20 Ptype planets orbiting both stars where the binary separation is below 1 AU (mostly discovered with the Kepler satellite, Borucki et al., 2010; Doyle et al., 2011; Welsh et al., 2012; Leung & Lee, 2013; Kostov et al., 2013). However, only a handful of Stype planet candidates in moderately close binary systems ( 30 AU) are known in the literature and they all were discovered using the radial velocity (RV) method.
A famous example is the Cephei binary system, which consists of a K giant primary of = 1.6 and a secondary star with a minimum mass of = 0.44 , separated by 19 AU. This system has a Jovian planet with a minimum mass of 1.7 (Campbell et al., 1988; Hatzes et al., 2003) orbiting on a stable orbit around the primary star at 2.0 AU (Haghighipour, 2006). A planet candidate on an Stype orbit is also evident in the RV data taken for the HD 196885 binary system (Correia et al., 2008). This system consists of an F8V primary of = 1.3 and a secondary star with a minimum mass of = 0.45 , orbital semimajor axis = 21 AU and eccentricity = 0.42. The double Keplerian bestfit for the RV data of HD 196885 reveals an Stype planet around the primary with 2.6 AU, 0.48 and a minimum mass of 3.0 . Chauvin et al. (2011) have carried out dynamical simulations which show that the planet’s orbit is more stable in a highly inclined configuration near the equilibrium points of the LidovKozai regime (Lidov, 1962; Kozai, 1962). Later, Giuppone et al. (2012) have confirmed the stability of a highly inclined configuration with a mutual inclination of 43 or 137, but they also found stable nearly coplanar configurations, where the planet’s orbit is either prograde or retrograde, with the retrograde orbits being less chaotic.
Another remarkable example is the Octantis binary (Ramm et al., 2009; Ramm, 2015; Ramm et al., 2016). This system consists of a = 1.6 K1 III giant primary and a lowmass secondary star separated only by = 2.5 AU, with moderate eccentricity of = 0.24. The binary inclination is well constrained at = 71, which yields a secondary mass of = 0.6 . A lower amplitude periodic RV variation is present in addition to the secondary star RV signal, and if these variations are due to an orbiting planet, then the Stype companion would have 1.2 AU, 0.1 and a minimum mass of about 2.0 . The planetary interpretation is problematic because the bestfit orbit (with semimajor axis ratio / 0.47) is located well outside the boundary for stability if one assumes a coplanar and prograde planet with respect to the binary’s orbit (Holman & Wiegert, 1999). However, Eberle & Cuntz (2010) and Goździewski et al. (2013) have shown that a nearly coplanar retrograde orbit is stable, even though the stable region is small due to nearby meanmotion resonances (MMR) at the 2:1, 3:1, and 5:2 period ratios.
The existence of prograde, retrograde, or even LidovKozai resonance Stype giant planets as a part of moderately compact systems remains a very challenging dynamical problem. Apart from the longterm stability problem, it is puzzling how planets can grow through coreaccretion or diskinstability mechanisms in such close binaries. These systems provide important clues on how planets could form and remain in stable orbits around a star under the strong gravitational influence of a close stellar companion.
In this work we study the HD 59686 singlelined binary system, which is composed of a 1.92 Kgiant and a lowmass star with a minimum mass of = 0.53 . This system was reported to have a massive ( 7.0 ) Jovian Stype planet orbiting at = 1.09 AU around the primary star (Ortiz et al., 2016). The binary itself, however, is very eccentric ( = 0.73), which challenges the planet’s orbital stability. We carry out an extensive statistical and dynamical analysis to the available RV data to demonstrate that this system has stable configurations and to further constrain its orbital parameters.
This paper is organized as follows: In Section 2 we review the physical configuration of the HD 59686 system. Section 3 describes the methodology of our dynamical fitting and longterm stability analysis. In Section 4 we introduce the bestfit results from our tests and we reveal the possible Stype planet configurations. In Sections 5 and 6 we present dynamical and stability results around the best fits based on our bootstrap and systematic parameter grid search analysis. Finally, in Section 7 we present conclusions based on our results and discuss the possible Stype planet configurations.
2 The HD 59686 binaryplanet system
2.1 System configuration
HD 59686 (= HR 2877, HIP 36616) is a bright, photometrically stable (V = 5.45 mag., van Leeuwen, 2007) horizontal branch (HB) red giant star with an estimated mass of = 1.92 0.21 , radius of = 13.2 0.3 (Reffert et al., 2015) and metal abundance of [Fe/H] = 0.15 0.1 (Hekker & Meléndez, 2007). With luminosity = 73.3 3.3 and effective temperature = 4658 24 K, HD 59686 is a typical K2 III giant star. More physical parameters of HD 59686 can be found in Reffert et al. (2015).
Based on 88 precise (5 – 8 m s) RV observations of HD 59686 taken at Lick Observatory between November 1999 and December 2011, Ortiz et al. (2016) reported that HD 59686 is actually part of a much more complicated threebody system. The Keplerian orbital solution for HD 59686 given in Ortiz et al. (2016) shows that the RV data has a large amplitude variation of = 4014.12 5.84 m s, whose characteristic RV shape reveals a stellar companion with = 0.53 on a highly eccentric orbit ( = 0.73). The large eccentricity and the fact that the binary recently passed through its periastron ( February 2008) allowed the orbital period to be well determined as = 11679.94 192.92 days, even though the time span of the observations does not cover a full period of the binary orbit. In addition, the RV data yielded a loweramplitude signal of = 136.92 3.31 m s with a derived period of = 299.36 0.28 days. This signal is due to a Jovian planet with a minimum mass of = 6.92 on a nearly circular ( = 0.05 0.02) Stype orbit around the primary K giant star. The planetary signal remained coherent over many periods and was further confirmed by Trifonov et al. (2015) using followup RV measurements in the nearinfrared taken with ESO’s VLT spectrograph CRIRES (Kaeufl et al., 2004).
2.2 Constraints on the Planetary Companion from the Hipparcos Intermediate Astrometric Data
HD 59686 was a target of the Hipparcos mission (HIP 36616). We analyze the Hipparcos Intermediate Astrometric Data of HD 59686 based on the rereduction by van Leeuwen (2007) in the same way as described in Reffert & Quirrenbach (2011). We ignore the stellar companion, since its period is much longer than the Hipparcos mission duration, and fit only the astrometric orbit of the planetary companion to the abscissa residuals, simultaneously allowing for adjustments in the standard astrometric parameters (position, proper motion, parallax) and keeping the spectroscopic parameters fixed. The best fit occurs at an inclination = 2.9 and longitude of the ascending node = 266.7. The joint 3 confidence region extends from = 1.6 to 12.8 and from = 206.8 to 308.2. Thus, a significant parameter range in and can formally be rejected as a possible solution, and in Reffert & Quirrenbach (2011) we argued that this is the best indicator for an actual detection of the astrometric orbit.
However, there are several concerns with the Hipparcos data of HD 59686.
(1) The reduced value in the van Leeuwen version of the Hipparcos Catalog is 0.71, already quite small. It indicates that either the errors are overestimated, or that the solution is in fact already quite satisfactory, with no need for a better model for the measurements. The reduced value after fitting for the astrometric orbit is 0.59, which is uncomfortably small.
(2) There is a clear correlation between the time of the year and the scan direction (not only for HD 59686, but also for other targets). This is particularly a problem for periods close to one year, which is the case for the planet orbiting HD 59686 (best fit spectroscopic period 299 days). On top of that, many Hipparcos measurements have been obtained during the same season, i.e. with the same scan direction. Thus, the Hipparcos measurements for HD 59686 are poorly constrained in the perpendicular direction (roughly coinciding with the right ascension direction), and can float freely to fit any astrometric orbit.
We believe that, as a result, the Hipparcos data for HD 59686 should be treated with caution. In fact, we will show later on that any solution with for the inner companion is highly improbable, so the detection of an astrometric orbit could not be brought in line with the observed radial velocity data. We conclude that most likely the astrometric orbit of the inner companion has not been detected in the Hipparcos data.
2.3 Dynamical considerations
As discussed in Ortiz et al. (2016), the stellar companion of HD 59686 must be either a lowmass (and lowluminosity) star such as a K dwarf or a white dwarf remnant. These scenarios are particularly important to trace the possible origin of the Stype planet. However, the question of whether the secondary is a K dwarf or white dwarf is of little importance for the goal of this paper, which is to study the current permitted (stable) orbital configuration.
At first look, it is unclear how the planet could remain stable in such configuration. The binary semimajor axis is = 13.6 AU, but the pericenter distance is only = 3.67 AU, leading to strong interactions with the planet, which has = 1.09 AU and = 1.03 AU. Assuming a minimum mass of = 0.53 the Hill radius of the secondary star can be approximated as:
(1) 
which would cover the Stype planet orbit entirely during the binary periastron passage. Due to the large eccentricity of the binary, however, one can define the Hill radius at the pericenter distance instead of (see Hamilton & Burns, 1992), which leads to a smaller value of 1.66 AU. This suggests that the planetsecondary separation close to the binary periastron would be 1.5 , making the survival of the planet still challenging.
A quick check using the empirical stability criterion of Holman & Wiegert (1999) reveals that the critical (upper limit) semimajor axis for the Stype planet is 1.03 AU. Considering the binaryplanet orbital uncertainties, we find that the Stype planet is most likely unstable, with an orbit slightly outside the stability region. Using similar empirical stability criteria from Eggleton & Kiseleva (1995), we find that the planet is most likely stable, while the criterion of Mardling & Aarseth (2001) suggests that the planet is unstable. In any case, these stability criteria agree that the planet is close to the stability border. Therefore, in this paper we aim to inspect the threedimensional orbital architecture of the HD 59686 system and study its longterm stability and dynamics.
3 Methodology
Our orbital analysis for HD 59686 is based on the multidimensional body modeling scheme, which was previously applied to the 2:1 MMR exoplanet pairs around HD 82943 (Tan et al., 2013), HD 73526 (Wittenmyer et al., 2014) and Ceti (Trifonov et al., 2014). Briefly, we model the RV data using a LevenbergMarquardt (LM) minimization scheme, which performs an body fit by integrating the equations of motion using the GraggBulirschStoer integration method (see Press et al., 1992). The output parameters from our fitting code are the planetary and secondary star RV semiamplitude (), orbital period (), eccentricity (), argument of periastron (), mean anomaly (), inclination () relative to the sky plane, and ascending node (), as well as the RV offset (RV). All orbital parameters are the osculating ones in the Jacobi frame (e.g., Lee & Peale, 2003) at the first RV observational epoch, which is JD = 2451482.024. Each fit comes with a reduced value (), the residual value, and the 1 uncertainties of the adjusted parameters obtained from the covariance matrix.
For HD 59686, we adopt a stellar mass of 1.92 and a stellar velocity jitter amplitude of 20 m s. Ortiz et al. (2016) have shown that the Lick data of HD 59686 are consistent with additional stellar radial velocity jitter of about 20 m s. The most likely reason for the notable RV noise in early K giants like HD 59686 are solarlike mode oscillations (Barban et al., 2004; Zechmeister et al., 2008), which have typical periods much shorter than the typical time sampling of our Lick data, and thus appear as scatter. Using the stellar parameters for HD 59686 from Reffert et al. (2015) and the scaling relation from Kjeldsen & Bedding (2011), we estimated a jitter amplitude of 16.1 2.9 m s, which agrees well with the observed jitter of other K2 III giants in the Lick survey (Frink et al., 2001; Hekker et al., 2006; Trifonov et al., 2014; Reffert et al., 2015). Therefore, for our dynamical modeling we adopt a uniform a priori stellar jitter value of 20 m s, which we quadratically add^{2}^{2}2 Alternatively, the RV jitter could be fitted as a free parameter of the RV model (e.g. Baluev, 2009). The best doubleKeplerian fit of HD 59686 optimized with an additional jitter term to the Lick data yields a jitter value of m s, which is consistent with the uniform jitter value of 20 m s adopted in this work. to the total RV data error budget.
All dynamical fits in our study are further tested for longterm dynamical stability. We integrate the orbits using a custom version of the WisdomHolman algorithm (Wisdom & Holman, 1991), modified to handle the evolution of hierarchical systems consisting of massive bodies (Lee & Peale, 2003). The bodies are assumed to be point masses, and mutual collisions between them are not considered in defining system stability. We also neglect General Relativity and companionstar tidal effects during the simulations. We integrate the individual fits for a maximum of 10 Myr, by adopting an integration time step equal to 1 day. Our integration setup corresponds to more than full binary orbits with about 300 steps per complete planetary orbit. We find that this setup is sufficient to resolve the planet’s orbit with high resolution and study the system’s longterm stability.
We define the HD 59686 system as stable if during the integration the companion bodies remain in orbits which do not deviate significantly from their initial bestfit configuration. The system’s stability depends primarily on the survival of the Stype planet. In most cases when the planet inclination is , the planet has relatively low mass to perturb the binary orbit significantly, and can well be approximated as a test particle in a twobody system. However, we also test fits with , where the mass of the Stype body becomes quite large as gets smaller, and thus can significantly influence the binary orbit during the orbital evolution. A simulation is terminated and the system is considered as unstable if at some point of the integration the semimajor axis or changes by more than 60% from their initial values, or if . Particularly, when and AU, the planet periastron distance to the central star, , is well within the physical radius of the K giant ( 0.06 AU), and the planet would collide with the star. Although these criteria for instability are somewhat arbitrary, our simulations show that even small chaotic deviations in and quickly accumulate, and there are no cases where the orbits change significantly without exceeding these criteria.
4.0cm



Coplanar edgeon prograde  


Parameter  HD 59686 Ab  HD 59686 B 


[m s]  137.0 (3.4)  4012.6 (20.6) 
[days]  299.1 (0.3)  11696.4 (196.4) 
0.05 (0.02)  0.730 (0.003)  
[deg]  121.1 (28.7)  149.4 (0.1) 
[deg]  299.5 (28.4)  259.2 (1.7) 
RV [m s]  248.6 (12.5)  
[deg]  90.0^{a}^{a}Fixed parameters.  90.0^{a}^{a}Fixed parameters. 
[deg]  0.0^{a}^{a}Fixed parameters.  0.0^{a}^{a}Fixed parameters. 
[deg]  0.0  
[AU]  1.089  13.611 
[]  6.97  558.41 
[m s]  19.59  
76.61  
0.995  


Coplanar edgeon retrograde  


Parameter  HD 59686 Ab  HD 59686 B 
[m s]  136.7 (3.3)  4013.7 (20.5) 
[days]  299.0 (0.3)  11669.3 (194.7) 
0.05 (0.02)  0.729 (0.003)  
[deg]  126.8 (28.3)  149.4 (0.1) 
[deg]  293.5 (28.4)  258.9 (1.7) 
RV [m s]  247.1 (12.4)  
[deg]  90.0^{a}^{a}Fixed parameters.  90.0^{a}^{a}Fixed parameters. 
[deg]  180.0^{a}^{a}Fixed parameters.  0.0^{a}^{a}Fixed parameters. 
[deg]  180.0  
[AU]  1.089  13.591 
[]  6.96  558.46 
[m s]  19.46  
75.63  
0.982  


Mutually inclined  


Parameter  HD 59686 Ab  HD 59686 B 
[m s]  130.1 (26.2)  4020.0 (182.1) 
[days]  300.5 (0.5)  11398.3 (1244.3) 
0.08 (0.02)  0.725 (0.003)  
[deg]  145.2 (19.7)  149.8 (4.7) 
[deg]  280.4 (20.0)  256.4 (21.0) 
RV [m s]  239.5 (12.7)  
[deg]  178.8 (0.3)  86.4 (3.3) 
[deg]  316.5 (10.3)  0.0^{a}^{a}Fixed parameters. 
[deg]  92.73  
[AU]  1.15  14.06 
[]  359.22  618.36 
[m s]  16.62  
55.50  
0.750  

Note. – bootstrap uncertainties, (0.0) covariance matrix uncertainties.
4 Best fits
4.1 Edgeon prograde and retrograde fits
The best coplanar and edgeon dynamical fit is generally consistent with the Keplerian fit shown in Ortiz et al. (2016). Our dynamical fit is close to a double Keplerian, since any significant gravitational perturbations on the planetary orbit (and thus on the induced RVs) are expected to be detected only after a few binary cycles, while the RV data currently cover only 40% of one full binary orbit. We first keep the orbital inclinations fixed at = = 90 and the difference between the lines of node = – = 0, which defines a planar and prograde configuration. The best fit in this orbital configuration has = 0.995 and leads to orbital elements of = 299.1 0.30 days, = 0.05 0.02, = 1.09 AU, and mass of = 6.97 for the planet, and orbital elements of = 11696.4 196.4 days, = 0.73 0.003, = 13.61 AU, and secondary star mass of = 558 for the binary. The full set of orbital elements and their bootstrap (see §5) and covariance matrix estimated uncertainties are given in Table 1, while the actual fit to the data (black curve in the upper panel) and its residuals are illustrated in Figure 1. The longterm evolution of the orbital semimajor axes and eccentricities are shown in Figure 3 (left panel). According to our stability criteria, this fit is stable only for about 42 kyr, before the planet collides with the star. A close examination of this fit indicates that the orbits show large variations in and small, but chaotic variations in . Eventually the secondary companion excites the planet eccentricity above , which interrupts our integration.
Since the best coplanar and prograde fit is unstable, we test how the fit quality and stability change if we allow noncoplanar orbits. We simplify this test by keeping the binary on an edgeon orbit with fixed = 90 and = 0. For the planet we also fix the inclination at = 90, but we systematically vary between 0 and 359 with a step of 1. Thus, in this test we keep the companion masses at their minimum, while the mutual inclination comes only from the difference between the longitudes of the ascending nodes = – following the expression:
(2) 
Figure 2 shows the results from this test. We plot the quality of the mutually inclined fit in terms of () as a function of (). With horizontal dashed lines are shown the 1, 2 and 3 confidence levels according to . The bestfit in Figure 2 appears at = 180, which is again a coplanar, but retrograde planet orbit. The best coplanar prograde fit has = 76.61, while the best coplanar retrograde fit has = 75.63, resulting in = 0.98. This difference is slightly below the 1 limit, and thus the retrograde fit does not represent a significant improvement to our model. In Figure 2 most of the edgeon fits with 145 are above 1 from the best fit and are unstable (red dots), while all fits with between 145 and 180 are within 1 and are stable (blue thick line) for at least 10 Myr.
Figure 3 (middle panel) shows a 50 kyr time span of the orbital evolution for the best coplanar retrograde fit. The semimajor axes and are nearly constant during the stability test. The planet eccentricity oscillates with a large amplitude between 0 and 0.35, but the system remains stable, with the bodies well separated from each other. Interestingly, the mean period ratio of this stable retrograde fit is 39, but the system is not in 39:1 MMR, as none of the resonance angles associated with the 39:1 MMR are librating. For the MMR, the resonance angles are
(3) 
where is positive for prograde motion and negative for retrograde motion, are the longitudes of periastron and are the mean longitudes. All fits with between 145 and 180 have similar behavior for , , , and , while oscillates with small amplitude around the initial fitted value. None of them seems to be locked in a MMR.
The most likely reason for the wider stable region for the retrograde orbits is that the individual MMR are of higher order for retrograde than prograde orbits. As demonstrated in Morais & Giuppone (2012), at the same meanmotion ratio, the MMR is of order for prograde versus order for retrograde orbits, which for the latter results in much narrower MMR libration widths and thus smaller phasespace overlap of neighboring MMR where the planet would be likely unstable. The findings of Morais & Giuppone (2012), however, were restricted to the dynamics of Stype planets in circular binary systems, while the dynamics of the planet in the highly eccentric HD 59686 binary is far more complex. A more detailed analysis of resonance width and overlap for prograde and retrograde orbits using the formalism developed by Mardling (2008) in the context of the HD 59686 system will be presented in a future paper (Wong & Lee, in preparation).
4.2 Inclined coplanar fits – constraining
Both prograde and retrograde edgeon bestfits suggest a coplanar configuration. Therefore, as a next step we test how the fit quality and stability for both configurations depends on the inclination (measured from the plane of the sky). For the prograde geometry we fix = 0 and = . We systematically vary and from 90 to 5 with a decreasing step of 1, and thus we gradually increase the companion masses by a factor of approximately . The same test is done for the retrograde fits, which are constructed by keeping = 180, = 180 . and varying from 90 to 5 with a step of 1. According to Equation (2), the prograde fits have a mutual inclination of = 0 and the retrograde fits = 180.
The results from this systematic test are illustrated in Figure 4, which shows a comparison between prograde and retrograde dynamical fits as a function of . The minimum for both prograde and retrograde cases is at = 1, which corresponds to the same best fits presented in Table 1 and Figure 2. The confidence levels (1, 2 and 3) in Figure 4 are measured from the best retrograde fit, and represent the same confidence levels as shown in Figure 2. Clearly, our body fits can only weakly constrain the orbital inclination from the RV data. Overall, the retrograde fits have better values, but in both configurations the fits are gradually becoming worse for lower , and thus higher planet and secondary star masses. For retrograde fits down to 2 the orbital inclination can be between 10 and 90, while for prograde fits the inclination can be between and 90. In both configurations, however, the inclination is unlikely less than , as the secondary star would then be at least a Gtype mainsequence star with about twice the minimum mass. Such a stellar companion should have been detected by Ortiz et al. (2016) via LBT angular differential imaging, but since it was not, we assume = 0.5 as a lower limit.
When it comes to stability, all retrograde coplanar fits in Figure 4 are stable for 10 Myr, including those at very low inclinations, while none of the prograde fits is longterm stable.
4.3 Mutually inclined fits – the global minimum
Finally, in our dynamical modeling of HD 59686’s RV data, we allow nonedgeon mutually inclined orbits by fitting independently and in the range between 0 and 180 and between 0 and 360. In this way, we allow our fits to adopt a large range of companion masses and we cover all possible orbital alignments.
Our mutually inclined bestfit has a strong minimum with = 0.75, yielding for the inner companion = 300.5 0.5 days, = 0.08 0.02, = 1.15 AU, and for the outer companion = 11398.3 1244.3 days, = 0.725 0.003, = 14.06 AU (see Table 1). Remarkably, this fit suggests that the inner companion has nearly faceon orbit with well constrained = 178.8 0.3. This means that the inner companion is no longer a planet but a stellar mass companion with = 359 () forming an inner binary pair with the K giant. The outer companion has nearly edgeon orbit with = 86.4 3.3 and mass = 618 (). The difference between the ascending nodes is = 316.5 10.3, and according to Equation (2) this leads to = 92.7 3.3. We achieve practically the same fit (within errors) with = 0.75 at = 1.15 0.6, = 93.4 2.5 and = 43.7 12.3, which is a mirror image of the above orbital configuration.
The value for this mutually inclined fit is 55.5, which is much lower than the best coplanar retrograde fit with = 75.6. This fit, however, has three additional fitting parameters compared to the edgeon coplanar fits, which must be taken into account when testing for significance. Following the approach, we assume that , and are systematically adjusted, while the rest of the orbital parameters are fitted by our body model. The confidence intervals in this case obey the distribution for 3 degrees of freedom. The difference between the fits is = = 20.1, suggesting that the best coplanar retrograde fit is between 3 and 4 worse than the mutually inclined fit, and thus the latter represents a significant model improvement.
Since the coplanar model with = 11 fitting parameters is “nested” within the mutually inclined model with = 14 parameters, another way to test the significance is the use of the test and determine the value following (Bevington & Robinson, 2003):
(4) 
where = is the number of additional parameters being tested, = is the number of degrees of freedom for the best mutually inclined model, with the number of data points. For , the probability for model improvement is = 0.000039, which is much lower than our adopted cutoff value of = 0.01, meaning that the null hypothesis is successfully rejected. Thus we conclude that the mutually inclined fit is indeed better when compared to the coplanar edgeon model.
This significant model improvement is intriguing and deserves a closer look. First, it should be noted that the and tests only work well for Gaussian errors and models that are linear in the parameters (or could be linearized in the uncertainty region of the parameters due to large enough sample size) (Press et al., 1992), which we do not have when we apply an body dynamical fit to the existing RV data for HD 59686.
Dynamical fitting of RV data consistent with two or more companions can be in principle sensitive to the true companion masses, but this has been proven to be very challenging even for the most extensively studied multiplanet systems (see Bean & Seifahrt, 2009; Correia et al., 2010; Nelson et al., 2016). The critical requirements to measure successfully mutual inclinations are: (1) high RV precision, (2) low velocity jitter, typically on the order of at most a few m s, (3) large set of RV data covering many orbital cycles, and (4) the signal discrepancy between the minimum mass coplanar fit and the mutually inclined fit must be larger than the RV noise. In this context, we note that the available Lick data for HD 59686 do not satisfy these criteria, with only 88 RVs (with precision of 5–9 m s) distributed over 11 years covering only 40 % of the outer binary orbit, and RV jitter of 20 m s. Thus it is unlikely that we would be able to tightly constrain the true companion masses (through ).
In the bottom panel of Figure 1, the RV residuals to the coplanar prograde fit are shown, and overplotted (red curve) is the difference of the mutually inclined model from the coplanar prograde one. Clearly, some of the outliers present in the coplanar prograde case are well modeled by the mutually inclined model with three additional fitting parameters. These outliers, however, lie in the relatively sparsely sampled epochs between JD = 2454200 and 2455200, which unfortunately coincides with the outer binary periastron passage, when the RV signal changes rapidly. Perhaps, the signal would be validated if we had more RV data following the mutually inclined fit prediction at that orbital phase, but currently it is fair to conclude that we could be fitting noise rather than a true signal.
In addition, this mutually inclined best fit is extremely unstable. The inner companion has a nearly circular and highly inclined (polarlike) orbit with respect to the outer binary plane. Such orbits are potentially unstable due to the LidovKozai effect (Lidov, 1962; Kozai, 1962), which leads to periodic exchanges between and . Following the analytic expression given in Takeda & Rasio (2005), we can estimate the inner companion maximum eccentricity that can be reached through the LidovKozai cycle:
(5) 
With = 92.7 obtained from the best fit, 0.998, which exceeds our stability criterion. As expected, our direct numerical integration (right panel of Figure 3) shows that the best mutually inclined fit does not survive even 600 years. The inner companion is quickly excited to and it eventually collides with the K giant. Obviously, such a highly inclined stellar triple is very unstable, which seems to be a good argument against this solution.
5 Bootstrap statistics
Our goal in this section is to obtain parameter estimates and confidence regions based on the empirical distribution of constructed orbital parameters using bootstrap resampling. We analyze the distribution of the adjustable fitting parameters around the bestfit by randomly drawing RV data points with replacement (e.g., Efron, 1979; Press et al., 1992; Tan et al., 2013) and perform a dynamical fit to each RV data set obtained in this way.
We create a total of = 5000 bootstrapped data sets and we fit each sample with strictly coplanar edgeon (prograde and retrograde) and mutually inclined configurations. All fits to the bootstrapped data sets are integrated for a maximum of 10 Myr and their stability is examined. The 1 confidence levels from the distributions are used to estimate the asymmetrical bestfit parameter uncertainties. These estimates are listed in Table 1, along with the bestfit covariance matrix errors. For the mutually inclined fits, a bootstrapped sample is rejected if the fit suggests or (i.e., ). This ensures that our samples are consistent with a reasonable secondary companion mass in the range to about 1.1, and thus are consistent with the observational constraints given in Ortiz et al. (2016). The total fraction of fits with highly inclined () outer companion is 6% of the total constructed bootstrapped RV data sets. In addition, 1.5 % of the fits are unable to converge when we try to fit a mutually inclined configuration, meaning that in order to get 5000 fits we have to create a total of 5400 bootstrapped RV data sets.
Figure 5 shows the results from the bootstrap analysis for the edgeon () prograde () and retrograde () configurations. In each panel we illustrate the distribution of planet versus binary orbital elements (, , , , ), semimajor axes , and dynamical masses . Solid contours show the 1, 2 and 3 (68.27%, 95.45%, 99.73%) confidence levels from the twodimensional parameter distributions. In all panels, the green dot represents the position of the bestfit values from the prograde or retrograde dynamical fit to the original data set, while the green error bars are the estimated uncertainties from the covariance matrix (see Table 1). With blue dots we mark the unstable configurations, while with red points we show the configurations which survive for 10 Myr.
Clearly, the distributions of orbital elements for the prograde and retrograde configurations are very similar. For both configurations, the covariance matrix errors estimated from the original data and the 1 (68.3%) confidence region from the bootstrap analysis are roughly consistent with each other. The main difference between the configurations comes from the stability results. We find that 97% of the retrograde configurations are stable for 10 Myr, while for the prograde case this number is only 4%. We note, however, that the stable prograde configurations fall mostly within the 1 confidence region and from their distribution we can identify potentially stable regions of the parameter space. As can be seen in Figure 5, the stable fits seem to cluster around 145, 257 and 259, 0.05 to 0.07, and a few discrete values in which avoid initial integer period ratios of 39 and 40. This result implies that a small, but robust set of stable prograde fits does exist. Therefore, we cannot eliminate the possibility that the HD 59686 system is in a stable prograde configuration.
The results of the bootstrap analysis for the mutually inclined configuration are shown in Figure 6. In particular, we aim to quantitatively estimate the inclination distribution () and see how often polarlike Stype companion orbits will occur in the resampled data. Therefore, in Figure 6 we introduce three additional panels:  ,  and a comparison between the mutually inclined and prograde coplanar values. The last of these panels also shows the test probability for significant improvement when three additional parameters are added. We find that 86.1% of the fits lead to significant improvement, applying our chosen threshold of = 0.01.
The orbital parameter distributions are wider than in the coplanar cases but seem to agree with the bestfit errors. The binary orbit covariance matrix errors are sometimes larger than the bootstrap 1 contours, but in general the distribution peak is consistent with the bestfit estimate. The inner companion inclination is found mostly at lower values, driving towards brown dwarf and star like masses. On the other hand, within 3, stays between 75 and 105 (i.e., between 1 and 0.96), leading to similar to those of the edgeon cases. The distribution of clusters around the bestfit and is well constrained around 90, leading to nearly perpendicular triplestar orbits. We suspect, however, that the distribution in might be the result of a model degeneracy. By randomly scrambling the data with repetition, the possibly problematic outliers with sparse cadence are not always removed (see Figure 1 and discussion in Section §4.3). Even worse, for some bootstrapped data sets, we repeat these points, while removing data points with lower residuals. We find that none of the mutually inclined fits based on bootstrap analysis and shown in Figure 6 is longterm stable. The average bootstrap survival time is only a few hundred years before the inner binary pair collides, which is consistent with the orbital evolution of the best mutually inclined fit to the original data.
6 Parameter grid search
A detailed picture of the dynamical properties of the HD 59686 system can also be assessed using the parameter gridsearch technique (e.g., Lee et al., 2006). We systematically vary a pair of orbital parameters, which we then keep fixed while the rest of the parameters in the model are adjusted to minimize . This method allows us to systematically inspect the multidimensional parameter space around the best fit and study the properties and longterm stability of nearby fits.
6.1 Coplanar prograde grids
Our edgeon prograde bootstrap analysis reveals that stable fits are clustered around initial , = 0.05 to 0.07, and / 39 or 40 (within 1 from the best fit). These three orbital parameters are also the least constrained parameters from our fitting, especially when compared to , and . Therefore, any grid combination including , and yields an adequate search for prograde coplanar stable fits. Figure 7 shows stability results for an edgeon coplanar grid in the /  space, where we fix and at their bestfit values of 299.1 days and 149.4, respectively. Since from the bootstrap analysis we know that stability appears when 0 we keep = fixed, leading to an initially aligned configuration with = 0. In Figure 7, stability within the 1 confidence region from the best fit is achieved when 0.04 to 0.07 and / 39 or 40, which confirms the results from the bootstrap analysis.
As a next step, we test for stability in eight  grids, each with = 0.02 to 0.09 in steps of 0.01. and are varied around the best coplanar fit in the range of 11000 to 12500 days and 60 to 210, respectively, while the remaining parameters in the dynamical model are adjusted. Figure 8 shows the survival time resulting from this test together with the confidence levels. Since these grids are constructed with three systematically varied parameters (, and ), we consider Figure 8 as a threedimensional parameter cube, where each grid is a separate  slice placed on a lowerresolution “z”axis constructed for different . Thus, the significance levels (black contours) shown in Figure 8 are calculated for three degrees of freedom. Clearly, most of the fits are highly unstable, except those located between integer period ratios, with nearly aligned orbits () and between 0.05 and 0.07. We find that these stable islands cover the largest area when = 0.06.
Knowing that and are critical stability parameters, in Figure 9 we show results for  grids, where we fix = 0.06 and = = 149.4, and we systematically adopt inclination of = 90, 75, 60, 45 and 30. In this way we study the stability of fits in the  parameter space for initially aligned orbits and increasing companion masses. In these grids, the stable regions now extend through a large range of and cross inside the confidence regions. Configurations with period ratios near an integer initially are highly unstable, while stable configurations can be found between initial integer period ratios. The stable regions are evident in the = 90 to 60 grids, but they become smaller with decreasing or increasing masses. The = 45 grid contains only a few marginally stable fits, which are likely unstable beyond 10 Myr, while all configurations for the = 30 grid (not shown in Figure 9) are unstable. These results identify a lower limit for the inclination of 45 for stable prograde coplanar configurations. This inclination limit happens to coincide with the secondary star mass constraints discussed in Ortiz et al. (2016).
Table 2 gives the orbital parameters and corresponding errors for the best stable fit among these grids, which has an initial / 39.3. Figure 10 shows the orbital evolution of this stable fit. The evolution of the semimajor axes and in the upper left panel shows that this configuration remains longterm stable with well separated orbits. The binary eccentricity has very small amplitude variations around 0.73, while the planet remains nearly circular, with varying between 0.04 and 0.11 (upper right panel). The bottom two panels of Fig. 10 show the evolution of the secular apsidal angle = , which exhibits a clear libration around 0 with a semiamplitude of 37, while the mean period ratio during the integration is 39.4, close to the initial /. We examine this configuration for librating resonance angles (see Eq. 3) associated with the nearest 39:1 and 40:1 meanmotion commensurabilities, and confirm that this stable prograde fit is not involved in a MMR. Such an orbital evolution is characteristic for secular apsidal alignment (e.g. Lee & Peale, 2003; Michtchenko & Malhotra, 2004). The libration of around 0 is critical for the stability of our system, since it helps the lowermass Stype object to retain small eccentricities while being significantly perturbed by the secondary star. We have investigated all stable prograde islands shown in Figures 7, 8 and 9, and they all exhibit similar evolution, with librating around 0, circulating MMR angles , small and a noninteger mean period ratio /. Thus, we conclude that if the HD 59686 system is indeed prograde, then it must be locked in secular apsidal alignment to stabilize the orbits.
Two additional remarks on the secular apsidal alignment and nonMMR nature of the stable islands in Figures 7, 8 and 9 are in order. First, one may be concerned that the noninteger initial and mean period ratios may not represent the true period ratio due to the large mass of the secondary star. However, the Hamiltonian in Jacobi coordinates in Equation (11) of Lee & Peale (2003) shows that the perturbations to the Keplerian motions from the interactions between the secondary star and the planet remain small and the semimajor axes and period ratio should be nearly constant, throughout most of the binary orbit (including the initial epoch when the secondary is AU from the primary). Even when the secondary is at periastron, we can estimate from the lowest order term in the perturbations to the Kepler motions in Equation (11) of Lee & Peale (2003) that the full amplitude of the variation in the period ratio should be if is small. These results are consistent with the evolution of shown in Figure 10, where is near the initial value most of the time and shows scatter of .
Second, Wong & Lee (in preparation) have systematically studied the stability of circumprimary planetary orbits in the HD 59686 system, with initial conditions in grids of and for several values of and mean anomalies. For the coplanar prograde case, they confirm the existence of islands stabilized by secular apsidal alignment. They also find islands that are stabilized by MMR, but these are at higher and do not fit the observed planet.
6.2 Coplanar retrograde grids
We repeat the grid analysis for retrograde coplanar configurations by fixing = 180 for each fit. We find that all fits within 3 from the best fit are stable for at least 10 Myr. Particularly, all fits in the  grids are stable despite the large companion masses for = 30 and even = 15. The stability of the  grids in the range = 90 to 15 is in agreement with the results presented in Figure 4 and Section §4.2, where all the retrograde coplanar inclined best fits are stable down to = 5. We conclude that the best retrograde coplanar fit is well within a large stable phase space region, not necessary involved in a MMR. Therefore, no meaningful stability constraints can be obtained from the retrograde coplanar grids, except that the retrograde orbits yield very strong candidates for the HD 59686 system configuration.
6.3 Mutually inclined grids
We construct twelve separate  grids, where for each grid we adopt to 330, with a step size of 30. In these grids and are varied from 5 to 175, with a step size of , corresponding to 50 different values. Fits with , below 5 and above 175 are not considered, since these inclinations lead to nearly faceon orbits and the dynamical masses of the companions will be very large. In fact, as discussed in Section §4.2, other observations constrain to values greater than 0.5, but since the dynamical model does not reject more massive secondary stellar companion, we construct symmetrical grids with the same range of and . Thus, we study mutually inclined configurations covering almost all possible system geometries. The results from this test are shown in Fig. 11.
5cm



Parameter  HD 59686 Ab  HD 59686 B 


[m s]  137.0 3.3  4011.6 3.7 
[days]  299.2^{a}^{a}Fixed parameters.  11772.7^{a}^{a}Fixed parameters. 
0.06^{a}^{a}Fixed parameters.  0.731 0.005  
[deg]  149.3 ^{a}^{a}Fixed parameters.  149.4 0.1 
[deg]  272.0 1.4  259.8 0.4 
RV [m s]  243.2 4.0  
[deg]  90.0^{a}^{a}Fixed parameters.  90.0^{a}^{a}Fixed parameters. 
[deg]  0.0^{a}^{a}Fixed parameters.  0.0^{a}^{a}Fixed parameters. 
[deg]  0.0  
[AU]  1.09  13.67 
[]  6.97  558.5 
[m s]  19.84  
78.57  
0.970  



The only stable region we find in these grids corresponds to nearly coplanar retrograde geometries with 145. The orbital evolution of these stable fits is very similar to that of the best retrograde coplanar fit discussed in §4.1 and shown in the middle panel of Fig. 3. For nearly coplanar retrograde orbits, the amplitude of the variations in the inner companion eccentricity is rather large. It is about 0.35 for strictly coplanar orbits and increases with increasing mass of the secondary star (i.e., with decreasing ) to about 0.42 at maximum, when and 145. Meanwhile, the and oscillation frequency is highest in the coplanar configuration and decreases with decreasing . For fits near the stability boundary of 145, another shorter term eccentricity variation is visible on top of the main secular eccentricity oscillation, which has a period of a few hundred years. In this stable region, also exhibits small variations around the initial best fit value, except for the coplanar case, where remains constant at 180. We do not find evidence that any of the stable retrograde fits are in MMR.
The small retrograde corners of Fig. 11 are intriguingly stable, suggesting that the inner companion’s dynamical mass could be as high as 10 times its minimum mass (i.e. 0.1), converting the planet to a massive brown dwarf or at the extreme even a very lowmass Mdwarf star. A massive, highly inclined and retrograde inner companion may be consistent with the Hipparcos astrometry and dynamical modeling, but to preserve the system’s stability the outer binary companion must also have a small and that is not supported by observations. On the one hand, if we assume that the outer companion is a main sequence star, then we can limit its inclination to (marked by a blue dashed lines in Fig. 11), beyond which it would be a solar mass star and should have been detected in Ortiz et al. (2016). All stable solutions with secondary support a planet mass object for the inner companion. On the other hand, the LBT observations obtained in Ortiz et al. (2016) would not be sensitive to a white dwarf secondary with 1 . Thus, in principle, both companions of HD 59686 can have very small , making the system a hierarchical retrograde triple of a Kgiant, Mdwarf and a white dwarf. Such an exotic system would be stable in a retrograde orbit, but then in all cases the excited innerbinary eccentricity would oscillate with much larger values than the one currently observed, i.e., our observations would have caught the system at a very special time. Therefore, we conclude that the inner companion is most likely of planetary origin.
Overall, we identify a large and confident stable region for the HD 59686 system, which turns out to be at high mutual inclinations of 145, corresponding to nearly coplanar, but retrograde orbits. Mutual inclinations with between 30 and 145 lead to instability on very short time scales due to LidovKozai effects, and thus such orbital configurations are very unlikely. Nearly coplanar and prograde best fits with are also unstable.
7 Summary and Conclusions
HD 59686 is without any doubt a very interesting threebody system. It consists of a singlelined spectroscopic binary with a K giant primary with = 1.92 , a lowmass secondary star with at least = 0.53 and at most 1 , and an additional Stype planet with at least 7 . This system has a challenging architecture, since the secondary star orbits beyond the planet orbit ( = 1.09 AU) on a relatively close ( = 13.6 AU) and very eccentric ( = 0.73) orbit. As a result, the binary periastron distance is only = 3.7 AU, leading to strong interactions with the planet, and thus challenging the system’s longterm stability.
In this paper, we performed a detailed orbital and stability analysis of the HD 59686 system via dynamical modeling of RV data and longterm body integrations. We aimed to refine the orbital parameters by stability constraints, which can provide clues on the formation history. This is important since only a handful of Stype planetary candidates in compact binary systems are known in the literature, and the HD 59686 system illustrates how planets could form and remain stable in an Stype orbit around a star under the strong gravitational influence from a close stellar secondary.
Our global best fit with the lowest suggests a triple star system with nearly polar orbits ( = 92.7 3.3), instead of a binary with an Stype planet. We have shown, however, that such orbits quickly lead to instability due to the LidovKozai effect. Within only 600 yrs the LidovKozai effect excites the eccentricity of the inner object to a value which leads to collision with the primary K giant star. Orbital fits with parameters similar to the nearpolar global best fit are very unstable and experience a similar fate, although they all have lower values compared to the fits corresponding to a coplanar configuration. We conclude that the near polar configurations cannot represent the true system configuration, and that their small values are most likely a result of model degeneracies, which come from the limited number and accuracy of RV data points. These conclusions are supported by our bootstrap statistical analysis.
We find that HD 59686’s planet can survive only on nearly coplanar, most likely retrograde orbits. We find that when the system’s mutual inclination is = 180 (i.e., coplanar and retrograde), the system is fully stable for a large set of orbital solutions and companion masses. Longterm stability is also preserved for nearly coplanar retrograde configurations with 145 180.
Although most of the coplanar prograde fits consistent with HD 59686’s RV data are unstable, we have demonstrated that stable prograde fits in fact do exist. Our bootstrap and grid search analysis shows that a fraction of prograde fits (mostly within from the best fit) are stable for at least 10 Myr. These fits are located in narrow strips of the orbital period space where the initial period ratio / is not an integer number, with the best chances of stability near / 38.4 and 39.4. Therefore, the system can survive only between highorder MMR, while the MMR themselves have a destabilizing effect on the Stype planet. However, we find that the planetary and are also very important parameters which control the planet stability. The stability results indicate that the planet must be initially on a nearly aligned orbit with the binary () and have 0.06. In such a configuration, longterm stability is ensured by secular apsidal alignment between the binary and planetary orbits, with librating around with relatively small amplitude.
The stable islands shrink when we assume lower inclinations and more massive companions, while keeping coplanarity. Below , all prograde coplanar configurations are unstable, which suggests that if the system in indeed prograde and coplanar, then 90 45, with most stable fits at = 90. The orbital dynamics of these stable prograde fits with larger masses similarly shows secular apsidal alignment where librates around 0.
As a final discussion point, we note that there are arguments in favor of and against both prograde and retrograde configurations. For example, the retrograde stable region is very large and it can explain the RV data with great confidence, but forming a retrograde planet requires some exotic scenarios (see discussion in Ortiz et al., 2016). On top of that, looking at the retrograde planet eccentricity evolution (middle panel of Figure 3), we estimate that 18% of the time (within 3 of the bestfit value), and only 8% of the time (within 1). If the system is indeed in a retrograde configuration, then the Lick RVs must have been obtained in a phase which has rather low probability when is as low as . On the other hand, small is not a problem for the prograde stable fits where most of the time the planet eccentricity is in the range (Figure 10). However, the stable prograde islands are very narrow, and whether it is possible to form a massive planet in a narrow stable region in secular apsidal alignment with the eccentric close binary is a problem that deserves a closer look in the future.
References
 Baluev (2009) Baluev, R. V. 2009, MNRAS, 393, 969
 Barban et al. (2004) Barban, C., De Ridder, J., Mazumdar, A., et al. 2004, in ESA Special Publication, Vol. 559, SOHO 14 Helio and Asteroseismology: Towards a Golden Future, ed. D. Danesy, 113
 Bean & Seifahrt (2009) Bean, J. L., & Seifahrt, A. 2009, A&A, 496, 249
 Bevington & Robinson (2003) Bevington, P. R., & Robinson, D. K. 2003, Data reduction and error analysis for the physical sciences
 Borucki et al. (2010) Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977
 Campbell et al. (1988) Campbell, B., Walker, G. A. H., & Yang, S. 1988, ApJ, 331, 902
 Chauvin et al. (2011) Chauvin, G., Beust, H., Lagrange, A.M., & Eggenberger, A. 2011, A&A, 528, A8
 Correia et al. (2008) Correia, A. C. M., Udry, S., Mayor, M., et al. 2008, A&A, 479, 271
 Correia et al. (2010) Correia, A. C. M., Couetdic, J., Laskar, J., et al. 2010, A&A, 511, A21
 Doyle et al. (2011) Doyle, L. R., Carter, J. A., Fabrycky, D. C., et al. 2011, Science, 333, 1602
 Dvorak (1986) Dvorak, R. 1986, A&A, 167, 379
 Eberle & Cuntz (2010) Eberle, J., & Cuntz, M. 2010, ApJ, 721, L168
 Efron (1979) Efron, B. 1979, Ann. Statist., 7, 1
 Eggenberger & Udry (2010) Eggenberger, A., & Udry, S. 2010, in EAS Publications Series, Vol. 41, EAS Publications Series, ed. T. Montmerle, D. Ehrenreich, & A.M. Lagrange, 27–75
 Eggleton & Kiseleva (1995) Eggleton, P., & Kiseleva, L. 1995, ApJ, 455, 640
 Frink et al. (2001) Frink, S., Quirrenbach, A., Fischer, D., Röser, S., & Schilbach, E. 2001, PASP, 113, 173
 Giuppone et al. (2012) Giuppone, C. A., Morais, M. H. M., Boué, G., & Correia, A. C. M. 2012, A&A, 541, A151
 Goździewski et al. (2013) Goździewski, K., Słonina, M., Migaszewski, C., & Rozenkiewicz, A. 2013, MNRAS, 430, 533
 Haghighipour (2006) Haghighipour, N. 2006, ApJ, 644, 543
 Hamilton & Burns (1992) Hamilton, D. P., & Burns, J. A. 1992, Icarus, 96, 43
 Hatzes et al. (2003) Hatzes, A. P., Cochran, W. D., Endl, M., et al. 2003, ApJ, 599, 1383
 Hekker & Meléndez (2007) Hekker, S., & Meléndez, J. 2007, A&A, 475, 1003
 Hekker et al. (2006) Hekker, S., Reffert, S., Quirrenbach, A., et al. 2006, A&A, 454, 943
 Holman & Wiegert (1999) Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621
 Kaeufl et al. (2004) Kaeufl, H.U., Ballester, P., Biereichel, P., et al. 2004, in Society of PhotoOptical Instrumentation Engineers (SPIE) Conference Series, Vol. 5492, Society of PhotoOptical Instrumentation Engineers (SPIE) Conference Series, ed. A. F. M. Moorwood & M. Iye, 1218–1227
 Kjeldsen & Bedding (2011) Kjeldsen, H., & Bedding, T. R. 2011, A&A, 529, L8
 Kostov et al. (2013) Kostov, V. B., McCullough, P. R., Hinse, T. C., et al. 2013, ApJ, 770, 52
 Kozai (1962) Kozai, Y. 1962, AJ, 67, 591
 Lee et al. (2006) Lee, M. H., Butler, R. P., Fischer, D. A., Marcy, G. W., & Vogt, S. S. 2006, ApJ, 641, 1178
 Lee & Peale (2003) Lee, M. H., & Peale, S. J. 2003, ApJ, 592, 1201
 Leung & Lee (2013) Leung, G. C. K., & Lee, M. H. 2013, ApJ, 763, 107
 Lidov (1962) Lidov, M. L. 1962, Planet. Space Sci., 9, 719
 Mardling (2008) Mardling, R. A. 2008, in Lecture Notes in Physics, Berlin Springer Verlag, Vol. 760, The Cambridge NBody Lectures, ed. S. J. Aarseth, C. A. Tout, & R. A. Mardling, 59
 Mardling & Aarseth (2001) Mardling, R. A., & Aarseth, S. J. 2001, MNRAS, 321, 398
 Michtchenko & Malhotra (2004) Michtchenko, T. A., & Malhotra, R. 2004, Icarus, 168, 237
 Morais & Giuppone (2012) Morais, M. H. M., & Giuppone, C. A. 2012, MNRAS, 424, 52
 Nelson et al. (2016) Nelson, B. E., Robertson, P. M., Payne, M. J., et al. 2016, MNRAS, 455, 2484
 Ortiz et al. (2016) Ortiz, M., Reffert, S., Trifonov, T., et al. 2016, A&A, 595, A55
 Press et al. (1992) Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 1992, Numerical recipes in FORTRAN. The art of scientific computing
 Ramm (2015) Ramm, D. J. 2015, MNRAS, 449, 4428
 Ramm et al. (2009) Ramm, D. J., Pourbaix, D., Hearnshaw, J. B., & Komonjinda, S. 2009, MNRAS, 394, 1695
 Ramm et al. (2016) Ramm, D. J., Nelson, B. E., Endl, M., et al. 2016, ArXiv eprints, arXiv:1605.06720
 Reffert et al. (2015) Reffert, S., Bergmann, C., Quirrenbach, A., Trifonov, T., & Künstler, A. 2015, A&A, 574, A116
 Reffert & Quirrenbach (2011) Reffert, S., & Quirrenbach, A. 2011, A&A, 527, A140
 Roell et al. (2012) Roell, T., Neuhäuser, R., Seifahrt, A., & Mugrauer, M. 2012, A&A, 542, A92
 Takeda & Rasio (2005) Takeda, G., & Rasio, F. A. 2005, ApJ, 627, 1001
 Tan et al. (2013) Tan, X., Payne, M. J., Lee, M. H., et al. 2013, ApJ, 777, 101
 Thebault & Haghighipour (2014) Thebault, P., & Haghighipour, N. 2014, ArXiv eprints, arXiv:1406.1357
 Trifonov et al. (2014) Trifonov, T., Reffert, S., Tan, X., Lee, M. H., & Quirrenbach, A. 2014, A&A, 568, A64
 Trifonov et al. (2015) Trifonov, T., Reffert, S., Zechmeister, M., Reiners, A., & Quirrenbach, A. 2015, A&A, 582, A54
 van Leeuwen (2007) van Leeuwen, F. 2007, A&A, 474, 653
 Welsh et al. (2012) Welsh, W. F., Orosz, J. A., Carter, J. A., et al. 2012, Nature, 481, 475
 Wisdom & Holman (1991) Wisdom, J., & Holman, M. 1991, AJ, 102, 1528
 Wittenmyer et al. (2014) Wittenmyer, R. A., Tan, X., Lee, M. H., et al. 2014, ApJ, 780, 140
 Zechmeister et al. (2008) Zechmeister, M., Reffert, S., Hatzes, A. P., Endl, M., & Quirrenbach, A. 2008, A&A, 491, 531