CHARACTERIZATION OF THE KOI-94 SYSTEM WITH TRANSIT TIMING VARIATION ANALYSIS: IMPLICATION FOR THE PLANET–PLANET ECLIPSE

2.1.1. Data Processing

2.1.2. Transit Parameters

Before analyzing the TTV signals, we revise the transit parameters of KOI-94c, KOI-94d, and KOI-94e obtained by the Kepler team (Table 2) so that they are consistent with the light curves obtained in Section 2.1.1. Here we first use the parameters publicized by the Kepler team to phase-fold the observed transit light curves, and then refit those phase curves to obtain the revised transit parameters.

Table 2. Parameter Values Estimated by Other Authors

Parameter	KOI-94c	KOI-94d	KOI-94e
Transit parameters determined by the Kepler teama
t0 (BJD − 2454833)	138.00718 ± 0.00093	132.74047 ± 0.00019	161.23998 ± 0.00079
P (days)	10.423707 ± 0.000026	22.343001 ± 0.000011	54.31993 ± 0.00012
a/R⋆	15.70 ± 0.37	26.10 ± 0.62	47.2 ± 1.1
Rp/R⋆	0.02544 ± 0.00012	0.06856 ± 0.00012	0.04058 ± 0.00013
b	0.019 ± 0.048	0.305 ± 0.014	0.387 ± 0.014
Parameters determined by Hirano et al. (2012)
Ω (deg)b	$-6_{-11}^{+13}$	⋅⋅⋅	$-5_{-11}^{+13}$
u1	0.10 ± 0.06
u2	0.61 ± 0.08
Parameters determined by Weiss et al. (2013)
m (M⊕)	$15.6_{-15.6}^{+5.7}$	106 ± 11	$35_{-28}^{+18}$
e	0.43 ± 0.23	0.022 ± 0.038	0.019 ± 0.23

Notes. ^aData from MAST archive http://archive.stsci.edu/kepler/. ^bIn this paper, we define the reference direction so that the spin–orbit angle λ measured by the Rossiter–McLaughlin effect be equal to Ω for the orbital inclination in the range [0, π/2] (Fabrycky & Winn 2009). With this choice, the reference line points to the ascending node of a virtual circular orbit whose angular momentum is parallel to the stellar spin vector.

Download table as:  ASCIITypeset image

In the first step, we fit each of the detrended light curve centered at the transit (for ∼1.7 times its duration) to obtain the times of transit centers t_c, using a Markov chain Monte Carlo (MCMC) algorithm. Here we use a light curve model by Ohta et al. (2009), and fix a/R_⋆, R_p/R_⋆, and b to the values obtained by the Kepler team, assuming e = 0. We model the limb-darkening using a quadratic law in Equation (7) and adopt the limb-darkening coefficients u₁ and u₂ obtained by Hirano et al. (2012) (all these parameters are summarized in Table 2). Since the detrend procedure in Section 2.1.1 can remove only the out-of-transit outliers, we also exclude in-transit 5σ outliers of this fit, if any, and fit the light curve again. Using the series of t_c obtained in this way, we construct the phase-folded transit light curve for each planet.

As the second step, we fit the resulting phase-folded transit light curve for a/R_⋆, R_p/R_⋆, b, u₁, and u₂ using the same light curve model as above. In this way, we obtain the revised values of the set of parameters shown in Table 3 and the corresponding best-fit light curves (Figures 2–4).⁷ In this fit, all the parameters converge well in the case of KOI-94d. In contrast, a/R_⋆ and b of KOI-94c and KOI-94e do not converge well moving back and forth between several local minima in a strongly correlated fashion, because they show smaller transit depths and the ingresses/egresses of their transits are less clear. For this reason, we impose an additional constraint that all the planets share the same host star: we convert the well-constrained a/R_⋆ for KOI-94d into stellar density ρ_⋆ via ρ_⋆ ≈ (3π/GP²)(a/R_⋆)³ (Seager & Mallén-Ornelas 2003), and calculate the corresponding values and uncertainties of a/R_⋆ for KOI-94c and KOI-94e. The phase curves of KOI-94c and KOI-94e are fitted with prior constraints centered on these values and with Gaussian widths of their uncertainties. With this prescription, all the parameters of KOI-94c and KOI-94e converge well. Therefore, it does not make sense here to discuss the consistency of ρ_⋆ calculated from the transit parameters to check the possible false positives. It is important to note, however, that the limb-darkening coefficients for each planet obtained individually are consistent within their 1σ error bars (those obtained by Hirano et al. (2012) are different from our values because they fixed smaller R_p/R_⋆ for KOI-94d; see Table 2). This supports the notion that these three planets are indeed revolving around the same host star.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3. Revised Transit Parameters Obtained from Phase-folded Light Curves

	KOI-94c	KOI-94d	KOI-94e
a/R⋆	$15.7798_{-0.0016}^{+0.0016}$	$26.24_{-0.19}^{+0.20}$	$47.4400_{-0.0093}^{+0.0094}$
Rp/R⋆	$0.025673_{-0.000075}^{+0.000074}$	$0.07029_{-0.00015}^{+0.00014}$	$0.04137_{-0.00013}^{+0.00012}$
b	$0.021_{-0.015}^{+0.020}$	$0.299_{-0.025}^{+0.022}$	$0.3840_{-0.0055}^{+0.0048}$
u1	0.40 ± 0.06	0.40 ± 0.02	0.36 ± 0.07
u2	$0.10_{-0.09}^{+0.08}$	0.14 ± 0.03	$0.19_{-0.10}^{+0.11}$
$\chi ^2/\mathrm{{\rm dof}}$	25843/23611	15473/13766	6822/5934

Notes. The quoted error bars denote 1σ confidence intervals obtained from the posterior distributions. The values of a/R_⋆ and b for KOI-94c and KOI-94e are determined with the prior information about the stellar density based on the result for KOI-94d (see the main text).

Download table as:  ASCIITypeset image

2.1.3. TTV Signals

