Kepler-90: Giant Transit-timing Variations Reveal a Super-puff

The Kepler mission has detected thousands of exoplanets via the transit method, in which the brightness of stars is measured precisely and continuously, seeking evidence for periodic dimming events due to eclipsing planets (Borucki et al. 2010). A small percentage of the discovered planets exhibit transit-timing variations (TTVs), which can be an indicator of gravitational interaction with other (possibly nontransiting) objects, and have been proposed as a powerful probe of planet masses and eccentricities (Agol et al. 2005; Holman & Murray 2005). TTVs have already been used to detect nontransiting planets (Nesvorný et al. 2012; Fox & Wiegert 2018) and to determine planet masses and eccentricities (Lithwick et al. 2012; Hadden & Lithwick 2017). Most of the existing analyses have been on shorter-period planets (under 100 days) with small-amplitude TTVs, exhibiting a characteristic sinusoidal structure: these are easy to identify in the data but require many transits. A sinusoidal signal can determine only a limited number of parameters and cannot break the degeneracy between the mass and eccentricity (Lithwick et al. 2012). In the purely sinusoidal case, only three independent measurements can be extracted from the observed TTVs, while there are seven parameters to be determined for each planet.

While exoplanet mass and eccentricity statistics are relatively abundant, characterization of the full orbital parameters is rare, and is even rarer for the long-period planets. The radial velocity (RV) is insensitive to the mutual inclinations, and direct imaging can only place a loose bound on the dynamical parameters (Wang et al. 2018). Combining RV with TTVs can break some, but not all, of the degeneracies inherent to RV or small TTV alone (Petigura et al. 2018). However, orbital parameters can be extremely valuable as they provide important information about the exoplanet formation, migration, and evolution that is not available otherwise (Winn & Fabrycky 2015; Petigura et al. 2018).

Super-puffs are planets with small masses and large radii, corresponding to a very low density. Most super-puffs are found close to their host star (0.1∼0.3 au) and their masses are typically measured from TTVs (Jontof-Hutter et al. 2014; Masuda 2014; Ofir et al. 2014). Two possible explanations of the low inferred densities are proposed: the absorption from a dusty atmosphere (e.g., Wang & Dai 2019) and the effect of exo-rings (see, e.g., Piro & Vissapragada 2020 for a discussion). The atmosphere accretion and migration mechanisms are discussed by Choksi & Chiang (2020). In this paper we will define a super-puff as a planet with an apparent density less than 0.3 g cm⁻³, regardless of the origin.

The Kepler-90 system has more known transiting planets than almost any other system. Among its seven confirmed planets, the inner five have orbital periods ranging from 7 to 125 days, while the outer two are 90g and 90h with orbital periods of 210.5 days and 331.6 days, respectively (Cabrera et al. 2013). Based on the six recorded transits of 90g and 3 transits of 90h during the four years of Kepler observations, both planets show significant TTVs. Most notably, the sixth observed transit of 90g exhibits a 25 hr anomaly, the largest TTV found to date.

We analyzed the TTVs using a likelihood analysis pipeline developed by Robnik & Seljak (2020): this analysis models the non-Gaussian noise of Kepler data, and performs a Fourier-based Gaussian process analysis, fitting star variability parameters with the optimal prior power spectrum jointly with the planet parameters. The resulting outputs are TTVs and their associated uncertainties. We use the long-cadence data (30 minute sampling), supplementing with the short-cadence data when available (transits 4–6 of Kepler-90g and transit 3 of Kepler-90h). We use the same analysis to output the timing duration variability (TDV) of each event and the associated error, for a total of 18 data points, 9 TTVs, and 9 TDVs.

The TTV and TDV observations, the errors, together with the model prediction as a function of the orbital parameters, allow us to form a data likelihood. Given the data likelihood of observations, as well as the priors on the orbital parameters, the goal is to determine their posterior distribution.

