MNRAS A New Planet in the Kepler-159 System From Transit Timing Variations a

A New Planet in the Kepler-159 System From Transit Timing Variations

Chris Fox, Paul Wiegert
The University of Western Ontario, Department of Physics & Astronomy, London, Ontario, Canada
Contact e-mail:
   Fox, Wiegert
Last updated 2018 October 2
Accepted for publication by the Monthly Notices of the Astronomical Society

The Kepler Space Telescope has discovered thousands of planets via the transit method. The transit timing variations of these planets allows us not only to infer the existence of other planets, transiting or not, but to characterize a number of parameters of the system. Using the transit timing variations of the planets Kepler-159b and 159c, the transit simulator TTVFast, and the Bayesian Inference tool MultiNest, we predict a new non-transiting planet, Kepler-159d, in a resonant 2:1 orbit with Kepler-159c. This configuration is dynamically stable on at least 10 Myr time scales, though we note that other less stable, higher-order resonances could also produce similar TTVs during the three-year window Kepler was in operation.

planets and satellites: detection, methods: numerical, techniques: photometric, stars: individual: KIC 5640085
pubyear: 2018pagerange: A New Planet in the Kepler-159 System From Transit Timing Variations A New Planet in the Kepler-159 System From Transit Timing Variations

1 Introduction

While the Kepler Space Telescope primarily found planets via the transit method (Borucki et al., 2010), additional planets have been discovered by examining the changes in the transit timings. These transit timing variations (TTVs) result from gravitational perturbations due to other planets, which can hurry or delay the time that a planet transits its parent star (Agol et al., 2005; Holman & Murray, 2005). Here we look at the case of the Kepler-159 system (also known as KIC 5640085, KOI-448, 2MASS J19481684+4052076, WISE J194816.85+405207.5). This star has two confirmed transiting planets, Kepler-159b and Kepler-159c (hereafter named 159b and 159c). The latter shows strong transit timing variations, as high as 5 hours, suggesting a significant gravitational influence acting upon it.

To assess the possible existence of a third planet, we used two publicly available pieces of software. The first is TTVFast (Deck & Agol, et al., 2014). This program simulates the orbits of planets around a star and outputs calculated TTVs. The second piece of software used was MultiNest (Feroz et al., 2009), a Bayesian Inference tool.

Under the assumption that the observed TTVs are produced by an unseen planet, these two programs allow us to determine the best set of planetary parameters that reproduce the observed TTVs. We find a degenerate set of parameters, though with our new planet consistently in (or near) an orbital resonance with 159c. We will refer to this inferred planet as Kepler-159d.

2 Target Details

The relevant stellar information of Kepler-159 comes from Muirhead et al. (2012) and Mathur et al. (2017). Transit Timing data is taken from Holczer & Mazeh et al. (2016). Planetary radii were taken from Rowe et al. (2015). The host star is a M0V star with an estimated mass of 0.52 , radius of 0.50 , and effective temperature of 389380 K (Muirhead et al., 2012; Mathur et al., 2017). There are not any direct measurements of the mass, but we used initial mass estimates from Chen & Kipping (2017) (using their 2 values) and the mass-radius relation for sub-Neptunes from Wolfgang et al. (2016).

The two transiting planets have estimated radii based on the observed transit depths. Their radii suggest both are likely gaseous (Rogers, 2015). Figure 1 shows the transit data for both planets. The inner planet, 159b, has a period of 10.14 days. All TTVs are within 30 minutes of the expected time, and the majority are consistent with zero. Kepler 159b is periodic to within measurement error. However 159c has an average period of 43.58 days and shows significant TTVs over nearly all of its 30 transits, all well outside the measurement error. This planet’s behaviour is the primary reason for predicting the existence of a third body.

Figure 1: Transit Timing Variations of Kepler-159b and Kepler-159c, created using data from Holczer & Mazeh et al. (2016)
Parameter Kepler-159b Kepler-159c
Period (days) 10.1396240.000022 43.5825880.000044
Semi-major axis (AU) 0.07370.0028 0.19490.0075
Radius () 1.8700.230 2.6400.320
Mass (Chen)
Mass (Wolfgang)

Semi-major axis computed using stellar mass and planetary periods from Holczer & Mazeh et al. (2016). Planetary radii from Rowe et al. (2015). Mass estimates are from Wolfgang et al. (2016) and Chen & Kipping (2017).

