PLANET HUNTERS. VII. DISCOVERY OF A NEW LOW-MASS, LOW-DENSITY PLANET (PH3 C) ORBITING KEPLER-289 WITH MASS MEASUREMENTS OF TWO ADDITIONAL PLANETS (PH3 B AND D)^*

Since the first planet discoveries in the 1990s (Wolszczan & Frail 1992; Mayor & Queloz 1995), more than 1400 planets have been discovered, according to the Exoplanet Archive (Wright et al. 2011). Kepler was launched in 2009 and obtained precise photometric measurements for ∼160, 000 stars with nearly continuous coverage for four yr. More than 3500 candidates have been discovered via the planet transit method (Borucki et al. 2011; Batalha et al. 2013; Burke et al. 2014), with nearly 1000 of them now confirmed (e.g., Lissauer et al. 2014; Rowe et al. 2014). The transit planet search (TPS) algorithm (Jenkins et al. 2002, 2010) has been used to search Kepler stellar light curves for the characteristic dips in flux indicative of a planet. Those that meet a certain significance threshold are placed on the Threshold Crossing Event (TCE) list (Tenenbaum et al. 2013, 2014). The TCEs are further examined and can be upgraded to Kepler Object of Interest (KOI) candidate status or downgraded as false positives (FPs), which may be erroneous signals or real, astrophysical phenomena like eclipsing binaries (Prša et al. 2011; Slawson et al. 2011), either on the target star or in the background.

The photometric transit technique can determine the radius of a planet, but generally not the mass and hence does not immediately indicate if a transit signal is due to a planet or a binary star system. In certain cases, these transiting planet candidates can be confirmed as planets through a statistical elimination of astrophysical FPs, such as eclipsing binaries (Fressin et al. 2012). Validation by multiplicity has been used on hundreds of candidates in multiple planet candidate systems (Lissauer et al. 2014; Rowe et al. 2014). This method uses the fact that there are very few FPs in multiple transiting planet candidate systems due to the rare occurrence of two independent, unlikely events (either multiple FPs or a FP and a true planet) both occurring for the same star. However, neither of these methods can determine the planetary masses or orbital parameters, which is necessary for further characterization of the planets. To date, the only way to determine masses purely from a photometric light curve is through the detection and inversion of transit timing variations (TTVs), deviations from a perfectly linear ephemeris due to gravitational interactions of neighboring planets (see Section 1.3). In this paper, strong TTVs of the planets orbiting PH3 allow for the determination of planetary masses and thus their confirmation as planets.

1.1. Planet Hunters

1.2. KOI-1353

1.3. Transit Timing Variations and Mean Motion Resonances

Since the transit technique provides a planet's radius, a mass measurement can both confirm the planetary nature of the candidate as well as measure the planetary density. In some systems with multiple planets, TTVs can determine or provide constraints on the planet masses (Miralda-Escudé 2002; Agol et al. 2005; Holman & Murray 2005). TTVs of transiting planets can even be used to indicate the presence of non-transiting planets (Ballard et al. 2011; Nesvorný et al. 2012), co-orbital (Trojan) planets (Ford & Holman 2007), or exomoons (Kipping 2009). TTVs were first used to confirm a pair of planet candidates orbiting Kepler-9 (Holman et al. 2010). Since then, numerous other Kepler candidates have been confirmed using TTVs (e.g., Lissauer et al. 2011; Steffen et al. 2012; Xie 2014).

One potential disadvantage of measuring planet masses with TTVs is that, in some circumstances, the masses become degenerate with the eccentricity (Lithwick et al. 2012). This occurs when: (1) the planet system is close to, but not within, a first-order MMR; and (2) the measurement errors are significant enough to only detect the sinusoidal TTV caused by the nearby resonance. In this case, the amplitude of the TTV depends both on the planet masses and on the free eccentricity, z_free, which is the component of eccentricity that is not driven by the resonant forcing. If the z_free is known, the mass is also known. Lithwick et al. (2012) examines the degeneracy closely for systems near MMR. In a z_free = 0 system, the TTVs of pairs of interacting planets will be anti-correlated. If the two sets of TTVs are not anti-correlated, this implies that significant free eccentricities exist. However, the absence of a phase shift does not necessarily mean that free eccentricities are zero (Lithwick et al. 2012). This could instead be an unlikely moment in time when free eccentricities do exist, but their phases cancel out. If z_free is small, then one can calculate an upper limit on the mass by assuming z_free = 0 (Lithwick et al. 2012). Determining z_free allows one to calculate the true mass.

The mass–z_free degeneracy can also be broken with more precise data. An O−C TTV signal for a pair of planets near MMR results in a periodic signal (sinusoidal for cases of low z_free). This period is equal to