Solving this inverse problem is challenging because there are 13 parameters, relatively few data points (9 TTVs and 9 TDVs), and the expensive-to-evaluate and highly nonlinear dynamical model. TTVs can yield useful constraints even with such few transits because TTVs can vary by orders of magnitude even with slight changes in orbital parameters (Rodriguez et al. 2018). In contrast, TDVs are mainly sensitive to secular effects and provide more limited information about orbital parameters.

The main challenge is finding the solution: most of the prior space is very far from the solution and the current full dynamical N-body integration models are expensive and difficult to use in inverse problem optimization. Analytic gradients were not available for TTVs until very recently (Agol et al. 2020) and due to the high nonlinearity of the model, the local gradient direction may not be very useful even if it were available. A dynamical fit to the sinusoidal small TTV signals usually requires many transits to determine a unique solution. As a consequence, a full dynamical analysis has predominantly been performed on short-period inner planets (P < 100 days) with many transits (Hadden & Lithwick 2017).

Since orbit fitting using the full N-body dynamics is computationally expensive, many works turn to analytic TTV predictions (Lithwick et al. 2012), which often suffice for the small TTVs. Analytic methods are however not very useful for our TTV analysis: these methods are based on first-order perturbative expansion in terms of eccentricity, characterizing TTVs as sinusoidal signals, while the observed TTV signals are nonsinusoidal (Figure 1). We will argue that this is a result of the strong planetary interactions between Kepler-90g and 90h: 90g has a 25 hr TTV anomaly, 0.5% of the orbital period, while typical observed planetary TTVs are seconds to minutes. The small errors of individual TTVs, typical of long-period transits, can make it easier to observe the departures from the sinusoidal model despite the fewer transits. As a result, the analytic TTV calculations are not accurate enough for the orbit fitting. However, while the full dynamical analysis is difficult, when applied to large TTVs it also provides a unique opportunity to go beyond the sinusoidal analysis and measure the angular parameters that cannot be otherwise constrained.

Zoom In Zoom Out Reset image size

While any existing body in the system may perturb 90g and induce TTVs, a 25 hr excursion is too large to be caused solely by a small object like an asteroid or a moon, or distant objects such as the inner planets (Cabrera et al. 2013). The most likely source of this anomaly is the neighboring gas giant, Kepler-90h. Thus, the simplest hypothesis is to model the TTVs with interactions between 90g and 90h, and we first model the TTVs with three bodies: the host star, 90g, and 90h. We verified through numerical experiments that Kepler-90g and 90h do not induce observable TTVs on the interior planets, consistent with the theoretical prediction (Agol et al. 2005). Similarly, the inner planets do not produce significant TTVs on 90g and 90h (see Appendix A). It is also recognized in Cabrera et al. (2013) that Kepler 90e/f and 90g/h are dynamically decoupled, so it is safe to ignore the inner planets when studying the 90g/h pair.

To obtain the model prediction, we use a nonrelativistic N-body integration code TTVFast (Deck et al. 2014), as relativistic effects are negligible for such long periods. Given one set of orbital parameters, this program integrates over 1500 days in 0.01 day step size and detects all the transiting events. Then, we compare the modeled TTVs with the observed signal to evaluate the likelihood of one particular model, and we run an optimization routine to find the set of planetary parameters that maximize the posterior. We describe the details of the optimization procedure in Appendix A.

We obtain the best-fit solution that greatly reduces the TTV residuals from up to 25 hr to less than 1 minute (Figure 1). We searched extensively for additional solutions, by scanning over the parameter space on fine grids (see Appendix A), and by running the dynesty sampler on an uninformative prior (same as the one reported in Table 1; Speagle 2020). Both methods converged to the same best solution, and the second best one has an unacceptable χ² > 100 for 5 degrees of freedom (18 data points and 13 free parameters) so to the best of our knowledge the solution is unique.

Table 1. Summary of Parameter Priors and Posteriors