Table 1: Properties of Kepler-159 Transiting Planets

The inclinations of 159b and 159c must be near 90 for a transit to occur. Note that when we use the term inclination, we are referring exclusively to the value as measured from the plane of the sky, not the mutual inclination of the planetary orbits. The impact parameters of -b and -c are and respectively (Rowe et al., 2015). These would give best-fit inclinations of and .

Transit simulations were set to start at the first transit of the inner planet, providing the starting position for 159b. This corresponds to a BJD of 2454970.890506. The timing of the first transit of 159c then gives its relative location at the start of the simulation.

We adopt the TTV measurements of Holczer & Mazeh et al. (2016). Of the 30 observed data points for 159c, there are 4 that are considered outliers by Holczer & Mazeh et al. (2016). These data points are transit numbers 9, 19, 26 and 27. Similarly, 159b had 10 points excluded. All such outliers were removed from the runs reported on here. We also ran simulations which included these points, and while they have an effect on the reported log-Evidence values, the conclusions are unchanged if the outliers are included.

3 Methods and Setup

3.1 Parameters

Each planet has 6 orbital elements plus a mass for a total of 7 parameters. Of the 14 total parameters of the two transiting planets, only six are known to high precision. These are the periods, inclinations, and relative starting positions (mean anomalies).

Periods of both 159b and 159c were kept fixed for all simulations. It was expected that slight deviations from 90 inclination (with respect to plane of the sky) would have little impact on the observed TTVs (Agol et al., 2005). Some test simulations allowing inclination to vary by 2 from 90 (chosen to maintain the transit condition) were performed and confirmed this expectation. Thus, to reduce the parameter space, for all simulations the inclination was set to 90 for both known transiting planets.

Because the mean anomaly is measured from the argument of periastron, the mean anomaly for each planet was not fixed but rather predetermined from the eccentricity, argument of periastron, and timing of the first transit. As such, the mean anomaly of 159c was different for each simulation but not a free parameter.

This leaves 8 unknown parameters, plus 7 for the presumed third planet. This is a 15-dimensional parameter space to search. To further reduce computational requirements, we reduced this to 11 parameters. This reduction was accomplished by presuming 159b was in a circular orbit, using its orbital plane as the reference for the system, and using a constant estimate of its mass. Thus, 159b’s eccentricity, argument of periastron, longitude of the ascending node, and mass are removed from the priors. This is discussed in more detail in section 3.4.2.

3.2 TTVFast and MultiNest

We used TTVFast (Deck & Agol, et al., 2014), a symplectic numerical integrator, to compute the transit times for each hypothesized set of orbital parameters. These times were converted to TTVs, then compared to the observed TTV data from Holczer & Mazeh et al. (2016). Another set of parameters would then be chosen by MultiNest, and the values compared again. To search the parameter space efficiently, we used the Bayesian Inference algorithm MultiNest (Feroz et al., 2009, 2013). This iterative process can find solutions far more quickly than searching the entire parameter space. We used the Python interface for MultiNest, PyMultiNest by Buchner et al. (2014) for this project.

Due to the size and complexity of the parameter space, finding the best fits required an iterative process and dozens of sets of runs. Each "run" produces a single result from MultiNest, estimated from millions of samples of the parameter space and as many simulations from TTVFast. Repeated runs, even with the same priors, would often return different best-fits. This is indicative of multiple maxima in our parameter space, particularly associated with different resonances between 159d and 159c. It is known that resonances produce strong TTVs and that there can be degeneracy between resonances in this regard (Boué et al., 2012). To examine this effect, we performed additional runs with periods of 159d constrained to be near the resonances (periods within 10 days of resonance) of 159c. This narrowing of ranges allowed us to determine the best-fits associated with each resonance independent of the others. We used the default parameters of PyMultiNest with exception of the evidence tolerance, which we set to 0.15 (from default of 0.5). This reduced tolerance value allows for a more fine-tuned search across the peaks of this rapidly varying landscape.

3.3 Likelihood

The Likelihood, we used is based on the usual statistic and is given by:


where is the calculated TTV, is the observed TTV, and is the error in the observed TTV for point , respectively (Ivezic et al., 2014). The log-likelihood is the value we aim to maximize in the simulation.


