VERY LOW DENSITY PLANETS AROUND KEPLER-51 REVEALED WITH TRANSIT TIMING VARIATIONS AND AN ANOMALY SIMILAR TO A PLANET–PLANET ECLIPSE EVENT

Since its launch in 2009, NASA's Kepler telescope has discovered more than 3000 transiting planetary candidates, more than one third of which are the members of multi-transiting systems (Batalha et al. 2013). With longer observation spans, the detection limits have been extended to smaller and longer-period planets, and even a compact solar-system analog has been recently discovered (Cabrera et al. 2014). Characterization of such multi-planetary systems (not necessarily transiting) is especially valuable because reproducing the architecture of multiple planets is theoretically more demanding than for a single planet, and hence the resulting constraints on the formation theories become much tighter.

Fortunately, multi-transiting systems exhibit various features advantageous for their confirmation and characterization (e.g., Ragozzine & Holman 2010). These include transit timing variations (TTVs; Miralda-Escudé 2002; Holman & Murray 2005; Agol et al. 2005), deviations of transit times from the strict periodicity due to the mutual gravitational interaction among the planets. By modeling the TTVs of several transiting planets consistently, we can precisely estimate their masses, which are usually inaccessible with the photometric observations alone. The TTVs become especially prominent for planets in near mean-motion resonances (Agol et al. 2005), which turned out to amount to several percent of the Kepler sample (Mazeh et al. 2013; Xie et al. 2013). Indeed, planetary masses obtained from TTVs are all the more valuable for systems for which radial velocities are difficult to obtain, as are the cases for most of the Kepler target stars.

A wealth of the Kepler data has also revealed a new phenomenon called a planet–planet eclipse (PPE), which was first and only observed in the Kepler-89 (KOI-94) system (Hirano et al. 2012a; Masuda et al. 2013). In this event, two planets transit the host star simultaneously and even overlap with each other on the stellar disk. In addition to its rareness and astronomical interest, this phenomenon tightly constrains the relative angular momentum of the two planets involved, which could give some clues to unveil the history of their orbital evolutions through the presence or absence of dynamical planet–planet interaction.

This paper focuses on the Kepler-51 (KOI-620) system, one of the multi-transiting systems found by Kepler. This system hosts three transiting planet candidates, two of which were confirmed by Steffen et al. (2013). They made sure that these two planets, Kepler-51b (KOI-620.01) and Kepler-51c (KOI-620.03), are revolving around the same star by confirming that their TTVs are anti-correlated, and showed that they are indeed planetary by giving mass upper limits based on the long-term stability of the system. In this paper, we perform a numerical analysis of their TTVs to more fully characterize the system and to confirm that the mass of KOI-620.02 is also in the planetary range. We also discuss a very intriguing light curve of Kepler-51 recently made public on the Space Telescope Science Institute MAST archive, which shows a feature similar to a PPE event.

We adopt the stellar properties in Table 1 taken from the NASA Exoplanet Archive.¹ As an initial guess for the limb-darkening coefficients for the quadratic law, we adopt (u₁, u₂) = (0.36, 0.28), the values for (T_eff, log g, Z, ξ) = (6000 K, 4.5 dex, 0.0, 0.0 km s⁻¹) in the grid of Claret & Bloemen (2011). Linear ephemerides and transit parameters are retrieved from the MAST archive (Table 2) as a starting point for the iterative determination of these parameters in Section 3.

3.1. Data Processing

3.2. Transit Times and Transit Parameters

From the transit light curves obtained in Section 3.1, we determine the transit times and transit parameters of the three planets by iterative fit using a Markov chain Monte Carlo (MCMC) algorithm. We repeat the following two steps: (1) we fit each transit for the time of the transit center t_c using the light curve model by Ohta et al. (2009). Here we assume e = 0 and fix the values of P, R_p/R_⋆, b (of each planet), u₁, u₂, and ρ_⋆. From the series of transit times, period P and time of a transit center t₀ are extracted by linear regression. (2) Using the transit times obtained in step (1), we phase fold the transits of each planet and fit the three phase curves simultaneously for R_p/R_⋆, b (of each planet), u₁, u₂, and ρ_⋆. Here the values of P are fixed at those in step (1) and e = 0 is assumed for all the planets.

