Revisiting the Architecture of the KOI-89 System

We analyzed the long-cadence, Pre-search Data Conditioning (PDC) light curves downloaded from the Mikulski Archive for Space Telescopes.⁴⁴4https://archive.stsci.edu The transits of each planet were fitted iteratively as described in Masuda (2015), using the analytic light curve model for the quadratic limb-darkening law (Mandel & Agol, 2002) as implemented in batman (Kreidberg, 2015), multiplied by a quadratic polynomial function of time that accounts for longer-term trends. We fit both transit times and mid-durations letting impact parameters to vary for each transit, and the fitted transit times were used to stack transits without binning; then the stacked transits were fitted to update transit-shape parameters for the next iteration. After 10 iterations, the light curve around each transit was divided by the best-fit polynomial. We use these normalized and detrended transit light curves for photodynamical modeling. An overlapping transit around $\mathrm{BJD}_{\mathrm{TDB}}=2454833+912.7$ was omitted from this analysis, but was later used for a consistency check.

We fit osculating Jacobi elements defined at $\mathrm{BJD}_{\mathrm{TBD}}=2454833+149$ and mass ratios for each planet, as well as quadratic limb-darkening coefficients (as parameterized in Kipping, 2013) and mean density of the star.⁵⁵5Because the light-curve model alone does not constrain masses and lengths separately, we fix the stellar mass to be $1\,M_{\odot}$ and fit only the density of the star as well as ratios of masses and radii of the star and planets without loss of generality. We used TTVFast (Deck et al., 2014) to compute model mid-transit time, sky-plane velocity, and impact parameter during each transit of each planet. They were used to calculate planet locations on the stellar surface assuming a linear function of time,⁶⁶6Ignoring the change in the sky-projected velocity during transits has a negligible effect for such long-period planets. and the locations, along with planet-to-star radius ratios and limb-darkening coefficients, were used to calculate the relative flux losses. We again adopted the Mandel & Agol (2002) model for a quadratically limb-darkened star, as implemented in batman (Kreidberg, 2015), and computed the long-cadence fluxes $\bm{m}$ with a supersampling factor of 11. The residuals between the data $\bm{f}$ and the transit model $\bm{m}$ were modeled as a Gaussian process whose covariance matrix consists of a Matérn-3/2 covariance and a white-noise term. Thus the log-likelihood of the model $\ln\mathcal{L}$ is

where

Here $t_{i}$ and $\sigma_{i}$ are the time and PDC flux error of the $i$th data point, respectively. This likelihood was evaluated using celerite (Foreman-Mackey et al., 2017).

We assumed a Gaussian prior for the mean stellar density with a central value of $0.53\,\mathrm{g\,cm^{-3}}$ and a width of $0.06\,\mathrm{g\,cm^{-3}}$ based on the results in Section 2. For the other parameters, we adopted uniform or log-uniform priors on the parameters as specified in Table 2. The posterior samples were obtained using the nested-sampling algorithm MultiNest (Feroz & Hobson, 2008; Feroz et al., 2009, 2013) and its python interface PyMultiNest (Buchner et al., 2014). We typically used 4000 live points and a sampling efficiency of 0.8, and kept updating the live points until an evidence tolerance of 0.5 is achieved. We allowed for the detection of multiple posterior modes.

As the simplest possible explanation, we first adopt a two-planet model. We first fit the data assuming log-uniform priors for the planet-to-star mass ratios $m_{1}/M_{\star}$ and $m_{2}/M_{\star}$, and a truncated uniform prior $e_{1}<0.8$ and $e_{2}<0.45$ for the orbital eccentricities, where the subscripts 1 and 2 denote KOI-89.01 and KOI-89.02, respectively. These prior limits on $e_{1}$ and $e_{2}$ ensure that the orbits of the two planets do not cross, i.e., the periastron distance of the outer orbit is larger than the apastron distance of the inner orbit. Crossing orbital configurations lead to rapid instabilities and are thus short-lived.