Parameter	90g Prior	90h Prior	90g Mean ±1σ	90h Mean ±1σ
m(M⊕)	(5, 50)	(5, 500)	${15.0}_{-0.8}^{+0.9}$	${203}_{-5}^{+5}$
P(day)	(209, 212)	(330, 333)	${210.48}_{-0.05}^{+0.05}$	${331.657}_{-0.008}^{+0.009}$
e	(0, 0.15)	(−0.075, 0.075) a	${0.049}_{-0.007}^{+0.011}$	$-{0.011}_{-0.003}^{+0.002}$
i(°)	(89.8, 90.0)	(89.8, 90.0)	${89.92}_{-0.01}^{+0.03}$	${89.927}_{-0.007}^{+0.011}$
Ω(°)	0 b	(−180, 180) a	0 b	$-{2}_{-6}^{+6}$
ω(°)	(0, 180)	(−180, 180) a	${90}_{-20}^{+20}$	$-{3}_{-2}^{+3}$
ω+M(°)	(0, 360)	(0, 360)	${198}_{-2}^{+2}$	${297}_{-1}^{+1}$

Notes. The orbital parameters used in the analysis: m is the planetary mass in Earth mass unit; P is the orbital period in days; e is eccentricity; i is orbital inclination; Ω is the mutual inclination; ω is the argument of the periapsis; M is mean anomaly. All the angles are measured in degrees. Details about the choice of priors can be found in Appendix A. The middle columns show the prior range of Kepler-90g and 90h planetary parameters, respectively. The right columns contain the mean and 1σ errors from the posterior distribution obtained through MCMC sampling.

^a The marked 90h parameters are measured from the respective parameter of 90g to improve MCMC sampling if strongly correlated. ^b Ω for 90g is set to 0, defining the x-axis.

Download table as:  ASCIITypeset image

From the Markov Chain Monte Carlo (MCMC) posterior analysis we find the best-fit masses are ${15.0}_{-0.8}^{+0.9}$ M_⊕ and ${203}_{-5}^{+5}{M}_{\oplus }$ for 90g and 90h, respectively, consistent with the TTVs in 90g being significantly larger than in 90h.

The total χ² is 9.19 (reduced χ² = 1.8), and is a good fit to the data. We perform a posterior analysis of the two-planet model, which is numerically challenging due to the highly correlated orbital parameters. We reparameterized the orbital parameters to break the degeneracy and speed up the convergence. Details are described in Appendix A. The posterior means and variances were obtained using the MCMC emcee (Foreman-Mackey et al. 2013) package and are shown in Table 1. The posterior plots reveal some remaining strong correlations between the parameters. However, because of the extreme sensitivity of TTVs and TDVs to these parameters, the marginalized posteriors are narrow relative to the prior for all of them.

The near-perfect fit we obtained suggests that there is no need to seek further solutions with a third body. Nevertheless, we describe such a search in Appendix A, confirming that no solutions exist that do not require mutual perturbations of the two planets as the main origin of the TTVs.

The long-term stability of Kepler-90g and 90h system with a large TTV can be qualitatively understood in terms of the so-called Hill stability, which quantifies how close the two planets can be while still orbiting a star. A sufficient condition for a coplanar two-planet system with initially circular orbits to be stable is ${\rm{\Delta }}\gt 2.4{({\mu }_{1}+{\mu }_{2})}^{1/3}$, where μ₁ and μ₂ are the planet to star mass ratios of the inner and outer planet, respectively, and Δ = a₂/a₁ − 1, where a₁ and a₂ refer to the semimajor axis (Gladman 1993). In the Kepler-90g and 90h system, Δ = 0.35, and the right-hand side gives 0.19, so the system is Hill stable.

