KEPLER 453 b—THE 10th KEPLER TRANSITING CIRCUMBINARY PLANET^*

Identification of Kepler-453 b as a circumbinary planet candidate was made by visual inspection of a subset of the Kepler eclipsing binary star light curves (Prsa et al. 2011; Slawson et al. 2011), in particular, those with orbital periods greater than ∼1 day that show both primary and secondary eclipses. Once identified as a circumbinary planet candidate, the system was given the Kepler Object of Interest (KOI) number 3151, though this same system had previously been named KOI-1451 and rejected as a false-positive; thus, KOI-1451 is formally the correct KOI number. Because we identified the binary as hosting a planet, the target was placed on Short Cadence (∼1 min integration sampling) starting in Quarter 13 (2012 March 29). The NExScI Exoplanet Archive lists for KOI-1451 a Kepler magnitude of 13.552, temperature of ${T}_{\mathrm{eff}}=5618$ K, and a surface gravity of $\mathrm{log}g=4.586$, while the Mikulski Archive for Space Telescopes (MAST) reports ${T}_{\mathrm{eff}}=5425$ K and $\mathrm{log}g=4.803$. Note that in their study on the rate of occurrence of circumbinary planets, Armstrong et al. (2014) independently identified KIC 9632895 as a circumbinary planet candidate.

The orbital period of the binary is 27.322 days and the depth of the primary eclipse is ∼8%. The secondary eclipse is shallow, only $\sim 0.25\%$, and flat-bottomed, indicating a total eclipse. Since the primary eclipse is likely to be total (i.e., annular), the relative depth of the eclipses suggests that the secondary star contributes only a small fraction to the total luminosity and is likely a low-mass small star.

The upper panel of Figure 1 shows a one-year long section of the light curve, after normalizing each Quarter with a simple cubic polynomial. The light curve exhibits quasi-periodic variations of ∼0.5% on a timescale of tens of days, with the largest peak-to-peak variation being ∼1.5%. The rms of the mildly detrended PDC-MAP light curve, after removing all of the eclipses, is 0.22%. We interpret these modulations as being caused by starspots on the primary star.

Zoom In Zoom Out Reset image size

The lower panels of Figure 1 show a sample of eclipse events, all in Short Cadence, along with our best-fit model curve. The leftmost is a secondary eclipse, followed by a planetary transit, then two primary eclipses. The timescale in each panel is identical and nicely demonstrates the much longer duration of the transit compared with that of the eclipses. The two primary eclipses are consecutive in time, yet the first shows a notable deviation in its residuals. This is the signature of a starspot on the primary star being occulted by the secondary star. Such events are common, but are not readily visible in Long Cadence eclipse profiles.

Figure 2 shows the phase-folded eclipse profiles, based on the Long Cadence Kepler observations that span 1470.5 days (from times −46.5 to 1424.0 in BJD–2,455,000). The preliminary model fit to the binary star data only (no planet included) is generally excellent, though there is a marked increase in the scatter of the residuals for the primary eclipse. These are due to the secondary star covering starspots (and possibly other features) on the primary star. Our ELC photodynamical model, described in Section 3, does not include starspots, thus the minimization of the residuals leads to a relatively large scatter above and below the best-fit model. If we were to use this preliminary fit, we would need to correct for a bias in the model—the starspots should induce upward-only residuals, and as this is not the case, the model is too shallow compared to the true eclipse depth. It is important to note that the starspot occultations lead not only to residuals in flux, but also to shifts in the measurements of the timing of the eclipses. We return to this point in Section 3 when describing our photodynamical modeling methodology.

Zoom In Zoom Out Reset image size

The center-of-eclipse times for the 49 primary and 50 secondary eclipses were measured in two ways: (i) by using the best-fit model eclipse profile (described in Section 3) as a template and sliding it to best match the individual eclipses; and (ii) using the technique described in Steffen et al. (2011) and Welsh et al. (2012), which uses a template created from a polynomial fit to the folded eclipse profile. The results were very comparable, so we adopt the former method because it should be less-sensitive to starspot-induced deviations in the mean eclipse profile. The mean uncertainty is 6.8 s for the primary eclipses (7.6 s for the 37 Long Cadence eclipses, 4.3 s for the 12 Short Cadence eclipses). The shallowness of the secondary eclipses made their timing measurements much more difficult: the mean uncertainty is 117.3 s (123.5 s for the 37 Long Cadence eclipses, 99.7 s for the 13 Short Cadence eclipses).

A linear ephemeris was derived from the eclipse times and a time series of Observed minus Computed ($O-C$) values was made. When compared with the local slope of the light curve surrounding each primary eclipse, the $O-C$ times show a significant anti-correlation. This anti-correlation indicates that the spin of the star is prograde with respect to the binary orbit (Mazeh et al. 2015 see also Sanchis-Ojeda et al. 2012). Such behavior was also seen in Kepler-47 (Orosz et al. 2012a). Using a linear fit to the primary star $O-C$ values versus the local slopes, the eclipse times were then statistically corrected for the starspot-induced timing shifts; the uncertainties in the times were increased (by ∼10 s in quadrature) to account for noise in this correction. The eclipse measurements made after date BJD–2,455,000 = 1015 used Short Cadence observations, though because of the boosting of the uncertainties in the starspot correction, the error bars are similar to those of the Long Cadence $O-C$. The secondary eclipse $O-C$ values were not correlated with their local slope, so no correction was made. A common linear ephemeris was then fit to both sets of eclipse times, and Figure 3 shows the resulting $O-C$ diagram. The rms of the $O-C$ residuals is 10.1 s for the primary and 100.0 s for the secondary, and both appear consistent with noise. This indicates that the circumbinary object has no measurable gravitational effect on the binary, at least on a timescale of a few years. Consequently, the object is of relatively low mass and only an upper limit on its mass can be robustly determined.