Fixing a/R_⋆, R_p/R_⋆, b, u₁, and u₂ at the values in Table 3, we refit the transit light curves of each planet to find the values of t_c given in Tables 8–10 in Appendix A, with the values of reduced χ² and 1σ upper/lower limits obtained from the posterior. The column labeled as O−C tabulates the residuals of a linear fit to t_c versus transit number, in which linear ephemerides in Table 4 are extracted. The values of reduced χ² of the linear fits in this table indicate the significant deviations of the transit times from the linear ephemerides (i.e., TTVs) for all the three planets, as shown in Figure 5. Note that the TTV of KOI-94c shows the modulation with the period of (1/P_c − 2/P_d)⁻¹ ≃ 155 days, which clearly comes from near 2:1 resonance of KOI-94c and KOI-94d as pointed out by Xie et al. (2013) (see also Appendix B).

Zoom In Zoom Out Reset image size

Table 4. Linear Ephemerides of KOI-94c, KOI-94d, and KOI-94e

Parameter	KOI-94c	KOI-94d	KOI-94e
t0 (BJD − 2454833)	138.00826 ± 0.00038	132.74103 ± 0.00012	161.23888 ± 0.00046
P (days)	10.4236888 ± 0.0000053	22.3429698 ± 0.0000036	54.319849 ± 0.000035
χ2/dof	10.4	13.7	29.2

Download table as:  ASCIITypeset image

We numerically analyze the TTV signals in Figure 5 to constrain the masses, eccentricities, and longitudes of periastrons of KOI-94c, KOI-94d, and KOI-94e (nine parameters in total). In this section, we fix the mass of KOI-94d at m_d = 106 M_⊕, the best-fit RV value obtained by Weiss et al. (2013). Since this value is only marginally consistent with m_d = 73 ± 25 M_⊕ obtained by Hirano et al. (2012) based on the out-of-transit RVs taken by the Subaru telescope, we also investigate the case of m_d = 73 M_⊕ here.

2.2.1. Calculation of the Simulated TTV

Orbits of the planets are integrated using the fourth-order Hermite scheme with the shared time step (Kokubo & Makino 2004). The time step (typically ∼0.025 days) is chosen so that the fractional energy change due to the integrator during ∼1000 days integration should always be smaller than 10⁻⁹. All the simulations presented in this section integrate planetary orbits beginning at the same epoch T₀(BJD) = 2455189.1703 (the first transit time of KOI-94d), until BJD = 2456133.0 (approximately the last transit time of KOI-94d we analyzed).

For each planet, the initial value of the orbital inclination i is fixed at the value obtained from a/R_⋆ and b (see Table 3), assuming that e = 0. Here, all the i values are chosen in the range [0, π/2].⁸ The values of Ω_d and Ω_e are fixed at those in Table 2. Since we have no information on the nodal angle of KOI-94c, we assume Ω_c = 0°. Initial semi-major axes are calculated via Kepler's third law with M_⋆ = 1.25 M_☉ (Hirano et al. 2012), orbital periods in Table 4, and planetary masses adopted in each simulation. The phases of the planets are determined from the transit ephemerides: for each planet, we convert the transit time closest to T₀ into the sum of the argument of periastron and mean anomaly ω + f, taking account of the non-zero eccentricity if any, and then move it backward in time to T₀, assuming a Keplerian orbit.

The mid-transit times of each planet are determined by minimizing the sky-plane distance D between the star and the planet, where the roots of the time derivative of D are found by the Newton–Raphson method (Fabrycky 2010). Then these transit times are fitted with a straight line and thereby the TTVs (= residuals of the linear fit), as well as the linear ephemeris (P and t₀), are extracted. We compute the chi squares of the simulated TTVs obtained in this way as

$$\begin{equation} \chi ^2_{j} = \sum _{{{i:{\rm observed}\atop{\rm transits}}}} \left[ \frac{{\rm TTV}^{(j)}_{\rm sim}(i) - {\rm TTV}^{(j)}_{\rm obs}(i)}{\sigma ^{(j)}_{\rm obs}(i)} \right]^2, \quad (j = {\rm c}, {\rm d}, {\rm e}), \end{equation} \tag{ 1 }$$

where ${\rm TTV}^{(j)}_{\rm sim}(i)$ and ${\rm TTV}^{(j)}_{\rm obs}(i)$ are the i-th values of simulated and observed TTVs of planet j, respectively, and $\sigma ^{(j)}_{\rm obs}(i)$ is the observational uncertainty of the i-th transit time of planet j.

Note that we do not fit the transit times directly but only the deviations from the periodicity in our analysis, assuming that they provide sufficient information on the gravitational interaction among the planets. Indeed, although the initial values of semi-major axes are chosen to match the observed periods, periods derived from the simulations are different typically by ∼0.01 days. This is because strong gravitational interaction among massive, closely packed planets in this system causes the oscillations of their semi-major axes with amplitudes dependent on the parameters of the planets adopted in each run. We will show that this simplified method still yields reasonable results in the last part of Section 2.4.

2.2.2. Estimates for the TTV Amplitudes

Before directly fitting the observed TTV signals, we evaluate the contribution from each planet to the TTVs of KOI-94c, KOI-94d, and KOI-94e. We divide the four planets into six pairs and integrate circular orbits for each pair using the best-fit masses by Weiss et al. (2013) listed in Table 2. Semi-amplitudes of the resulting TTVs of KOI-94c, KOI-94d, and KOI-94e are shown in Table 5. Considering the uncertainties of transit times listed in Tables 8–10 (typically $1.4\,\mathrm{min{\rm utes}}$, $0.3\,\mathrm{min{\rm utes}}$, and $0.9\,\mathrm{min{\rm utes}}$ for KOI-94c, KOI-94d, and KOI-94e, respectively), this result indicates that KOI-94b has negligible contribution for the TTVs of the other three. In the following analysis, therefore, we integrate the orbits of the other three planets (KOI-94c, KOI-94d, and KOI-94e) only. We also find that the TTVs of KOI-94c and KOI-94e are mainly determined by the perturbation from the neighboring planet KOI-94d, while that of KOI-94d depends on both of its neighbors. Such dependence is naturally understood from the architecture of this system (see Figure 1). Consequently, each planet's TTV mainly depends on the parameters listed in the rightmost column of Table 5, where we define $\boldsymbol {e}_j = (e_j \cos \varpi _j, e_j \sin \varpi _j)$ (j = c, d, e). Note that the TTV of each planet is insensitive to its own mass. This is why m_j (j = c, d, e) is not included in the row for planet j.