Figures 2 and 3 compare the posterior models with the data for KOI-89.01 and KOI-89.02, respectively. The constraints on the parameters are summarized in Table 2. Here the mass and radius ratios and stellar density derived from the dynamical modeling were converted into physical masses and radii by sampling the stellar mass from a Gaussian with the central value $1.59\,M_{\odot}$ and width $0.08\,M_{\odot}$ (Section 2). In Figure 4, mid-transit times and durations computed for posterior samples of the model parameters are compared with the values measured by fitting individual transits (see Section A.2). This comparison is to illustrate that the variations in transit times and durations are well captured in our photodynamical modeling, in which we model the whole transit light curves as shown in Figures 2 and 3. We note that the transit times and durations derived by fitting the simultaneous transits of the two planets (shown as diamonds in this figure), which were not included in the photodynamical modeling, are in reasonable agreement with the posterior predictions conditioned on the other transits. This serves as a sanity check of the modeling.

We checked for long-term stability of our solution by integrating orbital configurations randomly drawn from our posterior samples. The integrations were performed using the IAS15 integrator (Rein & Spiegel, 2015) as implemented in the REBOUND package (Rein & Liu, 2012). The stepsize is automatically controlled to maintain numerical errors near machine precision. All the solutions turned out to be dynamically stable at least for $10^{7}$ inner orbits (i.e., 2.3 Myr) and no sign of instability was found.

From this analysis, we found lower and upper limits on the masses of KOI-89.01 and KOI-89.02, respectively. The former comes from strong TTVs of KOI-89.02 (see top-right panel of Figure 4), and the latter is due to the weak, if any, TTVs of KOI-89.01 (top-left panel of Figure 4). The orbital eccentricities are constrained to be low, mainly by the TTV signals and partly by a combination of the transit shape and prior constraint on the mean stellar density.

The orbital misalignment in the sky plane ($\Omega_{1}$ in Table 2, approximately equivalent to the mutual inclination given that both planets transit) is found to be small from the lack of significant transit duration variations (TDVs), which would otherwise be induced through the nodal orbital precession caused by the planets’ mutual gravitational interactions (Miralda-Escudé, 2002). We note that the high ($\sim 0.9$) impact parameter of KOI-89.02 makes its transit duration particularly sensitive to the perturbations perpendicular to its orbit, and that a weak variation is apparent in Figure 4 even for our solution with a small mutual inclination. See Appendix B for more detailed discussion on this mutual inclination constraint.

The weak constraints on mass ratios from this analysis hint that the result can be sensitive to the adopted prior. Thus, we repeated the analysis with uniform priors on mass ratios whose ranges were selected to be $[0,200]\,M_{\oplus}/M_{\odot}$ and $[0,50]\,M_{\oplus}/M_{\odot}$ for KOI-89.01 and KOI-89.02 respectively, based on the above result and the radii ($\approx 3\,R_{\oplus}$) of the planets. The range of $\Omega_{1}$ was also narrowed down to $[-30^{\circ},30^{\circ}]$ to speed up the convergence. The results are summarized in the right part of Table 2. As expected, the marginal posteriors for the mass ratios changed slightly, although we again recover firm lower and upper mass limits for KOI-89.01 and KOI-89.02 respectively. Additionally, the solution still points to near-circular, well-aligned orbits. We adopt this latter result assuming uniform priors on the mass ratios as the canonical one.

We note that the above solution is consistent with a previous dynamical analysis performed by Ahlers et al. (2015). Although their modeling of the light curve shape favored $e_{2}\gtrsim 0.5$, they found that such solutions are not dynamically stable for $10^{8}\,\mathrm{yr}$ or longer. Instead, they found that the system can be long-term stable if $e_{2}\lesssim 0.35$ and the orbits are co-planar. Here we presented a solution that satisfies this requirement and also explains the transit shapes and TTVs. We will show in Section 4 that this remains the case even if we take into account the possible effect of rapid stellar rotation on the transit shapes.