$$\begin{equation} P_{\rm {TTV}}=\frac{1}{j/P_j - (j-1)/P_{j-1}}, \end{equation} \tag{ 1 }$$

where the outer and inner planets are in a near j: j − 1 resonance, respectively, and j is a positive integer (Agol et al. 2005; Lithwick et al. 2012). This period is called the super-period or TTV period. After one fits the transit midpoints with a linear period and the TTV period, a combination of high quality data and large amplitude TTVs can show the proportionally smaller residual chopping signal (Agol et al. 2005; Fabrycky et al. 2012). This chopping signal is best seen for the inner planet of a near-resonant system and is not seen as strongly in the outer planet, as the position of the inner planet at the time of the outer planet's transit changes very slowly when near a j: j − 1 MMR. The chopping signal scales with the mass of the outer planet (Holman et al. 2010) and with the eccentricity (Nesvorny & Vokrouhlicky 2014). This additional constraint allows one to break the mass–z_free degeneracy to calculate the mass of each planet.

2.1. Orbital Fit

Using quarters 1–16 of the Kepler data, we extracted and flattened each transit using the IDL AutoKep program (Gazak et al. 2012). For the high signal-to-noise transits of the outer planet, we used short cadence data where available. We then used a new, modified version of the IDL program TAP (Carter & Winn 2009; Gazak et al. 2012; Eastman et al. 2013) to fit for the orbital parameters of each planet: impact parameter (b), duration (T), the ratio of planet radius to stellar radius (R_p/R_*), the midtransit times, linear limb darkening, quadratic limb darkening (Kipping 2013), and white and red noise. The ratio of semi-major axis to the radius of the star (a/R_*) and the inclination (i) can be derived from these parameters. For the purpose of transit fitting, circular orbits were assumed, which we found to be a good assumption (see Section 3). See Figures 1–3 for the transit fits of PH3 b, c, and d, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The linear and quadratic limb darkening is very poorly constrained for PH3 b and c when analyzing each planet individually. Therefore, we used the posterior distribution for the linear and quadratic limb darkening from the high signal-to-noise transits (depth ≈12500 ppm) of PH3 d's fit as a prior for the inner and middle transit fits. Starspot anomalies are seen in a few transits of PH3 d. We masked out obvious spots, but we did not account for smaller scale starspot anomalies. This causes a slight residual in the light curve of PH3 d, primarily due to the ninth transit (see Figure 3). For PH3 b, TAP was unable to simultaneously fit for both midtransit times and the orbital properties. Therefore, we fit the phase-folded model first for the orbital parameters. Then, holding the orbital parameters fixed, we fit for the midtransit times.

The parameter a/R_* from the TAP fit is poorly constrained for PH3 b and PH3 c, and their best-fit values are inconsistent with PH3 d's best fit. For example, PH3 c's a/R_* is greater than PH3 d's, but it is on an interior orbit (see Table 2). Therefore, we revise the a/R_* for each planet using Kepler's Third Law with the planet's period and the stellar parameters M_* and R_* (derived from stellar density; see Section 2.2) to get (a/R_*)_rev. The two values agree within errors, the unrevised a/R_* approximately one σ higher than (a/R_*)_rev. We use the revised a/R_* to calculate the planet's incident flux (S) and semi-major axis (a). The best-fit transit light curves for PH3 b, c, and d are shown in Figures 1–3, respectively. The orbital and planetary properties are shown in Table 2.

Table 2. Orbital and Planetary Parameters of PH3 b, c, and d

	PH3 b	PH3 c	PH3 d
P (days)	34.5450 ± 0.0005	66.0634 ± 0.0114	125.8518 ± 0.0076
T0 (JD−2454000)	965.6404 ± 0.0040	975.6436 ± 0.0068	1069.6528 ± 0.0077
MP ( M⊕)	7.3 ± 6.8	4.0 ± 0.9	132 ± 17
RP ( R⊕)	2.15 ± 0.10	2.68 ± 0.17	11.59 ± 0.19
ρ (g cm−3)	4.1 ± 3.9	1.2 ± 0.3	0.47 ± 0.06
arev (AU)	0.21 ± 0.01	0.33 ± 0.02	0.51 ± 0.03
S (S⊕)	24.8 ± 4.4	10.7 ± 1.8	4.4 ± 0.8
tdur (hr)	$3.178^{+0.055}_{-0.053}$	$3.557^{+0.072}_{-0.065}$	$8.067^{+0.026}_{-0.024}$
a/R*	$71.1^{+10}_{-20}$	$117.8^{+21}_{-42}$	109.5 ± 1.2
(a/R*)rev	45.9 ± 0.5	70.7 ± 0.7	108.6 ± 1.1
b	$0.04^{+0.66}_{-0.68}$	$0.05^{+0.70}_{-0.76}$	$0.394^{+0.026}_{-0.029}$
i (deg)	$89.59^{+0.30}_{-0.48}$	$89.73^{+0.20}_{-0.38}$	$89.794^{+0.017}_{-0.016}$
RP/R*	$0.0197^{+0.0011}_{-0.0006}$	$0.0246^{+0.0022}_{-0.0009}$	$0.10620^{+0.00049}_{-0.00050}$
MP/M* (× 10−5)	2.0 ± 1.9	1.1 ± 0.2	36.43 ± 4.66
ecos (ω)	−0.0215 ± 0.0255	−0.0035 ± 0.0022	0.0032 ± 0.0066
esin (ω)	−0.0113 ± 0.0239	−0.0108 ± 0.0122	0.0033 ± 0.0086

