Discovery of a Third Transiting Planet in the Kepler-47 Circumbinary System

Kepler-47 (KIC 10020423, KOI 3154, 2MASS J19411149+4655136) was observed in long-cadence mode (≈30 minute samples) from Kepler Quarter 1 (BJD 2,454,964.51 or 2009 May 13) through the end of the nominal spacecraft operation in the second month of Quarter 17 (BJD 2,456,424.00 or 2013 May 11). There were a total of 65,428 long-cadence observations during this span (including flagged data), for a duty cycle of 91.6%. Starting in Quarter 14 (BJD 2,456,107.13 or 2012 June 28), Kepler-47 was on the list of short-cadence targets (≈1 minute samples), where it remained until the cessation of normal spacecraft operations. There were a total of 409,046 short-cadence observations, including flagged data. We downloaded the FITS tables from the Mikulski Archive for Space Telescopes (MAST). These files were processed using the SOC release 9.0.3, and they include a ≈1 minute correction to the time stamps that was not available at the time of original work reported in Orosz et al. (2012b).

The Kepler light curves of Kepler-47 required detrending, owing to instrumental trends and modulations due to star spots. We arrived at our final detrended light curves using a two-step process. First, an interactive technique described in Orosz et al. (2012b) and Bass et al. (2012) was used. In this step, the light curve was broken up into segments using data gaps and sudden jumps in the flux as end points. For each segment, data in the eclipses and transits were masked out, and cubic splines were fit to the remaining data. The segments were normalized by the spline fits, and pieced together to produce a detrended light curve. While effective, this technique is not entirely reproducible; it is also difficult to use with the short-cadence data, owing to limitations in our implementation of the method. Thus, an automated and reproducible detrending algorithm was devised as follows: The light curve that was detrended manually was modeled as described in Section 5. Once a good fit was found, the model was used to precisely determine the times of the eclipses and transits and the durations of these events. The times and durations were then used to determine which data points occurred during an eclipse or a transit event, which allowed us to assign those points zero weight during the normalization.

In our original work, only long-cadence Kepler data were available to us. As noted above, Kepler-47 was on the list of short-cadence targets for Kepler Quarters Q14 through Q17. In our view, using short-cadence data (when available) instead of long-cadence data is preferable because the ingress and egress phases of eclipses and transits are better-resolved. Thus, a "raw" light curve using long-cadence data from Quarter 1 through Quarter 13 and short-cadence data from Quarter 14 through Quarter 17 was made, and data with quality flags with values greater than 16 were eliminated (in our experience, obviously poor cadences almost always have data quality flags greater than 16). The overall duty cycle of the retained data was 91.6%. Segments centered at the mid-eclipse and transit times with widths of three times the duration of the event in question were extracted from the raw light curve. In a few cases, these segments included both an eclipse event and a transit event; also in a few cases, the segments were truncated owing to gaps in the data. A fifth-order Legendre polynomial was fit to each segment, where data that occurred during an eclipse or a transit event were given zero weight. The normalized segments were assembled to give the final normalized and trimmed light curve. We also produced a normalized light curve consisting entirely of long-cadence observations for the purposes of plotting the relatively weak transits.

Owing to the relative faintness of the secondary star (it is about 0.5% as bright as the primary in the optical), Kepler-47 is a single-lined binary, so only radial velocity measurements for the primary star are available. We used 11 radial velocity measurements of the primary in the original work (Orosz et al. 2012b). Four of these came from the HRS spectrograph on the Hobby–Eberly Telescope (HET), six of these came from the Tull Coudé spectrograph on the Harlan J. Smith 2.7 m telescope (HJST) at McDonald Observatory, and one measurement came from the HIRES spectrograph on the Keck I telescope. Since that time, Kostov et al. (2013) have provided an additional four measurements with the SOPHIE spectrograph on the 1.93 m telescope at Haute-Provence Observatory. These measurements are summarized in Table 1 and displayed in Figure 1. Small differences were found in the velocity zero-points between the different observatories, so these were fitted for and removed. To arrive at our final radial velocity curve, we subtracted systemic offsets of 4.5680 km s⁻¹ for the HJST radial velocities, 4.5999 km s⁻¹ for the HET radial velocities, and 4.2685 km s⁻¹ for the Kostov et al. (2013) velocities. Given that we have only one measurement from Keck, we cannot determine the systematic offset and we do not use that measurement in our analysis. The uncertainties in the individual measurements in each set were scaled to give χ² = N, where N = 14 is the number of radial velocity measurements. The offset velocities with the scaled uncertainties are also given in Table 1.

Zoom In Zoom Out Reset image size

The numerical integrator that ELC uses is a symplectic, 12th-order Gaussian Runge–Kutta (GRK) routine based on methods and codes devised by Hairer et al. (2006). For situations where only the mutual gravitational forces between the five bodies are considered, the ordinary differential equations (ODEs) are second-order and can be written such that only the positions of the bodies appear explicitly. GRK is a collocation method in which a system of equations based on the Gaussian quadrature nodes must be solved at each time step. In particular, for the 12th-order method, there are six nodes that determine the coefficients of the collocation polynomial. These coefficients are solved for iteratively. With a suitable initial guess, the convergence of this iteration is quadratic—and hence, very fast. If needed, the equations of motion can be modified to account for extra effects such as precession due to General Relativity (GR) and/or precession due to tidal bulges on the stars (Eggleton et al. 1998; Mardling & Lin 2002). In the case of the GR correction, the velocities of the bodies appear explicitly in the equations of motion, and the first-order ODEs must be solved using a different iteration scheme. This iteration scheme is only to first order, so the convergence is somewhat slower than the second-order scheme.