Zoom In Zoom Out Reset image size

The Kepler Simple Aperture Photometry "SAP" data were used throughout this paper, with the exception of measuring the starspot modulation rms amplitude. The light curve was detrended and normalized using the method described in J. A. Orosz et al. (2015, in preparation): windows around each eclipse event were kept and the rest of the light curve discarded. The duration of the windows were three times the width of the eclipse. For each window, a fifth-order Legendre polynomial was fit to the out-of-eclipse portion of the light curve. The eclipse was then restored and the data normalized by dividing by the polynomial. Any points with a Kepler pipeline Data Quality flag greater than 16 (indicative of some anomaly with the observation) were omitted prior to the detrending. The mean Combined Differential Photometric Precision (CDPP) over a 3 hr baseline as reported on MAST is 121 ± 34 ppm. Prior to the eclipse and transit fitting, the light curve was heavily detrended (to minimize starspot effects) and normalized. The rms scatter outside of the eclipse and transit events was measured to be 0.0159% (159 ppm), consistent with the CDPP value, and was a factor of 1.17 times larger than the mean SAP error bar. This difference suggests that either the SAP error bars are slightly underestimated or there are additional high-frequency variations that the detrending did not remove.

Kepler-453 was observed from the McDonald Observatory with the High Resolution Spectrograph (Tull 1998) on the Hobby–Eberly Telescope (HET), and with the Tull Coude Spectrograph (Tull et al. 1995) on the Harlan J. Smith 2.7 m Telescope (HJST). A total of 11 spectra were obtained in 2013, spanning 47 nights. The HRS spectra cover a wavelength range from 4780 to 6800 Å at a resolving power of R = 30,000 and a typical continuum signal-to-noise ratio (S/N) ratio of 40:1 at 5500 Å. The data taken with the Tull spectrograph span the entire optical spectrum and have a resolving power of R = 60,000 with a typical continuum S/N ratio around 20:1 at 5500 Å. We determined the RVs by cross-correlating the spectra with the RV standard star HD 182488 and eleven radial velocities were extracted. Note that a 0.181 $\mathrm{km}\ {{\rm{s}}}^{-1}$ zero-point offset was found between the two spectrographs and subtracted from the HJST Tull velocities. As expected from the light curve, the secondary star is sufficiently faint that it was not detected in the spectra. By injecting a synthetic spectrum signal into the observed spectra and attempting to recover that signal, an upper limit of ∼4% was found for any non-primary star light. This upper limit is over an order of magnitude higher than the secondary star's light contribution estimated via the photodynamical modeling described in Section 3, and therefore this is not valuable in determining the secondary star's flux contribution. However, it is interesting with regard to the presence of any contaminating third light. The radial velocities for the primary are listed in Table 1 and the ELC model fit (Section 3.1) is shown in Figure 4. Note that the velocities in Table 1 do not include the zero-point velocity offset correction. The uncertainties listed are from the data calibration, but when used in the fitting process the uncertainties were boosted by a factor of 1.6 to achieve a reduced ${\chi }^{2}$ near 1.0.

Several independent spectral analyses were carried out using the higher signal-to-noise HET spectra, including an Stellar Parameter Classification (SPC) analysis (Buchhave et al. 2012), an analysis similar to the one for KOI-126 (Carter et al. 2011), a MOOG-based analysis, and by fitting a theoretical template to the spectrum (Tal-Or et al. 2013). Each gave similar results, with the best fit ${T}_{\mathrm{eff}}$ ranging from 5480 to 5620 K, and metallicity $[{\rm{m}}/{\rm{H}}]$ ranging from −0.10 to +0.09. The spectra have relatively low S/N ratio for this type of spectral analysis and the SPC method has worked reliably on noisy data, hence we prefer its results: ${T}_{\mathrm{eff}}=5527\pm 50$ K, $[{\rm{m}}/{\rm{H}}]\;=\;+0.09\pm 0.08$, and $\mathrm{log}g=4.56\pm 0.10$. This surface gravity measurement was guided, but not constrained by, the photodynamical modeling value of $\mathrm{log}g=4.57$. The agreement in metallicity between methods was not as tight as for ${T}_{\mathrm{eff}}$, thus we conservatively inflate the uncertainty in both metallicity and ${T}_{\mathrm{eff}}$ to ±0.1 dex and ±100 K, respectively. In addition, the projected stellar rotation velocity, ${V}_{\mathrm{rot}}\mathrm{sin}i$, was measured to be 1.9 $\mathrm{km}\ {{\rm{s}}}^{-1}$.

Based on the values retrieved from MAST, the seasonal mean of the contamination for Kepler-453 is 7.0% ± 1%. Such contamination will dilute the eclipse and transit depths and could significantly bias the inferred radii. Therefore we undertook several direct imaging investigations to determine the amount of contamination.