Starting from the values in Table 2 (and in Section 2 for u₁ and u₂), all the parameters converge sufficiently well after five iterations. The resulting transit parameters, ephemerides, and transit times are summarized in Tables 3–6. The quoted best-fit parameters (a/R_⋆ to ρ_⋆ and t_c) denote the median values of their posteriors, and uncertainties exclude 15.87% of values at upper and lower extremes. The corresponding best-fit transit models with the phase-folded transits are shown in Figure 1. As reported in previous analyses (Steffen et al. 2013; Mazeh et al. 2013), we find significant TTVs for all three planets as shown in Figure 2. Note that the TTV amplitude of KOI-620.02 in our analysis is about twice as large as that first reported by Mazeh et al. (2013), who analyzed the first twelve quarters of Kepler data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

Parameter	Kepler-51b	Kepler-51c	KOI-620.02
t0 (BJD − 2454833)	159.10653 ± 0.00033	295.3131 ± 0.0018	212.03246 ± 0.00039
P (days)	45.155314 ± 0.000019	85.31644 ± 0.00022	130.178058 ± 0.000071
a/R⋆	$61.5_{-1.2}^{+1.5}$	$94.1_{-1.9}^{+2.2}$	$124.7_{-2.5}^{+3.0}$
Rp/R⋆	$0.07414_{-0.00061}^{+0.00059}$	$0.094_{-0.017}^{+0.028}$	$0.10141_{-0.00085}^{+0.00084}$
b	$0.251_{-0.138}^{+0.073}$	$1.017_{-0.023}^{+0.034}$	$0.250_{-0.141}^{+0.075}$
u1	$0.375_{-0.036}^{+0.040}$
u2	$0.311_{-0.087}^{+0.083}$
ρ⋆ (g cm−3)	$2.16_{-0.13}^{+0.15}$
χ2/dof	12681/12417

Download table as:  ASCIITypeset image

Table 4. Transit Times of Kepler-51b (KOI-620.01)

Transit	tc	1σlower	1σupper	χ2/dof	O − C
Number	(BJD − 2454833)				(days)
0	159.10975	0.00072	0.00072	2.14	0.00323
1	204.26437	0.00078	0.00076	1.86	0.00253
2	249.41453	0.00120	0.00152	3.24	−0.00262
3	294.57446	0.00251	0.00159	2.12	0.00199
4	339.72399	0.00083	0.00088	2.32	−0.00379
5	384.87799	0.00078	0.00079	4.04	−0.00510
6	430.03405	0.00076	0.00076	1.78	−0.00436
8	520.34240	0.00151	0.00168	0.80	−0.00663
9	565.49926	0.00106	0.00148	3.29	−0.00509
10	610.65682	0.00087	0.00095	1.00	−0.00285
11	655.81302	0.00080	0.00084	1.38	−0.00196
12	700.97595	0.00204	0.00156	2.19	0.00566
13	746.12646	0.00082	0.00086	1.10	0.00085
14	791.28654	0.00102	0.00129	1.79	0.00562
15	836.43982	0.00074	0.00074	2.24	0.00358
16	881.59882	0.00072	0.00071	0.91	0.00727
17	926.75475	0.00083	0.00078	1.42	0.00789
18	971.90566	0.00181	0.00262	1.95	0.00348
19	1017.05878	0.00083	0.00088	1.62	0.00129
20	1062.21217	0.00075	0.00075	2.50	−0.00064
21	1107.36887	0.00095	0.00097	0.94	0.00075
22	1152.52090	0.00088	0.00088	0.96	−0.00253
23	1197.67687	0.00097	0.00097	0.87	−0.00188
24	1242.83059	0.00087	0.00087	0.99	−0.00347
25	1287.98482	0.00086	0.00088	0.92	−0.00456
26	1333.14289	0.00091	0.00090	0.95	−0.00179
27	1378.29779	0.00088	0.00088	0.86	−0.00220
28	1423.45442	0.00091	0.00090	1.00	−0.00089
29	1468.61324	0.00089	0.00089	0.97	0.00261

Download table as:  ASCIITypeset image

Table 5. Transit Times of Kepler-51c (KOI-620.03)