In this case, summation includes transit data from both 159b and 159c.

3.4 Priors

The priors used for the search were constrained based on physical grounds. First we tested the two-planet case, then we moved to the three-planet scenario. Several priors differed between these two cases. In this section, we describe the ranges of parameter space explored in our scenarios.

3.4.1 Testing the Two-Planet case

First, simulations were performed to ensure that the observed TTVs were not mutually induced. 159b was used as the reference planet, with longitude of the ascending node of 0°. The timing of the first transit of 159b serves as the start time for all simulations. This time along with the timing of the first transit of 159c thus determines the position of 159c at the start of each simulation. Masses of both planets were allowed to vary from Earth-mass to 10 . Most orbital parameters were allowed to vary freely so long as the initial conditions (initial transit times) were maintained. The periods were known and thus fixed, and the eccentricity was constrained such that the two orbits would not cross. The stellar mass was allowed to vary from 0.3 to 1.0 .

3.4.2 Adding the Third Planet

After it was determined that the two planets could not mutually induce the observed TTVs, we looked at adding a third planet. This required new priors for 159d, as well as modifications to the priors of 159b and 159c.

Kepler 159b: The inner planet, 159b shows no trends in its TTVs. It has a smaller radius from 159c suggesting it is likely a less massive planet. We use it as a constraint for the system; any good solution must maintain its lack of TTVs. 159b has a radius estimate of 1.87 (Table 1). Based on the work of Chen & Kipping (2017), and supported by the radius-mass relationship for sub-Neptune-radius planets established by Wolfgang et al. (2016), we estimate the mass of 159b to be 4.0-6.0 , though both of these sources allow for a wide range of mass values. Because we did not expect 159b to be a major player dynamically, we kept the mass of 159b to a fixed value of 4.0 . We set the longitude of the ascending node of 159b to 0 as the orbital reference plane for the system. The time of first transit is used as the starting point of each simulation, and thus not computed by TTVFast. We expect this planet to be in a circular orbit due to its proximity to the host star, so eccentricity was fixed to 0.

Kepler 159c: To encompass the mass possibilities for the outer transiting planet, 159c, its mass was allowed to vary from 2.0 to 35 , covering the likely ranges for sub-Neptune planets estimated by Wolfgang et al. (2016) for the planetary radius of 2.64 (see Table 1). Eccentricity was limited to a maximum of 0.66 to prevent crossing with the orbit of 159b. Initially, the full range of values of the longitude of the ascending node were explored. This resulted in some solutions with ascending node values of 159c and 159d near 180 from each other. The inclination values for the planets were always within 15 of 90. This combination of out of phase ascending nodes and near-90 inclinations correspond to the two planets orbiting in opposite directions. Retrograde planets appear to be valid solutions, but we deem it more likely that the planets all orbit in the same sense. So, we then restricted the ascending node prior to -45 to +45 for the solutions reported here. This ensured the orbital axes of the planets were not anti-aligned.

Kepler 159d: The parameters for the hypothesized planet, 159d, were completely unknown. Like 159c, the entire range of ascending node values was initially explored, but ultimately limited to -45 to +45 to prevent retrograde motion. Similarly, inclination was constrained to 45 to 135. We limited the eccentricity to no higher than 0.7 to prevent orbital crossing with 159c for periods as high as 250 days. Its period was explored from 50 days (just outside 159c’s orbit) to 300 days, and its mass allowed to range up to 30 Jupiter masses. We found that 159d’s period always settled near a resonance with 159c, and ultimately we did individual runs with the period prior limited to 10 days around each resonance.

The stellar mass was set to the value from Muirhead et al. (2012) and Rowe et al. (2015) of 0.520 for all reported runs.

In total we varied 11 parameters. The mass of 159b was static and its orbit assumed circular. The three static values of 159c were the period, inclination, and (effectively) mean anomaly. The new planet had no static parameters. We used uniform priors in all tests, summarized in Table 2.

