KEPLER-18b, c, AND d: A SYSTEM OF THREE PLANETS CONFIRMED BY TRANSIT TIMING VARIATIONS, LIGHT CURVE VALIDATION, WARM-SPITZER PHOTOMETRY, AND RADIAL VELOCITY MEASUREMENTS^*

One of the objects identified with possible transiting planets is the Kp 13.549 mag (where Kp is the magnitude in the Kepler passband) star KIC 8644288 (2MASS J19521906+4444467, K00137). After a possible transiting planet has been detected, the Kepler data are subjected to a set of statistical tests to search for possible astrophysical false-positive origin of the observed signal. These data validation tests for the first five Kepler planet discoveries are described by Batalha et al. (2010b). Additional tests, including measurement of the image centroid motion during a transit (Wu et al. 2010) all gave a high probability that the signals seen were real. The application of these techniques to Kepler-10b is described in detail by Batalha et al. (2011).

Two separate transiting objects were immediately obvious in the LC data. K00137.01 has a transit ephemeris of T₀[BJD] = (2455167.0883 ± 0.0023) + N*(7.64159 ± 0.00003) days and a transit depth of 2287 ± 9 ppm. (All transit times and ephemerides in this paper are based on UTC.) K00137.02 has an ephemeris of T₀[BJD] = (2455169.1776 ± 0.0013) + N*(14.85888 ± 0.00004) days and a transit depth of 3265 ± 12 ppm. It was noted that the orbits of these two objects were very near a period ratio of 2:1. After filtering these transits of K00137.01 and K00137.02 from the light curve, we searched again for transiting objects and found a third planet candidate in the system, K00137.03, which has a shorter orbital period than the other two transiting planets. K00137.03 has the ephemeris T₀[BJD] = (2454966.5068 ± 0.0021) + N*(3.504725 ± 0.000028) days, and a transit depth of only 254 ± 8 ppm. The light curve for K00137 is shown in Figure 1. The upper panel shows the raw uncorrected "photometric analysis" (PA) light curves that come out of the data processing pipeline, and the lower panel shows the corrected PDC light curves. The two transit events that look significantly deeper than the others, near BJD 2454976.0 and BJD 2455243.5 are simultaneous transits of K00137.01 and K00137.02.

Zoom In Zoom Out Reset image size

The folded light curves for each of the three candidate planets are shown in Figure 2. This figure shows the light curves folded on the ephemeris given above for each KOI. The significant width of the ingress and egress for Kepler-18c and Kepler-18d is a result of the transit time variations discussed in Section 6.

Zoom In Zoom Out Reset image size

The first FOP step is to obtain high spectral resolution, low signal-to-noise ratio (S/N) spectroscopy in order to verify the Kepler Input Catalog (KIC) stellar classification, and to search for evidence of stellar multiplicity in the spectrum. Spectra from the Hamilton echelle spectrograph on the Lick Observatory 3 m Shane Telescope were obtained on the nights of 2009 August 8 and 9, and on 2009 September 1 UT. These spectra showed no convincing evidence for RV variability at the 0.5 km s⁻¹ level, and no hints of any contaminating spectra. These spectra were cross-correlated against a library of synthetic stellar spectra as described by Batalha et al. (2011), in order to derive basic stellar parameters to compare with the KIC values. These spectra yielded T_eff = 5250 ± 125 K, log g = 4.0 ± 0.25, and Vsin i = 2 ± 2 km s⁻¹. The height of the cross-correlation peaks ranged from 0.82 to 0.89, indicating a very good match with the library spectra.

High-resolution imaging of the surroundings of a KOI is an important step to identify possible sources of false-positive signals. We need to ensure that the detected transit signal is indeed originating on the selected target star, and not on a background star that was unresolved in the original KIC imaging of the Kepler field.

A seeing limited image was obtained at Lick Observatory's 1 m Nickel telescope using the Direct Imaging Camera. This image was a single one-minute exposure taken in the I band, with seeing of approximately 15. The only nearby star is a faint object (approximately 5 mag fainter than K00137) about 55 to the north. A J-band image of the 1' × 1' field of view around K00137 was taken as part of a complete J-band survey of the Kepler field of view using the wide-field camera on the United Kingdom Infrared Telescope. These images have a typical spatial resolution of 08–09 and an image depth of J = 19.6 (Vega system). This image also shows only a faint source about 55 to the north.

