TRANSIT TIMING OBSERVATIONS FROM KEPLER. II. CONFIRMATION OF TWO MULTIPLANET SYSTEMS VIA A NON-PARAMETRIC CORRELATION ANALYSIS

Visual inspection of TT data sets revealed KOIs with MTPCs including two planet candidates which appear to have TTVs that are anti-correlated with each other. In this and two companion papers (Fabrycky et al. 2012; Steffen et al. 2012), we develop complementary methods to establish the statistical significance of apparently correlated TTVs measured in systems of MTPCs. Physically, TTVs (relative to a linear ephemeris) can only be measured at the times of transit. Since TTs can only be measured at discrete times and transits of two planets are rarely coincident with each other, one cannot calculate the standard correlation coefficient, since that would require data sampled either continuously or at common times. In this section, we develop tools to quantify the extent of the correlations between TTV curves and the statistical significance of the TTVs in MTPC systems.

Even if we restrict our attention to systems with only two planet candidates, there is an astounding variety of potential TTV signatures (e.g., Veras et al. 2011). For some systems (e.g., in 1:2 mean motion resonance (MMR) with small libration amplitude), the TTV signature may be relatively simple (e.g., nearly sinusoidal) for the time spans and timing precision of interest. In these cases, one might be well served by fitting a parametric model to the observations (e.g., polynomial, sinusoid). However, for other systems (e.g., slightly offset from resonance, modest eccentricity, more than two planets), the TTV signature can be quite complex. Often the shape and amplitude of the TTVs changes from year to year (or even longer timescales; e.g., Figure 2 of Ford & Holman 2007). Given the diversity and potential complexity of TTV signatures, it is necessary to consider a broad range of functional forms and a large number of model parameters. Normally, this would raise concerns about exploring the high-dimensional parameter space and potentially over-fitting data.

On a simplistic level, one could address this problem by smoothing the TTV observations to obtain a continuous "TTV curve" for each planet candidate and calculate the standard correlation coefficient between the two smoothed TTV curves. Of course, the results will depend on the choice of smoothing algorithm. We overcome the challenges of working with a discrete data set with potentially complex structure by the application of Gaussian Processes (GPs). While one could think of the GPs as a fancy smoothing algorithm, there are several attractive features of GPs that make them well suited for this application. In particular, GPs are infinite-dimensional objects, providing them with enormous flexibility for modeling the data. At the same time, their mathematical properties make it practical to marginalize over the infinite-dimensional parameter space to calculate, not just the predictions of a GP for the TTV curve, but also the full posterior predictive probability distribution for the TTV curve at any finite number of times (see Section 3.3). Thus, we are able to naturally account for the uncertainties in the model TTV curve in a fully Bayesian manner.

Once we obtain a continuous function for the "TTV curve" via our GP model, we calculate a Pearson correlation coefficient between the GP models for each pair of neighboring planets. To establish the statistical significance of the apparently correlated TTVs, we calculate the distribution of correlation coefficients calculated similarly, but applied to synthetic data sets. Each synthetic data set is a random permutation of the actual data. If the correlation coefficient calculated from the actual data is more extreme than the 99.9th percentile of the correlation coefficients calculated from the simulated data, then we conclude that TTVs are sufficiently statistically significant to consider the planet pair confirmed.

Readers who are not interested in the details of the statistical methods used to establish the significance of TTVs in MTPC systems may choose to skip Sections 3.2–3.4. GPs have been extensively studied and applied in the statistics and machine learning communities (Rasmussen & Williams 2006). GPs have also been applied in a variety of astronomical contexts (e.g., Rybicki & Press 1992; Gibson et al. 2012). A GP defines a distribution over functions. "Process" refers to a method for generating a time series of data. "Gaussian" refers to the defining property of a GP that the joint distribution of any finite number of measurements (at discrete times) has a multi-variate Gaussian distribution, even when conditioned on any combination of observations. While this assumption may seem to be restrictive, GPs are extremely versatile. The assumption of Gaussianity is essential for making it practical to perform computations on the infinite-dimensional function space. In particular, for a GP that describes functions of time only (f(t)), a particular GP ($\mathcal {GP}[m(t),k(t,t^{\prime })]$) is specified by a mean value, m(t), and a covariance function, k(t, t'). Once we adopt a form for the covariance function (i.e., for particular values for hyperparameters; see Section 3.3), it becomes practical to perform Bayesian inference and calculate the posterior predictive probability distribution for the values of f(t) at any number of times.

For our application, we can intuitively think of TTVs as measurements of the value of a function (f(t)) that is related to how much a planet is ahead (or behind) schedule in its orbit.²⁴ Of course, Kepler can only make measurements of this function at times of transit (or occultation) as viewed by Kepler. Using a Bayesian framework, we can calculate the posterior probability distribution for $f(\mathbf {t}^*)$ at hypothetical observation times (t*_i), conditioned on the actual measurements of $f(\mathbf {t})$ (i.e., TTVs), and hyperparameters ($\mathbf {\theta }$) that specify the form of the covariance function. We perform this procedure for each planet candidate in a system. Then, we can calculate the correlation coefficient of the GPs for a pair of transiting planet candidates. We focus our attention on neighboring pairs of transiting planet candidates, meaning we do not search for significant correlations between planets if there is an additional transiting planet candidate with an intermediate period. In order to establish the statistical significance of the correlation coefficient, we perform the same procedure on synthetic data sets, generated by permuting the order of the TTVs (along with their measurement uncertainties). We compare the correlation coefficient for the actual TTVs to the distribution of correlation coefficients for the synthetic data sets to determine a false alarm probability (FAP) for the existence of TTVs.

We model the TTVs of each KOI as an independent GP. Following Rasmussen & Williams (2006), the prior probability distribution for the values at times ($\mathbf {t^*}$) of a GP with zero mean is

$$\begin{equation} \vspace*{2.5pt}f_* \equiv f(t^*) \sim \mathcal {N}(\mathbf {0},K(t^*,t^*)),\vspace*{2.5pt} \end{equation} \tag{ 1 }$$

where K(t*, t*) is the correlation matrix. We assume that each observed TT, y_i, is normally distributed about the true TT (f(t_i)) and that independent measurement errors have a variance of σ²_{obs, i}. Then the joint distribution for the actual observations and a set of hypothetical measurements is

$$\begin{equation} \vspace*{2.5pt}\left[\begin{array}{@{}c@{}}\mathbf {y} \\ \mathbf {f^*}\end{array} \right] \sim \mathcal {N}\left(\mathbf {0}, \left[\begin{array}{@{}cc@{}}K\left(\mathbf {t},\mathbf {t}\right)+ {{\mathbf{\sigma}}^2 _{\rm obs}} & K\left(\mathbf {t},\mathbf {t^*}\right) \\ \ms K\left(\mathbf {t^*},\mathbf {t}\right) & K\left(\mathbf {t^*},\mathbf {t^*}\right) \end{array} \right] \right),\vspace*{2.5pt} \end{equation} \tag{ 2 }$$

where $\mathbf {\sigma }^2_{{\rm obs}}$ is a diagonal matrix with entries of σ²_{obs, i}. We can calculate the posterior predictive distribution by conditioning the joint prior distribution on the actual observations, $\mathbf {y}$. The standard results for the mean ($\bar{\mathbf {f}}^*$) and covariance ($\mathrm{cov}(\mathbf {f^*}))$ of the posterior predictive distribution are

$$\begin{eqnarray} \vspace*{2.5pt}\bar{\mathbf {f}}^* & = & K(\mathbf {t^*},\mathbf {t})\left[K(\mathbf {t},\mathbf {t})+\mathbf {\sigma }^2_{{\rm obs}}\right]^{-1} \mathbf {y}\vspace*{2.5pt} \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} \mathrm{cov}(\mathbf {f^*}) & = & K(\mathbf {t^*},\mathbf {t^*}) - K(\mathbf {t^*},\mathbf {t}) \left[K(\mathbf {t},\mathbf {t})+\mathbf {\sigma }^2_{{\rm obs}}\right]^{-1} K(\mathbf {t},\mathbf {t^*}).\quad\quad \end{eqnarray} \tag{ 4 }$$