Eighteen short exposures (10 s each) of the target were obtained in the SDSS-r filter using the LCOGT (Brown et al. 2013) FTN 2 m robotic telescope in Haleakala, Hawaii over the span of 4 nights in 2012 March. The data were median-combined and provided an image with a FWHM of 1.70 arcsec. According to the KIC there should be a ${\rm{\Delta }}\mathrm{Kp}=2.8$ mag fainter star located 4.5 arcsec south of the target (KID 9632896; g-r = 4.4 mag). This star was not detected, although fainter stars further away were clearly detected. The closest star to the target in the image is about 15 arcsec away and has r = 18.3 mag according to the KIC. That star is also the closest star to the target detected by the UBV survey of Everett et al. (2012).

To better understand any photometric contamination due to non-target light captured in the Kepler aperture, we observed Kepler-453 through J, H, and Ks filters on the WIYN High-Resolution Infrared Camera (WHIRC—Meixner et al. 2010) at the Kitt Peak National Observatory on 2013 October 19 UT. We employed a standard five-point dithering pattern and 30 s exposure times in each filter. Unfortunately, the conditions were not photometric and the seeing was ∼0.9 arcsec. We estimated detection limits in each filter following the procedure of Adams et al. (2012). We used the IDL aper routine to measure the contribution from the target's PSF in non-overlapping, concentric annuli centered on the star. We define a detection limit as a 5σ signal above the measured stellar PSF background in each annulus. The innermost annulus in each filter is defined as the measured FWHM of the stellar PSF. The detection limits are presented in Table 2. The expected nearby star KIC 9632896 was again not detected. A star ∼8.5 arcsec away toward the SW was detected in the J band (and weakly in the H band and barely in the K band), with ∼0.5% the brightness of Kepler-453.

Table 2.  WIYN/WHIRC Detection Limits

			${\rm{\Delta }}(\mathrm{mag})$
Filter	FWHM ['']	15	20	25	50
J	0.87	2.6	3.8	4.1	4.3
H	0.97	2.1	3.3	3.5	3.6
Ks	0.87	2.8	3.2	3.3	3.4

Note. Tabulated values give the detection limits in relative magnitudes at four radial distances from the source.

Download table as:  ASCIITypeset image

Finally, an even deeper J-band observation was taken with UKIRT as part of the J-band survey of the Kepler field for the UKIDDS Survey. The J-band detections were converted into expected Kepler magnitudes via the transformations from Howell et al. (2012). The missing KIC star was again not detected. The star 8.5 arcsec away was detected (expected Kepmag = 19.12, based on its J-band magnitude of 17.64), along with a very faint star 7.26 arcsec due south of Kepler-453, with an expected Kepmag of 19.85. In summary, the star reported as KIC 9632896 was not detected, and no stars were detected that could contribute any significant contamination within the Kepler aperture.

An upgraded "photodynamical" version of the ELC code (Orosz & Hauschildt 2000) was developed. The upgrade allows for dynamics, instead of Keplerian kinematics, by integrating the Newtonian equations of motion under the assumption of point-masses and Newtonian gravity with general relativistic corrections (Mardling & Lin 2002; see also Hilditch 2001 and Ragozzine & Wolf 2009). A 12th order Gaussian Runge–Kutta symplectic integrator, based on the code of Hairer & Hairer (2003) was used.

Transits and eclipses are modeled using the prescription of Mandel & Agol (2002), replacing the Gimenez method (Gimenez 2006) formerly used in ELC. Quadratic limb darkening is used, following the method of Kipping (2013), which more naturally handles correlation and limits on the coefficients. A total of 26 parameters are specified in the model: the five standard Keplerian orbital parameters for each orbit ($P,{T}_{c},i,e,\omega$), the three masses, the three radii, the stellar temperature ratio, two quadratic limb darkening coefficients for each star, the longitudinal nodal angle Ω of the planet's orbit, and four seasonal contamination levels. Since the Keplerian parameters evolve rapidly with time, their values presented in Table 3 are instantaneous "osculating" values valid at the reference epoch. In particular, the time of conjunction, T_conj, is the time of conjunction of the body and the barycenter; for the binary this is approximately a primary mid-eclipse time, but for the planet it need not be close to an actual transit. Furthermore, while this fiducial time serves to set the position of the planet in its orbit at the reference epoch, it does not accurately define the planet's position when the model is integrated away from the reference time. To faithfully reproduce our solution the instantaneous positions and velocities given in Table 4 should be used. Although the nodal longitude of the binary is not measurable with our data, it must be specified; we set the angle to be zero. In the actual fitting, ratios and other combinations of parameters are often more robustly constrained by the data, e.g., mass ratios, semimajor axis/radius, temperature ratio, $e\mathrm{cos}\omega$, and $e\mathrm{sin}\omega$, and these are therefore used. Both a genetic algorithm and a Markov Chain Monte Carlo optimization were used to explore parameter space and find the least-squares best solution along with the uncertainties (defined as the interval for which the ${\chi }^{2}$ is less than the minimum ${\chi }^{2}+1$ for each marginalized parameter).

Table 3.  Kepler-453 System Parameters