We perform speckle observations at the WIYN 3.5 m telescope, using the Differential Speckle Survey Instrument (DSSI; Horch et al. 2011; Howell et al. 2011). DSSI provides simultaneous observations in two filters by employing a dichroic beam splitter and two identical EMCCDs as the imagers. We generally observe simultaneously in V and R bandpasses, where V has a central wavelength of 5620 Å, and R has a central wavelength of 6920 Å, and each filter has an FWHM of 400 Å. The speckle observations of K00137 were obtained on 2010 June 22 UT and consisted of three sets of 1000, 80 ms individual speckle images. Along with a nearly identical V-band reconstructed image, the speckle results reveal no companion star near K00137 within the annulus from 005 to 18 to a limit (5σ) of 4.2 mag fainter in both V and R than the Kp 13.549 target star.

Near-infrared adaptive optics (AO) imaging of K00137 was obtained on the night of 2009 September 8 UT with the Palomar Hale 5 m telescope and the Palomar High Angular Resolution Observer (PHARO) near-infrared camera (Hayward et al. 2001) behind the Palomar AO system (Troy et al. 2000). PHARO, a 1024 × 20124 HgCdTe infrared array, was utilized in 25.1 mas pixel⁻¹ mode, yielding a field of view of 25''. Observations were performed in the J filter. The data were collected in a standard five-point quincunx dither pattern of 5'' steps interlaced with an off-source (60'' east) sky dither pattern. Data were taken with integration times per frame of 30 s with a total on-source integration time of 7.5 minutes. The AO system guided on the primary target itself and produced a central core width of FWHM = 0075. Figure 3 shows this Palomar image of K00137. There are two additional sources within 15'' of the primary target. The first source, located 56 to the north of K00137, is 4 mag fainter at J; the second, located 150 to the southeast, is 5.8 mag fainter. To produce the observed transit depths in the blended Kepler photometry, the δ_J = 4 mag star would need to have eclipses that are 0.09 and 0.5 mag deep, and the δ_J = 5.8 mag star would need to have eclipses that are 0.15 and 1.3 mag deep. Such deep eclipses from stars separated by more than 4'' (>1 Kepler pixel) would easily be detected in the centroid motion analysis, but no centroid motion was detected between in and out of transit, indicating that these two stars were not responsible for the observed events. Six additional sources were detected at J within 15'' of the primary target. But all of these sources are δ_J > 8 mag fainter than K00137 and could not produce the observed transit events (i.e., even if the stars dimmed by 100%, the resulting transit in the blended photometry would not be deep enough to match the observed transit depths).

Zoom In Zoom Out Reset image size

Source detection completeness was evaluated by randomly inserted fake sources of various magnitudes in steps of 0.5 mag and at varying distances in steps of 1.0 FWHM from the primary target. Identification of sources was performed both automatically with the IDL version of DAOPhot and by eye. Magnitude detection limits were set when a source was not detected by the automated FIND routine or was not detected by eye. Within a distance of 1–2 FWHM, the automated finding routine often failed even though the eye could discern two sources, but beyond that distance the two methods agreed well. A summary of the detection efficiency as a function of distance from the primary star is given in Table 1.

A major source of false-positive planet indication in the Kepler data is background eclipsing binary stars within the photometric aperture of Kepler-18, which, when diluted by Kepler-18 itself, can produce a planetary-size transit signal. We perform a direct measurement of the source location via difference images. Difference image analysis takes the difference between average in-transit pixel images and average out-of-transit images. Barring pixel-level systematics, the pixels with the highest flux in the difference image will form a star image at the location of the transiting object, with an amplitude equal to the depth of the transit. Performing a fit of the Kepler pixel response function (PRF; Bryson et al. 2010) to both the difference and out-of-transit images quantifies the offset of the transit source from Kepler-18. Difference image analysis is vulnerable to various systematics due to crowding and PRF errors which will bias the result (Bryson et al. 2010). These types of biases will vary from quarter to quarter. We ameliorate these biases by computing the uncertainty-weighted robust average of the source locations over available quarters. Table 2 gives the offsets of the transit signal source from Kepler-18 averaged over quarters 1 through 8 for all three-planet candidates. We see that the average offsets are within 1σ of Kepler-18, with Kepler-18b being just over 2σ. From all of these lines of evidence, we conclude that there are no other objects within 4 mag near K00137 from 005 out to 15''.