Many of the symplectic integrators used in the field of solar system dynamics (e.g., Wisdom & Holman 1991) rely on a factorization (or splitting) of the Hamiltonian in a Lie algebra sense. These integrators are very fast, but are best applied to problems where there is one dominant central mass. One advantage of GRK integrators is that they do not approach the issue as a Lie algebra factorization and hence do not have the restriction of just one large mass. The trade-off is that the step size is determined by the shortest period, which results in longer run-times (we typically use a time step h that is about 400 times smaller than the smallest period in the system). However, we note that, even with the relatively low speed compared to other integrators, it takes much less than one second to solve the coupled equations of motion for the five bodies in Kepler-47 over a 1500 day time span using the GRK method.

The GRK scheme does have another advantage, in that it is possible to find positions and velocities at intermediate times to high accuracy with relatively few computations. In our implementation, the solutions at the six internal nodal points used for the collocation for each time step are saved. Using the collocation points with divided differences, it is easy to determine the positions and velocities of each body at intermediate times to the same order of accuracy (namely, 12th) as at the integration times (e.g., t_start, t_start + h, t_start + 2h, etc.). The number of operations needed to achieve the 12th-order accuracy is comparable to what one would need for a spline interpolation scheme, which would give the intermediate values with a much lower-order accuracy.

The integrator works in Cartesian coordinates relative to the system's barycenter. The unit of mass is the solar mass, the unit of distance is the astronomical unit, the unit of time is the day, and the unit of velocity is astronomical units per day. The adopted value of the Gaussian gravitational constant is k = 0.01720209895 (Clemence 1965). Newton's gravitational constant, expressed in ${(\mathrm{au})}^{2}{{M}_{\odot }}^{-1}{(\mathrm{day})}^{-2}$, is then G = k². For convenience, instantaneous Keplerian parameters valid at some reference epoch for each body are used to specify the initial conditions (our adopted reference time is T_ref = BJD 2,454,965.000). These parameters are the orbital period P, the time of barycentric transit (which is the same as the time of inferior conjunction) T_conj, the inclination i, the eccentricity e, the argument of periastron ω, and the nodal angle Ω. The coordinate system is Jacobian, so the orbital parameters of the secondary star are given relative to the primary, the orbital parameters of the first planet are given relative to the binary's center of mass, and so on. Because of this, the specified conjunction time for the binary occurs very near a primary eclipse, whereas nothing observable has to happen at the conjunction times of the planets (the actual transits can occur either earlier or later than the barycentric conjunction times). The Keplerian parameters are converted to Cartesian coordinates using the algorithms given in Murray & Dermott (1999).

After the ODE is solved, the plane-of-sky positions of each body are corrected for light travel time following the method outlined in Carter et al. (2011). The times of the primary and secondary eclipses, and the times of transits of each planet across the primary and secondary are then computed by finding the times when the projected separation between two given bodies in the plane of the sky is minimized. For computational convenience, we minimize the square of the plane-of-sky distance vector, which is just the dot product of that distance vector with itself. That dot product is minimized when its derivative is zero. That derivative is two times the dot product of the distance vector with its velocity vector (which itself is readily found using the appropriate components of the velocities that are output from the ODE solver), and is zero when the distance vector itself is zero or when the distance vector is perpendicular to its velocity vector. Several iterations of the secant method are used to find the zeros of the derivative. By comparing the results found using the first-order ODE solver with the results found using the second-order ODE solver, we estimate the times of conjunction to be accurate to ≈1 μs or better. An eclipse or a transit will occur at a conjunction if the projected minimum separation on the sky is less than the sum of the radii of the two bodies.

For the specific case of Kepler-47, the corrections due to GR result in an apsidal advance of about 000017 per orbital period. The rate of apsidal advance due to tidal bulges on each star is between 10 and 50 times smaller, depending on the values of the tidal Love numbers used (we considered 0.0 ≤ k₂ ≤ 0.01 for the primary and 0.0 ≤ k₂ ≤ 0.2 for the secondary). In what follows, we have included the GR correction in the ODEs and neglected the corrections for the apsidal advance due to tidal bulges on the stars.

For situations where the bodies are spherical with linear or quadratic limb darkening laws, ELC can use the algorithm of Mandel & Agol (2002) or the algorithm of Giménez (2006a) to compute the light curves during eclipses or transits. In the present work, we have used the Mandel & Agol (2002) routine with a quadratic limb darkening law, because it is much faster than the Giménez (2006a) routine.

Our model as applied to Kepler-47 has 42 free parameters in total. Each orbit needs six Keplerian parameters for a total of 24. However, the nodal angle of the binary is fixed at zero, so there are really only 23 free parameters to describe the orbits. For each orbit, we used the combinations $e\cos \omega$ and $e\sin \omega$ rather than e and ω separately, because the former pair of variables usually are less correlated with each other compared to the latter pair of variables. Because there are five bodies, five parameters describing the masses of the bodies and five parameters describing the sizes of the bodies are needed. For the binary, we used the primary mass M₁ and the binary mass ratio Q_2,1 ≡ M₂/M₁. For the two stars, we used the fractional radii R₁/a and R₂/a to parameterize their sizes. To parameterize the planetary radii, we used the ratio R₁/R_p, which is the ratio of the primary star's radius to the particular planet radius in question. There are two radiating bodies in the system, so a total of six parameters are required to specify the radiative properties: the two temperatures and two limb darkening coefficients for each star. In this case, the temperature of the primary was fixed at 5636 K (Orosz et al. 2012b), so only five free parameters were used. For convenience, the temperature ratio T₂/T₁ was used instead of directly using T₂. ELC has a table of model atmosphere specific intensities (for solar metallicity) derived from the NextGen models (Allard et al. 1997; Hauschildt et al. 1999), and these were used to compute the baseline fluxes of the primary and secondary star (in the Kepler bandpass) given their temperatures and gravities. The standard quadratic limb darkening law given by

$$\begin{eqnarray}&&I(\mu )/{I}_{0}=1-{u}_{1}(1-\mu )-{u}_{2}{\left(1-\mu \right)}^{2}\end{eqnarray} \tag{ 2 }$$