Table 5. Semi-amplitude of the Simulated TTV (in Units of Minutes) for Each Planet Pair

${\backslashbox[2pc]{{\rm TTV}}{{\rm Pair}}}$	KOI-94b	KOI-94c	KOI-94d	KOI-94e	Major Parameters for TTV
KOI-94c	≲ 0.05	⋅⋅⋅	11	≲ 0.05	md, $\boldsymbol {e}_{\rm c}$, $\boldsymbol {e}_{\rm d}$
KOI-94d	≲ 0.05	0.47	⋅⋅⋅	2.1	mc, me, $\boldsymbol {e}_{\rm c}$, $\boldsymbol {e}_{\rm d}$, $\boldsymbol {e}_{\rm e}$
KOI-94e	≲ 0.05	0.15	0.83	⋅⋅⋅	md, $\boldsymbol {e}_{\rm d}$, $\boldsymbol {e}_{\rm e}$, (mc, $\boldsymbol {e}_{\rm c}$)

Download table as:  ASCIITypeset image

2.2.3. Results

TTV of KOI-94e. Since the TTV of KOI-94e is mainly determined by $\boldsymbol {e}_{\rm d}$ and $\boldsymbol {e}_{\rm e}$ (and m_d, of course, which we fix at the RV value), we fit it first so as to constrain these parameters. We calculate $\chi ^2_{\rm e}$ for |e_dcos ϖ_d|, |e_dsin ϖ_d| ⩽ 0.06 and |e_ecos ϖ_e|, |e_esin ϖ_e| ⩽ 0.25 (which well cover the 1σ regions for these parameters obtained from the RVs) at the grid-spacing of 0.01, fixing e_c = 0 and planetary masses at the best-fit values from the RVs. However, we cannot fit the observed TTV well in both m_d = 106 M_⊕ and 73 M_⊕ cases. The best fit for the former case, which gives $\chi ^2_{\rm e} = 122$ for 4 dof, is shown in Figure 6. For this reason, in addition to the fact that we have only eight transits observed for KOI-94e, we decide not to use the TTV of this planet to constrain the system parameters, but fit only the TTVs of KOI-94c and KOI-94d. The large discrepancy in the amplitudes of simulated and observed TTVs may suggest another source of perturbation which is not included in our model, such as a non-transiting planet or other minor bodies.

Zoom In Zoom Out Reset image size

Grid-search for an Initial Parameter Set. Then we perform the grid-search fit to the TTVs of KOI-94c and KOI-94d to find an appropriate initial parameter set for the following MCMC analysis. Based on the estimates given in Table 5, we fit these TTVs separately as follows. We first fit the TTV of KOI-94c varying $\boldsymbol {e}_{\rm c}$ and $\boldsymbol {e}_{\rm d}$ in |e_jcos ϖ_j|, |e_jsin ϖ_j| ⩽ 0.10 (j = c, d) at the grid-spacing of 0.01, and find one minimum of $\chi ^2_{\rm c}$ for both m_d = 106 M_⊕ and m_d = 73 M_⊕ cases. Next, for all the sets of (e_c, ϖ_c, e_d, ϖ_d) in 2σ (m_d = 106 M_⊕ case) or 1σ (m_d = 73 M_⊕ case) confidence regions around the minimum, we run integrations varying e_ecos ϖ_e and e_esin ϖ_e from −0.1 to 0.1 at the grid spacing of 0.01, m_c from 0 to 24 M_⊕ at the grid spacing of 6 M_⊕, and m_e from 7 to 57 M_⊕ at the grid spacing of 10 M_⊕ (all of these cover the 1σ intervals from RV), to find the set of eight parameters that best fits the TTV of KOI-94d.

MCMC Fit to the TTVs of KOI-94c and KOI-94d. Choosing the above set as initial parameters, we then simultaneously fit the TTVs of KOI-94c and KOI-94d using an MCMC algorithm. In this fit, we use $\chi ^2_{\rm c} + \chi ^2_{\rm d}$ as the χ² statistic. The resulting best-fit parameters and their 1σ uncertainties are summarized in Table 6 for the two choices of m_d (the second and third columns). The best-fit simulated TTVs are plotted in Figures 7 and 8 for KOI-94c and KOI-94d, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 6. Best-fit Parameters Obtained from TTV Analysis

Parameter	Value (md = 106 M⊕)	Value (md = 73 M⊕)	Value (TTV Only)
KOI-94c
mc (M⊕)	$11.8_{-1.5}^{+1.6}$	$13.9_{-2.7}^{+2.7}$	$9.4_{-2.1}^{+2.4}$
eccos ϖc	$0.0329_{-0.0055}^{+0.0047}$	$0.0092_{-0.0050}^{+0.0264}$	$0.0143_{-0.0059}^{+0.0080}$
ecsin ϖc	$-0.0104_{-0.0042}^{+0.0038}$	$-0.0031_{-0.0061}^{+0.0067}$	$0.0045_{-0.0079}^{+0.0091}$
$\chi ^2_{\rm c}$	84	62	56
KOI-94d
md (M⊕)	106 (fixed)	73 (fixed)	$52.1_{-7.1}^{+6.9}$
edcos ϖd	$0.055_{-0.014}^{+0.011}$	$-0.016_{-0.011}^{+0.064}$	$-0.022_{-0.011}^{+0.014}$
edsin ϖd	$0.012_{-0.012}^{+0.011}$	$0.009_{-0.018}^{+0.018}$	$0.008_{-0.018}^{+0.021}$
$\chi ^2_{\rm d}$	66	48	43
KOI-94e
me (M⊕)	$15.9_{-2.2}^{+2.4}$	$12.9_{-2.3}^{+3.0}$	$13.0_{-2.1}^{+2.5}$
eecos ϖe	$0.067_{-0.019}^{+0.014}$	$-0.069_{-0.018}^{+0.120}$	$-0.078_{-0.014}^{+0.021}$
eesin ϖe	$0.042_{-0.017}^{+0.012}$	$-0.022_{-0.016}^{+0.032}$	$-0.025_{-0.014}^{+0.017}$
$(\chi _{\rm c}^2 + \chi _{\rm d}^2) /{\rm dof}$	150/57	110/57	99/56