After completion of the reconnaissance spectroscopy and high-resolution imaging, K00137 showed no evidence that might refute the planetary nature of the transit event in the Kepler light curve and the target was approved for high precision RV observations. We obtained 14 relative velocities with the Keck I High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) between 2009 September 1 and 2010 August 28 UT. We used the same spectrograph configuration as is normally used for precise Doppler measurements of solar-type stars (cf. Marcy et al. 2008; Cochran et al. 2002). Decker B5 (087 × 35) was used for the first five RV observations taken in 2009, and decker C2 (087 × 14'') was used for all of the later spectra. Exposures taken with the B5 decker entrance to the HIRES spectrometer suffer RV errors of up to ±15 m s⁻¹ when the moon is full and the relative Doppler shift of the star and solar spectrum is less than 10 km s⁻¹ (i.e., within a line width). This 15 m s⁻¹ error was determined from 10 stars for which observations were taken with and without moonlight subtraction. The first five RV measurements of Kepler-18 were made with the B5 decker and hence suffer from such errors. Decker C2 is long enough to permit us to do accurate subtraction of any moonlit sky spectrum. Sky subtraction was not possible with the spectra taken with decker B5. The iodine absorption cell was used as the velocity metric. This HIRES configuration can give a velocity precision as good as 1.0 m s⁻¹, depending on the stellar spectral type, the stellar rotation velocity Vsin i, and on the S/N of the observation. The exposure times for the K00137 spectra ranged from 2300 to 2700 s, and the final S/Ns ranged from 63 to 80 pixel⁻¹, or 127 to 162 per 4.1 pixel resolution element. The Doppler RV analysis algorithm (Butler et al. 1996; Johnson et al. 2009) computes an uncertainty in each data point from the variance about the mean of the individual spectral chunks into which the spectrum is divided. Main-sequence stars having T_eff near 5400 K have been measured for precise RVs for roughly 100 stars using the same HIRES instrument. Such stars typically show an additional noise of 2 m s⁻¹, caused by surface velocity fields and instrumental effects. We have modeled this additional intrinsic short-term stellar variability by adding that 2.0 m s⁻¹ "jitter" in quadrature to the uncertainties of each RV data point computed by the RV code. The measured HIRES relative velocities are given in Table 3. These velocity measurements are shown as the black points in Figure 4. The solid blue line in this figure is the model fit from the joint Markov Chain Monte Carlo (MCMC) solution of the RV data and the light curve, presented in Section 4. The rms scatter of the RV observations around this model fit is 4.3 m s⁻¹, whereas the rms scatter of the raw RVs is 9.6 m s⁻¹. The inferred planetary masses from this solution are given in the first row of Table 8. Also shown as the dashed red line in Figure 4 is a two-planet RV solution that includes only Kepler-18c (K00137.01) and Kepler-18d (K00137.02). The black line in Figure 4 shows the velocities from the full multi-body dynamical solution discussed in Section 6. The masses from this dynamical solution are adopted in Section 6 as our best determination of the planet masses.

Zoom In Zoom Out Reset image size

We determined stellar parameters using the local thermodynamic equilibrium line analysis and spectral synthesis code MOOG²⁷ (Sneden 1973), together with a grid of Kurucz²⁸ ATLAS9 model atmospheres. The method used is virtually identical to that described in Brugamyer et al. (2011). We analyzed a spectrum of Kepler-18 obtained with the Keck I HIRES spectrograph on 2010 August 27 UT. We first measured the equivalent widths of a carefully selected list of 53 neutral iron lines and 13 singly ionized iron lines in a spectrum of the Jovian satellite Ganymede, taken using the same instrumental setup and configuration as that used for Kepler-18. MOOG force-fits abundances to match these measured equivalent widths, using declared atomic line parameters. By assuming excitation equilibrium, we constrained the stellar temperature by eliminating any trends with excitation potential; assuming ionization equilibrium, we constrained the stellar surface gravity by forcing the derived iron abundance using neutral lines to match that of singly ionized lines. The microturbulent velocity ξ was constrained by eliminating any trend with reduced equivalent width (W_λ/λ). Our derived stellar parameters for the Sun (using our Ganymede spectrum) are as follows: T_eff = 5785 ± 70 K, log g = 4.54 ± 0.09 dex, microturbulent velocity ξ = 1.17 ± 0.06 km s⁻¹, and Fe abundance log () = 7.54 ± 0.05 dex.