was used (where $\mu =\cos \theta$ is the projected distance from the center of the stellar disk), but with the "triangular" sampling technique of Kipping (2013) with coefficients given by ${q}_{1}={({u}_{1}+{u}_{2})}^{2}$ and ${q}_{2}=0.5{u}_{1}{({u}_{1}+{u}_{2})}^{-1}$. Finally, to account for light from other sources in Kepler's aperture, four seasonal contamination parameters were used.

As was shown in Section 4, many of the primary eclipse profiles show clear evidence of spot crossing events. Because our light curve model does not include star spots, fitting the primary eclipses where there are clear spot crossing events might lead to biased results. To avoid this possible bias in the fitting, we fit the statistically corrected primary eclipse times along with only a small subset of the light curve. We chose ten primary eclipse profiles (seven in long cadence and three in short cadence) that do not have spot crossing events to fit, as well as ten secondary eclipses that are close in time to the ten selected primary eclipses. We also fit the 25 transits of the inner planet, the six transits of the middle planet, and the four transits of the outer planet. Finally, we also included the radial velocities of the primary star in the fit.

We used four algorithms to explore parameter space: an adaptation of a simple grid-search algorithm (Bevington 1969), a genetic algorithm (Charbonneau 1995); a simple search based the MCMC algorithm outlined by Tegmark et al. (2004), and a Differential Evolution MCMC (DE-MCMC) algorithm of Ter Braak (2006; see also Nelson et al. 2014). All four of these algorithms are bounded, meaning the fitting parameters never wander outside predefined ranges. In addition, the grid-search code, the genetic code, and the DE-MCMC code can be run in parallel on multiple CPU cores. We assumed the likelihood function follows a χ² distribution, where χ² is based on the sum of squares of differences between observed quantities (a photometric measurement, a radial velocity measurement, or an eclipse or transit time) and predicted quantities normalized by the measurement uncertainty for each quantity. Each of these algorithms needs to have fixed ranges for each free parameter. Some parameters, like some of the orbital parameters for the binary, are constrained reasonably well. We also had some good starting solutions for the parameters of the inner planet from our previous work (Orosz et al. 2012b), whereas some of the orbital parameters (for example, the eccentricity parameters) for the outer planet were constrained much less well. Apart from the rough value of the orbital period and the fact that the inclination of the orbit must be near 90° in order for transits to occur, the orbital parameters for the middle planet were initially not known.

Because the parameter space is vast, we looked for optimal solutions in three main stages: (i) updating the model from the previous work by adding the middle planet to the model and finding some initial fits; (ii) exploring parameter space to find something close to the "best" model, get some rough idea of the uncertainties in the fitted and derived parameters, and also look for possible correlations between various parameters; and (iii) finding realistic uncertainties on the fitted and derived parameters. The three stages are described below.

Stage i: the inner planet has seven additional transits and the outer planet has one additional transit since the previous work, so the orbital parameters for those two planets were revised first using the DE-MCMC code. Next, the middle planet was added to the model, and initial orbital parameters were found "by hand." Educated guesses for the time of barycentric conjunction, the inclination, the nodal angle, etc. were used to initialize the simple MCMC code, and the model fits were inspected by eye after running a few iterations of that optimizer. It was not unduly difficult to find a model that produced transits of the middle planet at roughly the correct times.

Stage ii: after a few decent models were found, longer runs using the genetic algorithm and the DE-MCMC were done. Both of these codes require an initial "population" of models (we typically used population sizes between 100 and 200), and we have various ways of generating the initial population: (a) totally random values between the specified lower and upper parameter bounds for each free parameter; (b) one or more "elite" models (e.g., ones that are already a good match to the data) along with randomly generated models; and (c) one or more "elite" models along with "mutated" copies of the elite models where a few randomly chosen parameters are "tweaked" by adding or subtracting small offsets. In practice, options (b) and (c) work the best for cases where one has a few dozen or more free parameters.

For this intermediate stage of fitting, we adopted ranges for the free parameters given in Table 4. We ran the genetic algorithm or the DE-MCMC code for a few pilot runs of several thousand generations to confirm that the prior ranges included support for the entire range with nontrivial likelihood. After each pilot run, plots of parameter values versus the generation number and plots of the χ² versus each parameter were inspected to find instances where the allowed range of a parameter was either too large or too small. We verified that the likelihood falls to extremely small values by the time model parameters reach any of these boundaries (with the possible exception of hard physical boundaries).

Table 4.  Priors for Fitting Parameters

Orbit	Parameter	Lower	Upper
Binary and Stellar	Period (days)	7.4480	7.44899
	Tconj (BJD 2,455,000+)	−29.316	−29.300
	$e\cos \omega$	−0.1	0.0
	$e\sin \omega$	−0.05	0.01
	Inclination (deg)	89.00	90.25
	M1 (M⊙)	0.80	1.0
	Q	0.340	0.373
	R1/a	0.046	0.061
	R2/a	0.0186	0.0202
	T2/T1	0.39	0.83
	q1 (primary)	0.0	1.0
	q2 (primary)	0.0	1.0
	q1 (secondary)	0.0	1.0
	q2 (secondary)	0.0	1.0
Inner	Period (days)	49.41	49.52
	Tconj (BJD 2,455,000+)	−31.6	−31.14
	$e\cos \omega$	−0.08	0.10
	$e\sin \omega$	10.114	0.069
	Inclination (deg)	89.03	90.32
	Nodal angle (deg)	−0.390	0.315
	M3 (M⊕)	0.01	100.0
Middle	Period (days)	187.15	187.60
	Tconj (BJD 2,455,000+)	45.10	45.75
	$e\cos \omega$	−0.39	0.35
	$e\sin \omega$	−0.1	0.2
	Inclination (deg)	90.18	90.62
	Nodal angle (deg)	−2.8	0.5
	M3 (M⊕)	0.01	100.0
Outer	Period (days)	303.00	303.42
	Tconj (BJD 2,455,000+)	−55.8	−55.4
	$e\cos \omega$	−0.19	0.23
	$e\sin \omega$	−0.21	0.11
	Inclination (deg)	90.10	90.28
	Nodal angle (deg)	−3.0	1.2
	M3 (M⊕)	0.01	100.0