Transit	tc	1σlower	1σupper	χ2/dof	O − C
Number	(BJD − 2454833)				(days)
0	295.31257	0.00378	0.00384	0.98	−0.00057
1	380.64295	0.00358	0.00354	0.97	0.01337
2	465.95289	0.00287	0.00283	1.41	0.00687
3	551.26161	0.00319	0.00304	0.99	−0.00086
4	636.56677	0.00324	0.00325	2.04	−0.01214
7	892.51469	0.00384	0.00393	1.90	−0.01355
8	977.84149	0.00360	0.00364	1.16	−0.00319
10	1148.45861	0.00327	0.00327	1.00	−0.01896
11	1233.80785	0.00322	0.00324	0.89	0.01385
12	1319.11072	0.00331	0.00342	0.95	0.00027
14	1489.75414	0.00337	0.00340	0.88	0.01080

Download table as:  ASCIITypeset image

Table 6. Transit Times of KOI-620.02

Transit	tc	1σlower	1σupper	χ2/dof	O − C
Number	(BJD − 2454833)				(days)
0	212.02417	0.00066	0.00066	2.67	−0.00829
1	342.20715	0.00063	0.00062	2.28	−0.00337
2	472.39116	0.00064	0.00064	2.08	0.00258
3	602.57341	0.00063	0.00063	2.17	0.00678
5	862.93196	0.00076	0.00070	3.88	0.00921
6	993.10424	0.00064	0.00065	2.35	0.00343
7	1123.28307	0.00065	0.00066	1.12	0.00420
8	1253.44963	0.00062	0.00063	0.89	−0.00730
9	1383.62994	0.00064	0.00064	0.99	−0.00505

Download table as:  ASCIITypeset image

Several comments should be added to our revised values of transit parameters in comparison to those in Table 2. With the longer baselines, we refine the orbital periods of the three planets with better precision than the previous values. We also find the larger values for R_p/R_⋆, albeit with relatively large uncertainties. This is because the slight variations of transit depths we identify in the archived light curves, probably due to the starspot activities (see also the discussion in Section 5). In fact, our analysis completely neglects such spot effects, and so the values of R_p/R_⋆ we determined may be overestimated. The constraint on R_p/R_⋆ is especially poor for Kepler-51c, whose grazing transit causes the strong correlation between its planetary radius and impact parameter. The values of impact parameters we determine are marginally consistent with those in Table 2, corresponding to the slightly different value of stellar density; these parameters would be determined more precisely with spectroscopic constraints on the stellar mass and radius.

In addition, we find that the difference between the impact parameters of Kepler-51b and KOI-620.02 is tightly constrained due to the strong correlation, in spite of their relatively large uncertainties: the MCMC posteriors for the two parameters yield $b\,({\rm KOI\mathchar`-620.02}) - b\,({\rm Kepler\mathchar`-51b}) = -0.001 \pm 0.01$. Since this difference is closely related to the minimum separation during a simultaneous transit of the two planets, it has an important role in assessing the occurrence of the PPE, as will be discussed in Section 5.

In this section, we perform a numerical analysis of the TTVs for the planetary parameters (especially their masses) to confirm KOI-620.02 as a planet and to more fully characterize the system. For simplicity, we assume the coplanar orbits for the three planets and fix the stellar mass at M_⋆ = 1.04 M_☉. We make no attempt to model transit parameters other than the transit time.

We define transit centers as the minima of the star–planet distance in the plane of the sky (D), and calculate the simulated transit times in the following way. We integrate the planetary orbits using the fourth-order Hermite scheme with the shared time step (Kokubo & Makino 2004). From the position and velocity of each planet, we calculate the time derivative of D and search for its root applying the Newton–Raphson method (Fabrycky 2010). All the simulations presented in this section are performed between BJD = 2454980 and BJD = 2456345, beginning at the same epoch T₀(BJD) = 2455720 (close to the center of the observation time).

We fit the three planets' TTVs simultaneously for the mass M_p, transit time closest to the epoch T_c, orbital period P, eccentricity e, and argument of periastron ω (measured from the sky plane) of each planet. Here P, e, and ω are the osculating orbital elements defined at the epoch T₀. Since we assume the coplanar orbits, we fix the initial values of the orbital inclinations i = π/2 and longitudes of the ascending nodes Ω = 0. Chi-squares are computed from the simulated transit times as

$$\begin{eqnarray} &&\chi ^2 = \sum _{j:\,{\rm planets}}\ \sum _{\begin{array}{@{}c@{}}{\scriptstyle i:}\,{\scriptstyle \rm observed}\\ {\scriptstyle \rm transits}\end{array}} \left[ \frac{t_{c,j}(i) - t^{\rm sim}_{c,j}(i)}{\sigma _j (i)} \right]^2, \end{eqnarray} \tag{ 1 }$$