Notes. Best-fit parameters for the orbital and planetary properties. The period is the mean period given over the length of observations. The best-fit a/R_* from TAP is poorly constrained for PH3 b and PH3 c and are obviously inconsistent with PH3 d. Therefore, we revise it using the Newton's version of Kepler's Third Law with M_*, R_*, and P to get (a/R_*)_rev.

Download table as:  ASCIITypeset image

2.2. Stellar Fit

We derived the stellar density, ρ_*, from the Markov Chain orbital fit analysis of the outer planet, converting the observed transit characteristics into the density of the star at each link in the chain to derive the posterior distribution of the stellar density (Seager & Mallén-Ornelas 2003). In computing the density, we calculate the sky velocity of the planet during transit; this in turn depends on the eccentricity, e, and longitude of periastron, ω, or equivalently the eccentricity vector, e = (ecos ω, esin ω). Each component of the eccentricity vector of the planet was constrained by the standard deviation of its posterior distribution derived from the transit timing analysis, which was found to have a negligible effect on stellar density (see Section 3). This resulted in a stellar density of ρ_* = 1.577 ± 0.066 g cm⁻³.

We obtained a spectrum of KIC 7303287 with the HIRES instrument on the Keck I telescope (Vogt et al. 1994) and performed Spectroscopy Made Easy (SME) on the data (Valenti & Piskunov 1996; Valenti & Fischer 2005). To further constrain the stellar mass (M_*), radius (R_*), and age (t_*), we fit the Padova PARSEC (v1.1) isochrone models (Bressan et al. 2012) to the observed properties of the PH3 host star: (1) the effective temperature (T_eff) and metallicity ([Fe/H]) derived from the SME analysis; (2) the stellar density, ρ_* (rather than log g, as log g is poorly constrained and weakly dependent on the mass of the star on the main sequence); (3) the spectral energy distribution derived from measured stellar magnitudes in three available databases: Sloan Digital Sky Survey g, r, i, and z (Abazajian et al. 2009), Two Micron Sky Survey J, H, and Ks (Cutri et al. 2003), and WISE W1 and W2 (Cutri 2012). For each band we used the reported uncertainty, σ_i, but also accounted for systematic uncertainty in the calibration of the flux by adding in quadrature a magnitude uncertainty, σ₀, that we allowed to float, and multiplied the likelihood by $\Pi _i(\sigma _i^2+\sigma _0^2)^{-1/2}$ (Ford 2006). We assumed a uniform prior on the extinction, A_V, and chose an R_V = 3.1 Milky Way extinction law. We also allowed the distance, D, and t_* to vary.

There are six stellar model parameters: (M_*, [Fe/H], t_*, D, A_V, σ₀). We computed a Markov Chain using the affine-invariant population Markov chain method (Goodman & Weare 2010; Foreman-Mackey et al. 2013) with eighteen chains (three times the number of free parameters). At each point in the chain, the likelihood was computed from the Gaussian probability of the agreement with the nine observed magnitudes, the stellar density, the effective temperature, and the metallicity. To compute these properties, the stellar isochrones were linearly (or log-linearly) interpolated from a grid of stellar models. The stellar isochrones assumed scaled solar abundances ($\rm {[}\alpha \rm {/Fe]}=0$) and sampled metallicity in intervals of 0.1 dex, while we sampled age in intervals of 0.05 dex and mass in intervals that vary with age and metallicity. Table 3 gives the results of this analysis, which show that this is a young, solar-type star. These stellar parameters are consistent with those derived by Rowe et al. (2014) and are moderately consistent with Huber et al. (2014), which are also shown for comparison. The small uncertainties on the stellar parameters are due to the precise determination of the density from the light curve and transit timing analysis and from the precise temperature and metallicity measurement from SME. However, these uncertainties do not account for possible systematic errors in the analysis. Hence, we re-ran our isochrone analysis using the Dartmouth isochrones (Dotter et al. 2008) and found results that were consistent within 1σ in the mean and yielded standard deviations of similar magnitude. Also, see Torres et al. 2012) for a discussion of the systematic error and biases of SME and a comparison to stellar parameter classification technique and MOOG.