We repeated the process described above for the spectrum of Kepler-18. We then took the difference, on a line-by-line basis, of the derived iron abundance from each line. Our resulting iron abundance is therefore differential with respect to the Sun. To estimate the rotational velocity of Kepler-18, we synthesized three 5 Å wide spectral regions in the range 5640–5690 Å and adjusted the Gaussian and rotational broadening parameters until the best fit (by eye) was found to the observed spectrum. The results of our analysis yield the following stellar parameters for Kepler-18: T_eff = 5345 ± 100 K, log g = 4.31 ± 0.12, ξ = 1.09 ± 0.08 km s⁻¹, [Fe/H] = +0.19 ± 0.06, and Vsin i < 4 km s⁻¹. The sky-projected rotational velocity of the star is very small. For such small Vsin i values, disentangling line-broadening effects due to the instrument, macroturbulence, and rotation is difficult, at best, and requires higher signal to noise and resolution than our spectrum offers. In our MOOG analysis, we have therefore chosen to quote an upper limit for Vsin i, which we estimate by assuming that all broadening is due to stellar rotation.

A completely independent analysis of the same spectrum was done using the stellar spectral synthesis package SME (Valenti & Piskunov 1996; Valenti & Fischer 2005). This analysis of Kepler-18 gives T_eff = 5383 ± 44 K, log g = 4.41 ± 0.10, [Fe/H] = +0.20 ± 0.04, and Vsin i = 0.4 ± 0.5 km s⁻¹. The agreement between the MOOG and the SME analyses of this star is excellent. We have arbitrarily selected the MOOG parameters as our adopted values.

We used a transit light curve model multiplied by instrumental decorrelation functions to measure the transit parameters and their uncertainties from the Spitzer data as described in Désert et al. (2011b). We compute the transit light curves with the IDL transit routine OCCULTSMALL from Mandel & Agol (2002). In the present case, this function depends on one parameter: the planet-to-star radius ratio R_p/R_⋆. The orbital semimajor axis to stellar radius ratio (system scale) a/R_⋆, the impact parameter b, and the time of mid-transit T₀ are fixed to the values derived from the Kepler light curves. The limb-darkening coefficients are set to zero since these Spitzer light curves do not have enough photometric precision to detect the curvature in the transit light curve that would be produced by limb darkening.

The Spitzer/IRAC photometry is known to be systematically affected by the so-called pixel-phase effect (see, e.g., Charbonneau et al. 2005; Knutson et al. 2008). This effect is seen as oscillations in the measured fluxes with a period of approximately 70 minutes (period of the telescope pointing jitter) and an amplitude of approximately 2% peak to peak. We decorrelated our signal in each channel using a linear function of time for the baseline (two parameters) and a quadratic function of the PSF position (four parameters) to correct the data for each channel. We checked that adding parameters to the correction function of the PSF position (as in Désert et al. 2009) does not improve the fit significantly. We performed a simultaneous Levenberg–Marquardt least-squares fit (Markwardt 2009) to the data to determine the transit and instrumental model parameters (seven in total). The errors on each photometric point were assumed to be identical and were set to the rms of the residuals of the initial best fit obtained. To obtain an estimate of the correlated and systematic errors (Pont et al. 2006) in our measurements, we use the residual permutation bootstrap, or "Prayer Bead," method as described in Désert et al. (2009). In this method, the residuals of the initial fit are shifted systematically and sequentially by one frame, and then added to the transit light curve model before fitting again. We allow asymmetric error bars spanning 34% of the points above and below the median of the distributions to derive the 1σ uncertainties for each parameter as described in Désert et al. (2011a).