Download table as:  ASCIITypeset image

Note that uncertainties of e_ccos ϖ_c, e_dcos ϖ_d, and e_ecos ϖ_e are relatively large for m_d = 73 M_⊕ case. This is because the posterior distributions of these parameters have two peaks, the smaller of which lies close to the best-fit value for m_d = 106 M_⊕ case. Considering this fact, the two results are roughly consistent with each other. Nevertheless, a total χ² in m_d = 73 M_⊕ case is smaller by 40 for 57 dof than in m_d = 106 M_⊕ case. This suggests that the TTV alone favors m_d smaller than the RV best-fit value, as will be confirmed in the next subsection.

In order to obtain a solution independent of RVs, we perform the same MCMC analysis of TTVs of KOI-94c and KOI-94d, this time also allowing m_d to float. Since the above analyses suggest that the eccentricities of all the planets are small, we choose circular orbits with the best-fit RV masses as an initial parameter set. The resulting best-fit parameters are summarized in Table 6, and the corresponding best-fit simulated TTVs are shown in Figures 9 and 10. As expected, we find a solution with small eccentricities and with m_d smaller than the RV best-fit value. This solution is similar to that for m_d = 73 M_⊕ case, except that m_d is even smaller.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

While the values of e_d and e_e obtained in our TTV analysis are consistent with the RV values in Table 2, the best-fit e_c obtained from the TTV is ∼1.8σ smaller than the RV best fit (e_c = 0.43 ± 0.23). Considering the marginal detection of KOI-94c's RV and the dynamical stability of the system, however, the TTV value seems to be preferred. In fact, this value is robustly constrained by the clear TTV signal of KOI-94c; in the grid-search performed in Section 2.2.3, we searched the region where e_c ≲ 0.14 to fit the TTV of this planet, but the resulting $\chi ^2_{\rm c}$ strongly disfavored large e_c regions in both m_d = 106 M_⊕ and m_d = 73 M_⊕ cases.

The TTV values of m_c and m_e are consistent with the RV results, but m_e is constrained to a rather lower value than the RV best fit. Using this value, along with the photometric values of R_p/R_⋆ and ρ_⋆, and spectroscopic value of M_⋆, the density of KOI-94e is given by ρ_e ∼ 0.3 g cm⁻³. This implies that it is one of the lowest-density planets ever discovered.

The largest discrepancy arises in the value of m_d, mass of KOI-94d, for which the TTV best-fit value is smaller than the RV value by ∼4σ. The worse fit in the case of m_d = 106 M_⊕ is mainly due to the fact that the observed TTV amplitude of KOI-94c is smaller than expected from this value of m_d. As shown in Table 5, the TTV of this planet is completely dominated by the perturbation from KOI-94d, and so the parameters relevant to this TTV are m_d, $\boldsymbol {e}_{\rm c}$, and $\boldsymbol {e}_{\rm d}$. Table 5 also shows that m_d = 106 M_⊕ leads to the TTV semi-amplitude of ∼11 minutes for e_c = e_d = 0, which is much larger than the observed TTV amplitude of KOI-94c (see Figure 5). As a result, the values of $\boldsymbol {e}_{\rm c}$ and $\boldsymbol {e}_{\rm d}$ are fine-tuned to fit the observed signal, resulting in strict constraints on these parameters. The problem is that these values of $\boldsymbol {e}_{\rm c}$ and $\boldsymbol {e}_{\rm d}$ do not fit the TTV of KOI-94d well. In the grid-search, we first fit the TTV of KOI-94c alone, and the best fit gives $\chi ^2_{\rm c} = 66$ for m_d = 106 M_⊕. However, this value is largely increased in the simultaneous MCMC fit to the TTVs of KOI-94c and KOI-94d (Table 6), which means that these two TTVs cannot be explained with the same set of $\boldsymbol {e}_{\rm c}$ and $\boldsymbol {e}_{\rm d}$. On the other hand, this tension does not exist in m_d = 73 M_⊕ case, in which both of the grid-search and MCMC return the same values for the best-fit $\chi ^2_{\rm c}$. In fact, the analysis using the analytic formulae of TTVs from two coplanar planets lying near j: j − 1 resonance (Lithwick et al. 2012) also supports the above reasoning, suggesting that m_d expected from the TTV of KOI-94c is rather small (see Appendix B).

For these reasons, it is clear that the TTV favors the solution with m_d smaller than the RV best-fit value. It is also true, however, that the RV of KOI-94d is detected with high significance, in contrast to those of the other planets. Indeed, we calculate the RVs using the best-fit TTV parameters and find that the resulting amplitude is much smaller than that observed by Weiss et al. (2013). Since its origin is not yet clear, this discrepancy motivates further investigation of KOI-94 including additional RV/TTV observations.

Finally, we note again that in the above analysis we just fit the TTVs rather than the transit times of KOI-94c and KOI-94d. For this reason, our solution corresponds to the planetary orbits whose periods are slightly different from the actually observed values. One may argue that such an approximate method leads to an incorrect solution. In order to assure that this is not the case, we fit the transit times of the two planets, choosing m, a, ecos ϖ, esin ϖ, and ω + f of KOI-94c, KOI-94d, and KOI-94e at time T₀ as free parameters (fifteen parameters in total). We run two MCMC chains starting from (1) circular orbits with RV best-fit masses and (2) a local minimum reached by the Levenberg–Marquardt method (Markwardt 2009) starting from the TTV best-fit parameters (rightmost column of Table 6). We find that the best-fit values of (m, e, ϖ) obtained in these ways show similar trends as the TTV best fit in Section 2.3, giving comparable reduced χ² values. Namely, the mass of KOI-94d is much smaller than the RV best fit and eccentricities of the three planets are close to zero.