Parameter	Best Fit	Uncertainty	Units
Binary
M1	0.944	±0.010	${M}_{\odot }$
M2	0.1951	±0.0020	${M}_{\odot }$
R1	0.833	±0.011	${R}_{\odot }$
R2	0.2150	±0.0014	${R}_{\odot }$
Orbital period	27.322037	±0.000017	days
${T}_{\mathrm{conj}}$	−34.574012	±0.000060	BJDa
Inclination	90.266	±0.052	degrees
$e\ \mathrm{sin}\omega$	−0.0520	±0.0037	⋯
$e\ \mathrm{cos}\omega$	−0.006339	±0.000016	⋯
Eccentricity	0.0524	±0.0037	⋯
Arg. periastron	263.05	±0.48	degrees
Semimajor axis	0.18539	±0.00066	AU
Limb darkeningbq11	0.41	±0.06	⋯
Primary q21	0.07	±0.11	⋯
Secondary q12	0.33	±0.11	⋯
Secondary q22	0.07	±0.07	⋯
T1	5527	±100	Kc
${T}_{2}/{T}_{1}$	0.5837	±0.0150	K
T2	3226	±100	K
$\mathrm{log}{g}_{1}$	4.571	±0.015	cgs
$\mathrm{log}{g}_{2}$	5.0630	±0.0050	cgs
Contamination s0	0.021	±0.027	⋯
Contamination s1	0.024	±0.027	⋯
Contamination s2	0.034	±0.045	⋯
Contamination s3	0.049	±0.031	⋯
Planet
Mp	0.2	±16.0	${M}_{\oplus }$
Rp	6.204	±0.039	${R}_{\oplus }$
Orbital period	240.503	±0.053	days
${T}_{\mathrm{conj}}$	69.020	±0.054	BJDa
Inclination	89.4429	±0.0091	degrees
$e\ \mathrm{sin}\omega$	−0.00322	±0.00023	⋯
$e\ \mathrm{cos}\omega$	−0.03575	±0.0088	⋯
Eccentricity	0.0359	±0.0088	⋯
Arg. periastron	185.1	±3.7	degrees
Semimajor axis	0.7903	±0.0028	AU
Nodal longitude Ω	2.103	±0.055	degrees
Mutual inclination	2.258	±0.039	degrees

Notes. The reported Keplerian parameters are the instantaneous (osculating) values at the reference epoch Tref = 2,454,964 BJD.

^aDate with respect to BJD–2,455,000. ^bThese are the quadratic limb darkening coefficients for the primary and secondary stars, following Kipping (2013). ^cThe primary star ${T}_{\mathrm{eff}}$ is spectroscopically determined.

Download table as:  ASCIITypeset image

Table 4.  Barycentric Initial Dynamical Parameters in Cartesian Coordinates

Parameter	Primary Star	Secondary Star	Planet b
$M\;\;({M}_{\odot })$	$9.436056276889058\times {10}^{-1}$	$1.950599079365290\times {10}^{-1}$	$4.358779010359746\times {10}^{-7}$
x (AU)	$-9.735136288118057\times {10}^{-3}$	$4.709297415736234\times {10}^{-2}$	$4.087938246546758\times {10}^{-1}$
y (AU)	$1.482084767646292\times {10}^{-4}$	$-7.169795822100196\times {10}^{-4}$	$8.301755937288406\times {10}^{-3}$
z (AU)	$-3.194184618234016\times {10}^{-2}$	$1.545207660656327\times {10}^{-1}$	$-6.896119623633258\times {10}^{-1}$
vx (AU day−1)	$6.615322808726673\times {10}^{-3}$	$-3.200177661512062\times {10}^{-2}$	$1.782393075783825\times {10}^{-2}$
vy (AU day−1)	$9.677067331539996\times {10}^{-6}$	$-4.681465367384790\times {10}^{-5}$	$7.498479158542154\times {10}^{-4}$
vz (AU day−1)	$-2.085636760223096\times {10}^{-3}$	$1.008928146684059\times {10}^{-2}$	$9.796765804214942\times {10}^{-3}$

Note. Valid at the reference epoch of 2,454,964.00 BJD (or time = −36.00 if using BJD–2,455,000). The (x, y) axes are in the plane of the sky, and the z direction points toward the observer.

Download table as:  ASCIITypeset image

We initially fit the normalized and detrended Kepler Long Cadence light curve plus the 11 radial velocity values. Because the normalized flat sections of the light curve far from any eclipse or transit contain no information on the system parameters, they were omitted from the fitting process. The uncertainties in the Kepler data were increased to account for additional "noise" that occurs when a starspot is eclipsed (see Figure 2). An increase of 1.44 was chosen to give a reduced ${\chi }^{2}\equiv 1.0$, based on a preliminary model fit to the eclipses. Our initial fit was acceptable, yielding a ${\chi }^{2}$ of 11,114.6 for 11,139 data points. However, we found the results unsatisfactory for the following reason: the best fit estimate for the mass of the planet was tightly constrained to zero mass. When the non-negative mass constraint was lifted, the preferred solution fell to roughly $-90\ {M}_{\oplus }$, at an unnerving 3σ confidence.