A unique and accurate solution to all planetary parameters enables a quantitative long-term analysis of the orbits. We ran an N-body integration from today's epoch using the WHFast integrator in the rebound package (Rein & Liu 2012). The time span is chosen to be 2 billion years, the same as the age of the star. We started from the best-fit solution and integrated using a 2 day step size, 1% of the shortest period in the system. As demonstrated in Figure 2, the two-planet orbit remains stable during the course of integration. Due to the strong planetary interaction, some of the orbital parameters are not constant. We explored variations of the orbital parameters by sampling from their posterior and found that they do not significantly affect the long-term orbit integration solution.

Zoom In Zoom Out Reset image size

We monitor the evolution of the semimajor axis, eccentricity, apsidal angle (Δω = ω_h − ω_g ), and resonant angle of Kepler-90g and 90h during the 2 billion years of integration (Figure 2). The resonant angle is defined as ϕ_g,h = 3λ_h − 2λ_g − ϖ_g,h , where λ stands for the mean longitude, ϖ is the longitude of the periapsis, and the subscripts denote planet 90g or 90h. If ϕ oscillates around zero, it is librating and the system experiences orbital lock. In our case ϕ circulates, suggesting that Kepler-90g and 90h are not locked in resonance. ⁴ We also performed a relatively short integration of 1 million years on 100 nearby solutions drawn from the MCMC posterior. All the solutions show circulating ϕ_g and ϕ_h throughout the integration.

Kepler observations show there is an excess of planet pairs observed just wide of resonances (Lithwick & Wu 2012; Fabrycky et al. 2014; Choksi & Chiang 2020). Kepler-90g and 90h also have a period ratio slightly wide of 2:3. This may be explained by the resonant repulsion theory: when the planets are formed in disk, eccentricity damping drives planet pairs wide of resonance, while they remain in resonance lock (Choksi & Chiang 2020). This process concludes after the first 10 million years and what happens afterward is unknown. In the Kepler-90g and 90h system, we observe a larger eccentricity than the 1% upper limit predicted by this theory, and the two planets are not in a resonance lock (Figure 2). One possible explanation is that planetary interactions raise the eccentricities and break the resonance lock after the disk dissipates.

The orbital parameter solution also reveals several features that may help preserve the long-term stability, withstanding strong perturbations. The eccentricities are small, 0.045 and 0.035, which means that the orbits are near-circular. However, 0.045 is significantly larger than 0.03, which was assumed to be an upper limit for the long-term stability (Cabrera et al. 2013).

The best solution is close to coplanar (mutual inclination is −2° ± 6°), and the arguments of periapsis are aligned (Δω librates over ±20°). These factors bring the two planet's perigees closer to each other and enhance the gravitational pull when they encounter. Dynamical studies show that in coplanar configuration, apsidal resonance can be produced by secular interactions among resonant planets, which also excite the eccentricities (Chiang et al. 2001; Chiang & Murray 2002). The apsidal alignment and relatively high eccentricities of the Kepler-90g and 90h system may be explained by the same mechanism.

We monitor the TTV evolution for 1000 years: Figure 3 shows that the TTV is very nonsinusoidal and dominated by secular variations. The orbital stability of Kepler-90g and 90h is not caused by orbital locks, but by the low eccentricity, small mutual inclination, compact orbits, and aligned periapsis.

Zoom In Zoom Out Reset image size

It is worth exploring if alternative solutions can be found without resorting to mutual perturbations between 90g and 90h. We tested alternative hypotheses that introduce a new seen or unseen body orbiting Kepler-90, such as a planet, asteroid, or a Trojan planet, and searched over the full seven-dimensional parameter space. The best planet solution gives a total χ² = 20 for 18 data points and 20 model parameters. It is a worse fit to the data and does not eliminate the mutual perturbations between 90g and 90h as the main origin of TTVs. As a consequence, hypothesis testing with Bayes factor shows a large Occam's razor penalty for this solution, due to the large prior in the seven-dimensional parameter space of the new body.