We measured the transit depths at 4.5 μm of 1521⁺²⁷⁷ _{− 245} ppm for Kepler-18c and 4083⁺²⁹⁸ _{− 285} for Kepler-18d. The values we measure for the transit depths in the Spitzer bandpass are in agreement at the 2σ level compared to the Kepler bandpass. This indicates that the transit depths of Kepler-18c and Kepler-18d are only weakly dependent on wavelength, if at all. This is in agreement with expectations for a dark planetary object, and indicates that there is no significant contamination from any nearby unresolved star of significantly different color. As discussed in Section 7, this is important evidence that the transit signals arise from planets and not from eclipsing stars blended with additional light.

These Spitzer observations provide a useful constraint on the kinds of false positives (blends) that may be mimicking the transit signal, such as eclipsing binaries blended with the target. If Kepler-18 were blended with an unresolved eclipsing binary of later spectral type that manages to reproduce the transit depth in the Kepler passband, the predicted depth at 4.5 μm would be expected to be larger because of the higher flux of the contaminant at longer wavelengths compared to Kepler-18. Since the transit depth we measure in the near-infrared is about the same as in the optical, this argues against blends composed of stars of much later spectral type. Based on model isochrones, the properties of the target star, and the transit depths measured with Spitzer at the 3σ level, we determine a lower limit to the blend masses of 0.86 and 0.79 M_☉ for Kepler-18c and Kepler-18d, respectively.

Kepler-18 was observed at the 29.4244 minute LC rate for the first three observing periods Q0–Q2. After the transits of K00137.01 and K00137.02 were detected, Kepler-18 was observed at the 58.85 s SC in order to facilitate the measurement of possible TTVs. These SC data were used for Q3 through Q7. Transit times and their errors for each member of the Kepler-18 system were determined through an iterative procedure; a single step of this procedure is described as follows.

First, the detrended photometric data within four transit durations at each epoch were shifted by the current best-fit mid-transit times, and the light curves were folded on these transit times to form a transit template. This template was then fit with a transit light curve model (Mandel & Agol 2002). Next, at each individual epoch, this light curve template was shifted in time and compared to the data by computing the standard χ² statistic. This statistic was computed for a dense, uniform sample of mid-transit times centered on the current best-fit time and spanning approximately five current best-fit timing errors. The time t₀, corresponding to the minimum χ², was recorded as the new best-fit mid-transit time. The χ² data were then fit with a quadratic function of time, χ²(t) = C(t − t₀)² + χ²₀. By choosing this functional form, we are assuming that the posterior likelihood is well described by a Gaussian function of the mid-transit time. The fidelity of the quadratic approximation was verified visually at each epoch. The timing error, σ_t, was found by solving for the time, t = t₀ + σ_t, at which the quadratic model indicated a Δχ² = 1. This corresponds to σ_t = C^−1/2.

This iterative procedure of template generation followed by timing estimate generally converged to the final values in two to three steps. The measured transit times for each observed transit of Kepler-18c are given in Table 5, and the measured transit times for Kepler-18d are given in Table 6.

A drift in the orbital inclination was measured for each planet by augmenting the nominal transit model with an epoch-dependent linear inclination (viz., i(E) = i(0) + (Δi/ΔE) × E), and then fitting for the linear coefficient. In this fit, transit times were fixed to their best-fit values while the remaining transit parameters were allowed to vary. The best-fitting solution was found by minimizing the standard χ² metric. The uncertainty was estimated by fitting a multivariate Gaussian to a sampling of the posterior parameter distribution. The drifts (Δi/ΔE) were found to be [10 ± 9, 4 ± 22, −3 ± 10] × 10⁻⁴ degrees per epoch for planets b, c, and d, respectively. Thus, we detect no secular inclination drift.

To model the transit times and radial velocities, we applied the numerical routines that were previously used for Kepler-9 and -11 (Fabrycky 2010; Holman et al. 2010; Lissauer et al. 2011a). That is, a Levenberg–Marquardt χ²-minimization routine drove three-planet dynamical integrations, which calculated the transit times resulting from given orbital parameters. As in Lissauer et al. (2011a), we chose the parameters (m_p, P, T₀, ecos ω, esin ω) for each planet: mass, orbital period, transit phase, and eccentricity vector components. These parameters are osculating Jacobian orbital elements at the epoch 2455168.0 [BJD]. As in Holman et al. (2010), we also used the integrations to calculate radial velocities of the star at each of the observed times. We assumed M_⋆ = 0.95 M_☉, assumed the orbits are edge-on and coplanar, and neglected light-travel time effects in these numerical calculations.