where $t^{\rm sim}_{c,j}(i)$, t_{c, j}(i), and σ_j(i) are the simulated central time, observed central time, and uncertainty of the ith transit of planet j, respectively. For simplicity, we adopt averages of 1σ upper and lower limits of transit times as σ_j(i).²

We first use the downhill simplex method by Nelder and Mead (Press et al. 1992) to find the minimum in the above χ², and then perform an MCMC analysis (Ford 2005, 2006) around this minimum. The median values of the MCMC posteriors, their 1σ uncertainties, and minimum value of χ² are shown in Table 7, and the corresponding best-fit simulated TTVs are plotted in Figure 2. We perform the same procedures also floating M_⋆ with the Gaussian prior based on M_⋆ = 1.04 ± 0.12 M_☉, and obtain the consistent results with no better constraint on M_⋆. Analysis taking account of the apparent non-coplanarity of Kepler-51c (i ∼ 894) does not alter the result, either.

Table 7. Best-fit Parameters Obtained from TTVs

Parameter	Kepler-51b	Kepler-51c	KOI-620.02
Mp (M⊕)	$2.1_{-0.8}^{+1.5}$	4.0 ± 0.4	7.6 ± 1.1
Tc (BJD − 2454833)	881.5977 ± 0.0004	892.509 ± 0.003	862.9323 ± 0.0004
P (days)	45.1540 ± 0.0002	$85.312_{-0.002}^{+0.003}$	$130.194_{-0.002}^{+0.005}$
ecos ω	−0.016 ± 0.006	$0.010_{-0.008}^{+0.013}$	$0.005_{-0.006}^{+0.011}$
esin ω	−0.04 ± 0.01	$-0.009_{-0.013}^{+0.009}$	$-0.006_{-0.010}^{+0.008}$
χ2	51 (29 transits)	21 (11 transits)	11 (9 transits)
χ2/dof	83/34

Note. P in this table is one of the osculating orbital elements at the simulation epoch T₀, and so its value is different from the average period obtained from the transits in Table 3.

Download table as:  ASCIITypeset image

In Figure 2, the sinusoidal TTVs of Kepler-51b and KOI-620.02 are well explained by their proximities to 2:1 and 3:2 resonances, respectively, with Kepler-51c: the periods of these two planets' TTVs inferred from the observed data (∼770 days and ∼2500 days) are in well agreement with the "super-period" P^j = 1/|j/P_outer − (j − 1)/P_inner| for a j: j − 1 resonance defined by Lithwick et al. (2012). In addition, the best-fit masses of all three planets, including the one of KOI-620.02, fall into the planetary regime. These facts strongly indicate that KOI-620.02 is a planet belonging to the same system as the other two.

Remarkably, the best-fit masses of all three planets are less than that of Neptune, in spite of their relatively large values of R_p/R_⋆. At least for the mass of Kepler-51c, our value is also supported by another study: Hadden & Lithwick (2013) analyzed the Kepler data through Quarter 12 using analytic TTV formulae derived by Lithwick et al. (2012). In their formulae, the TTV amplitudes of a pair of coplanar planets in near a j: j − 1 resonance are analytically given as functions of their masses, eccentricities, and orbital phases. Since the orbital phases are already constrained from transit observations, the formulae allow us to constrain the planets' masses and eccentricities in a degenerate way. Assuming e = 0, they obtain two estimates for the mass of Kepler-51c, one using the inner pair (9.7 M_⊕) and one using the outer pair (3.1 M_⊕). Although the non-zero eccentricities easily alter these estimates by a factor of a few (Lithwick et al. 2012), these values are consistent with the mass of Kepler-51c that we obtain here.