In Section 3.3, we showed that the two-planet model well explains the transit light curves of the two planets including TTVs. This is the simplest model given the current data, and there is no strong need to add another planet. Nevertheless, here we explore three-planet models motivated by an unusually large mass ($\gtrsim 20\,M_{\oplus}$) estimated for the inner KOI-89.01 given its size ($\approx 3\,R_{\oplus}$), as we will discuss in more detail later (Section 5.2). This is a solid lower limit imposed by the large TTVs exhibited by KOI-89.02, whose orbit is not particularly close to that of KOI-89.01, nor to any low-order mean-motion resonances with it.

Here we fixed $\log_{10}(m_{1}/M_{\star})=0.77$ based on the Chen & Kipping (2017) mass–radius relation as a “reasonable” value for $m_{1}$, and searched for solutions where the third planet (instead of KOI-89.01) may explain the TTVs of KOI-89.02. Here we assume that the three planets have mostly aligned orbits by fixing the longitudes of ascending nodes to be the same value and $\cos i_{3}$ to be zero, in the spirit of searching for solutions more typical of Kepler multi-planet systems.

Even under this simplification, given the sparse TTV data and the several new degrees of freedom introduced by the third planet, we expect that the solution will be complicated and multimodal. We therefore divided the prior range for $P_{3}$ into $24$ equally-spaced log-intervals between $1.05\,P_{2}$ and $4\,P_{2}$, obtained posterior samples separately using the same MultiNest fitting, and combined them by weighing the evidence value for each interval. Mathematically, the normalized “partial” prior for the $i$th bin $p_{i}(P_{3})$ (i.e., log-flat in the $i$th bin and zero otherwise) and the “full” normalized prior are related by $p(P_{3})=\sum_{i}p_{i}(P_{3})/N$ with $N$ being the number of $P_{3}$ bins. Therefore,

where $\bm{\theta}$ is the set of parameters other than $P_{3}$, $p_{i}(\bm{\theta},P_{3}|\bm{d})$ is the posterior conditioned on the partial prior $p(\bm{\theta})p_{i}(P_{3})$, and $p_{i}(\bm{d})$ is the corresponding evidence. Although the reliability of this “combined” posterior sample may be limited due to the limited accuracy of evidence evaluation, we believe that the sample is sufficient for our current purpose to capture the plausible range of three-planet solutions.

We find that the three-planet model can explain the data equally well as the two-planet model. The maximum likelihood value is slightly higher, but not high enough to favor the three-planet hypothesis over the two-planet one, considering the increased model complexity (e.g., in terms of the Bayesian information criterion). The constraints for the selected parameters are summarized in Table 3, although we note that $m_{3}$ and $P_{3}$ have multimodal marginal posteriors and the statistics here do not fully capture the complex nature of the solutions. We found significantly non-zero $m_{3}$, which is consistent with our inference in the two-planet model that “normal” values of $m_{1}$ are insufficient to explain the strong TTVs of KOI-89.02. On the other hand, we found relatively well-constrained posteriors for the parameters of the inner two planets, including the mass of KOI-89.02 and eccentricities of both KOI-89.01 and KOI-89.02. In this model two planets both have $10$–$20\,M_{\oplus}$ (for $M_{\star}=1.59\,M_{\odot}$) and low eccentricities.

The findings from our photodynamical analysis are summarized as follows.

• •

  Significant TTVs of KOI-89.02 suggest either an unusually large mass ($\gtrsim 20\,M_{\oplus}$) for KOI-89.01 given its radius, or the presence of an undetected third planet.

• •

  In both of the two-planet and three-planet models, the mass of KOI-89.02 is in the planetary range, and the orbits of KOI-89.01 and KOI-89.02 have low eccentricities.⁷⁷7The conclusion should remain unchanged even if the TTVs of KOI-89.02 is caused by a third planet in between KOI-89.01 and KOI-89.02, in which case there is even smaller room for solutions involving larger masses and eccentricities. In the two-planet model, a low mutual inclination is also required by the data.