The main data constraining the orbital model are the transit times of Table 5 (75 data points) and Table 6 (37 data points). In some of the calculations reported below, we allow K00137.03 to interact dynamically, but we include only its first and last observed transits, at t = 2454955.99237 ± 0.00823 and 2455530.77178 ± 0.00174 (two data points), as a way of keeping its period and phase fixed at observed values. We have measured all of its transit times, but we do not report or analyze them here: their individual signal to noise is low, they are not apparently constant, and we have not found a consistent dynamical solution to date. Perhaps future work will show a fourth (non-transiting) planet is required to fit the transit times of Kepler-18b. In the meantime we proceed with fitting the three transiting planets with a focus on the transit timing constraints for Kepler-18c and Kepler-18d.

By fitting only the transit times, allowing P and T₀ for each of the three planets to vary (six free parameters), and setting dynamical interactions to zero, we find χ² = 750.6 for 114 transit time measurements. This is clearly an unacceptable fit, and the obvious timing patterns call for a dynamical model. The observed ("O") residuals to this calculated ("C") ephemeris are called the O − C values and are plotted in Figure 7 along with the preferred dynamical model described below. As in the case of Kepler-11b/c, we clearly see that the source of the variations of transit times for Kepler-18c and Kepler-18d is their nearly 2:1 resonance. The expected variation occurs on a timescale:

$$\begin{equation} P_{{\rm TTV}} = 1/(2/P_d-1/P_c) = 268{\rm \,days}, \end{equation} \tag{ 1 }$$

the time it takes the line of conjunctions to sweep around inertial space, through both the line of sight and the apsidal lines of these planets (on the approximation that precession can be ignored on this timescale); see Agol et al. (2005). In Figure 8, we plot the periodogram of the O − C values for Kepler-18c and Kepler-18d. The peaks occur at 1/(260.4 ± 3.3 days) and 1/(265.1 ± 4.7 days) for Kepler-18c and Kepler-18d, respectively, very close to the simple expectation given above, which uniquely identifies the dynamical mechanism for TTVs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Before moving on to a full solution, we tried fitting the radial velocities (14 data points, Table 3) with the above-determined periods and phases, with circular orbits, allowing only the planetary masses to vary. The solution, listed in the second row of Table 8, was [m_b, m_c, m_d] = [12 ± 5, 15 ± 5, 28 ± 7] M_⊕; the χ² = 9.7 for 10 degrees of freedom (14 RV data points, minus 3 K-amplitudes, minus a constant RV offset). These masses are well within the 1σ error bars of the MCMC solution to the light curve and RVs presented in Section 4. Therefore, this model is sufficient to explain the radial velocities, and each planet is detected, but only marginally so for planet b. In particular, if we hold the mass of b at zero, the masses of the others become [m_c, m_d] = [18 ± 5, 24 ± 7] M_⊕ and the χ² = 15.6 for 10 degrees of freedom. These values are within about 0.5σ of the masses determined from the joint MCMC and RV solution presented earlier in Section 4.

We may also attempt to constrain the masses and orbital elements using only the transit times. Naturally, this requires full dynamical integrations, in which the planets cannot remain on circular orbits. However, for planet b we assumed a circular orbit at the dynamical epoch, since we are not attempting to fit its transit time variations. All of the other orbital parameters and masses were free to vary. The resulting χ² is 88.5 for 101 degrees of freedom (114 transit times, minus 5 parameters for planets c and d, minus 3 for K00137.03), which is quite acceptable. To be compared with the RV solution, the solved-for masses were [m_b, m_c, m_d] = [18 ± 9, 17.3 ± 1.7, 15.8 ± 1.3] M_⊕, shown in the third row of Table 8. These solutions had extremely low eccentricities (e < 0.003) for Kepler-18c and Kepler-18d. As above, the innermost planet is only very marginally detected. In contrast to the RV solution, however, the masses of the interacting planets are very precisely pinned down by the large variations seen in Figure 7.