Table 3. Stellar Parameters of PH3 (KIC 7303287)

KIC 7303287	This Paper	Huber et al. (2014)	Rowe et al. (2014)
M* ( M☉)	1.08 ± 0.02	$1.16^{+0.31}_{-0.17}$	1.04a
R* ( R☉)	1.00 ± 0.02	$1.60^{+0.83}_{-0.48}$	1.16 ± 0.22
Teff (K)	5990 ± 38	$6279^{+171}_{-215}$	5930 ± 107
log g (cgs)	4.47 ± 0.01	$4.10^{+0.25}_{-0.27}$	4.33 ± 0.15
ρ* (g cm−3)	1.58 ± 0.07	$0.40^{+0.61}_{-0.26}$	0.94 ± 0.42
[Fe/H]	0.05 ± 0.04	$-0.08^{+0.24}_{-0.30}$	−0.06 ± 0.10
L* (L☉)	1.15 ± 0.06	3.58a	1.50a
t* (Gyr)	0.65 ± 0.44	⋅⋅⋅	⋅⋅⋅
t*, gyro (Gyr)	1.0 ± 0.3	⋅⋅⋅	⋅⋅⋅
D (pc)	700 ± 14	⋅⋅⋅	⋅⋅⋅
AV	0.06 ± 0.01	⋅⋅⋅	⋅⋅⋅

Notes. For comparison, we also show the stellar parameters from Huber et al. (2014) and Rowe et al. (2014). We note that Rowe et al. (2014) used a different Keck HIRES spectrum than us as their observational input for stellar parameters. ^aParameter was not explicitly given, so calculated from other parameters as appropriate for the sake of comparison. Due to the unknown posterior distributions of these values, we do not propagate its error.

Download table as:  ASCIITypeset image

PH3 shows strong activity clearly attributable to starspots (see Figure 4). The long, deep duration of PH3 d shows several strong starspot anomalies, with the starspots primarily appearing at phases of about 1.2 hr before and 3.5 hr after the transit midpoint (see Figure 3). This may allow for a future investigation of its spin-orbit obliquity. Starspots may also cause small, minute-scale TTVs themselves (Mazeh et al. 2013). However, the orbital period of the outer planet is long compared to the rotational period, making this investigation difficult and outside the scope of our study. The starspots allow for a determination of the rotational period, t_{*, rot} = 8.648 ± 0.0009 days (McQuillan et al. 2013). Gyrochronology allows one to calculate the age of the star (t_{*, gyro}) using corrected B−V colors and the stellar rotation period (Barnes 2007). Using the rotational period, the B and V magnitudes from the MAST¹⁶ catalog (B = 14.816, V = 14.255), the earlier derived A_V = 0.06, and the Milky Way extinction law (R_V = 3.1), we find t_{*, gyro} = 1.0 ± 0.3 Gyr, which is somewhat larger but consistent with the age derived from the Padova isochrones.

Zoom In Zoom Out Reset image size

Interactions between transiting planets can result in TTVs, which can provide constraints on the orbital parameters and masses of the transiting planets (Agol et al. 2005; Holman & Murray 2005). We modeled the times of transit with TTVFast, a newly written code for computing transit times of multi-planet transiting systems (Deck et al. 2014). A symplectic integrator was used to compute the positions of the planets as a function of time (Wisdom & Holman 1991; Wisdom 2006), while Keplerian interpolation was used to find the times of transit (Nesvorný et al. 2013). The three planets were assumed to be on plane-parallel orbits, and the likelihood was computed using a χ² fit to the times of transit, with the error bars of each planet inflated by a constant factor, f_i, which was added to the model. For each planet, five parameters were used to describe the orbit: the time of the first transit, t_i, the initial orbital period, P_i, the eccentricity vector, (e_icos ω_i,e_isin ω_i), and the mass ratio of the planet to the star, μ_i. With i = 1, 2, 3, there are a total of 18 model parameters. A prior on the likelihood of $f_i^{-N_i/2}$, with N_i being the number of transits of planet i, was added to penalize large values of f_i; this resulted in a median $\chi ^2_{\rm {red}}$ for each planet of order unity (Ford 2006). A prior of 1/e_i was added for each planet to prevent small eccentricities from being disfavored (Ford 2006). One transit for each PH3 b and PH3 c were excluded as outliers with their contribution to χ² ≫ 1.0; their midtransit times were approximately 2456209 and 2455900 JD, respectively. See Table 4 for a list of observed, missed, and predicted future transit times for all three planets.