Kepler	First season contaminationa	0.000	0.035
	Second season contaminationb	0.000	0.035
	Third season contaminationc	0.000	0.035
	Fourth season contaminationd	0.000	0.035

Notes.

^aKepler Quarters Q1, Q5, ⋯ ^bKepler Quarters Q2, Q6, ⋯ ^cKepler Quarters Q3, Q7, ⋯ ^dKepler Quarters Q4, Q8, ⋯

Download table as:  ASCIITypeset image

After this intermediate state of optimization, several million models have been computed. The genetic algorithm and the DE-MCMC code explore parameter space in very different ways; after several long runs of each, we can be reasonably sure the optimal region in parameter space has been found.

Stage iii: We used a brute-force "stepper" to generate initial "seed" models for the final run of the DE-MCMC code. The stepper works as follows. Take the best overall model and choose a key fitting parameter, such as the primary mass. Offset that parameter by a small amount (say, 0.5% of the parameter value) and hold it fixed while optimizing the other parameters using the simple parallel grid-search code, where the parameter values from the optimal model are used as the initial guess. After that optimization is done, offset the same key parameter again and repeat the optimization, using the previous solution as the initial guess. After the stepper is done, you have a set of reasonably optimal models with a range of different values of the chosen key parameter. We ran steppers on the mass parameters (primary mass, binary mass ratio, and the planet masses), the radius parameters (primary and secondary fractional radii and planet radius ratios), and the planet orbital periods.

For the last run of the DE-MCMC code, we used a population size of 1600 models per generation with 900 "seed" models found from the steppers. The other 700 models were mutated copies of randomly selected seed models with five parameters in each model randomly tweaked. The code was ran for 17,000 generations, giving about 27 million models in total. Judging from plots of the parameters versus the generation number, convergence was achieved after about 4000 generations. However, to be conservative, the first 5000 generations were discarded, leaving about 19 million models. Figure 1 shows the best-fitting model radial velocity curve. The best fits to the transits of the middle planet, inner planet, and the outer planet are shown in Figures 2–4, respectively. Figure 8 shows the best-fitting model with the "clean" primary eclipses and Figure 9 shows the best-fitting model with the corresponding secondary eclipses. The best-fitting input parameters are given in Table 5, and Table 6 summarizes several derived parameters of interest. In order to allow others to reproduce our best-fitting model using other codes, we give the initial Cartesian barycentric coordinates for the best-fitting model in Table 7 and the initial orbital elements (which are traditionally used in dynamical studies) for the best-fitting model in Table 8. These coordinates and orbital elements are valid for the reference time BJD 2,454,965.000.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5.  Fitted Parameters from the Photometric-dynamical Model^a

Parameter	Binary Orbit	Inner Orbit	Middle Orbit	Outer Orbit
Period (days)	${7.4483648}_{-0.0000270}^{+0.0000038}$	${49.4643}_{-0.0074}^{+0.0081}$	${187.366}_{0.051}^{+0.069}$	${303.227}_{-0.027}^{+0.062}$
Tconj (BJD + 2,455,000)	$-{29.305248}_{-0.000030}^{+0.000031}$	$-{31.3573}_{-0.0049}^{+0.0043}$	${45.428}_{-0.047}^{+0.044}$	$-{55.638}_{-0.046}^{+0.067}$
$e\cos \omega$	$-{0.019895}_{-0.000024}^{+0.000025}$	${0.0139}_{-0.0025}^{+0.0020}$	${0.024}_{-0.036}^{+0.038}$	${0.026}_{-0.034}^{+0.029}$
$e\sin \omega$	$-{0.0208}_{-0.0020}^{+0.0019}$	${0.0157}_{-0.0011}^{+0.0015}$	$-{0.004}_{-0.028}^{+0.023}$	$-{0.035}_{-0.025}^{+0.026}$
i (deg)	${89.613}_{-0.040}^{+0.045}$	${89.752}_{-0.045}^{+0.063}$	${90.395}_{-0.012}^{+0.009}$	${90.1925}_{-0.0042}^{+0.0055}$
Ω (deg)	0.0	$-{0.088}_{-0.053}^{+0.056}$	$-{0.87}_{-0.10}^{+0.11}$	$-{1.28}_{-0.25}^{+0.23}$
Q2,1 ≡ M2/M1	${0.3569}_{-0.0020}^{+0.0021}$	⋯	⋯	⋯
Mi (M⊕)	⋯	${2.07}_{-2.07}^{+23.70}$	${19.02}_{-11.67}^{+23.84}$	${3.17}_{-1.25}^{+2.18}$
Parameter	Primary	Secondary	Inner planet	Middle planet	Outer planet
R/ab	${0.05347}_{-0.00016}^{+0.00016}$	${0.019332}_{-0.000065}^{+0.000056}$	⋯	⋯	⋯
R1/Rc	⋯	⋯	${33.51}_{-0.48}^{+0.54}$	${14.55}_{-1.35}^{+1.11}$	${21.97}_{-0.35}^{+0.42}$
q1d	${0.502}_{-0.063}^{+0.067}$	<0.66	⋯	⋯	⋯
q2e	${0.188}_{-0.028}^{+0.020}$	${0.38}_{-0.12}^{+0.34}$	⋯	⋯	⋯
Other parameters
T2/T1	${0.613}_{-0.006}^{+0.012}$	⋯	⋯	⋯	⋯
s0 (Season 0 contamination)	<0.0069	⋯	⋯	⋯	⋯
s1 (Season 1 contamination)	${0.0219}_{-0.0037}^{+0.0045}$	⋯	⋯	⋯	⋯
s2 (Season 2 contamination)	${0.0044}_{-0.0038}^{+0.0043}$	⋯	⋯	⋯	⋯
s3 (Season 3 contamination)	${0.0232}_{-0.0029}^{+0.0047}$	⋯	⋯	⋯	⋯