Finally, we generated a joint solution to the transit times and radial velocities. Graphically, this solution is given by Figure 7. A χ² = 103.4 for 114 degrees of freedom is achieved, an excellent fit to the data. The orbital parameters and their formal errors (the output of the Levenberg–Marquardt algorithm) are given in Table 7, and the planetary masses are given in the fourth row of Table 8.

In the previous subsection we found no drift in the inclinations of the three planets. We can use the 3σ upper limits on |Δi/ΔE| of [37, 70, 33] ×  10⁻⁴ degrees per epoch for Kepler-18b, c, and d to place limits on their mutual inclination (Miralda-Escudé 2002). To do this, we measure the value of |Δi/ΔE| seen in numerical simulation, which depends nearly linearly on the difference in nodal angle on the sky of two planets (Ballard et al. 2010). Taking the masses and orbits from the best-fit TTV/RV solution above, we simulated Kepler-18c and d with 1° (and 10°) of mutual inclination, we find |Δi/ΔE|_c = 1.8 × 10⁻⁴ degree per epoch (16 × 10⁻⁴ degree per epoch) and |Δi/ΔE|_d = 2.9 × 10⁻⁴ degree per epoch (26 × 10⁻⁴ degree per epoch). The interaction between Kepler-18b and c is also potentially observable; simulating their orbits with 1° (10°) of mutual inclination, we find |Δi/ΔE|_b = 1.3(11) × 10⁻⁴ degree per epoch and |Δi/ΔE|_c = 0.82(7.4) × 10⁻⁴ degree per epoch. The interaction between Kepler-18b and d is considerably weaker, with a drift of ∼2 × 10⁻⁴ degree per epoch even for 10° mutual inclination; therefore we ignore it when setting mutual inclination limits. To find the 3σ upper limit to the nodal difference of c and d, we compare the 10° calculated value of |Δi/ΔE|_d to its observational 3σ upper limit, so we find |Ω_c − Ω_d| < 13°. The inclination difference in the complementary direction is only i_d − i_c = 040 ± 024, so the limit on the true mutual inclination is i_cd < 13°. The nodal difference of b and c is more poorly constrained, as the calculated value of |Δi/ΔE| for both planets in the 10° mutual inclination simulation is small compared to its observational 3σ upper limit. Thus, a moderate (∼20°) mutual inclination is permissible, which would be large enough to affect the timing fits. In that case, the model should self-consistently fit the radial velocities, transit times, and transit durations. However, we defer dynamical interpretation of the mutual inclination of the inner planet to the others until more data are gathered, which should tighten the limit.

The a priori frequency of stars in the background or foreground of the target that are orbited by a transiting planet and are capable of mimicking the K00137.03 signal may be estimated from the number density of stars in the vicinity of Kepler-18, the area around the target in which such stars would go undetected in our high-resolution imaging, and the frequency of transiting planets with the appropriate characteristics. To obtain the number density (stars per square degree) we make use of the Galactic structure models of Robin et al. (2003), and we perform this calculation in half-magnitude bins, as shown in Table 9. For each bin we further restrict the star counts using the constraints on the mass of the secondaries supplied by BLENDER (see Figure 10), and the eccentricity limit for transiting planets discussed above (e ⩽ 0.3). These mass ranges are listed in Column 3, and the resulting densities appear in Column 4. Bins with no entries correspond to brightness ranges excluded by BLENDER. The maximum angular separation (ρ_max) at which stars of each brightness would escape detection in our AO/speckle imaging is shown in Column 5 (see Table 1 and Section 3.2). The result for the number of stars in each magnitude bin is given in Column 6, in units of 10⁻⁶.

To estimate the frequency of transiting planets that might be expected to orbit these stars (and lead to a false positive), we rely on the results from Borucki et al. (2011b), who reported a total of 1235 planet candidates among the 156,453 Kepler targets observed during the first four months of the mission. These signals have not yet been confirmed to be caused by planets, and therefore remain candidates until they can be thoroughly followed up. However, the rate of false positives in this sample is expected to be quite small (∼10% or less; see Morton & Johnson 2011), so our results will not be significantly affected by the assumption that all of the candidates are planets. We further assume that the census of Borucki et al. (2011b) is largely complete. After accounting for the additional BLENDER constraint on the range of planet sizes for blends of this kind (tertiaries of 0.32–1.96 R_Jup), we find that the transiting planet frequency is f_planet = 374/156, 453 = 0.0024. Multiplying this frequency by the star counts in Column 6 of Table 9, we arrive at the blend frequencies listed in Column 7, which are added up in the "Totals" line of the table (12.7 × 10⁻⁶).