Using $\rho _{\star } = 2.16_{-0.13}^{+0.15}\,{\rm g\,cm^{-3}}$ obtained from the transit light curves (Table 3) and M_⋆ = 1.04 ± 0.12 M_☉, we obtain R_⋆ = 0.88 ± 0.04 R_☉, which is consistent with the value in Table 1. Note that e = 0 is assumed in determining the value of ρ_⋆, though the correction is of order esin ω and smaller than the uncertainty of ρ_⋆ in Table 3, according to the TTV results in Table 7. This value of R_⋆, along with the values of R_p/R_⋆ in Table 3 and M_p in Table 7, gives the radius and density of each planet listed in Table 8. In addition to the fact that they have relatively large uncertainties, the planetary radii could be slightly overestimated due to the starspots, as discussed in Section 3. Nevertheless, these results show that the bulk densities of all the planets in the Kepler-51 system are arguably among the lowest of the known planets, falling below that of the recently discovered sub-Saturn radius planet KOI-152d (Jontof-Hutter et al. 2013).

Table 8. Planet Properties Obtained from Transit Light Curves and TTVs

Parameter	Kepler-51b	Kepler-51c	KOI-620.02
Mp (M⊕)	$2.1_{-0.8}^{+1.5}$	4.0 ± 0.4	7.6 ± 1.1
Rp (R⊕)	7.1 ± 0.3	$9.0_{-1.7}^{+2.8}$	9.7 ± 0.5
ρp (g cm−3)	$0.03_{-0.01}^{+0.02}$	$0.03_{-0.03}^{+0.02}$	0.046 ± 0.009
a (AU)	0.2514 ± 0.0097	0.384 ± 0.015	0.509 ± 0.020
e	0.04 ± 0.01	$0.014_{-0.009}^{+0.013}$	$0.008_{-0.008}^{+0.011}$
Teq (K)	543 ± 11	439 ± 9	381 ± 8

Notes. We adopt M_⋆ = 1.04 ± 0.12 M_☉ in calculating the values of semi-major axes. Equilibrium temperatures are calculated from T_eff in Table 1, a/R_⋆ in Table 3, and e using $T_{\rm eq} = \sqrt{R_{\star }/2a} (1-e^2)^{-1/4}\,T_{\rm eff}$.

Download table as:  ASCIITypeset image

While their densities are much lower than those of the planets in the solar system, it is possible to form such planets with sufficiently rich gas envelopes. Lopez & Fortney (2013) calculated the radii of low-mass (1–20 M_⊕) planets for various values of envelope fraction (0.01%–60%), incident flux (0.1–1000 F_⊕), and age (10 Myr–10 Gyr). According to their Equation (3), the observed radii of the Kepler-51 planets can be explained if they have about 10% (Kepler-51b), 30% (Kepler-51c), and 40% (KOI-620.02) of their masses in their H/He envelopes, for the age of 0.3 Gyr. Note that the required H/He fractions increase with the system age, because the older planets tend to have smaller radii for fixed masses due to the gradual cooling of the gas envelopes. This may imply that the host star Kepler-51 is actually young, as suggested by the KIC classification (Table 1). This idea is also supported by Walkowicz & Basri (2013), who determined the age of Kepler-51 as 0.53 Gyr, though they note that the ages are highly uncertain for very young stars.

On the other hand, it seems more difficult to explain how they acquired the above fractions of H/He envelopes. Although the simulations by Rogers et al. (2011) (see, e.g., their Table 2) show that such planets could be formed by core-nucleated accretion beyond the snow line followed by the inward migration to T_eq ∼ 500 K, their results are based on the somewhat arbitrary assumption that the planet migrates after the sufficient growth of its core and envelope. In situ accretion (Ikoma & Hori 2012) is also unlikely to account for the predicted atmospheric fractions, unless their natal disk was relatively cool and dissipated slowly.

In the above analysis, we adopt the value of M_⋆ in Table 1, but not that of R_⋆, which is only poorly constrained. More precise determination of the stellar (and hence planetary) mass and radius requires the constraints on the stellar parameters with spectroscopic observations.

In the double transit of Kepler-51b and KOI-620.02 that occurred around BJD = 2456346.8,³ we identify an increase of the relative flux near the transit center. Considering the fact that KOI-620 shows ∼12 mmag (∼ 1%) variation associated with its rotation (McQuillan et al. 2013), this "bump" can be naturally explained by a spot-crossing event (e.g., Silva 2003; Silva-Valio 2008; Rabus et al. 2009). Indeed, we find several transits of KOI-620.02 showing brief brightenings of similar amplitudes (∼0.2 %) as seen in this double-transit light curve. In addition, this double transit occurred during a gradual increase of the stellar flux, which indicates that a large star spot (or a group of starspots) was moving on the visible side of the star at that time. However, as we mentioned in the last part of Section 3, the small difference of their impact parameters obtained from the transit light curves requires that the PPE should have occurred in this double transit, provided that (1) their orbital planes are nearly aligned and that (2) their cosine inclinations have the same signs. Since the inner planet Kepler-51b overtakes the outer one KOI-620.02 in this double transit, the minimum sky-plane separation becomes small enough under these conditions.