Parameter Kepler-159b Kepler-159c Kepler-159d
Mass 4.0 [2.0, 35.0] [0.05, 20.0]
Period 10.1396236 d* 43.58258676 d* [50.0, 300.0] d
Eccentricity 0.0 [0.0, 0.66] [0.0, 0.7]
Inclination 90.0°* 90.0°* [45.0, 135.0]°
Asc. Node 0.0° [-45.0,45.0]° [-45.0,45.0]°
Arg. of Peri. 0.0° [0.0, 360.0]° [0.0, 360.0]°
Mean Anom. 90.0°* ** [0.0, 360.0]°

Parameters for 159b were always static, assuming a simple circular orbit, and used as a simulation constraint.
*These parameters are known, and used as starting values for the start of all simulations.
**The mean anomaly of 159c for each simulation is dependent on its first transit time, eccentricity and argument of periastron.

Table 2: Planetary Parameter Priors for Testing

4 Results

4.1 Two Planet Scenario

The best fit scenario for the 2 planet case, where 159b and 159c mutually induce TTVs, is shown in Figure 2. The posterior results are in Table 3. The reported Log-Evidence value is -2508. Not only do the TTVs from 159c not match, TTVs are induced upon the inner planet, inconsistent with the observations from Kepler. Further, the best-fit eccentricity of 159c is 0.55, which puts its closest approach very close to the orbit of 159b. For these reasons, we conclude that the two planets do not mutually induce the observed TTVs and the existence of a third planet was then assumed.

Figure 2: Simulated TTVs of 2-Planet scenario. A blue + symbol is the ttv computed using Holczer & Mazeh et al. (2016) data, (uncertainty included), and a red x represents the simulated TTV.
Parameter Kepler-159b Kepler-159c
Mass 0.1049244768 6.2583441849
Eccentricity 0.0077796061 0.5542004366
Asc. Node 0.0° 91.9345810745°
Arg. of Peri. 0.3672973405° 255.727757573°
Table 3: Best Fit Inputs for 2-Planet Simulations

4.2 Three Planet Scenario

We found several configurations that successfully reproduced the observed TTVs with a small degree of error. The best configurations correspond to the 2:1, 3:1, 4:1, 5:1 and 6:1 resonances of 159d with 159c. Though most of these are only near and not actually in resonance, we will refer to each as N:1 case (or resonance) for the sake of convenience.

The simulated TTVs are plotted on top of the observed TTVs in Figure 3. Residuals for each case are shown in Figure 4 along with the fit value. From these results, we conclude that an additional planet, Kepler-159d, exists in this system.

In all simulated best-fit solutions, the newly discovered planet d does not transit. For the smallest orbit, the 2:1 case, the required inclination for a transit is 90, with even tighter constraints on the higher period solutions. However, our results in all cases put the inclination well outside this range. Our results are thus consistent with the lack of observations of such a transit from Kepler.

The simulated TTVs exerted on 159b were essentially non-existent in all cases reported here; never varying from the zero point by more than 0.001 days. Thus, our solutions are consistent with the observed behavior of 159b.

The marginal values (posteriors) from MultiNest provided a range of values and errors (1 (68%) values) for each parameter. The second set of values reported from MultiNest are the Best-Fit values, which are the set of parameters that provided the single highest Log-likelihood value sampled at any point in the algorithm. The Best-Fit values do not have errors reported since they represent exact values. The Best-Fit values all reside inside the errors of the posteriors.

The results from MultiNest were fed back into TTVFast manually as a data check. The Best-Fit values provided excellent correlation between the plots. These are the plots shown in Figures 3 and 4. However, one should use caution with the outputted posterior values. The nominal values from Table 4, when plugged into TTVFast, do not provide good matches with observations; the TTV curves maintain an overall arching shape, but progressively deviate from the observed data with each successive transit. This misfitting occurs due to the central posterior values being applied without consideration of the correlations that occur between parameters.