Next we address the frequency of hierarchical triples, that is, physically associated companions to the target that are orbited by a larger transiting planet able to mimic the signal. The rate of occurrence of this kind of false positive may be estimated by considering the overall frequency of binary stars (34% according to Raghavan et al. 2010) along with constraints on the mass range of such companions and how often they would be orbited by a transiting planet of the right size (0.43–0.53 R_Jup; see Figure 11). The mass constraints include not only those coming directly from BLENDER, but also take into consideration the color and brightness limits mentioned earlier. We performed this calculation in a Monte Carlo fashion, drawing the secondary stars from the mass ratio distribution reported by Raghavan et al. (2010), and the transiting planet eccentricities from the actual distribution of known transiting planets (http://exoplanet.eu/), with repetition. Planet frequencies in the appropriate radius range were taken as before from Borucki et al. (2011b). The result is a frequency of hierarchical triples of 4.4 × 10⁻⁶, which we list at the bottom of Table 9.

Combining this estimate with that of background/foreground star+planet pairs described previously, we arrive at a total blend frequency of BF = (12.7 + 4.4) × 10⁻⁶ ≈ 1.7 × 10⁻⁵, which represents the a priori likelihood of a false positive. From a Bayesian point of view analogous to that adopted to validate previous Kepler candidates, our confidence in the planetary nature of K00137.03 will depend on how this likelihood compares to the a priori likelihood of a true transiting planet, addressed below.

The BF of 1.7 × 10⁻⁵ corresponds to false-positive scenarios giving fits to the Kepler photometry that are within 3σ of the best planet fit. We use a similar criterion to estimate the a priori transiting planet frequency by counting the KOIs in the Borucki et al. (2011b) sample that have radii within 3σ of the value determined from the best fit to K00137.03 (R_p = 2.00 ± 0.10 R_⊕). We find 284 that are within this range, giving a planet frequency PF = 284/156, 453 = 1.8 × 10⁻³.

This estimate does not account for the fact that the geometric transit probability of a planet is significantly increased by the presence of additional planets in the system (Kepler-18c and Kepler-18d in this case), given that mutual inclination angles in systems with multiple transiting planets have been found be relatively small (typically 1°–4°; Lissauer et al. 2011b). Furthermore, a planet with the period of K00137.03 would be interior to the other two, further boosting the chances that it would transit. To incorporate this coplanarity effect, we have developed a Monte Carlo approach, described fully in the Appendix, in which we simulate randomly distributed reference planes and inclination dispersions around this plane, from which a weighted distribution of the inclination with respect to the line of sight for a third planet is calculated. Inclination angles relative to the random reference plane are assumed to follow a Rayleigh distribution (see Lissauer et al. 2011b). Although the known planets carry some information on the inclination dispersion, the probability of transit still depends somewhat on the allowed range of dispersion widths. When the assumed prior for the inclination dispersion is uniform up to 4°, following Lissauer et al. (2011b), we find that the flatness of the system results in a very significant increase in the transit probability for K00137.03 from 11.7% to 97%. To be conservative, we adopt a larger range of possible inclination dispersions from 0° to 10°, motivated by the upper limit from the similar Kepler-9 system (Holman et al. 2010). With this prior, the transit probability for K00137.03 becomes 84%, or an increase by a factor of ∼7 over the case of a single transiting planet.

Thus, the likelihood of a planet is more than 700 times greater (PF/BF = 0.013/1.7 × 10⁻⁵ ≈ 700) than that of a false positive, which we consider sufficient to validate K00137.03 as a true planet with a high degree of confidence. We designate this planet Kepler-18b. We note that our PF calculation assumes that the 1235 candidates cataloged by Borucki et al. (2011b) are all true planets. If we were to suppose conservatively that as many as 50% are false positives (an unlikely proposition that is also inconsistent with other evidence; see Borucki et al. 2011b; Howard et al. 2010), the planet likelihood would still be ∼350 times greater than the likelihood of a blend, implying a false alarm rate sufficiently small to validate the candidate.