An affine-invariant ensemble Markov Chain Monte Carlo (MCMC) code was used to compute the posterior distribution on the parameters (Goodman & Weare 2010; Foreman-Mackey et al. 2013). An ensemble of 31 chains was run for 10⁶ generations, reaching a high degree of convergence; the first 60,000 generations were discarded when computing the posterior parameters. Figure 5 shows the resulting TTVs of the three planets with the 1σ confidence intervals computed from the chain; the plotted errors are inflated by f = (1.73, 1.79, 1.35) for the inner to outer planets, respectively. There is evidence of starspot crossings in the transit light curves, which we attempted to address by masking these regions in the light curve modeling. However, this approach cannot take into account the effect of smaller spots crossed by the planets during a transit that lead to variation in the depth of transit that is commensurate with the photometric uncertainties. Such spot crossings do not affect the out-of-transit data, and hence are not accounted for in the red-noise model used by TAP. Spot crossings cause the transit shape to be asymmetric, and hence can skew the fit for the mid-transit time if not accounted for. Hence, we include f_i to account for additional systematic uncertainties on the errors in the times of transit that are not accounted for in our light curve fitting. Figure 6 shows a "lava plot" (or "river plot"), demonstrating the wavy nature of TTVs. The results show that the phase of the 1325 day TTV super period is not exactly anti-correlated, requiring a small free eccentricity in the system, but this has a negligible contribution to ρ_*. Table 2 gives the best-fit parameters for the planets' orbital parameters and masses. All three-planet masses are well within the planetary regime, confirming them as true planets.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The outer two planets in KOI-1353 are dynamically interacting, and this effect is large since they are close to a 2:1 period ratio. Their TTVs appear as anti-correlated sinusoidal variations with a period of P_TTV = 1/(2/P₃–1/P₂) = 1325 ± 5 days, with P₂ = 66.0634 ± 0.0114 days, and P₃ = 125.8518 ± 0.0076 days (see Equation (1) and Figure 5). These two planets are near commensurate, but not in resonance. The observed chopping signal breaks the mass–eccentricity degeneracy for this pair of planets. Since the middle planet has the largest amplitude TTV, the mass of the outer planet is constrained more precisely than the mass of the middle planet.

We tested the long-term stability of the system by numerically integrating a random set of 1000 initial conditions which produced TTVs consistent with the Kepler data. These 1000 initial conditions and corresponding planetary masses were drawn from the posterior distribution resulting from fitting the Kepler data using MCMC and hence form a statistically representative set of solutions. We integrated these solutions for 10⁸ orbits of the inner planet, or ∼10⁷ yr, using a Wisdom–Holman mapping (Wisdom & Holman 1991) in combination with a third-order symplectic corrector (Wisdom 2006). We employed a time step of one day which resulted in a maximum fractional energy error of ∼ × 10⁻¹⁰.

We then determined the maximum full oscillation amplitude for the semi-major axis and eccentricity of each of the planets. Larger variations in semi-major axis compared to the average value are a sign of instability, whereas large eccentricity oscillations may result simply from large free eccentricities and, in this case, do not indicate unstable behavior. During the 10⁸ orbits of the inner planet, we only found deviations of less than a percent in the semi-major axes. The median eccentricity oscillations of the inner and outermost planet during these integrations were on the order of 30% of the mean value, but the middle planet had large eccentricity oscillations of full amplitude ∼150% the mean value. Out of all the initial conditions, the maximum eccentricity oscillations were on the order of 200% of the mean value. These results are summarized in Table 5.

Table 5. Results for Semi-major Axis and Eccentricity Oscillations about the Mean ($\overline{a}$ and $\overline{e}$, Respectively) after 10⁸ Orbits of PH3 b

Planet	Max [$(\Delta a)/\bar{a}]$	Median [$(\Delta e)/\bar{e}]$	Max [$(\Delta e)/\bar{e}]$
	(%)	(%)	(%)
PH3 b	0.0443259	34.0028	219.367
PH3 c	0.302436	143.627	212.864
PH3 d	0.0144348	27.6211	202.432

Download table as:  ASCIITypeset image

The integrations indicated that the orbits were quiescent. We followed a subset of 100 of them for 10 times longer, for 10⁹ orbits of the inner planet or 10⁸ years. We find no significant deviations from the oscillation amplitudes reported in Table 5.