Figure 3: Best-Fit Simulated TTVs for 159c for each N:1 Resonant Case. A blue + symbol is the observed data with error bar, and a red x represents the simulated TTV.
Figure 4: Best-Fit Residuals for 159c for each N:1 Resonant Case
Parameter 2:1 3:1 4:1 5:1 6:1
d mass
d period days
d eccentricity
d inclination degrees
d longitude of ascending node degrees
d argument of periastron degrees
d mean anomaly degrees
c mass
c eccentricity
c longitude of ascending node degrees
c argument of periastron degrees
Nested Sampling Global Log-Evidence 475.70.02 487.80.02 485.00.02 482.00.04 493.40.03
Table 4: Marginalized Posterior Values with 1 Sigma Error for Each Resonant Case.
Parameter 2:1 3:1 4:1 5:1 6:1
d mass 0.3221631382 0.9823871944 3.1845969511 6.9502940711 14.9134726287
d period 88.6719024826 130.940442456 173.231659674 216.17353574 262.642480693 days
d eccentricity 0.0328965561 0.1771459360 0.0491852474 0.0783972053 0.2050824733
d inclination 91.0428510417 80.000250445 80.2976935675 96.6614963008 85.5464215937 degrees
d longitude of ascending node 23.9249822741 -19.7243818881 13.3203978877 18.1524418105 25.9271945662 degrees
d argument of periastron 2.7596764301 239.442586615 177.585357201 98.3775876453 93.5510409237 degrees
d mean anomaly 276.813762079 44.9256736424 81.2210236607 169.477617267 198.029737345 degrees
c mass 7.1875572623 4.4947762292 14.8442353074 16.1473341603 8.4575152214
c eccentricity 0.0100266745 0.2330428640 0.1178657898 0.097737198 0.0673859389
c longitude of ascending node 14.1797256542 0.792090619646 -4.41129668525 -9.5079755857 7.67530949204 degrees
c argument of periastron 223.705031533 252.966717978 287.938758322 260.167384522 278.817603374 degrees
Log-Likelihood 516.27 527.96 520.85 518.80 535.39
Table 5: Best Fit Parameters for Each Resonant Case
Figure 5: Long-term (10 Myr) simulation results for best-fit systems. Each row contains the results for each case: the 2:1, 4:1, and 5:1. Each column corresponds to a different orbital element: periastron/apastron, semi-major axis, eccentricity, inclination.

5 Discussion

5.1 Long-term stability

The best-fit solutions provided by MultiNest represent good matches to the observed TTVs during the few year lifespan of Kepler but do not necessarily represent systems with long-term (million to billion year) stability. We make the assumption that the Kepler-159 system is stable on these longer time scales to provide an additional constraint. The solutions provided by MultiNest are subsequently simulated for 10 Myr into the past and future to test for signs of instability. The simulation code is symplectic and based on the Wisdom-Holman algorithm (Wisdom & Holman, 1991) using a time-step of 0.5 days, (less than 1/20th of the period of the innermost planet: 0.507 days), and includes post-Newtonian general relativistic effects. Some solutions which provide good matches in the short-term will be seen to have longer-term instability that make them unlikely to reflect the actual properties of the Kepler-159 system.

Integrations of the best-fit systems for 10 Myr into the past and the future are presented in Fig. 5. The primary result is that only the 2:1, 4:1 and 5:1 cases are stable. The 3:1 and 6:1 cases become unstable on time scales of only a few tens of thousands of years in both directions, resulting in a violent rearrangement of the planets. We will take this as evidence that these particular cases do not represent the actual state of Kepler-159, which is unlikely to be in such a precarious state.

The 2:1, 4:1 and 5:1 states all remain stable over this time frame without any discernible trends in the orbital elements. The 4:1 case undergoes significant changes in the inclinations (20) of both 159b and 159c, but not in its semi-major axes and eccentricities. The 5:1 case shows similar changes in 159b, but even wider variation (30) in the inclination of 159c. Further, in the 5:1 case the variation of eccentricity of 159c is higher and varies rapidly between 0.1 and 0.25. The 2:1 case shows similar variations in the inclination of 159b, but less so in 159c (10), and all planets in the 2:1 case show less variation in the semi-major axis and eccentricities.

The 2:1 resonant case can be verified to be in resonance by calculating the resonant argument (see e.g. Murray & Dermott (1999) equation 8.8 though we recast it for an outer perturber, since 159d is more massive) where we find the resonant argument librates with an amplitude of 50. The other higher-order cases show circulation of their resonant argument, indicating that the systems may be near but are not fully in mean motion resonance. Though it is common for exoplanetary systems to be near but not quite in resonance (Lissauer et al., 2011; Fabrycky et al., 2014), capture into first-order resonances like the 2:1 is preferred in models of planet formation including migration (Snellgrove et al., 2011; Papaloizou & Szuszkiewicz, 2005).

5.2 Other Resonances