Motivated by this fact, we fit the observed double-transit light curve with the PPE model by Masuda et al. (2013) for ΔΩ, the longitude of the ascending node of KOI-620.02 relative to that of Kepler-51b. Note that here we choose the plane of the sky as a reference plane, and so ΔΩ corresponds to the mutual inclination of the two planetary orbits in this plane (see also the lower panel of Figure 3). Also note that, in the following, we consider a general case where ΔΩ can take any value from −180° to 180°, though we only discussed the aligned (ΔΩ ∼ 0) case above as a motivation for the PPE scenario. This is because, in general, the occurrence of a PPE event is not limited to the aligned case (Masuda et al. 2013). The other parameters R_p/R_⋆, b, and ρ_⋆ are also floated except for u₁ and u₂, which are fixed at the values in Table 3. While we restrict the impact parameter b of Kepler-51b to be positive, we allow b of KOI-620.02 to be either positive or negative, taking into account that the two planets can have different signs of cosine inclinations. Using an MCMC algorithm, we find that this model gives a reasonably good fit with χ²/dof = 0.94, and obtain $\Delta \Omega = -25.3_{-6.8}^{+6.2}\,{\rm deg}$ for the sky-plane mutual inclination of the two planets (Figure 3 and Table 9).

Zoom In Zoom Out Reset image size

Table 9. Resulting Parameters of the PPE Fit to the Double-transit Light Curve

Parameter	Value
Kepler-51b
a/R⋆	63.65 ± 0.33
Rp/R⋆	0.0741 ± 0.0017
b	$0.016_{-0.012}^{+0.024}$
tc(BJD − 2454833)	1513.76694 ± 0.00083
KOI-620.02
a/R⋆	128.93 ± 0.66
Rp/R⋆	0.1019 ± 0.0011
b	$0.039_{-0.040}^{+0.038}$
tc(BJD − 2454833)	1513.78988 ± 0.00070
ρ⋆	2.393 ± 0.037
ΔΩ	$-25.3_{-6.8}^{+6.2}$
χ2/dof	807/859

Note. In this fit, the impact parameter of KOI-620.02 is allowed to be either positive or negative, while that of Kepler-51b is fixed to be positive.

Download table as:  ASCIITypeset image

This value, if true, indicates that the orbital planes of Kepler-51b and KOI-620.02 are significantly misaligned, which means that either of their orbital axes are tilted with respect to the stellar spin axis. This result may be in contrast to the spin–orbit alignments observed in five multi-transiting systems so far (Sanchis-Ojeda et al. 2012; Hirano et al. 2012a, 2012b; Chaplin et al. 2013; Albrecht et al. 2013), but agrees with the recent discovery that the spin–orbit misalignment is not confined to hot-Jupiter systems (Huber et al. 2013).

However, if their orbits are really significantly misaligned, it follows that this multi-transiting system is a rare object. For example, the conditional probability p(02|b) that the outer KOI-620.02 transits when the inner Kepler-51b is known to transit is p(02|b) ≃ a₀₂sin ϕ/R_⋆ ∼ 1/60 for the mutual inclination of ϕ ∼ 30 deg and a₀₂/R_⋆ ∼ 130 (see, e.g., Ragozzine & Holman 2010). This value is smaller by the factor of a_bsin ϕ/R_⋆ ∼ 30 than in the aligned case, where p(02|b) ≃ a_b/a₀₂ ∼ 1/2.

In fact, a further examination of the PPE model reveals that it is consistent with the results of the phase-curve analysis in Section 3 only when both planets have |b| ∼ 0 as in Table 9, making this scenario all the less likely. This situation is illustrated in Figure 4, where we examine the validity of the PPE model for all the possible values of the two planets' impact parameters. In this plot, the color scale shows the minimum value of the χ² difference between the PPE and non-PPE models, found by incrementing ΔΩ by 1° from −180° to 180°. The grid scale of b is 0.01, and the other parameters are fixed at the best-fit values in Table 9. In the dark-blue region, the PPE model significantly improves the fit, while the dotted green lines correspond to the constraint from the phase-folded transits we mentioned in Section 3. Comparing these two regions, we find that only the values of b close to the PPE best fit (red point with error bars) are consistent with the PPE interpretation of the bump. Moreover, the impact parameters required by the PPE model are only marginally consistent with those obtained in Section 3 (Table 3), though they are close to the ones obtained by the Kepler team (Table 2).