Although there is no analytic criterion for the stability of three-planet systems, one can gain some insight by studying each pair of planets individually. Both the inner pair and outer pair satisfy the two-planet Hill criterion for stability (e.g., Gladman 1993) as well as the heuristic three-planet stability criterion used by (Fabrycky et al. 2014) for other Kepler systems. Furthermore, neither pair is near a low-order MMR. Although integrations of a length shorter than the age of the system cannot prove the system's stability, our tests strongly suggest that these orbits are long-lived.

With such low eccentricities, masses, and period ratios of 1.91, neither pair of planets is close enough to the 2:1 MMR for it to substantially affect the dynamics of the system. It is interesting to also establish how close the system is to a Laplace (or three-body) resonance. Note, however, that typically, three-body resonances are only important if one or both of the pairs is itself very close to a two-body resonance.

The equation of motion for the three-body resonance angle of the form ϕ = pλ₁ − (p + q)λ₂ + qλ₃, where p and q are integers and λ_i is the mean longitude of the planet i, resembles at lowest order the angle of a simple pendulum (Aksnes 1988; Quillen 2011). The Hamiltonian for this is written as

$$\begin{equation} H = \frac{1}{2} A_{p,q} p_\phi ^2 +\epsilon _{p,q} \cos {\phi }, \end{equation} \tag{ 2 }$$

where the coefficients A_{p, q} and _{p, q} depend on which (p, q) resonance is being considered. Therefore, one can approximate the full width of the resonance region (where the angle ϕ oscillates) simply as $\Delta p_\phi \sim 4 \sqrt{\epsilon _{p,q}/A_{p,q}}$. This can then be converted into a width in terms of the semi-major axis of one of the planets.

We follow the work of Quillen (2011), which focuses on the three-body resonances of circular orbits not near any two-body resonances and which shows that for equal mass planets the width of the resonance in terms of the semi-major axis of the inner planet a₁ roughly scales as

$$\begin{equation} \frac{\Delta a_1}{a_1} \sim \frac{m}{m_\star } \frac{|\log {\delta }|}{\delta } \exp {\left[-(p+q)\frac{\delta }{2} \right]}, \end{equation} \tag{ 3 }$$

with m the mass of one of the planets and for semi-major axis ratios of order α ∼ 1 − δ. In other words, for close pairs of planets the width of the resonance depends linearly on the mass of the planet, relative to the mass of the star, and is significantly larger for closer pairs of planets (smaller δ). Moreover, unless δ is very small, the widths of the resonances with large values of p, q will be exponentially small. This already suggests that the Laplace resonance will not be important for PH3.

We follow the exact formulas given in Quillen (2011), which do not assume equal mass planets, for _{p, q} and A_{p, q} to determine how far away the system PH3 is from each (p, q) three-body resonance. Given the orbital periods of the outer two planets and the masses of all of them, we determine the semi-major axis a_{1, r} that the innermost planet would need such that the system would be in exact resonance (satisfying $\dot{\phi } = 0$), and we also determine the width of the resonance in terms of the semi-major axis of the innermost planet Δa₁. We then computed the number (a₁ − a_{1, r})/Δa₁, the dimensionless number of how many resonance widths in semi-major axis of the inner planet the system was away from exact resonance, for each (p, q) pair, with 0 < p, q < 30.

As expected from Equation (3), the resonances with larger values of (p, q) are entirely negligible. Their widths are too small for the system to be near any of them. The PH3 system is closest to the three-body resonances of the form (p, q) = i(1, 2), with i being an integer greater or equal to unity. However, the system is already ∼60 widths away even from the (1, 2) resonance and is orders of magnitude further from the majority of the resonances looked at. From this, we conclude that three-body resonances in this regime (not near a two-body resonance and with circular orbits) are unimportant for this system.

To constrain the masses of the interior planets' gas envelopes, we consider a scenario wherein they formed by core-nucleated accretion inside the snow-line and consist of an Earth-like composition rocky core (32% by mass iron and 68% by mass silicate) surrounded by an H/He envelope. For a given point in planet mass–radius-incident flux parameter space, planet interior structure models are used to constrain the distribution of each planet's mass between the H/He envelope and the heavy element interior. We model each planet's interior structure following an approach similar to Rogers et al. (2011), but with updated opacities from Valencia et al. (2013). An envelope metallicity of 30 times solar and a Bond albedo of 0.2 are assumed. We sample the posterior distribution on the envelope mass fraction by (1) drawing 10⁵ samples from the planet mass–radius–flux posterior distribution returned by the orbital, stellar, and TTV MCMC fits, and (2) computing the envelope mass fraction for each sample from the planet interior structure models.

6.1. PH3 b