The above analysis also demonstrates an advantage of full photodynamical modeling over a two-step approach in which the TTVs and TDVs are derived first and then fitted with a dynamical model. The latter approach usually employs additional constraint on the mean impact parameter or the in-transit velocity (i.e., combination of the mean stellar density, eccentricity, and argument of periastron) derived from the transit shape, to solve the degeneracy between eccentricities and inclinations (e.g., Sanchis-Ojeda et al., 2012). This works well if the impact parameter is well constrained, but for KOI-89.02 the eccentricity and impact parameter are strongly degenerate as shown explicitly in Appendix A.1. In such a case, the result of the dynamical fit can depend sensitively on how those constraints from transit shapes are incorporated, and the issue could become more serious if the correlate noise in the light curve (which we modeled here) significantly affects those constraints. For example, if one adopts the mean ($\approx 0.6$) of the marginal posterior of the KOI-89.02’s transit impact parameter (Appendix A.1) and fits TTVs and TDVs, the modeling may erroneously favor high-eccentricity and mutually inclined solutions as we explored in Section 3.3.1. For KOI-89, we were able to exclude these solutions based on stability argument in any case, but that may not always be possible. See also Dawson (2020) for an in-depth study of related issues.

The inference of a large stellar obliquity in Ahlers et al. (2015) was largely based on an apparently V-shaped transit light curve for KOI-89.02 (their Figure 4), which was argued not to be well fit by a high impact parameter transit. This lead them to obtain a solution involving a rapidly rotating, nearly pole-on star. Indeed, the pole-to-equator flux contrast induced by strong gravity darkening can make transits V-shaped, and a pole-on configuration was also required to explain the observed $v\sin i_{\star}$. The solution also involved a moderate impact parameter so that the planet can travel across the regions on the stellar surface with significantly different surface brightness; this in turn led to a large eccentricity to explain the short observed transit duration.

The crux of the disagreement is therefore that we do not find the stacked transit light curve of KOI-89.02 to be V-shaped (Figure 5, bottom). The origin of the difference of the derived transit shapes is unclear, but possible explanations include incorrect TTV corrections and distortion due to the binning of the light curve. Our stacked light curve, as well as individual transits of KOI-89.02, are well modeled by a high-impact parameter orbit and obviate the need for a large stellar obliquity. In addition, we found that a large orbital eccentricity as required in the previous solution results in rapid orbital instability (as was also pointed out in Ahlers et al., 2015) and is not favored from a dynamical point of view (see Section 3).

We focused on the KOI-89 system because it is one of three multi-transiting systems with claimed significant spin–orbit misalignments. Such systems might imply that the inner regions of planetary systems are occasionally misaligned with the stellar spin, for reasons unrelated to dynamical evolution of planetary orbits (e.g., misaligned disk). As discussed in Section 1, the large stellar obliquity in the other two systems, Kepler-56 and HD 3167, could plausibly be explained by dynamical evolution of the planetary orbits after formation. This left KOI-89 as the best known candidate for an obliquity excitation mechanism that does not alter planeatry orbits.

Our reanalysis of the KOI-89 data shows that there is no evidence for a high stellar obliquity as previously claimed by Ahlers et al. (2015), although a high obliquity is not necessarily ruled out. We therefore conclude that there is currently no clear evidence for spin–orbit misalignments unrelated to dynamical orbital evolution from the available obliquity constraints on stars with multiple transiting planets.