We found that periods near any N:1 resonance out to 6:1 could produce a quality fit. Given the importance of the resonance in this system, we then looked at other resonances. We ran simulations near the 7:1, 3:2, 5:2, 7:2, 4:3 and 5:3 resonances. Of these, the 3:2 resonance is particularly important as it is a common configuration in exoplanetary systems (Fabrycky et al., 2014). We searched several resonances of these forms with a typical range of days of the resonance, so long as it didn’t overlap another case. The exception is the 4:3 case, whose range of periods also encompassed the 5:4 resonance.

Case log-Evidence
7:1 432.2 71.1
3:2 432.9 64.2
5:2 476.4 34.3
7:2 469.7 39.9
4:3 135.9 367.8
5:3 338.6 164.3
Table 6: Simulation Fit Quality for Other Resonances

The quality of fits for these runs is summarized in Table 6. The 7:1 and 3:2 simulations produced fits that did not as strongly match the earlier N:1 cases. The 4:3 and 5:3 resulted in very poor agreement with the observed TTVs with large deviations well outside the observational error. Only the 5:2 and 7:2 provided fits of a quality comparable to the N:1 cases. However, creating these cases required relatively high eccentricities (0.15-0.30), unlike the N:1 cases. Searching these regimes with forced low eccentricities resulted in far worse fits (log-Evidences near 100 and values near 400). Long-term simulations of the 5:2 and 7:2 showed them to be highly unstable, resulting in crossing orbits after only a few thousands of years. Thus, we conclude that the 7:1, 3:2, 5:2, 7:2, 4:3, and 5:3 cases do not represent the true nature of the system.

5.3 Favoured Configuration

We favour the 2:1 resonance as the most likely scenario for a number of reasons. First, we note that the quality of each N:1 fit shown here is essentially as good as any other. The log-Evidence values are all similar, and the Best Fit values are all of order of degrees of freedom; we have 26 data points and the worst value is 32. Second, smaller planets are more common than massive planets, particularly for lower-mass stars, and the 2:1 case requires the smallest mass for 159d. Third, planets with periods just above an exact 2:1 resonance are among the most common seen, while those above 3:1 resonance are quite rare (Fabrycky et al., 2014). Finally, our long-term simulations show the 2:1 to be a very stable case with orbital parameters not varying significantly over 10 million years. We conclude that the case where 159d is in 2:1 resonance with 159c is the most probable scenario.

5.4 Planetary Masses

The estimated mass of the new planet, 159d, is dependent on the case chosen, but has relatively low uncertainty within each case. The best fit masses are 0.32, 0.98, 3.18, 6.90, and 14.91 for the 2:1, 3:1, 4:1, 5:1, 6:1 cases respectively. In the 2:1 case, near-identical TTV curve fits can be produced with a mass of 159d as high as 0.49 and as low as 0.3 . Larger masses of d also required an increase in mass of 159c, seen in the correlation on Table 8. Our reported value of the mass of 0.32 was chosen due to this run having the highest log-Evidence value. Regardless of resonant case, the most massive object in the system is the new planet 159d. In all of our results, the lowest ratio of the mass of 159d to 159c in any run was 6. In our preferred case (best fit of 2:1), the ratio of the masses of 159d to 159c is approximately 14.

There is significant uncertainty in the mass of 159c. In addition to the reported error, multiple runs for the same resonance could produce very different Best Fit values for the mass of 159c. The highest value found was 30 with the lowest being 6.5 , with nearly as good results in both log-Evidence and values. There is a correlation (see Table 8) between the masses of 159c and 159d. Our best fit of 159c’s mass for the 2:1 case is 7.2 . This is within the initial mass estimates from Table 1, and would equate to a mean density about 20 percent higher than Neptune. However, the posterior values () have significant uncertainty. If the actual mass were on the upper end of this error, we obtain a density similar to that of Mars. Given the relative insensitivity of this parameter and that other runs with similar fits could have masses of around 30 , we cannot say anything definitive about the planet’s mass or density, nor whether it is gaseous or terrestrial.