For our calculations, we adopt a Gaussian covariance matrix

$$\begin{equation} K(t_i,t_j; \sigma _r, \tau) = \sigma _r^2 \exp \left(-\frac{(t_i-t_j)^2}{2\tau ^2}\right), \end{equation} \tag{ 5 }$$

where σ_r and τ are hyperparameters, describing the amplitude and timescale of correlations among data points. This choice of a kernel function ensures that the GP is a smooth function (i.e., continuously differentiable) and allows there to be a single characteristic timescale for TTVs. It can be shown that covariance matrices of this form correspond to a Bayesian linear regression model with an infinite number of Gaussian basis functions (Rasmussen & Williams 2006).

In modeling the TTVs as a GP, we are making use of a posterior predictive distribution for all functions that are consistent with the observations and the covariance matrix specified by a set of hyperparameters $\mathbf {\theta }$. The "trick" of GPs is that we can marginalize over all functions analytically using matrix algebra. Fortunately, the log marginal likelihood conditioned on the hyperparameters can also be readily calculated via matrix algebra

$$\begin{eqnarray} \log p(\mathbf {y}|\mathbf {t}, \mathbf {\theta }) &=& {-}\frac{1}{2}\mathbf {y}^T \left(K_\theta +\mathbf {\sigma }^2_{{\rm obs}}\right)^{-1}\mathbf {y}\nonumber\\ && -\frac{1}{2}\log \left|K_\theta +\mathbf {\sigma }^2_{{\rm obs}}\right| - \frac{n}{2} \log 2\pi, \end{eqnarray} \tag{ 6 }$$

where K_θ is the correlation matrix K evaluated at the observation times, $\mathbf {t}$, using a set of hyperparameters $\mathbf {\theta } = \left\lbrace \log \sigma _r, \log \tau \right\rbrace$ (Rasmussen & Williams 2006). We set σ²_r = median(y_i²), which is equivalent to normalizing the TTV observations by their median absolute deviation. We adopt a value of τ that maximizes the marginal likelihood. We expect this to be a good approximation for data sets of interest, since the likelihood is typically sharply peaked for data sets with a significant structure in the TTV curve. Technically, the marginal likelihood can be bimodal, with one mode having small τ (corresponding to functions that model the observations well) and a second mode with large τ (corresponding to functions that basically ignore the measure TTs, attributing them to measurement noise). In principle, one could explicitly compare the marginal likelihood of each mode or even use the weighted average of the GPs corresponding to the two local maxima. In practice, we found that this was not a problem for the data sets considered. After verifying that the results were not sensitive to the initial guess for τ, we adopted an initial guess for τ equal to the shorter orbital period of the two planets or 90 days (for long-period planets). We use the tnmin minimization package provided by Craig Markwardt²⁵ that is based on the Newton method for nonlinear minimization. We find that this algorithm is robust for the data sets considered.

We also experimented with Matérn class of covariance functions, which can yield GPs that are less smooth than the Gaussian covariance function. The Matérn class is parameterized by both a scale τ and a new hyperparameter, ν. In the limit that ν → ∞, the Matérn class approaches the Gaussian covariance function. For half-integer values of ν, the Matérn covariance functions can be written as a polynomial times an exponential. In particular, we tested ν = 5/2,

$$\begin{eqnarray} K\left(t_i,t_j; \sigma _r, \tau, \nu =\frac{5}{2}\right) &=& \sigma _r^2 \left(1+\frac{\sqrt{5} \Delta t}{\tau } + \frac{5\Delta t^2}{3\tau ^2}\right)\nonumber\\ && \times \exp \left(-\frac{\sqrt{5}\Delta t}{\tau }\right), \end{eqnarray} \tag{ 7 }$$

where Δt = t_i − t_j (Rasmussen & Williams 2006). While the choice of covariance function significantly affects the smoothness of the predictive distribution, the significance of the correlation between two TTV data sets did not appear to be sensitive to the choice of a Gaussian or Matén ν = 5/2 covariance function. Based on our experience analyzing real and simulated data sets, we found that choosing a Gaussian covariance function and the method of estimating τ described above resulted in a highly robust algorithm. This allowed us to automate our analysis, so that we could perform the Monte Carlo simulations necessary to establish the false alarm rates.

Using the GP model for two planet candidates' TTV curves, we calculate the mean ($f_p(\mathbf {t^*})$) and variance ($\sigma ^2_p(\mathbf {t^*})$) of the predictive distribution for each planet, indicated by the index p. For $\mathbf {t}^*$, we adopt the observed TTs for both planets combined into a single vector. We calculate a weighted correlation coefficient between the two GP models based on these two samples, using weights