With the large relative uncertainty in its measured mass, the bulk compositions that are consistent with the measured properties of PH3 b encompass a wide range of possibilities, including rocky planet, water-planet, and H/He-enshrouded rocky core scenarios. PH3 b could be a rocky planet with its transit radius defined by a rocky surface; integrating the posterior probability distribution on PH3b's mass, radius, and incident flux, we find a 16% posterior probability that PH 3b is more dense than a pure silicate sphere. Most of the posterior probability on PH3 b's properties, however, falls in a low-density regime in which the planet must have a volatile envelope comprised of some a combination of astrophysical ices (dominated by water), hydrogen, and helium. At the low-mass extreme, there is a 12% posterior probability that PH3 b requires a hydrogen-dominated envelope and is less dense than a pure H₂O sphere, but scenarios where the planet consists of a mixture of rock and water are viable at intermediate masses (and notably at masses near the best fit). In the scenario where PH3 b has a rocky Earth-like composition heavy element interior surrounded by H/He, the H/He mass fraction of PH3 b is determined to be $0.24^{+0.28}_{-0.22}\%$ (see Figure 7) or ΔR_env/R_p = 22% ± 14% ($\Delta R_{\rm {env}} = 0.47^{+0.32}_{-0.30} \,R_{\oplus }$).

Zoom In Zoom Out Reset image size

6.2. PH3 c

6.3. PH3 d

PH3 d lies in the Jovian planet regime; its composition is dominated by hydrogen and helium. Interpolating among model grids from 1 Gyr-old Jovian planets from Fortney et al. (2007), we compute the planet core mass for 10⁵ draws from the mass–radius–flux posterior distribution of PH3 d. In this way, we estimate a heavy-element core mass of 14 ± 4 M_⊕ (see Figure 8), corresponding to an envelope mass fraction of $89^{+0.03}_{-0.02}\%$. We note that Fortney et al. (2007) consider fully differentiated planet with cores that are 50% ice and 50% rock by mass. If the heavy elements are distributed through the planet envelope or are made up of a different mixture of ice and rock, the inferred planet core mass will be affected.

Zoom In Zoom Out Reset image size

PH3 c avoided proper detection by other transit search routines, very likely due to the bias against detecting planets with large amplitude TTVs (García-Melendo & López-Morales 2011). This may be the reason why the Kepler pipeline misidentified this planet as a TCE with a period five times too large. The Quasiperiodic Automated Search (QATS) algorithm (Carter & Agol 2013), designed for detecting planets with TTVs, originally missed PH3 c as well due to PH3's strong stellar variability. However, an upgraded version of QATS can now successfully identify this planet with improved detrending (E. Kruse 2014, private communication).

It is curious that both planet pairs (outer/middle and middle/inner) have very near the same period ratio (∼1.91), with |1 − (P_out/P_mid)/(P_mid/P_in)| = 0.385%. There are 127 stars with 3 + confirmed planets, in which there are 186 unique sets of three consecutive planets¹⁷ (Wright et al. 2011). Only four have the same ratio of period ratios to within 0.385%. The Laplace resonant GJ 876 is one of these, as is the Laplace resonant candidate system HR 8799, although it is a directly imaged planetary system with large error bars in period. The other two are Kepler-207 and Kepler-229. In the 174 KOIs with 3 or more confirmed candidates, there are 242 unique sets of three consecutive planets. Six of these have the same period ratio to within 0.385%: KOIs 665, 757, 869, 1151, 1358, and 2693.

For the half of the distribution of the relative difference in the pair-wise period ratios that is not mathematically truncated at 1.0 (see the red histogram Figure 9), a log-normal fit can be reasonably applied to the distribution. A variable with a log-normal distribution indicates that the variable is the product of many independent, random variables. However, in the blue population, which includes PH3, a log-normal does not fit the allowed region of parameter space (<1). There is an excess of planet triplets near zero. This implies that the two populations are affected by different processes. However, any further analysis of this is beyond the scope of this paper and is left for future studies.

Zoom In Zoom Out Reset image size

As of now, there is only one confirmed exoplanet system in a Laplace resonance (∼4:2:1 MMR): GJ 876 (Rivera et al. 2010; Martí et al. 2013). GJ 876 b has a period of 61.12 days (Marcy et al. 1998; Delfosse et al. 1998), GJ 876 c has a period of 30.09 days (Marcy et al. 2001), and GJ 876 e orbits every 124.26 days (Rivera et al. 2010), while the period ratios of outer/middle and middle/inner of these three planets are the same to three digits, 2.03. The Laplace resonance of GJ 876 was further explored by Martí et al. (2013), who concluded that GJ 876's Laplace resonant solutions are stable, but surrounded by extremely unstable regions. HD 8799 has three confirmed planets (Marois et al. 2008) and also has regions of parameter space that are stable to the 4:2:1 MMR (Reidemeister et al. 2009).