All simulations were performed assuming a 4.0 of the inner planet, 159b. The mass was not set as a free parameters because it shows no signs of TTVs, suggesting its coupling to the other planets is not large. To confirm this hypothesis, we re-ran each resonant case with a different (but still static) mass of 159b. Modifying the mass by a factor of 0.25 to 3.0 (for a range of 1-12 ) resulted in only minor differences to the results for 159c. However, we found in the 2:1 case TTVs were induced on 159b as its mass was increased above about 8 , showing a linear increasing trend. From these tests, we conclude that 159b has a minor impact on the reported results for 159c and 159d, and we can place an upper limit on the mass of 159b of about 8 Earth masses.

5.5 Orbital Parameters

The orbit of 159b is used as the orbital reference plane (longitude of the ascending node of 0°). The other planets are inclined relative to the orbital plane of 159b. Further, 159c and 159d are consistently inclined with each other by anywhere from 10 to 20. Correlation plots on Table 8 show a strong relationship between the two values. The inclination with respect to 159b does not seem to matter much to this relation; it is the mutual inclination between 159c and 159d that is most relevant to the dynamics.

While the argument of periastron varies, the physical starting location of 159d relative to the other planets is quite consistent. The starting location of 159d was found to be on the opposite side of the star (argument of periastron + true anomaly 270) in all our best case fits. That is, it starts near superior conjunction. We did find the occasional run in which 159d started on the near side of the star, near inferior conjunction. In the 2:1 inferior conjunction case, the masses of 159c and 159d were much higher; 29 and 0.5 respectively. At no point did we ever find a result (in any N:1 resonant case) where the starting location of 159d was far from superior or inferior conjunction.

For most resonant cases, we found equally good solutions with very different eccentricities. In most cases, the eccentricities for both planets was less than 0.1, but in some cases they were 0.15 and higher. Long term simulations of these higher eccentricity solutions proved to be quickly unstable, regardless of the resonance case examined. In some cases (most notably the 3:1 and 6:1), we did followup simulations forcing the eccentricity prior to less than 0.1. While we found equally good fits with low eccentricities, long term simulations showed these systems to still be unstable. This suggests that low eccentricity is required but not sufficient for stability in this system.

5.6 Habitable Zone

Using the estimated size and temperature of Kepler-159 (Rowe et al., 2015; Mathur et al., 2017), the habitable zone for this system can be computed (Kane & Gellino, 2012; Kopparapu et al., 2013). This system has a conservative zone (runaway greenhouse to maximum greenhouse boundaries) of 0.24 AU to 0.45 AU, while the optimistic (recent Venus to early Mars) estimate puts it at 0.18 AU to 0.47 AU.

Semi-major Stellar Equil.
Case axis (AU) Flux () Temp. (K)* HZ region
2:1 0.313 0.525 217 conservative
3:1 0.406 0.313 191 conservative
4:1 0.489 0.216 174 no (near edge)
5:1 0.573 0.157 161 no
6:1 0.646 0.124 151 no

*Equilibrium Temperatures assume a bond albedo of 0.3.

Table 7: Relative Flux from Star for Each Planet Location

159b is too close to the star and is well inside the habitable zone. 159c is inside the optimistic zone. Whether the new planet, 159d, exists in the habitable zone depends on the resonant case, as shown in Table 7. In the 2:1 and 3:1 case, 159d would be inside the conservative zone. 159d in the 4:1 case resides just outside the outer edge of the optimistic zone.

5.7 Stellar Mass

We used a mass of 0.52 in all of our runs, based on Muirhead et al. (2012), Rowe et al. (2015) and Mathur et al. (2017). NASA’s Exoplanet Archive lists a "standard" mass of 0.72 , based on Rowe et al. (2014). This value is often seen in other sources, such as For comparison we ran some simulations with this higher mass estimate of 0.72 .

We found that changing the mass primarily affected the reported planetary masses with only slight differences in orbital parameters. In effect, the masses of the planets scale with the mass of the star with each resonant case requiring a mass typically around 50 percent higher than for the 0.52 star. The reported periods were identical to within 1.

6 Conclusions

We report the existence of a new planet in the Kepler-159 system, likely in a 2:1 resonant orbit, period of 88.7 days with 159c and with a mass similar to that of Saturn’s. We cannot with complete certainty rule out the 4:1 and 5:1 cases (corresponding masses of 3.1 and 6.9 ), but these are less likely. We can rule out the other resonant cases as they all prove to be unstable on the order of a few ten-thousands of years or less.