$$\begin{equation} \vspace*{2.5pt}w_{i} = \left(\sigma ^2_p(t_i) + \sigma ^2_{q}(t_i)\right)^{-1},\vspace*{2.5pt} \end{equation} \tag{ 8 }$$

where σ²_p and σ²_q refer to the variances of the posterior predictive distribution of the GP (i.e., the diagonal elements of $\mathrm{cov}(\mathbf {f^*}$) in Equation (3)). After concatenating the weight vectors for the observation times of the two planets, the weighted average mean and variance of the predictive distributions are

$$\begin{equation} \vspace*{2.5pt}\langle f_p \rangle = \Big[ \sum _i w_i f_p(t_i) \Big] \Big/ \sum _i w_i\vspace*{2.5pt} \end{equation} \tag{ 9 }$$

and

$$\begin{equation} \vspace*{2.5pt}\big\langle \sigma _p^2 \big\rangle = \Big[ \sum _i w_i \sigma ^2_p(t_i) \Big] \Big/ \sum _i w_i.\vspace*{2.5pt} \end{equation} \tag{ 10 }$$

We calculate the covariance between the two GPs evaluated at the actual observation time according to

$$\begin{eqnarray} &&\mathrm{cov}_{p,q} = \Big[ \sum _i w_i \left(f_p(t_i) - \left<f_p\right> \right) \left(f_q(t_i) - \left<f_q\right> \Big) \right] \Big/ \sum _i w_i.\nonumber\\ \end{eqnarray} \tag{ 11 }$$

Thus, the correlation coefficient between the two GPs is estimated by

$$\begin{equation} C = \frac{\mathrm{cov}_{p,q}}{\sqrt{\big(\mathrm{cov}_{p,p}+ \big \langle \sigma _p^2 \big \rangle \big)\big(\mathrm{cov}_{q,q}+\big \langle \sigma _q^2\big \rangle\big)}}. \end{equation} \tag{ 12 }$$