Notes.

^aT_ref = BJD 2,454,965.000. ^bFractional stellar radii. ^cRatio of primary radius to planet radius. ^dKipping's first quad-law triangular limb darkening coefficient. ^eKipping's second quad-law triangular limb darkening coefficient.

Download table as:  ASCIITypeset image

Table 6.  Derived Parameters from the Photometric-dynamical Model

Parameter	Binary Orbit	Inner Orbit	Middle Orbit	Outer Orbit
a (au)	${0.08145}_{-0.00037}^{+0.00036}$	${0.2877}_{-0.0011}^{+0.0014}$	${0.6992}_{-0.0033}^{+0.0031}$	${0.9638}_{-0.0044}^{+0.0041}$
e	${0.0288}_{-0.0013}^{+0.0015}$	${0.0210}_{-0.0022}^{+0.0025}$	${0.024}_{-0.017}^{+0.025}$	${0.044}_{-0.019}^{+0.029}$
ω (deg)	${226.3}_{-2.6}^{+2.8}$	${48.6}_{-2.3}^{+3.0}$	${352}_{-166}^{+96}$	${306}_{-67}^{+40}$
I1,i (deg)a	⋯	${0.166}_{-0.041}^{+0.037}$	${1.165}_{-0.093}^{+0.101}$	${1.38}_{-0.24}^{+0.20}$
I2,i (deg)b	⋯	⋯	${1.006}_{-0.050}^{+0.089}$	${1.247}_{-0.018}^{+0.024}$
I3,i (deg)c	⋯	⋯	⋯	${0.442}_{-0.102}^{+0.028}$
Parameter	Primary	Secondary	Inner planet	Middle planet	Outer planet
M (M⊙ or M⊕)	${0.957}_{-0.015}^{+0.013}$	${0.342}_{-0.003}^{+0.003}$	${2.07}_{-2.07}^{+23.70}$	${19.02}_{-11.67}^{+23.84}$	${3.17}_{-1.25}^{+2.18}$
R (R⊙ or R⊕)	${0.936}_{-0.005}^{+0.005}$	${0.338}_{-0.002}^{+0.002}$	${3.05}_{-0.04}^{+0.04}$	${7.04}_{-0.49}^{+0.66}$	${4.65}_{-0.07}^{+0.09}$
ρ (g cm−3)	${1.65}_{-0.01}^{+0.02}$	${12.42}_{-0.11}^{+0.11}$	${0.40}_{-0.40}^{+4.87}$	${0.30}_{-0.18}^{+0.38}$	${0.17}_{-0.07}^{+0.09}$

Notes.

^aMutual inclination between the orbital planes of the binary and the planets, given by $\cos I=\sin {i}_{p}\sin {i}_{b}\cos {\rm{\Delta }}{\rm{\Omega }}+\cos {i}_{p}\cos {i}_{b}$. ^bMutual inclination between the orbital planes of the inner planet and the other planets. ^cMutual inclination between the orbital planes of the middle planet and the outer planet.

Download table as:  ASCIITypeset image

Table 9.  Parameter Ranges Displayed in Figures 18–20

Orbit	Parameter	Lower	Upper
Binary	Period (days)	7.44820	7.44838
	Tconj (BJD 2,455,000+)	−29.30541	−29.30515
	$e\cos \omega$	−0.02000	−0.01979
	$e\sin \omega$	−0.028	−0.012
	Inclination (deg)	89.5	89.8
	M1 (M⊙)	0.8950	0.9905
	M2 (M⊙)	0.326	0.352
Inner	Period (days)	49.440	49.495
	Tconj (BJD 2,455,000+)	−31.375	−31.337
	$e\cos \omega$	0.006	0.024
	$e\sin \omega$	0.0122	0.0210
	Inclination (deg)	89.55	90.05
	Nodal angle (deg)	−0.32	0.12
	M3 (M⊕)	0.0	100.0
Middle	Period (days)	187.2	187.5
	Tconj (BJD 2,455,000+)	45.25	45.60
	$e\cos \omega$	−0.15	0.14
	$e\sin \omega$	−0.08	0.09
	Inclination (deg)	90.34	90.45
	Nodal angle (deg)	−0.4	−0.4
	M3 (M⊕)	0.0	100.0
Outer	Period (days)	303.2	303.37
	Tconj (BJD 2,455,000+)	−55.78	−55.49
	$e\cos \omega$	−0.10	0.13
	$e\sin \omega$	−0.110	0.055
	Inclination (deg)	90.17	90.21
	Nodal angle (deg)	−2.4	−0.4
	M3 (M⊕)	0.0	12.0

Download table as:  ASCIITypeset image

For each parameter's posterior distribution, we computed the mode and the lower and upper boundaries of the region that contains 68.3% of the area under the curve. We show plots of posterior distributions for several fitting and derived parameters in Figures 10–17, and plots of the two-parameter joint posterior distributions for the 28 major orbital parameters in Figures 18–20 (Table 9 gives the displayed parameters and their ranges).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When one examines the posterior distribution of a given parameter, the mode of the distribution might be taken as the "most likely" value of that parameter. We note, however, that a model constructed using the modes of each distribution will not necessarily provide a good fit to the data. Therefore, for the final adopted parameters, we take the parameters from the best-fitting model.

For several CBPs (Doyle et al. 2011; Welsh et al. 2012; Kostov et al. 2016), the gravitational pull of the planet is sufficiently large to measurably alter the positions and velocities of their host stars. These perturbations can be detected by measuring slight deviations from strict periodicity in the times of the stellar eclipses. From these eclipse-timing variations, the mass of the planet can be estimated. For Kepler-47, only upper limits on the masses could be made at the time of their discovery: <2.0 Jupiter masses (<635 M_⊕) for the inner planet and <28 Jupiter masses (<8900 M_⊕) for the outer planet (Orosz et al. 2012b).