Zoom In Zoom Out Reset image size

It should also be noted that, if the orbit of the outermost planet has such a large mutual inclination with respect to the inner two, their orbits will precess rapidly. In this case, the resulting transit duration variations (TDVs) would be fairly large for the middle grazing planet, Kepler-51c. Nevertheless, no significant TDVs are apparently seen in the transits of this planet (Figure 1). With no constraints on the nodal angle of Kepler-51c, it may still be possible that the orbit of this planet is also tilted with respect to that of Kepler-51b so that the effect of precession is canceled out, but it would require fine tuning of the parameters. Furthermore, we must happen to observe this system when all the planets have small eccentricities as indicated from TTVs, while their eccentricities would vary as the system evolves under such a large mutual inclination.

These arguments about the transit probability, impact parameters, and expected orbital precession, imply that a PPE event that is consistent with all of the observations would be exceedingly rare; in other words, the PPE interpretation for the observed anomaly is essentially refuted. The observed bump is, therefore, probably due to a star spot or just a correlated noise. If it is the spot crossing, the detailed starspot modeling may provide valuable information on the stellar obliquity (Dittmann et al. 2009; Silva-Valio et al. 2010; Sanchis-Ojeda et al. 2011, 2012, 2013; Nutzman et al. 2011), which is closely related to the orbital evolution history of the planets (e.g., Queloz et al. 2000; Winn et al. 2005). This is beyond the scope of this paper.

There are several possible approaches to strengthen the above interpretation of the anomaly. First, follow-up observation of the next double transit where a PPE event might occur is unreasonable, because we have to wait at least until 2092 even if the orbital planes of the two planets are completely aligned. Secondly, the more accurate determination of the impact parameters would be helpful. If it is confirmed that the impact parameters of the two planets differ from zero with better precision, the discussion based on Figure 4 is enough to exclude the PPE scenario. In contrast, if both planets have |b| ∼ 0 as suggested by the Kepler team (Table 2), we need to explain why the two planets did not overlap; if ΔΩ ∼ 0 in this case, we should have observed a large bump that is totally inconsistent with the observed one. In order to better determine the impact parameters, the better constraint on R_⋆ (or ρ_⋆ itself) would again be quite beneficial, because the prior knowledge on ρ_⋆ pins down their values, which are strongly correlated to that of ρ_⋆. Lastly, as mentioned above, a thorough analysis of the dynamical model taking account of mutual orbital inclinations and (if necessary) starspots would provide a more decisive conclusion on the origin of this anomaly. Yet another possibility is the analysis of the long-term dynamical stability, which may rule out the misaligned configuration.

We have discussed the two topics in this paper, characterization of the multi-transiting planetary system around Kepler-51 with TTV analysis (Sections 3 and 4) and interpretation of the light-curve feature similar to a PPE caused by the two planets in this system (Section 5). Here we briefly summarize each of the topics and give some additional comments.

1. Characterization of the Kepler-51 system. We analyzed the transit light curves and TTVs of the three planets in the Kepler-51 system, which lie close to a 1:2:3 resonance chain. Combining the planetary masses obtained from TTVs, and planet-to-star radius ratios and stellar density inferred from the transit light curves, we determined the properties of the three planets as follows: $M_{\rm b} = 2.1_{-0.8}^{+1.5} \,M_{\oplus }, R_{\rm b} = 7.1 \pm 0.3 \,R_{\oplus }, \rho _{\rm b} = 0.03_{-0.01}^{+0.02} \,{\rm g\,cm^{-3}}$ for Kepler-51b (KOI-620.01), $M_{\rm c} = 4.0 \pm 0.4 \,M_{\oplus }, R_{\rm c} = 9.0_{-1.7}^{+2.8} \,R_{\oplus }, \rho _{\rm c} = 0.03_{-0.03}^{+0.02}\,{\rm g\,cm^{-3}}$ for Kepler-51c (KOI-620.03), and M₀₂ = 7.6 ± 1.1 M_⊕, R₀₂ = 9.7 ± 0.5 R_⊕, ρ₀₂ = 0.046 ± 0.009 g cm⁻³ for KOI-620.02. From these results, as well as the sinusoidal modulation consistent with their proximities to the resonances, we confirmed KOI-620.02 as a planet in this system (Kepler-51d), which has an equilibrium temperature close to the inner edge of the habitable zone.