Obliquity measurements for a larger number of individual multi-transiting systems will serve as a benchmark for untangling the origin of high obliquities. This is particularly important for hot stars like KOI-89, given the statistical evidence that high obliquities may be the norm among stars hotter than $\approx 6000\,\mathrm{K}$, not only for those with hot Jupiters (Schlaufman, 2010; Winn et al., 2010a; Albrecht et al., 2012) but for other Kepler stars with smaller and/or longer-period planets (Mazeh et al., 2015); although Winn et al. (2017) found otherwise and the reason of the apparent discrepancy remains unclear. Given that the majority of the samples in these works are stars with single transiting planets, which may be a part of multi-planetary systems with misaligned orbits (Zhu et al., 2018; He et al., 2019), these statistical results would be less sensitive to multi-transiting systems around hot stars. Obliquity constraints in multi-transiting systems would thus play a key role in uncovering the connection between the spin–orbit misalignment in generic Kepler systems around hot stars and dynamical excitation of planetary orbits.

KOI-89.02 exhibits significant TTVs (Figure 4). The simplest explanation is that they are caused by KOI-89.01, and we showed that this model works well if KOI-89.01 is more massive than $\sim 20\,M_{\oplus}$. The value is unusually high given its size ($\approx 3\,R_{\oplus}$): such a planet is physically allowed, but we do not (yet) know of such planets (Figure 7, top). This motivated us to test the scenario in which the TTVs of KOI-89.02 are caused by another undetected planet. We found that this model also explains TTVs well for more “natural” masses for the two transiting planets (Figure 7, bottom) — at the cost of increasing the model complexity to a level that is not justified by the data. It is in principle possible to quantify which model is favored from a Bayesian perspective, but the results would hardly be convincing given the lack of our current knowledge on the possible mass–radius relation, and on the plausible properties of the third planet conditioned on the properties of the star and two transiting planets.

Nevertheless, we feel that the two-planet solution with a massive KOI-89.01 is worth some discussion in light of recent observations that found planets with similar properties (i.e., prior knowledge that assigns more plausibility to a model involving a dense sub-Neptune). For example, Armstrong et al. (2020) reported the discovery of TOI-849b from the Transiting Exoplanet Survey Satellite (TESS) data. It is an ultra-short-period (0.765524 days) planet with a mass of $39.1^{+2.7}_{-2.6}\,M_{\oplus}$ and a radius of $3.44^{+0.16}_{-0.12}\,R_{\oplus}$. The mass and radius indicate that the gaseous envelope, if any, should be very thin and Armstrong et al. (2020) argued that KOI-849b may be the remnant core of a former giant planet. As shown in the top panel of Figure 7, the mass and radius of KOI-89.01 as inferred from the two-planet model is similar to TOI-849b, and imply that mass fraction of the gaseous envelope should be a few percent or less. If the two-planet model is correct, KOI-89.01, with an orbital period of $85\,\mathrm{days}$, would suggest that such a massive core without a thick envelope can form further away from the star where the effect of photoevaporation is less significant. This poses a theoretical challenge similar to the one posed by super Earths: why would such a planet not grow into a gas giant? This property is even more puzzling in the presence of the outer less massive KOI-89.02 with a thicker atmosphere (as simply implied by its larger radius), whose core would have had a longer formation timescale than the inner KOI-89.01, all else being equal. Such a large contrast in the mean densities has been observed in the Kepler-36 system (Carter et al., 2012) and has been explained by photoevaporation (e.g., Lopez & Fortney, 2013), but the planets of KOI-89 are much further away from the star where this explanation would not work. This might then suggest that accretion onto planetary cores significantly depends on the surrounding environment.

Another implication of the two-planet solution is that such dense, (sub-)Neptune-sized planets may be abundant on orbits wider than have been probed with Doppler surveys to date. If so, Doppler follow-up observations of relatively long-period transiting Neptune-sized planets orbiting bright stars from K2 or TESS would reveal more such planets. HD 95338 b (Díaz et al., 2020), a planet with $42.4^{+2.2}_{-2.1}\,M_{\oplus}$ and $3.89^{+0.19}_{-0.20}\,R_{\oplus}$ (Figure 7, top) on a 55-day period orbit around a bright K dwarf observed by TESS, may be such an example.