We also explored an exomoon orbiting 90g, assuming circular orbits aligned with 90g orbit, which only require three extra parameters: mass, period, and phase. We assign the best-fit 90g parameters to the 90g-moon barycenter and require the moon mass to be lower than that of 90g. The lower bound of the moon period is determined by the Roche limit of 90g, which is 0.09 days, and the upper bound is 30 days, at 0.4 Hill radius. The moon phase is allowed to vary freely. Marginalizing over 16 parameters, we obtain a total χ² = 9.51, which is a good fit to the data, but the parameters for the two planets are nearly unchanged relative to the solution without the exomoon, and again the Bayes factor strongly disfavors this solution. We thus conclude that the two-planet solution we found is strongly favored and there is no need to add a third body, and there are no found solutions that do not involve the mutual perturbations between the two planets as the main origin of observed TTVs.

While Kepler-90h has mass and density typical of giant gas planets, the best-fit mass of Kepler-90g is 15 ± 0.9M_⊕, and the radius is 8.1 ± 0.8R_⊕ (Santerne et al. 2016). Combining the two we obtain the apparent 90g density of 0.15 ± 0.05 g cm⁻³, which suggests it is among the lowest-density planets found to date, most of which have been found with TTVs in the shorter-period planets (Jontof-Hutter et al. 2014; Masuda 2014; Ofir et al. 2014). One possible explanation is that it has a dusty atmosphere that inflates the observed radius (Wang & Dai 2019). Super-puffs typically reside close to the host star, less than 0.3 au, but in Lee & Chiang (2016) it is argued that super-puffs likely acquire their atmosphere through nebular accretion at around 1 au, where the disk gas is cooler and less dusty, and then migrate to the inner orbits where they are found. Kepler-90g is at 0.7 au, relatively close to the expected formation site, with the resonant repulsion and subsequent stable orbit responsible for not migrating inward. An alternative explanation is that these planets have large optically thick rings, which could give an appearance of a larger planet. We have searched for a signature of a tilted ring in the transit data (Heising et al. 2015), but found no evidence of it. This however does not mean that this explanation is ruled out, as the data do not have sufficient sensitivity to distinguish between most of the ring versus no-ring solutions.

Our analysis suggests that even with just nine transits between the two planets, a full dynamical analysis of large TTVs can constrain all of the orbital parameters for the two planets, and break the degeneracies present in low TTVs. While the low number of transits requires a difficult high dimensional optimization, the solution we found is unique, leading to an unprecedented precision of mass and orbital parameters determination, and enabling dynamical analyses such as the long-term evolution of the system. A systematic search for TTVs in Kepler data has revealed several candidates with periods above 100 days (Holczer et al. 2016; Ofir et al. 2018), but to date, there has been no systematic search of large TTVs. If more systems with large TTVs can be found it would open up a new window into the study of exoplanet formation and evolution models, and a new way to characterize exoplanet demographics.

We thank Eugene Chiang, Joshua Winn, and Daniel Tamayo for constructive discussions and comments. This material is based upon work supported by the National Science Foundation under grant Nos. 1814370 and NSF 1839217, and by NASA under grant No. 80NSSC18K1274. We acknowledge Ad futura Slovenia for supporting J.R.'s MSc study at ETH Zürich.

Software: ttvfast (Deck et al. 2014), rebound (Rein & Liu 2012), emcee (Foreman-Mackey et al. 2013).

A.1. Data Preparation and Analysis

We first describe how the TTV and TDV data and its errors are extracted from the Kepler lightcurves. We apply the method developed in Robnik & Seljak (2020). We use PDCSAP flux of the Kepler data ⁵ processed through the Pre-search data conditioning module (Jenkins et al. 2017), meaning that long-term trends and most of the systematics have already been removed. We use long-cadence and, where available, short-cadence lightcurves with 29.4 and 1 minute spacings, respectively.