The presence of the third planet dramatically changes this situation. The two outer planets gravitationally perturb one another, resulting in dynamical transit-timing variations. Due to their proximity, the planets perturb each other more than they affect the stars. In the discovery paper (Orosz et al. 2012b), the planets were assumed to be massless in the modeling process, as the interactions between the planets and the binary—and between each other—was determined to be small. In the present model, all five bodies have mass and the mutual interactions between each body are fully accounted for. We find that we can place meaningful constraints on the masses of the middle and outer planets. The 1σ range of the mass of the inner planet is between zero and 25.8 M_⊕, whereas the mass of the middle planet is different from zero at the ≈5σ limit and the mass of the outer planet is different from zero at the ≈2σ limit.

We computed the planet transit times of the best model with a zero-mass middle planet and the best model with a zero-mass outer planet, and compared those times to the transit times from the overall best-fitting model. The difference in the respective transit times are shown in Figure 21. In either case, the transit times of the inner planet do not change by more than ≈1 minute. On the other hand, a massless middle planet results in changes in the outer planet's transit times by up to ≈20 minutes. Likewise, a massless outer planet results in changes in the middle planet's transit times by up to ≈10 minutes. These changes in the transit times can be detected at the significance of a few σ, leading to the mass constraints that we have on the middle and outer planets.

Zoom In Zoom Out Reset image size