Incidentally, as in the case of the fit to TTVs, the transit times of KOI-94d are less well reproduced than those of KOI-94c. In fact, the fit to transit times gives a better $\chi ^2_{\rm c}$ than that to TTVs because the small linear trend apparent in the lower panel of Figure 9 is removed; on the other hand, $\chi ^2_{\rm d}$ values are not so different in both cases. The difficulty in completely reproducing the orbit of KOI-94d may also indicate the presence of the additional perturber discussed in Section 2.2.3.

A PPE is observed as an increase in the relative flux of a star (or a "bump" in the light curve) during the double transit of two planets. Assuming that the two overlapping planets are spherical and neglecting the effect of limb-darkening, the increase in the relative flux S is given by the area of overlapping region of the two planets divided by that of the stellar disk (Figure 11):

$$\begin{equation} S = \left\lbrace \begin{array}{@{}ll@{}}0 \quad {\rm for} \quad R_{p1} + R_{p2} < d,\\ \frac{1}{2\pi } R_{p1}^2 (2 \alpha - \sin 2 \alpha) + \frac{1}{2\pi } R_{p2}^2 (2 \beta - \sin 2 \beta)\\ \quad {\rm for} \quad |R_{p1} - R_{p2}| < d < R_{p1} + R_{p2},\\ (\min (R_{p1}, R_{p2}))^2 \quad {\rm for} \quad d < |R_{p1} - R_{p2}|, \end{array}\right. \end{equation} \tag{ 2 }$$

where d is the distance between the centers of the two planets in the plane of the sky, and the angles α and β are defined as

$$\begin{equation} \cos \alpha = \frac{R_{p1}^2 + d^2 - R_{p2}^2}{2 R_{p1} d}, \quad \cos \beta = \frac{R_{p2}^2 + d^2 - R_{p1}^2}{2 R_{p2} d}. \end{equation} \tag{ 3 }$$

Zoom In Zoom Out Reset image size

An alternative expression of S for |R_p1 − R_p2| < d < R_p1 + R_p2 can be obtained from the derivative of S with respect to d. Using dα/dd and dβ/dd obtained from Equation (3) and R_p1sin α = R_p2sin β, we obtain

$$\begin{eqnarray} &&\frac{{\rm d} S}{{\rm d} d} = - \frac{1}{\pi } \sqrt{ - d^2 + 2 \left(R_{p1}^2 + R_{p2}^2\right) - \left(\frac{R_{p1}^2 - R_{p2}^2}{d} \right)^2} \end{eqnarray} \tag{ 4 }$$

for |R_p1 − R_p2| < d < R_p1 + R_p2. Equation (4) can be integrated from R_p1 + R_p2 to d by changing the variable from d to

$$\begin{equation} \gamma \equiv \alpha + \beta = \pi - \arccos \left(\frac{R_{p1}^2 + R_{p2}^2 - d^2}{2 R_{p1} R_{p2}}\right). \end{equation} \tag{ 5 }$$

The result is

$$\begin{eqnarray} \nonumber S = \frac{1}{\pi } & \left[ \frac{R_{p1}^2 + R_{p2}^2}{2} \,\gamma - R_{p1} R_{p2} \sin \gamma \right. \\ &\left. - \left(R_{p1}^2 - R_{p2}^2\right) \,\arctan \left(\frac{R_{p1} - R_{p2}}{R_{p1} + R_{p2}} \tan \frac{\gamma }{2}\right) \right]. \end{eqnarray} \tag{ 6 }$$

Here, S is given as a function of a single angle γ. These two expressions of S show that the shape of a bump due to a PPE is solely determined by d as a function of time, which will be derived in the following subsection.

The effect of limb-darkening can be included in our model by multiplying S by a factor that corresponds to the limb-darkening at the position on the stellar disk over which a PPE occurs. We adopt the quadratic limb-darkening law

$$\begin{equation} I(\mu) / I(0) = 1 - u_1 (1 - \mu) - u_2 (1 - \mu)^2, \end{equation} \tag{ 7 }$$

where μ = (1 − r²)^1/2 and r is the radial coordinate on the stellar disk. For a PPE that occurs totally within the stellar disk, the approximation by Mandel & Agol (2002), which is valid for a small planet whose radius is less than about 0.1R_⋆, yields the modified relative increase S' as

$$\begin{equation} S^{\prime } = S \cdot \frac{I(r_*)}{\int _0^1 dr 2r I(r)} = S \cdot \frac{1 - u_1 (1 - \mu _*) - u_2 (1 - \mu _*)^2}{1 - u_1/3 - u_2/6}, \end{equation} \tag{ 8 }$$

where r_* is the distance to the overlapping region. During the whole PPE, r_* is given in terms of d, $\Delta R_{12}^2 \equiv R_{p1}^2 - R_{p2}^2$, and r_j (j = 1, 2), the radial coordinate of the planet j's center, as $r_* = \sqrt{ (r_1^2 + r_2^2)/2 - d^2/4 + (r_2^2 - r_1^2 + \Delta R_{12}^2 /2) \Delta R_{12}^2/2d^2 }$.

Hereafter, we adopt the Cartesian coordinate system (X, Y, Z) centered on the star, where the +Z-axis is chosen in the direction of our line of sight, and X- and Y-axes are in the plane of the sky, forming a right-handed triad. As stated in the note of Table 2, we align the +X axis with the ascending node of a virtual circular orbit whose angular momentum vector is parallel to the stellar spin vector, without loss of generality. Using R for the three-dimensional distance between the planet and the star, the position of a planet is expressed as

$$\begin{eqnarray} X &= R \,[ \cos \Omega \cos (\omega + f) - \sin \Omega \sin (\omega + f) \cos i], \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} &&Y = R \,[ \sin \Omega \cos (\omega + f) + \cos \Omega \sin (\omega + f) \cos i], \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} &&Z = R \sin (\omega + f) \sin i. \end{eqnarray} \tag{ 11 }$$