A lightcurve is modeled as a sum of the loss of light due to the transit U, stellar variability s, and noise n, composed of the Gaussian white noise and the non-Gaussian outliers. The planet j is parameterized by the amplitude of its transits A_j and for each transit, i by its duration D_ij and epoch T_ij . Each transit is calculated by integrating the quadratic limb-darkening model for the star's intensity profile over the planet's shadow and over the time spacing interval. Contributions of different transits are then added together.

We model stellar variability as a Fourier Gaussian process. A Fourier transformation introduces a convenient basis $s(\nu )\,={ \mathcal F }(s(t))(\nu )$. Two point correlations between the different Fourier components s(ν) vanish by stationarity, so it suffices to impose a hyperprior power spectrum P(ν) = 〈∣s(ν)∣²〉. The power spectrum is learned from the data in an iterative process where the largest planets are first found and removed and the power spectrum is updated to search smaller planets.

Noise is not correlated but is non-Gaussian: it is composed of the Gaussian central part and a small number of the outliers that deviate significantly from the average. A noise probability distribution can be modeled as a mixture model, whose parameters can be learned by examining residuals averaged over multiple stars,

$$\begin{eqnarray}&&n(t)=y(t)-s(t| s(\nu ))-U(t| T,D,A),\end{eqnarray} \tag{ A1 }$$

and are distributed according to the noise probability distribution p(n). Parameters can be extracted by minimizing the negative log-likelihood function

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal L }(T,D,A,s(\nu )) & = & -2\displaystyle \sum _{t}\mathrm{ln}p(n(t))\\ & & +\ \displaystyle \sum _{\nu }\left[\displaystyle \frac{| s(\nu ){| }^{2}}{P(\nu )}+\mathrm{ln}P(\nu )\right].\end{array}\end{eqnarray} \tag{ A2 }$$

The planet model for Kepler-90g and Kepler-90h depends on T, D, and A, which together contain 20 parameters, while s(ν) accounts for an additional thousands of parameters. Optimization is tractable if we iterate on minimization with respect to the stellar parameters at fixed planet's parameters and vice versa. The final result are the values of the TTVs, TDVs, and their errors, as well as the latent parameters. The procedure is near optimal, and has been shown to give smaller and more reliable errors relative to the other existing pipelines (Robnik & Seljak 2020). Posterior distribution of the parameters around their optimal values is somewhat non-Gaussian, but we summarize it by a single symmetrized 1σ error to make further analysis tractable.

A.2. Parameter Space Description and Choice of Priors

The simplest system we model is that of the hosting star and two planets, 90g and 90h. Cabrera et al. (2013) report the Kepler-90 star mass as 1.2 ± 0.1M_⊙. We fix the stellar mass to 1.2 solar mass in the optimization because TTVs are only sensitive to the mass ratio between the perturber and the hosting star.

In the two-planet scenario, there are seven parameters for each planet: planet mass, orbital period, inclination, the argument of ascending node, the argument of periapsis, and the mean anomaly at the epoch (HJD 2454833). We fix the argument of ascending node of Kepler-90g to zero, which effectively defines the x-axis, so we have 13 free parameters to optimize.

We determine the parameter priors based on the physical knowledge about this system. We anticipate the 90h mass to be much larger than that of 90g since the amplitude of the TTVs from 90g is 50 times those of 90h, which directly reflects the mass ratio between the two planets. In the fit, the 90g mass is allowed to vary from 5 to 50 Earth masses, and the 90h mass range goes from 5 to 500 Earth masses.