The even more remarkable implication of our analysis is that the densities all three planets in this system are among the lowest yet determined, though a more detailed study taking account of the presence of starspots might increase these values. In fact, such low-density planets are frequently seen in other compact multi-transiting planetary systems; these include the systems around Kepler-9 (Holman et al. 2010), Kepler-11 (Lissauer et al. 2011; Migaszewski et al. 2012; Lissauer et al. 2013), Kepler-18 (Cochran et al. 2011), Kepler-30 (Sanchis-Ojeda et al. 2012), Kepler-56 (Huber et al. 2013), Kepler-87 (Ofir et al. 2014), Kepler-89 (Weiss et al. 2013; Masuda et al. 2013), and KOI-152 (Jontof-Hutter et al. 2013), all of which have planets with sub-Saturn densities. Considering the fact that such systems frequently contain giant planets and planet pairs in near mean-motion resonances (and, of course, that they are the multi-transiting systems), they were probably formed via the convergent disk migration and subsequent resonance capture. As pointed out by Jontof-Hutter et al. (2013), the selection effects for planets suitable for the TTV analysis could be significant. Nevertheless, if the low-densities observed so far are intrinsic features of compact multi-transiting systems, they could be an important constraint on the formation mechanisms via disk migrations.

2. Anomaly similar to a planet–planet eclipse event. We also analyzed the double-transit light curve of Kepler-51b and KOI-620.02 around BJD = 2456346.8. The archived Kepler light curve shows a slight increase in the relative flux of Kepler-51, which could be explained by the PPE, the overlap of the two planets during their double-transit phase.

If the cosine inclinations of the two planets have the same signs, the impact parameters of the two planets strongly suggest that the PPE should have occurred in this double transit. Indeed, the PPE model well reproduces the observed anomaly for the sky-plane mutual inclination between the two planets of ∼25 deg, which implies that their orbital planes are misaligned. This result, if true, indicates that either of their orbital planes are tilted with respect to the stellar spin axis, and makes the Kepler-51 system another important piece of evidence that the spin–orbit misalignment is not confined to hot-Jupiter systems (Huber et al. 2013).

However, this interpretation of the anomaly seems unlikely for the following reasons. First, such a large mutual inclination significantly reduces the probability that both of Kepler-51b and KOI-620.02 transit. Second, the PPE model is consistent with the result of phase-curve analysis only for limited values of the two planets' impact parameters. Finally, the misaligned configuration would result in the rapid orbital precession, whose effect should have been readily detectable in the transit light curves of the middle grazing planet, Kepler-51c. Alternative interpretations of the anomaly include the correlated noise and the starspot crossing. If the latter is the case, it may provide us the information on the stellar obliquity (Dittmann et al. 2009; Silva-Valio et al. 2010; Sanchis-Ojeda et al. 2011, 2012, 2013; Nutzman et al. 2011), which is definitely valuable in unveiling the orbital evolution history of the planets in this system.

In any case, it is rewarding to explore the origin of this anomaly, because it serves as an example of the false positive of a PPE event. Compared with the case of the Kepler-89 (KOI-94) system, where a small light-curve modulation led to the clear detection of a PPE (Hirano et al. 2012a), the situation is less ideal for the Kepler-51 system analyzed in this paper. A detailed investigation of the possible phenomena (e.g., starspots) that could produce PPE-like features would help the future detection of this valuable event in such marginal conditions.

This work is based on the photometry of Kepler-51 provided by the Kepler mission, and the author gratefully acknowledges the Kepler team. The author thanks Masahiro Ikoma for a helpful discussion on the interpretation of the observed planetary densities. The author would also like to thank an anonymous referee, Dan Fabrycky, Jack Lissauer, Sean Mills, and Yasushi Suto for useful comments to improve the manuscript. This work was supported by the Program for Leading Graduate Schools, MEXT, Japan.