Uncertainties in the planet masses have improved by over an order of magnitude compared with the discovery paper (Orosz et al. 2012b). The mass estimates for the inner and outer planets are now <26 M_⊕ and ∼2–5 M_⊕, at the 1σ level. The middle planet has a mass ∼7–43 M_⊕ (1σ range). As usual, the radii are determined much better than the masses: 3.05, 7.0, and 4.7 R_⊕ for the inner, middle, and outer planets, respectively. The middle and outer planets have low bulk densities, <0.68 and <0.26 g cm⁻³ at the 1σ level (Saturn's density is 0.69 g cm⁻³). The densities of all currently known low-mass planets are shown in Figure 22; models indicate that such low-density planets must have substantial hydrogen and helium atmospheres (Lopez & Fortney 2014; Jontof-Hutter et al. 2016). There is also an apparent trend that highly irradiated planets tend to have high densities while planets with low incident fluxes have a range of densities. Kepler-47 c and Kepler-47 d fit this latter trend.

Zoom In Zoom Out Reset image size

Several studies have shown that the formation of CBPs at close distances to the binary is not expected to proceed efficiently, and that CBPs likely formed at large distances and migrated to their current orbits (Pierens & Nelson 2008, 2013; Meschiari 2012; Paardekooper et al. 2012; Pelupessy & Portegies Zwart 2013; Martin et al. 2013; Rafikov 2013; Kley & Haghighipour 2014, 2015; Lines et al. 2015b). Planets that survive the migration phase are expected to have orbits with small eccentricities and small mutual inclinations relative to the orbital plane of the binary stars. The planets in Kepler-47 have low masses compared with Jupiter or Saturn, and the orbits have low eccentricities: e < 0.030, 0.024, 0.05, 0.07 (at the +1σ level) for the binary and the inner, middle, and outer planets, respectively. All four orbits have mutual inclinations aligned to within 16 of one another (+1σ level). This nearly circular, co-planar, packed configuration is unlikely to have arisen as an outcome of strong gravitational scattering of the planets into their current orbits. Rather, the observations suggest that this planetary configuration is the result of relatively gentle migration in a circumbinary protoplanetary disk.

Numerical integrations have shown that the region around the binary where the orbit of a planet will become unstable extends to ≈0.18 au (Hinse et al. 2015). Thus, the semimajor axis of the innermost planet (0.2877 au) is much farther beyond the outer edge of the unstable region. The integrations of the observed system (using the nominal masses given in Table 6) also show evidence of stability, as the variations in the semimajor axis, eccentricity, ascending node, and inclination relative to the binary plane, i_rel, do not have appreciable secular changes and are remarkably flat in their maximum and minimum values over the 100 Myr integration timescale (see Figure 23). The angular momentum deficit (AMD) of the system is used to estimate the changes of secular perturbations that could induce instability on a long timescale. We scale this quantity relative to the present-day solar system value for the terrestrial planets. The Kepler-47 planets are more massive than our terrestrial system, thereby allowing for higher values in AMD. Finally, we show a periodogram for the cyclical variations in inclination within the system and find the precession timescales of the inner, middle, and outer planets to be ∼10 yr, 245 yr, and 738 yr, respectively (see Figure 23). Figure 24 shows the evolution of the impact parameters (for transits across the primary star) for the three planets. The horizontal dashed lines in the figure denote impact parameters where transits of the primary star could occur. As one might expect, the precessionary motion affects the duration of time for which a given planet can transit relative to our line of sight (Schneider 1994).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We examined the possibility of additional planets within the Kepler-47 system. Before exploring the possibility of the existence of an actual planet, we began by integrating the orbits of a large battery of test particles between the system's planets. The numerical scheme for these simulations used a modified version of the orbital integration package, mercury6, that is specifically designed to evaluate the orbits of CBPs efficiently (Chambers et al. 2002). We chose to investigate the test particles initially on coplanar (relative to the binary plane), circular orbits and varied the initial phase of the test particles between 0° and 360°, in increments of 2°. Additionally, we varied the starting semimajor axis of the test bodies from 0.25 to 1.05 au, such that we could evaluate whether non-transiting planets between either the inner-middle pair or the middle-outer pair could be stable for 10 Myr. Overall, we integrated test particles corresponding to a grid of 801 × 181 initial conditions. From this study, we found that virtually all of the test particles between the middle and outer planets became unstable (experiencing ejection from the system or a collision with another body), indicating that no stable (prograde) planetary orbits exist between these two planets. In other words, the region between the middle and outer planets is dynamically full (i.e., these planets are in a packed orbital configuration). However, there remained a broad region of parameter space where the test particles were stable between the inner and middle planets. Despite this stability, there is insufficient evidence to suggest that a massive body exists, because of the many dynamical interactions that would influence the observations of the transiting planets and thereby alert us to its presence. There are many regions where the test particles' eccentricities were significantly increased due to mean motion resonances with the middle planet and we expect these locations to be eroded over time in a similar manner to how the Kirkwood gaps in the solar system are.

Based on analytical studies, we used the estimate of Chambers et al. (1996), assuming that the binary can be approximated as a single star with a large oblateness, and found that the separation in mutual Hill radii of the inner-middle pair and middle-outer pair to be 32.90 ± 17.44 R_H,m and 12.34 ± 4.43 R_H,m, respectively. These values are well above the critical value of $2\sqrt{3}\approx 3.46$ for single-star systems (Gladman 1993) or the critical range of values (5–7) suggested by Kratter & Shannon (2014) for circumbinary systems. The mutual Hill radii of the middle-outer pair are quite close, even when considering the uncertainties, and this proximity does not allow for any additional planets to orbit between 0.7 and 0.96 au. Based on our numerical tests, assuming nominal mass estimates, the middle and outer planets are "packed." This conclusion is supported by other numerical studies (Smith & Lissauer 2009; Kratter & Shannon 2014; Quarles & Lissauer 2018).

We have calculated a dynamical Mean Exponential Growth factor of Nearby Orbits (MEGNO) map (Cincotta et al. 2003; Goździewski et al. 2001) of the middle planet using the best-fit osculating (Jacobian) initial conditions presented in this work. Initial values for the mean anomalies of the binary and each planet were set to zero. The MEGNO factor $\langle Y\rangle$ is a numerical technique to detect and differentiate between aperiodic (chaotic) and quasiperiodic (regular/stable) motion. The computation of the variational vector necessary to calculate MEGNO was done by solving the variational equations of motion in parallel to the five-body equations of motion. These equations of motion were solved using the accurate adaptive time-step ODEX integrator.¹⁹ A given integration was terminated once $\langle Y\rangle$ was greater than five. The result of computing MEGNO over a grid in semimajor axis versus eccentricity space is shown in Figure 25. The color coding follows the numerical value of MEGNO attained at the end of the integration. Yellow colors indicate chaotic and a blue ($| \langle Y\rangle -2.0| \lt 0.001$) color indicate quasiperiodic dynamics of the third planet. The binary and remaining two planets were started at their best-fit osculating Jacobian elements.

Zoom In Zoom Out Reset image size

Recently, the MEGNO technique was used to study possible positions of the middle planet in the framework of the five-body problem (Hinse et al. 2015). Compared to that study, the quantitative picture of the semimajor axis versus eccentricity phase-space structure has changed significantly in the present analysis. The main difference is that the eccentricity of the outer planet has decreased from 0.41 to almost zero in the newly found best-fit model. This demonstrates the overall destabilizing effect an eccentricity of 0.41 for the outer planet orbit has, as it would greatly limit the number of possible stable orbits of the middle planet. This is a fine example of the possibility of planet packing as demonstrated in Kratter & Shannon (2014). Finally, the location of mean-motion resonances (vertical structures in Figure 25) might change for different initial phases of the planets (when varied within their parameter uncertainties). However, the overall global topology structure of chaotic/quasiperiodic orbits of the third planet would not change dramatically.

While there is no dynamical evidence for an additional planet in the Kepler-47 system, another (low-mass) planet between the inner and middle planets could exist. Also, more planets could have stable orbits sufficiently far outside the orbit of the outer planet. Thus, it would be worthwhile to obtain additional high-precision photometry of Kepler-47 whenever possible, to look for additional transits. From Figure 24, we see that the precession timescales for the three known planets are relatively short, and that the middle and outer planets can transit the primary star for only a small fraction of their precession cycles (∼23% and ∼12% of the time for the middle and outer planets, respectively). If there are additional planets in the Kepler-47 system with similar precession cycles, they might precess into view at a later time.

Our new values for the masses of the stars in Kepler-47 are about 6% smaller than in the discovery paper, and the new radii are about 3% smaller than previously reported. We therefore update our analysis and discussion of the evolutionary state of the Kepler-47 binary. Figure 26 compares the new measurements against Dartmouth models (Dotter et al. 2008) for the nominal spectroscopic metallicity of [Fe/H] = −0.25 (Orosz et al. 2012b), in the same way as we did previously. There is still an age discrepancy for the primary star: the mass/radius combination yields an age of 3.5 Gyr (solid line in top panel), whereas the mass/temperature combination gives 11.5 Gyr. If the radii are given more weight than the temperatures, then the secondary star is now slightly larger than the theoretical predictions, as is seen in many other M dwarfs. As we noted in the discovery paper, the primary temperature derived from color indices points to a higher value than the spectroscopy does (which would actually be consistent with the 3.5 Gyr age). However, given the uncertainties in the reddening, the spectroscopy is perhaps more trustworthy.

Zoom In Zoom Out Reset image size

By way of comparison, Figure 27 shows the Dartmouth models for solar metallicity. The primary age from mass/radius is now 4.5 Gyr (solid line), and mass/temperature gives an age consistent with this within the error bar. The secondary radius would be perfectly consistent with the prediction for this metallicity, as would the secondary temperature. However, solar metallicity is formally ruled out by the spectroscopic analysis (Orosz et al. 2012b) at the 3σ level, so a scenario where the stars have solar metallicity is unlikely.

Zoom In Zoom Out Reset image size

We note also that, in its latest data release (DR2), the Gaia mission (Gaia Collaboration et al. 2016, 2018) has provided an accurate parallax for Kepler-47 (source ID 2080506523540902912) of π = 0.948 ± 0.029 mas, corresponding to a distance of ${1025}_{-29.7}^{+31.5}$ pc (Bailer-Jones et al. 2018), that enables a check on our dynamical radius measurement for the primary star. With a reddening estimate to the system of E(B − V) = 0.056 ± 0.020 from Green et al. (2015), the assumption that A_V = 3.1E(B − V), a bolometric correction from Flower (1996) based on the temperature, and the visual magnitude of the star (V = 15.395 ± 0.010, Henden et al. 2015), we obtain R = 0.928 ± 0.072 R_⊙. This is less precise than our dynamical value of R = 0.936 ± 0.005 R_⊙, but is in excellent agreement.

The habitable zone (HZ) is usually defined as the region around a star such that the stellar incident flux on the planet (i.e., the insolation) would allow an Earth-like planet with a CO₂/H₂O/N₂ atmosphere to permanently maintain liquid water on its solid surface. The situation is more complicated for a CBP system, in that: (i) there are two stars, typically of different temperatures and luminosities; and (ii) the orbital motion of the stars causes the HZ to "orbit" with the binary. Even for a planet on a purely circular orbit, the planet–star distances are always rapidly changing. The HZ is no longer static or spherically symmetric.

Therefore, for a CBP system, the definition of the HZ has to be slightly generalized. Here, we define a circumbinary HZ as the region in space where the time-averaged weighted insolation would allow conditions conducive to liquid water on the surface of an Earth-like planet. The weighted insolation is not simply the total flux incident at the top of the planet's atmosphere (though it reduces to this if one of the stars' flux contribution is negligible compared to the other, as is the case for Kepler-47). Adding the flux from each star is not appropriate because it is not simply a matter of the total energy; the spectral energy distribution is highly relevant. The planet's atmosphere acts to filter the energy received from each star before it reaches the planet's surface, and this filtering is strongly wavelength-dependent. The spectral energy distribution of the light emitted by the star can be reduced to a simple function of the star's effective temperature (T_eff) to a sufficient level of approximation. Thus, the T_eff-weighted sum of the fluxes from each star is used as the total flux received by the planet (Haghighipour & Kaltenegger 2013). To determine the inner and outer boundaries of the circumbinary HZ, we equate the weighted incident flux with the insolation corresponding to those boundaries for a single star as computed by Kopparapu et al. (2013, 2014).