To determine the prior on eccentricities we conduct the long-term integration using the rebound package (Rein & Liu 2012), to find the upper limit on eccentricities imposed by the dynamical stability. From a preliminary run, we observe a strong correlation between Kepler-90g and 90h eccentricities, so we assume that they are the same in the high e integration. We first draw a sample from the posterior distribution with e = 0.06 and ${\chi }_{{ttv}+{tdv}}^{2}=10.6$. Then, we manually change e and run long-term integration. For high-eccentricity orbits (e > 0.3), the system is unstable in 10⁵ years, much shorter than the age of the system. However, for low-eccentricity orbits (e < 0.15), the stability over millions of years is almost guaranteed, and does not help exclude regions in the parameter space (Figures 4 and 5). We thereby restrict the eccentricity of both planets to [0, 0.15].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The hard cutoff for inclinations is set to [89.8, 90]°. The prior range of the remaining angular arguments are determined in the following manner: the 90g argument of ascending node is fixed to 0° such that the initial 90g orbital plane defines the x-axis. The 90h argument of the ascending node, which represents the mutual inclination between the two orbital planes, is allowed to vary from −180° to 180°. The argument of periapsis of both planets is allowed to vary from −180° to 180° since we have no prior knowledge about this parameter. Finally, the initial guess of the mean anomaly at the epoch is obtained from the first transit time reported by Cabrera et al. (2013): Kepler-90g Epoch (HJD 2454833) = 147.0364; Kepler-90h Epoch (HJD 2454833) = 140.49631. In the fit we still allow it to vary from 0° to 360°.

A.3. Likelihood Function

The likelihood function we use in the optimization consists of two parts: the TTV term and TDV term,

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{log}L({\boldsymbol{p}})=\displaystyle \sum _{{\rm{j}}}\displaystyle \sum _{{\rm{i}}}-\displaystyle \frac{1}{2}\mathrm{log}\left(2\pi {\sigma }_{T,{ij}}^{2}\right)-\displaystyle \frac{{\left({\hat{T}}_{{ij}}-{T}_{{ij}}\right)}^{2}}{2{\sigma }_{T,{ij}}^{2}}\\ -\ \displaystyle \frac{1}{2}\mathrm{log}(2\pi {\sigma }_{D,{ij}}^{2})-\displaystyle \frac{{\left({\hat{D}}_{{ij}}-{D}_{{ij}}\right)}^{2}}{2{\sigma }_{D,{ij}}^{2}}.\end{array}\end{eqnarray} \tag{ A3 }$$

Here T_ij represents the observed transit times of the ith transit and the jth planet, ${\hat{T}}_{{ij}}$ is the corresponding value predicted by the N-body simulation, and σ_T,ij is the observational error of the transit times. Similarly, ${D}_{{ij}},{\hat{D}}_{{ij}}$, and σ_D,ij are the observed transit duration, predicted transit duration, and observed duration error of the ith transit and the jth planet. We have two planets so j = 1,2. i runs from 1 to 6 for Kepler-90g (j = 1) and 1 to 3 for Kepler-90h (j = 2).

We compute ${\hat{T}}_{{ij}},{\hat{D}}_{{ij}}$ given a set of initial parameters using the TTVFast package (Deck et al. 2014). While the recommended step size in the TTVFast is 1/20 of the shortest period, or 10 days in our case, we make a rather conservative choice of 0.01 day to ensure accuracy. We explore the consistency of TTVFast by running the integrator at different step sizes (Δt = 10⁻⁵, 10⁻⁴, 10⁻³... day) and monitor the resulting χ². We see no observable effect on χ² below Δt = 1 day.

A.4. Optimization Strategy

A.5. Point Search

A.6. Grid Search

A.7. Posterior Distribution

A.8. Inner Planets

The project was designed by Y.L, J.R., and U.S. J.R. reduced the TTV data; Y.L. performed the data analysis. The paper was written by Y.L., J.R., and U.S.

The Kepler stellar flux data are obtained from a public data repository https://exoplanetarchive.ipac.caltech.edu/bulk_data_download/. The reduced TTVs and TDVs can be found in data behind Figure 1.

The TTV data reduction code is publicly available in https://github.com/JakobRobnik/Kepler-data-analysis under an Apache-2.0 license, and a copy is preserved on Zenodo: doi:10.5281/zenodo.4540614.