Suppose that the two planets 1 and 2 are on Keplerian orbits whose semi-major axes are a₁ and a₂, respectively. Neglecting the corrections arising from the non-zero eccentricity, the two-dimensional position vectors $\boldsymbol {r}_j$ (j = 1, 2) of these planets in the plane of the sky can be written as

$$\begin{equation} \boldsymbol {r}_j = \left(\begin{array}{@{}l@{}}a_j \cos \Omega _j \cos (\omega _j + f_j) - b_j \sin \Omega _j \sin (\omega _j + f_j) \\ a_j \sin \Omega _j \cos (\omega _j + f_j)+ b_j\cos \Omega _j \sin (\omega _j + f_j) \end{array}\right). \end{equation} \tag{ 12 }$$

If the transits of the two planets are observed, a_j, b_j, R_pj, and the periods P_j are obtained as in Section 2. In this case, the relative motion of the planets is completely described except for the dependence on the relative nodal angle defined as

$$\begin{equation} \Omega _{21} = - \Omega _{12} = \Omega _2 - \Omega _1. \end{equation} \tag{ 13 }$$

Note that photometric surveys determine the absolute value of b, but not its sign. For a single transiting planet, b is conventionally defined to be positive (or equivalently, i is chosen to be in the range [0, π/2]), because the choice of its sign does not affect the transit signals. For multiple transiting planets, however, a different choice of the relative signs of b corresponds to a different orbital configuration. In this paper, we choose b₁ to be positive (i.e., 0 ⩽ i₁ ⩽ π/2), but allow b₂ to be either positive or negative (i.e., 0 ⩽ i₂ ⩽ π). The sign of b₂, as well as the relative nodal angle Ω₂₁, is determined from the observed data of a PPE event.

During a double transit by planets j = 1 and 2, their phases at time t are given by

$$\begin{equation} \omega _j + f_j = \frac{\pi }{2} + n_j \left(t - t_c^{(j)}\right), \end{equation} \tag{ 14 }$$

where n_j = 2π/P_j is the mean motion and $t_c^{(j)}$ is the central transit time of planet j in this double transit. We then expand Equation (12) to the first order of $n(t -t_c) \sim \mathcal {O}(a^{-1})$ to obtain

$$\begin{equation} \boldsymbol {r}_j = \boldsymbol {v}_j \left(t - t_c^{(j)}\right) + \boldsymbol {r}_0^{(j)}, \end{equation} \tag{ 15 }$$

where

$$\begin{equation} \boldsymbol {v}_j = - a_j n_j \left(\begin{array}{@{}l@{}}\cos \Omega _j \\ \sin \Omega _j \end{array}\right), \quad \boldsymbol {r}_0^{(j)} = -b_j \left(\begin{array}{@{}l@{}}\sin \Omega _j \\ {-} \cos \Omega _j \end{array}\right). \end{equation} \tag{ 16 }$$

These give the distance between the two planets d as⁹

$$\begin{equation} d^2 = | \boldsymbol {r}_2 - \boldsymbol {r}_1 |^2 = v^2 \left(t + \frac{r_0 \cos \theta _0}{v} \right)^2 + r_0^2 \sin ^2 \theta _0, \end{equation} \tag{ 17 }$$

where $\boldsymbol {v} \equiv \boldsymbol {v}_2 - \boldsymbol {v}_1$, $\boldsymbol {r}_0 \equiv \boldsymbol {r}_0^{(2)} - \boldsymbol {v}_2 t_c^{(2)} - \boldsymbol {r}_0^{(1)} + \boldsymbol {v}_1 t_c^{(1)}$, $v \equiv |\boldsymbol {v}|$, $r_0 \equiv |\boldsymbol {r}_0|$, and θ₀ is the angle between $\boldsymbol {v}$ and $\boldsymbol {r}_0$. Thus, the minimum value of d in this double transit

$$\begin{equation} d_{\mathrm{min}} \equiv r_0 \sin \theta _0, \end{equation} \tag{ 18 }$$

is reached at the time

$$\begin{equation} t = t_{\mathrm{min}} \equiv - \frac{r_0 \cos \theta _0}{v}. \end{equation} \tag{ 19 }$$

If d_min < R_p1 + R_p2, a PPE occurs during this double transit for a duration of

$$\begin{equation} \Delta t = \frac{2}{v} \sqrt{ \left(R_{p1} + R_{p2}\right )^2 - d_{\mathrm{min}}^2 }. \end{equation} \tag{ 20 }$$

Equations (18)–(20) can be readily understood by considering the geometry of the PPE shown in Figure 12. The explicit expressions of v, r₀, and cos θ₀ are

$$\begin{eqnarray} v^2 &= a_1^2 n_1^2 + a_2^2 n_2^2 -2 a_1 n_1 a_2 n_2 \cos \Omega _{21}, \end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray} \nonumber r_0^2 &=& \left(a_1 n_1 t_c^{(1)}\right)^2 + \left(a_2 n_2 t_c^{(2)}\right)^2 + b_1^2 + b_2^2\\ \nonumber && + 2\left(a_1 n_1 t_c^{(1)} b_2 - a_2 n_2 t_c^{(2)} b_1\right) \sin \Omega _{21}\\ && - 2\left(a_1 n_1 t_c^{(1)} a_2 n_2 t_c^{(2)} + b_1 b_2\right) \cos \Omega _{21}, \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} \nonumber \cos \theta _0 &=& \frac{1}{vr_0} \left[ -(a_1 n_1)^2 t_c^{(1)} - (a_2 n_2)^2 t_c^{(2)} \right.\\ && + (a_2 n_2 b_1 - a_1 n_1 b_2) \sin \Omega _{21}\nonumber \\ && +\left. a_1 n_1 a_2 n_2 \left(t_c^{(1)} + t_c^{(2)}\right) \cos \Omega _{21} \right]. \end{eqnarray} \tag{ 23 }$$

Zoom In Zoom Out Reset image size

The observed shape of a bump is characterized by its maximum height S_max, central time t_min, and duration Δt. If the bump is not saturated, i.e., d_min > |R_p1 − R_p2|, S_max is uniquely translated into d_min via Equation (2). In this case, we can use Equations (18)–(20), the expressions for the three observables of a bump (d_min, t_min, Δt), to calculate