Figure 28 shows a face-on view of the Kepler-47 system. The binary stars are shown in black, the planets in blue, the HZ in green, and the red circle marks the critical instability radius (Holman & Wiegert 1999). The darker green shaded area corresponds to the conservative HZ as defined by the runaway greenhouse and maximum greenhouse conditions. The lighter green regions show the optimistic HZ boundaries, as defined by the recent Venus and early Mars conditions (Kopparapu et al. 2013). The arrow shows the direction of the line of sight from the observer to Kepler-47. The sizes of the stars and planets are not to scale, and their locations are for the epoch of the orbital elements presented in Table 5. The outer planet is within the limits of the conservative HZ, and Planet D (the middle planet) skirts the inner boundary of the optimistic HZ.²⁰

Zoom In Zoom Out Reset image size

In Figure 29, we show the total (unweighted) insolation incident at the top of the atmosphere of the middle planet (in blue) and the outer planet (in black) as a function of time. The dotted lines show the limits of the optimistic HZ and the solid green lines mark the limits of the HZ defined by the moist greenhouse condition. The dashed line shows the boundary corresponding to the runaway greenhouse condition, which is also often used as the boundary for the inner edge of the conservative HZ. The three panels show three different timescales. The left panel spans two binary star periods, and the variations in insolation are dominated by the orbital motion of the primary star. Note the primary and secondary eclipses, as seen from the location of the planets. The middle panel shows one orbital period of the outer planet, and the insolation variations now show the effect of the eccentricities of the planet's orbit. The right-hand panel shows five orbits of the outer planet. The time-averaged total insolation for the outer planet is S = 0.865 in Sun–Earth units, with an rms variation of 7.2% about the mean. For the middle planet (d), the time-averaged insolation is S = 1.67, with an rms variation of 5.5%. For comparative purposes (only), if one assumes a Bond albedo of 0.34 (similar to that of Jupiter or Saturn; see de Pater & Lissauer (2015)), and complete redistribution of the incident energy, the equilibrium temperature is T_eq ≈ 241 K for the outer planet. In the course of this discussion, we must keep in mind that these are low-density planets, with thick H, He atmospheres, and therefore not likely to be habitable.

Zoom In Zoom Out Reset image size

Ground-based observations of eclipse and transit events would be useful in constraining the dynamical model of the Kepler-47 system. The times of primary and secondary eclipses are well-fit by the linear ephemerides given in Equation (1). Because the times of the planet transits do not follow a linear ephemeris, we computed (using ELC) the times, impact parameters, and depths of future transits (for the years 2019 through 2026) for each planet using the best-fitting model given in Table 5. The uncertainties on these quantities were estimated using sets of initial conditions drawn from the posterior samples of solutions found from the DE-MCMC run, and are given in Table 10. The formal uncertainties on the predicted times generally vary between about one to two hours for the inner planet, and about 0.75 to two hours for the middle planet and outer planet. With a depth of about 15%, the primary eclipses should be observable by modest-sized ground-based telescopes. On the other hand, owing to the relative faintness of the system (V = 15.395 ± 0.010, Henden et al. 2015), the 0.5% deep secondary eclipses will be more difficult to observe. With depths between about 0.1% and 0.4%, the transits of three planets will be more challenging still. According to our best-fitting model, the inner planet always transits during its ∼10 yr precession cycle, the middle planet transits about 23% of the time during its ∼245 yr precession cycle, and the outer planet transits about 12% of the time during its ∼738 yr precession cycle. We used the REBOUND software (Rein & Liu 2012) to integrate a set of 368 initial conditions drawn from the posterior sample to estimate the uncertainties on the predicted transit times near the times when the transits stop. We find that the uncertainties of the predicted transit times for the middle planet are about 1.4 days by the time this planet stops transiting in the year 2037.7 ± 3.4 (because the transits are grazing before they stop, the predicted stopping date has a much larger spread). Similarly, the uncertainties of the predicted transit times for the outer planet are about 2.1 days by the time this planet stops transiting in the year 2049.1 ± 6.8.