The cause of this non-physical result was traced to starspots on the primary star. As seen in Figures 1 and 2, starspots can distort the eclipse profile shape and induce spurious eclipse timing variations (ETVs). Thus starspot occultations "contaminate" the primary eclipses and compromise their reliability. Since the constraint on the mass of the planet arises from the perturbation the planet makes on the binary, correlated ETV noise variations can drive the model to prefer a negative mass planet. To mitigate this problem, we chose to fit only three primary eclipse profiles that showed very clean residuals. These Short Cadence data are shown in Figure 5. With the superb Kepler photometry, these three primary eclipses are sufficient to accurately derive the geometric binary-system parameters. Because the secondary eclipses do not show any evidence of starspot contamination, we include them all in the fitting. Since only three of the 49 primary eclipse profiles are fit, to capture the effect of the planet on the binary (the only way to measure the planet's mass in this single-planet system), we need to include as part of the data all of the remaining primary mid-eclipse times. The eclipse times need to be corrected for the starspot-induced variations, and this is done by de-correlating the times with the local light curve slope (see Mazeh et al. 2015 see also Holczer et al. 2015). The model is then penalized in a ${\chi }^{2}$ sense if its predicted eclipse times do not match the observed eclipse times. (The same 1.44 boost factor on the Kepler SAP error bars was maintained, for both the Long and Short Cadence.) As with the radial velocity observations, the eclipse times were given equal weight per point as the photometric data, i.e., no regularization of the different data sets was enforced. With this new strategy for modeling the data, the dynamical constraints in the data are preserved, even though not every individual eclipse profile is fit. The key advantage of this method is that it uses primary eclipse times that have been corrected for the effect of starspots. The individual eclipse profiles cannot be corrected for the effect of starspots without employing a prohibitively computationally expensive evolving-starspot model, requiring at least a dozen additional free parameters for even a rather simple two-starspot characterization (i.e., starspot latitudes, longitudes, radii, temperatures, onset times, lifetimes, etc.) Using this strategy, our photodynamical model now prefers a much larger and realistic mass for the planet. While the mass is still too small to measure, the uncertainty is ±16 ${M}_{\oplus }$, a much more plausible result.

Using the Short Cadence data when possible, there are 18,472 Kepler photometric data points, plus 11 radial velocities and 96 eclipse times. As evident in Figures 1 and 4–6, the ELC photodynamical model provides a good fit to the eclipse profiles, times of eclipse, radial velocities, and transits. Formally, the model fit yields a ${\chi }^{2}$ of 18725.9, or a reduced ${\chi }^{2}$ of 1.0093. In particular, the model matches the transit timing variations (TTVs) and the transit duration variations (TDVs). The time interval between the first two transits is 237.35 days, while the interval between the second and third transits is 235.56 days. The variation in transit duration is even more extreme, changing from 5.7 to 12.5 hr (see Figure 6). Table 5 lists the observed transit times, durations, and depths, and the predicted values from the photodynamical model. The large TTVs and the more than a factor of two difference in the duration of the transit are due to the "moving target" effect—the motion of the primary star about the barycenter. Thus the TTVs and TDVs follow well-defined curves as functions of the orbital phase of the binary and constitute an unambiguous signal of a circumbinary object, as no known false positive can reproduce this effect. This "smoking gun" signature will remain approximately true even in the presence of additional planets. However, if the transiting planet's orbit is significantly tilted with respect to the binary's orbit, the precession timescale may be sufficiently short that the changing impact parameter will alter the transit durations. In some cases, such as in Kepler-413, and now also in Kepler-453, the orbit precesses such that a large fraction of the time the planet does not transit the primary. We return to this topic in Section 4.3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Finally, the parameters from the photodynamical solution were input into the photodynam code of Josh Carter that was used to model the previous Kepler circumbinary planets (Carter et al. 2011; Pál 2012) and yielded good matches to the light curve. In addition, another completely independent photodynamical code developed by one of us (SMM), has been used to test and confirm the validity of the ELC code.

The secondary star has a mass ${M}_{2}=0.1951\pm 0.0020\ {M}_{\odot }$ and radius ${R}_{2}=0.2150\pm 0.0014\ {R}_{\odot }$, making it one of the lowest-mass stars with a precise dynamical mass and radius determination. It is just slightly less massive than Kepler-16 B (0.20255 ${M}_{\odot }$; Doyle et al. 2011; see also Winn et al. 2011 and Bender et al. 2012), KOI-126 C (0.2127 ${M}_{\odot }$; Carter et al. 2011), which in turn is very slightly smaller than CM Dra B,A (0.2141 and 0.2310 ${M}_{\odot }$; Morales et al. 2009). The ELC model provides a temperature ratio of ${T}_{2}{/T}_{1}$= 0.584, which, combined with the spectroscopically measured ${T}_{\mathrm{eff}}$ of 5527 K for the primary, gives a ${T}_{\mathrm{eff}}$ of 3226 K for the secondary star.

Dartmouth stellar evolution models (Dotter et al. 2008) were compared with the measured mass, radius, temperature, and metallicity of the stars ([Fe/H] = 0.09, and assuming no alpha-element enhancement). The results are shown in Figure 7 in the form of a radius versus mass diagram, with five isochrones spanning ages of 1–5 Gyr. For illustrative purposes, all currently known low-mass stars with masses and radii measured to better than 7% are shown in the figure. The best-match isochrone has an age of 1 Gyr, and was found by minimizing the distance between the observed values and the isochrone curves. "Distance" is defined as the quadrature sum of the uncertainty-normalized differences between points on the interpolated isochrone curve and the observed masses, radii, and temperatures of the stars. The distance is similar to the $\sqrt{{\chi }^{2}}$, but with the important distinction that this is not a fit, it is simply being used to compare isochrones.

Zoom In Zoom Out Reset image size

The sub-panels of Figure 7 show close-up views of the primary and secondary star, and for comparison an additional isochrone with [Fe/H] = 0.0 is shown as the dotted curve. The secondary star has a very mild preference for an age of over 5 Gyr but is consistent with a very large range of ages, while the primary star strongly prefers a young age, ∼1 Gyr or less. There is no formal discrepancy, but the tension between the two ages might be a manifestation of the common problem of low-mass stars being too large and too cool compared to stellar models (see, e.g., Torres 2013). However, we caution that for masses near the secondary star's mass, the isochrone models' points are relatively sparse and require large interpolation between the computed values, leading to possible systematic errors of ∼0.0015 of a solar radius. This is comparable to the size of the observational uncertainties, and so no conclusions beyond a general agreement should be drawn from the match between the secondary star and the model isochrones. Additional exploration of isochrones for a spread of metallicities showed that ages up to ∼1.75 Gyr are consistent with the mass, radius, and temperature of the primary star. We also compared Yonsei-Yale isochrones (Yi et al. 2001) with the primary star's observed mass, radius, and temperature, and obtained similar results.

The modulation in the light curve caused by starspots allows the rotation period of the primary star to be measured. Using the autocorrelation function (ACF) method of McQuillan et al. (2014), the rotation period is 20.31 ± 0.47 day. We also visually inspected and computed the ACF for individual Quarters in order to be certain that the true period was not half or twice this value. In comparison, the orbital period is 27.32 day, and the pseudosynchronous rotation period (Hut 1981, 1982) for the mild eccentricity of the orbit e = 0.051, is 26.90 day. Thus the star's rotation is not close to being synchronized with its orbit, though this is not unexpected given the relatively long orbital period. In general, the spin should synchronize much sooner than the orbit should circularize (at least in a 2-body system), so the system is likely to have been born with a low eccentricity. The predicted ${V}_{\mathrm{rot}}\mathrm{sin}i$ using the measured spin period and stellar radius, and assuming the spin inclination is aligned with the binary orbit inclination, is 2.08 $\mathrm{km}\ {{\rm{s}}}^{-1}$. This is in agreement with the observed low projected rotation velocity of the star 1.9 ± 1.0 $\mathrm{km}\ {{\rm{s}}}^{-1}$.

Walkowicz & Basri (2013) report an age of 2.31 Gyr for Kepler-453, based on their measured rotation period of 19.98 days which was determined via a Lomb–Scargle periodogram. Using our ACF-derived rotation period and the $B-V$ color index (corrected for reddening), we used five different gyrochronology relations and found age estimates spanning 2.11 to 2.87 Gyr (Barnes 2007, 2010; Mamajek & Hillenbrand 2008; Meibom et al. 2009 and Epstein & Pinsonneault 2014). These ages are older than the isochrone age, but they should be regarded more as upper limits because the star may have been influenced by tidal interactions in the binary, which are very slowly driving the spin period to match the binary orbital period of 27.32 day.

Finally, the observed ∼0.22% rms of the fluctuations in the light curve is slightly smaller than the solar value (∼0.3%), suggesting an older, rather than younger, age. The Ca ii H & K lines were not within the wavelength range of the HET spectra, and were too faint in the Tull spectra, so they could not be used as an activity or age indicator. The lack of any significant emission cores and the lack of the lithium 6708 Å absorption line rules out a very young age.

The ability to precisely measure the stellar masses and radii, especially for the low mass secondary star, and the mass ratio of ${M}_{2}{/M}_{1}=0.207\pm 0.003$, make the binary worthy of a more thorough investigation in its own right.

While the radius of Kepler-453 b is well determined at $6.20\pm 0.04\ {R}_{\oplus }$, we cannot reliably estimate its mass since it has no measurable effect on the binary during the course of our observations. If the planet were sufficiently massive, it would cause periodic changes in the binary eclipse times (ETVs) on roughly half the orbital period of the planet. And if the planet's orbit were not coplanar, it would cause a precession of the binary orbit (though measurable only if the binary orbit had nonzero eccentricity). This precession would be most readily detectable via a difference in the orbital periods of the primary and secondary stars, i.e., a divergence in the common-period $O-C$ diagram. Finally, a periodic change in the orbital inclination of the binary could lead to slight changes in the primary eclipse depth. As none of these are seen in the Kepler photometry, we are only able to estimate an upper limit on the mass. As mentioned earlier, the photodynamical model yields a 1σ upper value of ∼16 ${M}_{\oplus }$. We can also attempt to constrain the planet mass by matching the observed ETVs, both in their point-to-point fluctuations and in the long-term divergence between primary and secondary periods. However, there is little correlation between the observed primary eclipse-timings and the expected ETV signal that would be induced by a planet of any mass. Nor is there any significant difference in the orbital periods (mainly because the secondary star eclipse times are too noisy—see Figure 3). The observed 1σ upper limit is 1.7 s, which leads to an upper mass estimate of $\sim 30\ {M}_{\oplus }$, but with a large uncertainty of ∼±100 ${M}_{\oplus }$. We therefore adopt the photodynamical modeling result.

We can compare the above mass estimates/limits to the empirically expected mass, based on the radius of the planet. Several mass–radius relations have been published: Lissauer et al. (2011), Enoch et al. (2012), Kane & Gelino (2012), Weiss et al. (2013), Wu & Lithwick (2013), and Weiss & Marcy (2014). Assuming that these relations are applicable to circumbinary planets, the mass and radius pairs that fall within the range of validity of these relations span ∼17–43 ${M}_{\oplus }$. These masses are larger than what the photodynamical ELC model analysis prefers, though not inconsistent with the upper limits. Given the systematic uncertainties in the empirical relations, the agreement is reasonably good.

The photodynamical model predicts the planet's (sky-projected) inclination oscillates with a 103.4 year period, which is also the precession period of the planet's orbit. For the best-fit solution, the planet's inclination varies by 45. Conservation of angular momentum requires the binary's orbital plane to also oscillate, but with a miniscule $\sim 2\times {10}^{-5}$ deg amplitude—which is not surprising as the best-fit planet mass is near zero. Figure 8 illustrates the geometry of the binary and a spatial view of the rapid precession of the planet's orbit.

Zoom In Zoom Out Reset image size

As a consequence of this changing inclination of the planet's orbit, the impact parameter of the planet will vary by a large amount; over a fairly short timescale, the planet precesses from non-transiting to transiting, then will become non-transiting again. In these Kepler data the planet transits 50% of the time: three non-transits followed by three transits. However, this was very fortuitous—due to the relatively large mutual inclination (225) between the orbital planes, the majority of the time the planet does not transit the star at conjunction. Transits only occur when the impact parameter, defined as the minimum projected distance of the planet's limb from the center of the primary star during conjunction, relative to the radius of the primary, is less than unity. This occurs only 8.9% ± 0.2% of the time. In Figure 9 we show explicitly the sinusoidal oscillation of the planet's orbital inclination, along with the mutual inclination, and the evolution of the planet's impact parameter across the primary star. An interesting and important implication is that because of observational duration limitations, for every transiting circumbinary system like Kepler-453 detected, there are ∼11.5 similar systems that are not currently transiting. The importance of binary-induced planetary orbital precession in determining circumbinary transit probabilities was foreseen by the prescient Schneider (1994).

Zoom In Zoom Out Reset image size

In addition to the observed transit times, durations, and depths, Table 5 includes predicted near-future transit times. The near-future timescale is by necessity—we predict the transits will stop being visible after 2015 July and not return until 2066.

Following Dvorak (1986), Dvorak et al. (1989), and Holman & Wiegert (1999), we calculated the orbital period for which the planet is highly susceptible to a dynamical instability. With its orbital period 8.80 times that of the binary period, Kepler-453 b is comfortably above the critical period for circumbinary instability: ${P}_{\mathrm{planet}}{/P}_{\mathrm{crit}}\;=$ 2.40. Unlike previous Kepler circumbinary planets, Kepler-453 b is not skirting on the edge of the instability. Note that we have found an error in previously reported values for ${P}_{\mathrm{planet}}/{P}_{\mathrm{crit}}$ for Kepler 34, 35 (Welsh et al. 2012), and 38 (Orosz et al. 2012b). The reported values, 1.21, 1.24, and 1.42 should be 1.49, 1.34, and 1.30.

The above analysis suggests that the orbit of Kepler-453 b is dynamically stable. However, it is important to note that this critical stability limit was derived assuming the planet has negligible mass compared to the binary, and is initially on a circular and coplanar orbit. The perturbing effect of the planet on the dynamics of the binary is ignored. So to truly examine the long-term stability of the orbit of the planet, one has to integrate the equations of motion of the three-body system.

We first investigated the stability of the planet by applying the Mean Exponential Growth factor of Nearby Orbits (MEGNO) criterion (see Cincotta & Simo 1999, 2000; Goździewski et al. 2001; Goździewski & Maciejewski 2001 and Hinse et al. 2010). MEGNO is used to determine general regions of dynamical instability, chaotic zones, and the locations of orbital resonances. Figure 10 shows the result of computing the MEGNO²⁴ factor over a grid of the test planet's initial semimajor axis and eccentricity. All other orbital parameters are fixed and taken as the best-fit ELC values in Table 3. Each of the 240,000 pixels in the map was integrated for ∼18,500 binary periods (500,000 days). The color coding corresponds to the degree to which the orbit of the test planet is chaotic (yellow) or quasi-periodic (blue). The location of the planet's best-fit osculating Keplerian semimajor axis and eccentricity (from Table 3) is marked with the black dot. The MEGNO map shows that the planet orbit resides in a region of phase space that is not chaotic for at least the ∼1370 year duration of the MEGNO integration. Note that the yellow instability gaps dipping into the blue region of the map correspond to n:1 mean-motion resonances between the planet and the binary orbit. These resonances play an important role in establishing the final orbit of the planet during its migration from outer regions of the disk where it is presumably formed (Kley & Haghighipour 2014).

Zoom In Zoom Out Reset image size

To better determine the long-term stability of the planet, we integrated the orbits of the three-body system for 10 Myrs. The results are shown in Figure 11. As expected, the variations in the semimajor axis and eccentricity of the planet are negligibly small over the course of the integration. The long-term integrations also show that the mutual inclination between the planet and the binary is not caused by the gravitational interactions, i.e., it is a "free" inclination that is presumably primordial (tides are unable to significantly modify this inclination). Similarly, the eccentricity of the binary is also likely primordial. The planet's eccentricity is dominated by the forcing of the binary, but there is a free eccentricity component as well (see Leung & Lee 2013).

Zoom In Zoom Out Reset image size

Kepler-453 exhibits a dynamical interaction between the binary and the planet with a period of 103.4 years. The 10 Myr-long integration suggests that this interaction is unlikely to lead to an instability. The semimajor axes of both the binary and planetary orbits remain relatively constant for the entire duration of the simulation. Although the eccentricity of the planet varies on this timescale, the planetary periastron does not become comparable to the apastron of the binary orbit. The mutual inclination (i_rel) between the orbital planes exhibits a small variation of ∼0035; however precession of ${\omega }_{\mathrm{binary}}-{\omega }_{\mathrm{planet}}$ also influences the impact parameter which allows for only a limited window where transits can be observed. The ascending node of the planet deviates ∼2° from that of the binary which in turn limits the variation of the inclination vector. The lack of resonant phenomena or secular growth in the orbital elements implies that Kepler-453 lies within a stable 3-body configuration.

Unlike the HZ around a single star (e.g., see Kopparapu et al. 2013a), the HZ around a binary star is neither spherical nor fixed—it revolves with the binary. Nevertheless, time-averaged approximations are useful for a preliminary investigation of the orbit of the planet with respect to its location in the HZ. In this approximation, we use the time-averaged distance of bodies in elliptical Keplerian orbits (see Williams 2003). Using ${T}_{\mathrm{eff}}$ of 5527 K for the primary star and 3226 K for the secondary star, the secondary emits only 0.77% of the bolometric luminosity and ∼0.26% of the light in the Kepler bandpass. Thus we initially ignored the secondary star's flux contribution. Assuming a Bond albedo of 0.34, appropriate for both the gas giant planets in our Solar System and for the Earth itself, and assuming re-emission over a full sphere, we find the planet's equilibrium temperature ${T}_{\mathrm{eq}}$ to be 247 K and the time-averaged insolation S = 0.94 in Sun–Earth flux units. This is within the conservative HZ limits defined by Kopparapu et al. (2013a, 2013b), set by the runaway and maximum greenhouse criteria for the inner and outer boundaries. Even at the most extreme ranges of star–planet separation due to the eccentricity, ${T}_{\mathrm{eq}}$ only varies between 237 and 258 K.

We also calculated the properties of the HZ using the more detailed methodology presented by Haghighipour & Kaltenegger (2013) and the Multiple Star HZ calculator²⁵ (Müller & Haghighipour 2014). The results are presented in Figure 12, where the dark green region corresponds to the narrow (conservative) HZ and the light green corresponds to the nominal (extended) HZ as defined by Kopparapu et al. (2013a, 2013b, with coefficients updated on-line in 2014). The dashed circle represents the dynamical stability limit.

Zoom In Zoom Out Reset image size

We then computed the instantaneous insolation using our photodynamical model, which gave a mean insolation over one full precession cycle of $\langle S\rangle$ = 0.953. Figure 13 shows this varying insolation incident upon the planet from both stars as a function of time. The most obvious variation is due to the orbital motion of the primary star. Superposed on this is the variation due to the eccentric orbit of the planet. Finally, on a ∼103 year timescale, the effects of the oscillation of the mutual inclination are apparent. The overall fluctuation in insolation is 8.0% (rms) about the mean. Interestingly, from the perspective of the planet, not every stellar conjunction is seen as an eclipse. This is due to the mutual inclination of the orbits. Only when the planet crosses the plane of the binary (i.e., near the nodal points) will an eclipse be seen from the planet. Since the ratio of orbital periods is 8.8025, and there are two nodal crossings per orbit, the eclipses are visible every ∼4.4 binary periods, but this varies slightly with the oscillation of the mutual inclination.

Zoom In Zoom Out Reset image size

Although Kepler-453 b is a 6-${R}_{\oplus }$ planet and is unlikely to harbor conditions suitable for life, e.g., its density is <20% of the Earth's (at roughly the 3σ level), such a planet could host a large moon capable of sustaining life. Since no evidence for such a moon is apparent we do not pursue this further, except to mention that in general, the effect noted by Mason et al. (2013) should be considered: in a binary, tides can significantly change the spin evolution—and thus the activity—of the stars, rendering the system more, or less, habitable depending on the specifics of the binary.

Other than the outer planets of the Kepler-47 system, all transiting Kepler circumbinary planets fall within a factor of two of the critical radius for instability. As remarked upon in Orosz et al. (2012b), whether this is a selection effect or not, the observed close-to-critical orbits have an interesting consequence: the Kepler circumbinary planets tend to lie close to their HZ. Both Kepler-16 b and Kepler-47 c are in the HZ, and Kepler-453 b now joins this group of HZ circumbinary planets. Currently 3 out of 10 Kepler circumbinary systems reside in the HZ. Thus, although difficult to detect because of their large variations in transit timing and duration, the search for smaller Earth-size planets in circumbinary environments is a particularly exciting endeavor.