In order to establish the statistical significance of the putative TTVs for a pair of planet candidates, we apply the same methods described in Sections 3.3 and 3.4 to an ensemble of synthetic data sets. Each synthetic data set includes a random permutation of the TTVs for each of the planet candidates. The TTV and measurement uncertainty for a given transit remain paired after permutation. We subtract the best-fit linear ephemeris for each synthetic data set before generating a GP model and calculating the weighted correlation coefficients (C'_i) for the ith synthetic data set. It is important that we subtract the best-fit linear ephemeris before generating the GP model, so that any long-term trend in the TTVs is absorbed into the best-fit orbital period. We consider N_perm = 10⁴ permuted synthetic data sets. We estimate the false alarm probability (FAP_{TTV, C}) based on the fraction of synthetic data sets for which |C'| is greater than |C| for the actual pair of TTV curves. In cases where no synthetic data sets yield a |C'| as large as |C| for the actual observations, we report the FAP_{TTV, C} ⩽ 10⁻³.

In principle, the above process could result in identifying either a positive or a negative correlation coefficient. We expect that an isolated pair of interacting planets will have a negative correlation coefficient, reflecting that energy is exchanged between the orbits. While this is strictly true for two planet systems, one could conceive of a system with additional planets (e.g., a non-transiting planet that is perturbing both of the planets observed to transit) in which a pair of planets would have a positive correlation coefficient. Therefore, we conservatively calculate the FAP based on the distribution of the absolute value of the correlation coefficient, rather than only the fraction of synthetic data sets with C more negative than that measured for the actual observations.

Several arguments about multiplicity suggest that systems with MTPCs have a significantly lower false alarm rate than the overall sample of Kepler planets (Lissauer et al. 2011b). We conservatively adopt a threshold FAP for detecting TTVs of 10⁻³. Using such a threshold, we expect that the probability of claiming significant TTVs for a given pair of planets due to Gaussian measurement noise is less than 10⁻³. Since B11 reported 115, 45, 8, 1, and 1 systems with two, three, four, five, and six candidate transiting planets, there is a total of 238 pairs of neighboring planets which we could test (as opposed to 323 total pairs, including non-neighboring pairs). Thus, the expected number of statistical false alarms from the current sample of multiple planet candidate systems remains less than one, even if we were to consider every system with an FAP <10⁻³ as confirmed. In practice, we only claim to confirm those pairs for which the FAP was robust to outliers (and passed a series of additional tests).

Our estimates of the FAP assume that each TT measurement is independent and uncorrelated with other measurements. While there is correlated noise in the Kepler photometry, this does not directly translate into correlations among measured TTs, since the transits are measured at widely separated times. One possible mechanism for generating apparently anti-correlated TTs is measurement noise due to nearly contemporaneous transits of multiple planets. Therefore, our analysis excludes transits that are nearly coincident with the transit of another known planet candidate. Physically, it is extremely difficult for an alternative astrophysical process to cause the large TTVs observed in these systems. Aside from the gravitational perturbations of the other transiting planet, the most plausible mechanisms would be starspots. However, generating TT variations on timescales longer than a year would require an unusually long-lived spot complex. Furthermore, transit timing noise due to starspots would have an essentially random phase, except in rare cases where the planet orbital period were nearly coincident with a near integer multiple of the stellar rotation period (e.g., Désert et al. 2011). Such a coincidence for multiple planets in one system is even more unlikely. Finally, the observed TTV amplitude is much greater than what can be caused by starspots (Holman et al. 2010). As neither of the stars considered in Section 6 show large rotational modulation, any starspot-induced timing variations are negligible relative to TT measurement precision. Thus, starspots are not a viable explanation of the measured TTVs for either of the stars considered in Section 6. Indeed, if there were an autocorrelation of the TT residuals (TT observations relative to the GP model), the most likely cause would be actual timing variations due to gravitational perturbations that are not accurately described by our GP model. Since our GP model allows for only a single timescale, an orbital configuration that results in multiple TTV timescale would naturally lead to an autocorrelation of the TT residuals. We are optimistic that continued Kepler observations will allow us to detect TTVs on multiple timescales, providing more precise constraints on the masses and orbits of the planets in these systems (see Fabrycky et al. 2012).

To investigate the orbital stability of the system, we construct a nominal model based on the measured orbital periods with circular and coplanar orbits. In the nominal model, we adopt planet masses based on the adopted planet radii, and an empirical mass–radius relationship (M_p/M_⊕ = (R_p/R_⊕)^2.06; Lissauer et al. 2011a, 2011b). We verify that the nominal model is stable for at least 10⁷ years, including all transiting planet candidates in the system.

Second, we place upper limits on the masses of the systems with correlated TTVs by assuming that the real system is not dynamically unstable on a short timescale. For placing mass limits, we start from the nominal model, but include only the pair of planets with significant TTVs. (Including additional planet candidates would typically make the system even less likely to be stable.) We inflate the mass of each planet candidate by a common scale factor. We use the same scale factor for the mass of both planets, since the relative sizes are well determined from the transit light curves (with the possible exception of grazing transits). A misestimated stellar radius would cause both planets' sizes and hence nominal masses by a similar factor. As a misestimated stellar radius would not significantly change the planet–planet size ratio, the nominal mass ratio for our n-body simulations is insensitive to uncertainties in the stellar radius. Another possibility is that both planets' sizes could be significantly underestimated if the light from the target star were diluted by light from nearby stars. Again, this could significantly affect the planets' sizes, but not the planet–planet size ratio or the planet–planet mass ratio for our nominal model. A final possibility is that both planets with correlated TTVs are transiting a star other than the target star, which is diluted by the target star. In this case, the planet radii would be larger than estimated by B11, but the ratio of planet radii would still be accurately estimated (assuming neither is grazing). Importantly, in each of these cases dynamical stability would still provide upper limits for the masses that show the bodies to be planets, even if the planetary radii were significantly larger than estimated.

Basic information about transit parameters, stellar, and planetary properties was presented in B11. A few key parameters are reproduced or updated in Tables 1 and 3. In Table 4, we report the correlation coefficient (C) for several pairs of neighboring transiting planets, along with the false alarm probability (FAP_{TTV, C}) for the TTVs based on our correlation analysis and Monte Carlo simulations (see Section 3.5). The table also includes the ratio of transit durations normalized by orbital period to the one-third power (ξ), the 5th or 95th percentile of the distribution of ξ expected for a pair of planets around the target host star, the ratios of the rms and mean absolution deviation (MAD) of the TTs from the best-fit linear ephemeris for the two planets. In the final column, κ gives the ratio of the measured MAD of TTVs for the two planets to the ratio of the predicted TTVs for the same two planets, based on our nominal n-body models. A subset of these systems (Kepler-23 = KOI 168; Kepler-24 = KOI 1102) will be discussed individually in Section 6. For these systems, we show the measured TTs and the GP model for a subset of these systems to be discussed in more detail in Figures 1 and 2. To help illustrate how we calculate C_0.001 and FAP_{TTV, C}, we show the histogram of C' values from synthetic data in Figure 1 (bottom right panel) for Kepler-23b and c. The methods developed in this paper find significant and apparently correlated TTVs in several additional pairs of Kepler planet candidates that are to be discussed in Cochran et al. (2011), J.-M. Désert et al. (2011), Fabrycky et al. (2012), Lissauer et al. (2012), D. Ragozzine et al. (2012, in preparation), and Steffen et al. (2012).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. Statistics for Pairs of Neighboring Planet Candidates

KOIin	KOIout	Pout/Pin	$\frac{R_{p,{\rm in}}}{R_{p,{\rm out}}}$	C	FAPTTV, C	ξ	ξ5/95	$\frac{\mathrm{rms_{in}}}{\mathrm{rms_{out}}}$	$\frac{\mathrm{MAD_{in}}}{\mathrm{MAD_{out}}}$	κ
168.03	168.01	1.511	0.54	−0.863	<10−3	0.89	0.57	4.27	2.35	1.02
244.02	244.01	2.039	0.58	−0.774	<10−3	1.57	2.62	1.74	2.69	0.54
250.01	250.02	1.404	1.01	−0.825	<10−3	1.53	2.55	1.22	0.96	1.56
738.01	738.02	1.286	1.15	−0.692	0.0010	0.90	0.45	0.57	0.66	1.17
806.03	806.02	2.068	0.26	−0.841	<10−3	0.86	0.74	13.42	27.81	0.54
841.01	841.02	2.043	0.81	−0.803	<10−3	0.88	0.42	1.02	2.34	0.59
870.01	870.02	1.520	1.05	−0.223	0.2972	0.89	0.46	0.49	0.56	1.09
935.01	935.02	2.044	1.14	−0.237	0.3140	0.99	0.40	0.80	1.56	1.01
952.01	952.02	1.483	0.98	−0.728	<10−3	1.12	2.38	0.67	0.43	0.79
1102.02	1102.01	1.514	1.14	−0.905	<10−3	1.11	2.91	1.28	0.73	1.25

Download table as:  ASCIITypeset image

The detection of significant and anti-correlated TTVs provides strong evidence that the bodies responsible for the transits are in the same physical system. In principle, one might wonder whether the "transits" could actually be eclipses of stellar mass bodies. In order to account for the TTVs, the bodies still need to be in the same physical system. Given the similar orbital periods, any stellar companions would interact very strongly, raising serious doubts about the long-term dynamical stability of the system.

We report the maximum planet mass for which our n-body integrations did not result in at least one body being ejected from the system or colliding with the other body or the central star (Table 1; Figure 3). For each of the planets presented in Section 6, the maximum mass is less than 13 M_Jup, excluding a triple star system as a possible false positive. The observed TTVs suggest even lower maximum masses, but a complete TTV analysis will require a longer time series of observations. There are considerable uncertainties in the stellar masses, but the available observations preclude us from having underestimated the stellar mass by more than a factor of two or more, which would be necessary for the maximum stable mass to approach 13 M_Jup. Even ∼13 M_Jup is roughly half of the recently proposed criteria for exoplanets (∼25 M_Jup; Schneider et al. 2011). Therefore, regardless of the choice of definition, the masses are constrained to be in the planetary regime.

Zoom In Zoom Out Reset image size

The combination of TTVs and dynamical stability provides strong evidence that the transits are due to planets orbiting a common star. Thus, we promote these planet candidates to confirmed planets. KOI 168.03 and 168.01 become Kepler-23b and c, respectively. KOI 1102.02 and 1102.01 become Kepler-24b and c, respectively. In Section 4.3, we show that the probability of the planetary system orbiting a host star other than the original target star is very small for both of the cases considered in detail.

KOIs 168.02, 1102.03, and 1102.04 remain strong planet candidates. Given the low rate of false positives among the Kepler multi-planet candidates, it is quite unlikely that these KOIs are caused by a blend with a background object. The most likely form of a "false positive" would be another planetary system transiting a second physically associated and similar mass star in the Kepler aperture. Thus, we anticipate that further follow-up observations (such as high-resolution images) and analysis (such as BLENDER) could allow the remaining planet candidates to be validated as planets. Alternatively, continued Kepler TT observations may allow for dynamical confirmation of some of these candidates.

While the correlated TTVs and dynamical stability provide evidence for a planetary system, it is not yet obvious that the system must orbit the original target star. Here we consider the three alternate potential scenarios that could result in similar appearances: (1) both planets orbit a background star, (2) both planets orbit a significantly cooler star that is physically associated with the target star, and (3) both planets orbit one of two physically associated stars of similar mass.

The probability of the first case (planetary system around an unassociated star) can be quantified by considering the range of spectral type and magnitude differences (measured relative to the target star) which could result in a transit of the hypothetical background star mimicking the observed transit. Since dynamical stability precludes stellar masses and an object's radius is insensitive to its mass in the Jupiter-to-brown-dwarf-mass regime, there is a maximum size for the transiting body and the background star must be sufficiently bright that it could result in the observed transit depth (after accounting for dilution by the target star). The maximum difference in magnitude is ΔK_{p, max} = 5.3 mag for Kepler-23 and ΔK_{p, max} = 2.7 mag for Kepler-24. The potential locations for a background star are constrained by the observational limits on the centroid motion, i.e., the difference in the location of the flux centroid during transit and out of transit. We adopt maximum angular separation equal to the 3σ confusion radius, 03 for Kepler-23c and 09 for Kepler-24c. Using a Besançon galactic model with magnitude and position for each target star, we estimate the frequency of background blends that could match the observed transit depth to be ∼5 × 10⁻⁴ for Kepler-23 and ∼2 × 10⁻³ for Kepler-24. We have not included any constraints on the observed transit duration or shape. As Kepler continues to observe these systems, we expect that TTVs will eventually provide significant constraints on the orbital eccentricities. Incorporating such a constraint would be expected to rule out small host stars that would require an apocentric transit to match the observed duration. Even if we were to assume that large planets are as common as small planets, then these planetary systems are more than ∼10⁴ and ∼5 × 10² times more likely to orbit the target star than an unassociated background star. This excludes the vast majority of background blends, including those involving the reddest host stars.

Next, we consider the possibility that the target star might have an undetected stellar companion that could host the planetary system. We can exclude blends that would result in a color Δ(r − K) ⩾ 0.1 mag larger than expected based on a simple set of isochrones (Marigo et al. 2008). Combining photometry from the KIC and the Two Micron All Sky Survey, this typically rules out companions in the ∼0.6–0.8 M_☉ range. We will consider nearly equal-mass binaries later in this section. If the planets were to orbit a physically associated star other than the target, then the host star would be less massive than the target star and the corresponding upper mass limits for the planets would be further reduced, since the dynamical stability constraint is most closely related to the planet–star mass ratios. For Kepler-24 the constraints based on the transit depth and dynamical stability overlap, so there are no viable blend scenarios where the planetary system orbits a physically bound and significantly lower-mass secondary star. For Kepler-23, we compute the frequency of plausible blend scenarios using the observed frequency of stellar binaries (Raghavan et al. 2010) and the observed distribution of orbital periods and primary–secondary mass ratios (Duquennoy & Mayor 1991). We conservatively assume that large planets are as common as small planets and do not impose constraints based on the transit duration or shape. We find a blend frequency of 0.11, indicating that the Kepler-23 system is at least 9 times more likely to be hosted by the primary target than by a physically bound and significantly lower-mass secondary star.

Finally, we consider the potential for the target to be a binary star with two similar mass stars. In many cases spectroscopic follow-up observations would have detected a second set of spectra lines. However, we cannot totally exclude a long-period binary with two stars that happen to have the same radial velocity at the present epoch. Since the two stars would have similar properties, the planet properties are largely unaffected, aside from the ∼50% dilution causing the planet radii to increase by ∼40%. As we are not overly concerned about which of the two similar stars hosts the planet, this scenario essentially amounts to unseen dilution, a regular concern among faint transiting planets.

One of the common reasons for a KOI to have been identified as a likely false positive or to have received a vetting flag of 3 in B11 is a measurement of significant difference in the transit depth of odd and even numbered "transits." This can occur if the apparent transit is due to an EB, where the odd and even "transits" differ in which star is eclipsing and which is being eclipsed. (Typically, the EB must also be diluted and/or grazing in order for the depth to be consistent with a planet.) The Kepler pipeline provides an odd–even depth test statistic that can be used to identify KOIs potentially due to an EB (Steffen et al. 2010). Inspecting the odd–even test statistic is also advised to check that the inferred orbital period is not half the true orbital period. Such a misidentification can arise for low signal-to-noise candidates, such as KOI 730 (see Lissauer et al. 2011b; D. Fabrycky et al. 2012, in preparation). We have verified that the odd–even test statistic is less than 3 for each of the planets with significantly correlated TTVs that is discussed in Section 6, as well as the other planet candidates in these systems. In the course of our analysis, we noted that the originally reported period for KOI 168.02 was an artifact at one-third the period of the updated period for KOI 168.02.

Another of the common reasons for a KOI to receive a vetting flag of 3 in B11 was a measurement of significant difference in the location of the centroid of the target star during transit and out-of-transit. Centroid motion can be due to a background EB that is blended with the target star. On the other hand, small but statistically significant centroid motion does not necessarily imply that the KOI is not a planet. For example, local scene crowding can induce an apparent shift in the photometric centroid, as well as introduce biases in point-spread function (PSF)-fitted estimates of the location of the transiting object via difference images. Removing these spurious sources of centroid motion requires extensive analysis as described in S. T. Bryson et al. (2012, in preparation). Here, we limit ourselves to planet candidates where the above biases are negligible, so the apparent centroid motion is within the 3σ statistical error due to pixel-level noise for most of the quarters of data analyzed. None of the planet candidates with correlated TTVs that are discussed in Section 6 have statistically significant centroid motion.

For targets with MTPCs, the ratio of transit durations can be used as a diagnostic to reject blend scenarios (Holman et al. 2010; Batalha et al. 2011; Lissauer et al. 2012; R. Morehead et al. 2012, in preparation). Therefore, we perform a light curve fitting based on Q0–6 data primarily to measure transit durations. We construct folded light curves based on the measured TTs (see Table 2). Otherwise, we follow the fitting procedure of Moorhead et al. (2011). Note that the durations reported in Table 1 are based on when the center of the planet is coincident with the limb of the star and that this differs from the durations reported in B11 that were based on the time interval between first and fourth points of contact. The former definition of duration is less sensitive to the uncertainties in the planet size, impact parameter, and limb darkening (Colón & Ford 2009; Moorhead et al. 2011). Note that the planet candidates discussed in this paper have faint host stars, so there is often a large uncertainty in the impact parameters. In most cases where the impact parameter is not well constrained, we adopt a central transit, similar to both B11 and Moorhead et al. (2011).

We report the normalized transit duration ratio (ξ ≡ (T_{dur, in}/T_{dur, out})(P_out/P_in)^1/3) for each pair of neighboring planets in Table 4. We also calculate ξ₅ and ξ₉₅, the 5th and 95th percentile of ξ values obtained from Monte Carlo simulations of an ensemble of systems with two planets with the measured orbital periods and a distribution of impact parameters and eccentricities. Simulations are discarded if both planets do not transit or if transits of one planet would not have been detected (due to grazing transits that would result in a reduced S/N). In Table 4, we report ξ_5/95 which is simply ξ₅ for pairs with ξ < 1 and ξ₉₅ for pairs with ξ > 1. Typically, the distribution of ξ is not far from symmetric, so ξ_0.95 ∼ 1/ξ_0.05. For a few pairs where the two planets have substantially different radii, the simulated distribution of ξ is asymmetric due to the minimum signal-to-noise criterion. In all cases, the measured ξ is consistent with a pair of planets transiting a common host star.

The Kepler FOP has obtained high-resolution spectra of KOI host stars from the 10 m Keck I Observatory, the 3 m Shane Telescope at Lick Observatory, 2.7 m Harlan J. Smith Telescope at McDonald Observatory, or the 1.5 m Tillinghast Reflector at Fred Lawrence Whipple Observatory (FLWO). The choice of observatory and exposure time was tailored to produce the desired signal to noise. For some faint stars, an initial low-S/N reconnaissance spectra was used for initial vetting, before obtaining a second higher S/N spectrum that was used for analysis. For Kepler-23, stellar parameters are based on spectra from McDonald Observatory and the Nordic Optical Telescope. Spectra were analyzed by computing the correlation function between the observed spectrum and a library of theoretical spectra using the tools described in Buchave et al. (2012). This method provides measurements of effective temperature (T_eff), metallicity ([M/H]), and surface gravity (log (g)), as well as a means of recognizing binary companions that would produce a second set of spectral lines. In the case of Kepler-24, we adopt the stellar atmospheric parameters from the KIC and reported in Borucki et al. (2011), since a high-quality spectrum is not yet available.

We report the adopted stellar atmospheric parameters in Table 3. We also update the stellar mass and radius based on Bayesian comparison to Yonsei–Yale isochrones (Yi et al. 2001).

The Kepler mission follow-up observing program includes speckle observations obtained at the WIYN 3.5 m telescope located on Kitt Peak. Speckle observations of Kepler-23 were used to provide high spatial resolution views of the target star to look for previously unrecognized close companions that might contaminate the Kepler light curve. The speckle observations make use of the Differential Speckle Survey Instrument (DSSI), a recently upgraded speckle camera described in Horch et al. (2011) and Howell et al. (2011). The DSSI provides simultaneous observations in two filters by employing a dichroic beam splitter and two identical EMCCDs as the imagers. The details of how we obtain, reduce, and analyze the speckle results and specifics about how they are used eliminate false positives and aid in transit detection are described in Torres et al. (2011), Horch et al. (2011), and Howell et al. (2011). The latter paper also presents the speckle imaging results for the 2010 observing season.

Classical imaging systems provide complementary observations with a wider field of view. In particular, the Lick Observatory 1 m Nickel Telescope took an I-band image of Kepler-23 with a pixel scale of 0368 pixel⁻¹ and seeing of ∼15. For Kepler-24, the 2 m Faulkes Telescope North (FTN) provides SDSS-r' band images with a pixel scale of 0304 pixel⁻¹ in the default 2 × 2 pixel binning mode, and a typical seeing of ∼ 12. For each target a stacked image is generated by combining several images taken during the night, or on separate nights, when the target is at different positions on the sky. This is done in order to achieve increased sensitivity for faint stars without saturating the bright stars, and to average out the diffraction pattern of the spider vanes supporting the secondary mirror (since FTN has an alt-azimuth mount the diffraction pattern is different at different sky positions, relative to the positions of the stars).

Kepler-23 (KOI 168, KID 11512246, Kp = 13.4) hosts three small planet candidates (168.03, 168.01, and 168.02) with orbital periods of 7.10, 10.7, and 15.3 days and radii of ∼3.2, 1.9, and 2.2 R_⊕, respectively (B11; N. Batalha et al. 2012, in preparation). In the course of this analysis, we recognized that the originally reported period for KOI 168.02 was an artifact with one-third the true period. Initial vetting of KOI 168.01 (planet c, 10.7 days) resulted in a vetting flag of 2 and an estimated EB probability of ∼3.2 × 10⁻⁵ (B11), indicating that this candidate is very unlikely to be due to a background EB. KOI 168.02 and 168.03 (b) were not vetted in time for B11. Based on the analysis described in Section 5, we find no evidence pointing toward a false positive for any of the three planet candidates. An analysis of possible centroid motion resulted in 3σ radii of confusion of R_{c, 23c} ∼ 03 and R_{c, 23b} ∼ 5''. While the radius of confusion for Kepler-23b is relatively large, the mean offset of the centroid during transit is only ∼0.1σ, consistent with measurement uncertainties. Due to the faint stars and relatively low S/N of the transits, the constraint for Kepler-23b is relatively weak when compared to Kepler planet candidates transiting brighter stars. Nevertheless, the centroid measurements are still powerful results. Since the typical optimal aperture for photometry of Kepler-23 is 9 pixels in area (i.e., equivalent to a circular aperture with a radius of ∼67), excluding locations for a potential background EB to within 5'' eliminates 80% of the optimal aperture phase space for blends with background objects. The analysis in Section 4.3 shows that the probability of a background EB causing one of Kepler-23b or c is less than ∼10⁻⁴ and the probability of the planets being around a lower-mass and physically bound star is less than ∼11%.

An I-band image was taken with the Lick Observatory 1 m Nickel Telescope on 2010 July 10. There were no companions visible from ∼2'' to 5'' away from the target star down to 19th magnitude (see Figure 4). On 2011 June 11 and 12, the FOP obtained speckle images of Kepler-23 in a band centered on 880 nm with width 55 nm using DSSI at WIYN. Seven and five integrations, each consisting of one thousand 40 ms exposures were coadded. These observations reveal no secondary sources with a 4σ limiting delta magnitude of 2.94 at 02, increasing to 3.5 at 1'' and 3.57 at 18. As no new nearby stars were identified in either set of observations, we adopt ∼2.5% contamination of the Kepler aperture used by the Kepler pipeline based on pre-launch photometry, the optimal aperture used for photometry and the PSF.

Zoom In Zoom Out Reset image size

The Kepler FOP obtained two medium-resolution reconnaissance spectra of Kepler-23 from McDonald Observatory (HJD = 2455153.643746) and FIES (2455054.576815). The spectra result in stellar atmospheric parameters of T_eff = 5760 ± 124 K, log (g) = 4.0 ± 0.14, and [M/H] = −0.09 ± 0.14 (Buchave et al. 2012), consistent with the stellar properties in the KIC. We compare to the Yonsei–Yale isochrones and derive values for the stellar mass (1.11^+0.09_−0.12 M_☉) and radius (1.52^+0.24_−0.30 R_☉) that are slightly smaller than those from the KIC. We estimate a luminosity of ∼2.3 L_☉ and an age of ∼4–8 Gyr.

For Kepler-23b, the S/N of each transit is only ∼1.7, resulting in significant timing uncertainties (σ_TT = 47 minutes). All three planet candidates are clearly identified after folding the light curve at the best-fit orbital period (Figure 5). Ford et al. (2011) noted that Kepler-23c was a TTV candidate based on an offset in the transit epoch between the best-fit linear ephemerides based on Q0–2 and Q0–5. TTVs for Kepler-23b were not recognized based on Q0–2 observations alone.

Zoom In Zoom Out Reset image size

The correlation coefficient between the GP models for TTs of Kepler-23b and c is C = −0.86, well beyond the distribution of correlation coefficients calculated based on synthetic data sets with scrambled TTs (see Figure 1). The FAP for such an extreme correlation coefficient is <10⁻³. The correlation coefficient between the GP models for TTs of Kepler-23c and KOI 168.02 is −0.23 with an FAP ∼12%. This is well above our threshold for claiming a planet detection. We do not yet detect statistically significant TTVs for 168.02, as expected due to the smaller predicted TTV signal for KOI 168.02 and the sizable timing uncertainties. In combination with the analysis of Kepler light curves and FOP observations described above, the strongly anti-correlated TTVs of Kepler-23b and c provide a dynamical confirmation that the two objects orbit the same star.

Kepler-23b and c have short orbital periods and lie close to a 3:2 resonance: P_23b/P_23c = 1.511. Therefore, we expect a relatively short TTV timescale:

$$\begin{eqnarray} &&P_{\rm TTV} = 1 / (3/ 10.7421\;{\rm days} - 2 / 7.1073\;{\rm days}) = 470\;{\rm days}. \nonumber\\ \end{eqnarray} \tag{ 13 }$$

The observed timescale (∼455 days) and phase of TTVs for Kepler-23b and c are consistent with expectations based on n-body integrations using nominal planet mass estimates and circular, coplanar orbits (see Figure 6). The amplitude of the observed TTVs is ∼5 times larger than the nominal model. Such differences are not unexpected, as even modest eccentricities can significantly increase the TTV amplitude (e.g., Veras et al. 2011). The prediction of the nominal model for the ratio of the TTV amplitudes of the two confirmed planets is significantly more robust than the predictions of individual amplitudes. Indeed, the observed MAD TTV from a linear ephemeris is within 5% of that predicted by the nominal model. Our GP models for the TTs of Kepler-23b and c are non-parametric, so there is nothing in our model that would cause the GPs to match the ratio of TTV amplitudes, the TTV period, or the phase of TTV variations. The fact that our general non-parametric model naturally reproduces these features provides an even more compelling case for the confirmation of these planets via TTVs.

Zoom In Zoom Out Reset image size

While the currently available data are insufficient for robust determinations of the planet masses, we fit two n-body models to the TTV observations. If we assume initially circular and coplanar orbits, then the best-fit masses are ∼22 ± 6 and 12 ± 2 M_⊕. If we assume eccentric coplanar orbits, then the best-fit masses are ∼15 and 5 M_⊕, corresponding to orbital solutions with low eccentricities. Significantly smaller or larger masses are possible, but these require large eccentricities. Either model represents a significant improvement in the quality of the fit relative to a non-interacting model. Despite the sizable uncertainties in the estimates for the planet masses, the requirement of dynamical stability provides firm upper limits (2.7 and 0.8 M_Jup) on the planet masses (Figure 3).

Kepler-24 (KOI 1102, KID 3231341, Kp = 14.9) hosts four planet candidates (04, 02, 01, and 03) with orbital periods of 4.2, 8.1, 12.3, and 19.0 days and radii of 1.7, 2.4, 2.8, and 1.7 R_⊕, respectively (Figure 7; B11; N. Batalha et al. 2012, in preparation). KOI 1102.02 (b) and 1102.01 (c) were identified, but not vetted in time for B11. Based on the analysis described in Section 5, we find no evidence pointing toward a false positive for any of KOI 1102.01–1102.04. The centroids during transit for Kepler-24b and c differ from those out of transit by ∼1.4σ and 2.8σ, consistent with measurement uncertainties, particularly when considering the actual distribution of the apparent centroid motion during transit is highly non-Gaussian for other well-studied cases (Batalha et al. 2011). The ∼3σ radii of confusion are 12 and 09 for Kepler-24b and c, respectively. Based on a preliminary analysis of data through Q8, the centroid offsets for 1102.04 and 1102.03 are less than 3σ. Again, the limited centroid motion enables us to exclude the locations for potential background EBs to a small fraction of the optimal aperture phase space for blends with background objects. The analysis in Section 4.3 shows that the probability of a background EB causing one of KOI 1102.01 (c) or 1102.02 (b) is less than ∼2 × 10⁻³. The only acceptable blend scenario involves a pair of similar mass and physically bound stars.

Zoom In Zoom Out Reset image size

A UKIRT J-band image (Figure 4, left) reveals a faint companion that is not in the KIC ∼25 to the southeast, with brightness intermediate to the outlying faint stars KIC 3231329 (Kp = 19.0) and KIC 3231353 (Kp = 20.0). An FTN image also reveals the nearby star, ∼5 mag fainter than the target in the r' band (Figure 4, right). We estimate that it is ∼4–5 mag fainter than Kepler-24 in the Kepler band. As most of its light falls in the Kepler image and star is not in the KIC, it results in ∼1%–3% dilution in addition to the ∼8% dilution estimated by the pipeline for a total of 8.6%.

The Kepler FOP has not yet obtained a spectrum of Kepler-24. We estimate the stellar properties based on multi-band photometry from the KIC: T_eff = 5800 ± 200, log (g) = 4.34 ± 0.3, and [M/H] = −0.24 ± 0.40. We estimate uncertainties based on a comparison of KIC and spectroscopic parameters for other Kepler targets (Brown et al. 2011). We derive values for the stellar mass (1.03^+0.11_−0.14 M_☉), radius (1.07^+0.16_−0.23 R_☉), luminosity ∼1.16^+0.36_−0.60 L_☉, and age ⩽7.7 Gyr by comparing to the Yonsei–Yale isochrones. We caution that there are significant uncertainties associated with the stellar models and derived properties.

Kepler-24 is sufficiently faint that TTs for Kepler-24c and Kepler-24b have sizable uncertainties (σ_TT = 30 and 25 minutes). Ford et al. (2011) did not find significant evidence for TTVs in either of the planet candidates, but noted that Kepler-24b was likely to have significant TTVs (∼26 minutes), based on numerical integrations of two planet systems with nominal masses and circular orbits.

The correlation coefficient between the GP models for TTs of Kepler-24b and c is C = −0.905, well beyond the distribution of correlation coefficients calculated based on synthetic data sets with scrambled TTs (see Figure 2). The FAP for such an extreme correlation coefficient is <10⁻³. In combination with the analysis of Kepler light curves described above, this provides a dynamical confirmation that the two objects orbit the same star.

Kepler-24b and c have short orbital periods and lie near a 3:2 resonance: P_24c/P_24b = 1.514. Therefore, we expect a relatively short TTV timescale:

$$\begin{eqnarray} &&P_{{\rm ttv}} = 1 / (2 / 8.1451\;{\rm days} - 3/ 12.3336\;{\rm days}) = 421\;{\rm days}. \nonumber \\ \end{eqnarray} \tag{ 14 }$$

The observed timescale (∼400 days) and phase of TTVs for Kepler-24b and c are consistent with the expectations based on n-body integrations using nominal planet mass estimates and circular, coplanar orbits (see Figure 8). Our GP models for the TTs of Kepler-24b and c are non-parametric, so there is nothing in our model that would cause the GPs to match the period and phase of the TTV variations. The fact that our general non-parametric model naturally reproduces these features provides an even more compelling case for the confirmation of these planets via TTVs.

Zoom In Zoom Out Reset image size

The best-fit amplitudes of the observed TTVs are ∼7 and 9 times larger than the nominal model for Kepler-24b and c, respectively. Such differences are not unexpected, as even modest eccentricities can significantly increase the TTV amplitude (e.g., Veras et al. 2011). Like Kepler-23, the prediction of the nominal model for the ratio of the TTV amplitudes of the two confirmed planets is within 25% of the observed ratio (1.25). Still, there could be significant differences between the true masses and the nominal masses, e.g., if there are significant eccentricities and/or perturbations from other planets in the system. Assuming that KOI 1102.03 is indeed a planet in the same system, then Kepler-24c would be near a 3:2 MMR with both Kepler-24b (inside) and KOI 1102.03 (outside), leading to complex dynamical interactions.

While the currently available data are insufficient for robust determinations of the planet masses, we fit two n-body models to the TTV observations. If we assume initially circular and coplanar orbits, then the best-fit masses are ∼17 ± 4 and 5 ± 3 M_⊕. If we assume eccentric coplanar orbits, then the best-fit masses are ∼102 ± 21 and 56 ± 16 M_⊕, corresponding to orbital solutions with eccentricities of ∼0.2–0.4. Either model represents a significant improvement in the quality of the fit relative to a non-interacting model. While the best-fit eccentric model represents a significant improvement in the goodness-of-fit relative to the best-fit circular model, such large eccentricities may complicate long-term stability, particularly if all four of the planet candidates were confirmed. Despite the sizable uncertainties in the estimates for the planet masses, the requirement of dynamical stability provides a firm upper limit on the planet masses (1.6 and 1.6 M_Jup; Figure 3).

We do not detect a statistically significant correlation coefficient between TTVs of either KOI 1102.03 or 1102.04 and any of the other planet candidates. The scatter of the TTVs for 1102.04 is consistent with the measurement uncertainties. One can understand the lack of a TTV detection for KOI 1102.04 analytically, since the amplitude of a TTV signal scales with the orbital period, KOI 1102.04 has an orbital period roughly half that of Kepler-24b, and the TTV signal for Kepler-24b is near the limit of detectability. Furthermore, the transit of KOI 1102.04 is shallower than Kepler-24b, so the TTVs are less precise. If KOI 1102.04 is a planet in the same system, then the current innermost period ratio would be P_24b/P_1102.04 = 1.92, ∼5% less than 2, whereas it is more common for systems near a 1:2 MMR to have a period ratio ∼5% greater than 2 (Lissauer et al. 2011b). This could be related to the presence of four planet candidates with orbital period ratios near a 2:4:6:9 resonant chain. The current period ratios for the outer two pairs (P_24c/P_24b = 1.514 and P_1102.03/P_24c = 1.54) are slightly greater than 1.5, so none of the neighboring pairs are necessarily "in resonance." Such deviations from an exact period commensurability appear to be common among Kepler planet candidates (Lissauer et al. 2011b).

There may be early hints of excess scatter in the TTVs of 1102.03, but this is not statistically significant based on the current data set. If KOI 1102.03 is a planet in the same system, then Kepler-24c would feel significant perturbations from both Kepler-24b (on the inside) and 1102.03 (on the outside). This may help explain why the observed TTVs for Kepler-24b are significantly greater than those predicted by the nominal two-planet model that includes only planets b and c. We hope that this paper will motivate more detailed analysis of the TTVs and orbital dynamics of this fascinating planetary system.