$$\begin{eqnarray} v &= \frac{2}{\Delta t} \sqrt{(R_{p1} + R_{p2})^2 - d_{\mathrm{min}}^2}, \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} \nonumber r_0^2 &= d_{\mathrm{min}}^2 + (v t_{\mathrm{min}})^2\\ &= d_{\mathrm{min}}^2 + \left(\frac{2 t_{\mathrm{min}}}{\Delta t} \right)^2 \left[ (R_{p1} + R_{p2})^2 - d_{\mathrm{min}}^2 \right], \end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray} \tan \theta _0 &= -\frac{d_{\mathrm{min}}}{v t_{\mathrm{min}}} = -\frac{\Delta t}{2 t_{\mathrm{min}}} \frac{d_{\mathrm{min}}}{\sqrt{(R_{p1} + R_{p2})^2 - d_{\mathrm{min}}^2}}. \end{eqnarray} \tag{ 26 }$$

Furthermore, Equations (21) and (22), the explicit expressions for v and r₀, are rewritten as

$$\begin{eqnarray} \cos \Omega _{21} &= \frac{a_1^2 n_1^2 + a_2^2 n_2^2 - v^2}{2 a_1 a_2 n_1 n_2}, \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray} \sin \Omega _{21} &= \frac{r_0^2 - x_1^2 - x_2^2 + 2 (\boldsymbol {x}_1 \cdot \boldsymbol {x}_2) \cos \Omega _{21}}{2 | \boldsymbol {x}_1 \times \boldsymbol {x}_2 |}, \end{eqnarray} \tag{ 28 }$$

where we define $\boldsymbol {x}_1 =(a_1 n_1 t_c^{(1)}, b_1)$ and $\boldsymbol {x}_2 =(a_2 n_2 t_c^{(2)}, b_2).$ In this way, the relative nodal angle Ω₂₁ can be specified explicitly from (d_min, t_min, Δt) along with the photometrically obtained parameters $(a_j, b_j, n_j, t_c^{(j)}, R_{pj})$, provided that the sign of b₂ is determined. Although we do not present any general procedure to determine the sign of b here, it is possible to do so at least empirically, as described in Section 4. Equation (28) shows that sin Ω₂₁ is only weakly constrained when a₁n₁/b₁ ≃ a₂n₂/b₂, in which case the coefficients in front of sin Ω₂₁ in Equations (22) and (23) are close to zero because $t_c^{(1)} \simeq t_c^{(2)}$.

In the case of d_min < |R_p1 − R_p2|, only the upper limit on d_min can be obtained, and so the entire shape of the bump is required to determine Ω₂₁.

Note that the formulation in this section is also valid even in the presence of non-Keplerian effects. In such a case, the Keplerian orbital elements in this formulation should be interpreted as the osculating orbital elements around a certain double transit.

In order to determine the relative nodal angle crucial in predicting the next PPE, we refit the transit light curve around BJD = 2454211.5, in which the PPE was observed. We model the light curve as described in Section 4 including the effect of limb-darkening to obtain R_pd, R_pe, a_d, a_e, b_d, b_e, $t_c^{\rm (d)}$, $t_c^{\rm (e)}$, and Ω_ed. Here, the same prior as used in fitting the phase-folded light curve is assumed to fix the value of a_e. We adopt the limb-darkening coefficients obtained from the light curve of KOI-94d, u₁ = 0.40 and u₂ = 0.14, which are determined with the best precision among the three planets. The result is summarized in Table 7.

Table 7. Best-fit Parameters for the PPE Light Curve

Parameter	Best-fit Value
Rpd/R⋆	0.06954 ± 0.00043
Rpe/R⋆	0.04123 ± 0.00070
ad/R⋆	26.21 ± 0.082
ae/R⋆	47.44 ± 0.027
bd	0.2951 ± 0.0092
be	0.3693 ± 0.0095
$t_c^{\rm (d)}$ (BJD − 2454833)	378.51372 ± 0.00023
$t_c^{\rm (e)}$ (BJD − 2454833)	378.51785 ± 0.00060
Ωed (deg)	1.13 ± 0.52
$\chi ^2/{\mathrm{{\rm dof}}}$	746/687

Download table as:  ASCIITypeset image

In this fit, R_pd/R_⋆ is different from our revised mean value ($R_{p\rm {d}}/R_{\star } = 0.070288_{-0.00015}^{+0.00014}$) by 1.7σ. We suspect that this difference comes from the systematics introduced in the detrending procedure. When an artifact or some other astrophysical processes (e.g., star spots) accidentally increase the relative flux just before (or after) the transit, the result of detrending is biased toward such features in setting the baseline of the transit light curve. If we use a wider range of data points, the effect of such a small feature is averaged out and does not change the result so significantly. In contrast, if a narrow region around a transit is used, the baseline is somewhat distorted and the resulting detrended light curve becomes either deeper or shallower. Such systematics are averaged in the phase-folded light curves, but may be significant in an individual transit. In the case of the double transit light curve analyzed here, the relative flux begins to increase just before the ingress, and the depth of the transit is shallower in the first half of the transit than in the latter, making the light curve slightly asymmetric. This feature, along with the fact that the revised R_pd gives better χ² in all the other transit light curves than R_pd given by the Kepler team, suggests that the discrepancy in R_pd is caused by such an incidental brightening. To test this scenario, we repeat the analysis above changing the span of detrending from ∼1 day to ∼0.5 days, and find that the resulting mean transit parameters are consistent with those above, but the depth of the double-transit light curve becomes deeper, consistently with our revised parameters.

The occurrence of the PPE in the multi-body context can be assessed in a similar way as in Section 4; we compare R_pd + R_pe with d_min calculated from Equations (18), (22), and (23), but this time the variation of orbital elements must be taken into account. Specifically, we need to evaluate the variations of the parameters relevant to d_min, i.e., scaled semi-major axis a/R_⋆, mean motion n, transit center of the double transit t_c, impact parameter b, and nodal angle Ω (as long as the eccentricities are small). In order to give an estimate for these variations, we integrate the orbits of the three planets using the best-fit (m, e, ϖ) in the rightmost column of Table 6 up to BJD = 2461132.4, the double transit in which the PPE is expected from the two-body analysis in Section 4. The results are the following:

• 1.  
  The oscillation amplitudes of a_d and a_e are less than 5 × 10⁻⁵ AU (∼0.03%) and 4 × 10⁻⁴ AU (∼0.13%), respectively. Since these are much smaller than the observed uncertainties of a_d/R_⋆ and a_e/R_⋆ (∼1%), the multi-body effect can be neglected for these two parameters.
• 2.  
  Corresponding to the oscillations in semi-major axes above, n_d and n_e also show the modulations whose peak-to-peak amplitudes are ∼2π(0.01 day)⁻¹ and ∼2π(0.1 day)⁻¹, respectively. These are much larger than the uncertainties in n_d and n_e that come from those in P_d and P_e, and so the multi-body effect is important for these parameters.
• 3.  
  The differences of the periods calculated from the transit centers in the first ∼1000 days (the range in which we analyzed the TTVs) of integration and from the whole orbit are at most comparable to the observed uncertainties of these parameters in Table 4. Thus the uncertainties of $t_c^{\rm (d)}$ and $t_c^{\rm (e)}$ can be evaluated using those of P and t₀ in the table. Note that the effect of TTV is taken into account in obtaining these errors.
• 4.  
  Monotonic increase in i_d and decrease in i_e lead to at most ∼30% variations of b_d and b_e, larger than the observed errors. The multi-body effect is dominant for these parameters.
• 5.  
  Ω_d/Ω_e also monotonically increases/decreases, but only by <003. This means that the uncertainty in the relative nodal angle Ω_ed is completely dominated by the error in Table 7, and the multi-body effect is negligible.

The above results indicate that the multi-body effect is the most significant for b_d and b_e. In order to relate the values of these parameters to the occurrence of the PPE during the double transit around BJD = 2461132.4, we use Equation (18), (22), and (23) to calculate the maximum value of d_min during this double transit in terms of b_d and b_e, varying (1) a_d, a_e, $t_c^{\rm (d)}$, $t_c^{\rm (e)}$, and Ω_ed within 1σ intervals calculated from the photometric errors, and (2) n_d and n_e by the amplitudes of modulations estimated above. The region of (b_d, b_e) plane in which d_min < R_pd + R_pe, i.e., the PPE occurs in this double transit, is shown in Figure 14 with light-gray shade.¹² When we vary the parameters in the set (1) within their 2σ and 3σ intervals, the edges of the shaded region can be as narrow as black solid and dashed lines. In fact, the gray-shaded region is mainly determined by the difference between b_d and b_e, as seen from this figure, and the edge of this region is found to be most sensitive to the uncertainties in $t_c^{\rm (d)}$ and $t_c^{\rm (e)}$. The former fact originates from the well-aligned orbital planes of KOI-94d and KOI-94e: since their orbital planes are nearly parallel in the plane of the sky, the minimum separation during the double transit in which KOI-94d overtakes KOI-94e is nearly the same as the difference between b_d and b_e. However, if their transit times are too far away from each other, such closest encounter may occur out of the double transit. This explains the latter feature.

Zoom In Zoom Out Reset image size

The PPE occurrence for different choices of (m, e, ϖ) can be judged by evaluating the variation of b_d and b_e in this plot. Since the variations of a, n, t_c, and Ω_ed would be of the same order as long as the resulting TTVs are consistent with the observation, we fix the shaded area in Figure 14 determined by these parameters. We perform a similar MCMC calculation as in Section 2.3 to obtain the distribution of (b_d, b_e) in the double transit at issue; we fit the observed TTVs again using the same χ², but this time extend the orbit integration up to BJD = 2461132.4 and record the final values of b_d and b_e calculated via b = acos i/R_⋆.¹³ The resulting distribution of (b_d, b_e) is plotted with red points in Figure 14. We also repeat the same procedures fixing m_d = 106 M_⊕ and m_d = 73 M_⊕, and the distributions for these cases are plotted with blue points for comparison. In these calculations, we choose the initial values of i_d and i_e based on (b_d, b_e) = (0.2951, 0.3693) in Table 7, a red point with error bars in Figure 14, rather than the mean values obtained from the phase-folded light curves in Table 3. This is because the mean parameters do not take account of the actual occurrence of the PPE. The black, dark-gray, and light-gray points around the double-transit value in Figure 14 respectively show its 1σ, 2σ, and 3σ confidence regions based on the posterior distribution of the double-transit fit. Here, the difference between b_d and b_e is rather sharply constrained by the minimum separation between the planets, namely, the height of the bump caused by the PPE. Even considering this uncertainty in the initial (b_d, b_e), as well as the significant variation of impact parameters (or inclinations) due to the multi-body effect, (b_d, b_e) around BJD = 2461132.4 are well inside the region where the PPE occurs, at least within 1σ of transit and TTV parameters for all the three solutions.

For the best-fit (m, e, ϖ) obtained from the TTV alone, the expected height of the bump is much larger than in the last PPE (Figure 15), and so the detection of this PPE is highly feasible. In contrast, for (m, e, ϖ) based on the RV values of m_d, the bump height is comparable to the last PPE. This difference may be used to settle the difference of m_d values in RV and TTV analyses. For the RV-based solutions, the blue distributions in Figure 14 indicate that the PPE may not even happen when the variation of b_d is too large (corresponding to the large m_c values), and/or ≳ 2σ deviations of $t_c^{\rm (d)}$ and $t_c^{\rm (e)}$ from the linear ephemerides make the shaded region too narrow for the PPE to occur (see solid and dashed lines).

Zoom In Zoom Out Reset image size

We also check the other two double transits before the one discussed above (around BJD = 2457982 and BJD = 24583612), in case that the variations of b_d and b_e lead to the PPE which never happens without the interaction among the planets. In both of them, we find that the PPE does not happen for any possible values of b_d and b_e (from 0 to 1). Therefore, we can safely conclude that the next PPE will still occur during the same double transit as predicted by the two-body calculation, even when we include the mutual gravitational interaction among the planets.