Further, we have gained insights into the other aspects of the system. All planets have low (<0.1) eccentricities. There is a strong correlation in the orbital planes of 159c and 159d, differing by an average of 12. The mass of 159c and 159b are uncertain, but are less than the mass of 159d by at least a factor of 10.

The Kepler-159 system is an illustration of the degeneracy of TTVs produced by a non-transiting planet. Firming up certain parameters will require further observations, either from more transits, radial velocity measurements (which we estimate from reflexive motion should be approximately 20 m/s in the 2:1 case, and higher for the other cases), or possibly transit duration variations. We hope to explore this system further in the future.

7 Acknowledgements

We thank the anonymous referee for their work which much improved the paper. We also thank Stephen Kane for his useful feedback and advice. This work was funded in part by the Natural Sciences and Engineering Research Council of Canada.


  • Agol et al. (2005) Agol E. et al., 2005, MNRAS, 359, 567
  • Borucki et al. (2010) Borucki W. et al., 2010, Science, 327, 5968
  • Boué et al. (2012) Boué G. et al., 2012, MNRAS, 422, L57
  • Buchner et al. (2014) Buchner et al., 2014, A&A, 564, A125
  • Chen & Kipping (2017) Chen J. & Kipping D., 2017, MNRAS, 473, 2753
  • Deck & Agol, et al. (2014) Deck K. & Agol E., et al., 2014, ApJ, 787, 132
  • Fabrycky et al. (2014) Fabrycky, D. C., Lissauer, J. J., & Ragozzine, D. et al. 2014, ApJ, 790, 146
  • Feroz & Hobson (2007) Feroz F. and Hobson M., 2007, MNRAS, 384, 2
  • Feroz et al. (2009) Feroz F., et al. 2009, MNRAS, 398, 4
  • Feroz et al. (2013) Feroz F., et al., 2013, arXiv:1306.2144v2
  • Holczer & Mazeh et al. (2016) Holczer T. & Mazeh T.,, 2016, ApJS, 225, 9
  • Holman & Murray (2005) Holman M. & Murray N. (2005) Science, 307, 1288
  • Ivezic et al. (2014) Ivezic, Z, Connolly, A., VanderPlas J., Gray, A., 2014, Statistics, Data Mining and Machine Learning in Astronomy, Princeton University Press.
  • Kane & Gellino (2012) Kane S. & Gellino D., 2012, PASP, 124, 323
  • Kopparapu et al. (2013) Kopparapu R, et al., 2013, ApJ, 765, 131
  • Lissauer et al. (2011) Lissauer, J. J., Ragozzine, D., & Fabrycky, D. C. et al. 2011, ApJS, 197, 8
  • Mathur et al. (2017) Mathur S, et al., 2017, ApJS, 229, 30
  • Murray & Dermott (1999) Murray C., & Dermott, S., 1999, Solar System Dynamics, Cambridge University Press.
  • Muirhead et al. (2012) Muirhead P, et al., 2012, ApJ, 750, 2, L37
  • Papaloizou & Szuszkiewicz (2005) Papaloizou, J. C. B. & Szuszkiewicz, E. 2005, MNRAS, 363, 153
  • Peale (1976) Peale S.J. 1976, ARA&A, 14, 215
  • Peale & Lee (2002) Peale S. J. & Lee M. H. 2002, ASP Conf. Ser., 294, 197
  • Rogers (2015) Rogers L., 2015, ApJ, 801, 41
  • Rowe et al. (2014) Rowe J., et al. 2014, ApJ, 784, 45
  • Rowe et al. (2015) Rowe J., et al. 2015, ApJ, 217, 16
  • Snellgrove et al. (2011) Snellgrove, M. D., Papaloizou, J. C. B., & Nelson, R. P. 2001, A&A, 374, 1092
  • Wisdom & Holman (1991) Wisdom J. & Holman M., 1991, ȷ, 102, 1528
  • Wolfgang et al. (2016) Wolfgang A., Rogers L. & Ford E., 2016, ApJ, 825, 19
d mass d period d eccentricity d inclination d ascending node d periastron d mean anomaly c mass c eccentricity c ascending node
d period
d ecc
d inc
d asc node
d peri
d mean anom
c mass
c ecc
c asc
c arg peri
Table 8: Correlation Plots for the 2:1 Resonance Case
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description