In Figure 10, we show the mass–radius diagram of PH3's planets compared to other confirmed planets (those with masses, radii, and errors from http://exoplanets.org/; Wright et al. 2011), two recent notable results discussed below, and mass—radius-composition models provided by Zeng & Sasselov (2013).

Zoom In Zoom Out Reset image size

With a density of 1.2 g cm⁻³, PH3 c joins a growing population of low-mass (M < 10 M_⊕), low-density planets that require significant H/He gaseous atmospheres, e.g., Kepler-11 (ρ = 1.7, 1.28, 0.66, 0.58, and 0.69 g cm⁻³ for Kepler-11 b, c, d, e, and f, respectively, Lissauer et al. 2013), Kepler-87c (ρ = 0.152 g cm⁻³, Ofir et al. 2014), GJ 1214 b (ρ = 1.87 g cm⁻³, Charbonneau et al. 2009), and Kepler-36c (ρ = 0.89 g cm⁻³, Carter et al. 2012). In the Kepler-79 system (Xie 2014), Jontof-Hutter et al. (2014) calculate more precise masses of its four planets via TTVs and finds all four with low densities. The least dense, Kepler-79d, has ρ = 0.09 g cm⁻³ with a mass of just 6.0 M_⊕. Kepler-51 is even more extreme with two confirmed planets (Steffen et al. 2012) both with densities of 0.03 g cm⁻³ and masses of 2.1 M_⊕ and 4.0 M_⊕ (Masuda 2014). The third planet in the system confirmed by Masuda (2014) has ρ = 0.05 g cm⁻³ and M = 7.6 M_⊕.

The dominant method for discovering these low-mass, low-density planets is via TTV analysis. However, TTV discoveries may be biased to low-density results since, at constant mass, planets with larger radius and thus lower density will be detected more readily and with more precise midtransit times (Jontof-Hutter et al. 2014). The likelihood of transit detection approximately increases with $R^2_{\rm {P}}$. Therefore, although more massive planets would create larger, more detectable TTV signals, the balance is likely in favor of a bias toward discovering low-density planets.

As noted in Masuda (2014), a number of low-density planets are found near MMR. However, TTV detections are best-suited for planets near MMR as they provide the clearest and largest amplitude signals, so this feature may be a selection effect. More discoveries of low-mass, low-density planets and future theoretical work examining this new population of low-mass, low-density planets will be needed to address whether this is a truly a selection effect or a signature of planet formation or migration.

The data presented in this paper are the result of the efforts of the PH volunteers, without whom this work would not have been possible. Their contributions are individually acknowledged at http://www.planethunters.org/authors. The authors thank the PH volunteers who participated in identifying and analyzing the candidates presented in this paper.

E.A. acknowledges funding by NSF Career grant AST 0645416; NASA Astrobiology Institutes Virtual Planetary Laboratory, supported by NASA under cooperative agreement NNH05ZDA001C; and NASA Origins of Solar Systems grant 12-OSS12-0011. K.M.D. acknowledges the support of an NSF Graduate Research Fellowship. D.F. acknowledges funding support for PlanetHunters.org from Yale University and support from the NASA Supplemental Outreach Award, 10-OUTRCH.210-0001 and the NASA ADAP12-0172. L.A.R. gratefully acknowledges support provided by NASA through Hubble Fellowship grant #HF-51313 awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. K.S. gratefully acknowledges support from Swiss National Science Foundation Grant PP00P2_138979/1. M.E.S. is supported in part by an Academia Sinica postdoctoral fellowship.

This research has made use of the Exoplanet Orbit Database and the Exoplanet Data Explorer at exoplanets.org. The Zooniverse is supported by The Leverhulme Trust and by the Alfred P. Sloan foundation. PH is supported in part by NASA JPL's PlanetQuest program. The Talk system used by PH was built during work supported by the NSF under grant No. DRL-0941610. We gratefully acknowledge the dedication and achievements of Kepler science team and all those who contributed to the success of the mission. We acknowledge use of public release data served by the NASA/IPAC/NExScI Star and Exoplanet Database, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This research has made use of NASA's Astrophysics Data System Bibliographic Services. This paper includes data collected by the Kepler spacecraft, and we gratefully acknowledge the entire Kepler mission team's efforts in obtaining and providing the light curves used in this analysis. Funding for the Kepler mission is provided by the NASA Science Mission directorate. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX13AC07G and by other grants and contracts. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX13AC07G and by other grants and contracts. The data presented herein were partly obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation.