TWO TRANSITING LOW DENSITY SUB-SATURNS FROM K2

K2-24 was observed during K2 Campaign 2 with nearly continuous photometry from 2014 August 23 to November 13. The star is listed as EPIC- 203771098 in the Mikulski Archive for Space Telescopes (MAST). It was selected for K2 observations based on K2 Guest Observer proposal GO2104 (PI: Petigura). We list the star's identifying information, coordinates, and photometric properties in Table 1.

Table 1.  Stellar Parameters of K2-24

Parameter	Units	Value	Source
Identifying Information
EPIC ID	⋯	203771098	EPIC
2MASS ID	⋯	16101770–2459251	2MASS
α R.A.	h:m:s	16:10:17.69	EPIC
δ Decl.	d:m:s	−24:59:25.19	EPIC
Photometric Properties
Kp	mag	11.65	EPIC
B	mag	12.22 ± 0.20	APASS
V	mag	11.28 ± 0.10	APASS
${g}^{\prime }$	mag	13.99 ± 0.90	APASS
${r}^{\prime }$	mag	10.71 ± 0.22	APASS
${i}^{\prime }$	mag	10.64 ± 0.01	APASS
J	mag	9.63 ± 0.02	2MASS
H	mag	9.29 ± 0.02	2MASS
Ks	mag	9.18 ± 0.02	2MASS
Spectroscopic Properties
${T}_{\mathrm{eff}}$	K	5743 ± 60	SM, this paper
$\mathrm{log}g$	dex	4.29 ± 0.07	SM, this paper
[Fe/H]	dex	0.42 ± 0.04	SM, this paper
$v\mathrm{sin}i$	km s−1	$\lt 2$	SM, this paper
${S}_{{\rm{HK}}}$	⋯	0.128	this paper
$\mathrm{log}{R}_{\mathrm{HK}}^{\prime }$	dex	−5.26	this paper
Derived Properties
${\mu }_{\alpha }$	mas yr−1	−60.6 ± 2.5	Zacharias et al. (2012)
${\mu }_{\delta }$	mas yr−1	−65.4 ± 2.4	Zacharias et al. (2012)
${M}_{\star }$	${M}_{\odot }$	1.12 ± 0.05	SM, iso, this paper
${R}_{\star }$	${R}_{\odot }$	1.21 ± 0.11	SM, iso, this paper
${\rho }_{\star }$	g cm−3	0.89 ± 0.23	SM, iso, this paper
${L}_{\star }$	${L}_{\odot }$	1.44 ± 0.33	SM, iso, this paper
Distance	pc	181 ± 17	SM, iso, this paper
Age	Gyr	3.2–6.9	SM, iso, this paper

Note. SM: SpecMatch spectrum synthesis code (Petigura 2015). iso: isochrones interface to the Dartmouth suite of stellar isochrones (Dotter et al. 2008; Morton 2015).

Download table as:  ASCIITypeset image

Table 2.  Planet Parameters

Parameter	Units	b	c
Light Curve Fitting
T0	BJD ${}_{\mathrm{TDB}}-2454833$	2072.7948 ± 0.0007	2082.6251 ± 0.0004
P	days	20.8851 ± 0.0003	42.3633 ± 0.0006
i	deg	${89.25}_{-0.61}^{+0.49}$	${89.76}_{-0.21}^{+0.18}$
${R}_{P}/{R}_{*}$	%	${4.31}_{-0.08}^{+0.17}$	${5.94}_{-0.04}^{+0.10}$
T14	hr	${5.48}_{-0.04}^{+0.07}$	${6.47}_{-0.03}^{+0.04}$
T23	hr	${4.95}_{-0.11}^{+0.05}$	${5.70}_{-0.06}^{+0.03}$
${R}_{*}/a$	⋯	${0.035}_{-0.002}^{+0.005}$	${0.019}_{-0.000}^{+0.001}$
b	⋯	${0.37}_{-0.24}^{+0.22}$	${0.22}_{-0.16}^{+0.17}$
${\rho }_{\star ,\mathrm{circ}}$	g cm−3	${1.00}_{-0.33}^{+0.21}$	${1.47}_{-0.23}^{+0.31}$
a	AU	0.154 ± 0.002	0.247 ± 0.004
${S}_{\mathrm{inc}}$	${S}_{\oplus }$	60 ± 14	24 ± 5
${T}_{\mathrm{eq}}$	K	767 ± 177	606 ± 139
${R}_{P}$	${R}_{\oplus }$	5.68 ± 0.56	7.82 ± 0.72
Circular RV model (adopted)
K	m s−1	4.5 ± 1.1	4.6 ± 1.2
γ	m s−1	⋯	−2.5 ± 0.9
${dv}/{dt}$	m s−1 yr−1	⋯	−23.9 ± 9.7
${\sigma }_{\mathrm{jit}}$	m s−1	⋯	3.4 ± 0.7
${M}_{P}$	${M}_{\oplus }$	21.0 ± 5.4	27.0 ± 6.9
ρ	g cm−3	0.63 ± 0.25	0.31 ± 0.12
Eccentric RV model
K	m s−1	5.1 ± 1.2	5.3 ± 1.1
$e\mathrm{cos}{\omega }_{\star }$	⋯	0.20 ± 0.09	0.00 ± 0.09
$e\mathrm{sin}{\omega }_{\star }$	⋯	−0.06 ± 0.16	−0.02 ± 0.15
e	⋯	${0.24}_{-0.11}^{+0.11}$	<0.39 (95%)
γ	m s−1	⋯	−2.7 ± 1.0
${dv}/{dt}$	m s−1 yr−1	⋯	−22.5 ± 9.2
${\sigma }_{\mathrm{jit}}$	m s−1	⋯	2.9 ± 0.6
${M}_{P}$	${M}_{\oplus }$	23.2 ± 5.3	31.0 ± 6.4
ρ	g cm−3	0.70 ± 0.26	0.36 ± 0.12

Download table as:  ASCIITypeset image

We extracted the photometry of K2-24 from the Kepler pixel data, which we downloaded from the MAST. Our photometric extraction pipeline is described in Crossfield et al. (2015) and Petigura et al. (2015). In brief, during K2 operations the telescope is torqued by solar radiation pressure, causing it to roll around the boresight. This motion causes stars to drift across the CCD by ∼1 pixel every ∼6 hr. As stars sample different pixel-phases, inter-pixel sensitivity variations cause the apparent brightness of the star to change. We solve for the roll angle between each frame and an arbitrary reference frame. We model the time- and roll-dependent brightness variations using a Gaussian process. We also adjust the size of our circular extraction aperture to minimize the residual noise in the corrected light curve. This balances two competing effects: larger apertures yield smaller systematic errors while smaller apertures incur less background noise. The circular extraction aperture (r = 3 pixel) is shown in Figure 1. Figure 2 shows both the raw and corrected photometry for K2-24. Our calibrated photometry is available as an online supplement.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We searched through the calibrated photometry using the TERRA algorithm (Petigura et al. 2013). While two sets of transits are clearly visible by eye in the detrended K2 photometry for this star, we rely on TERRA to search through the photometry of 10,000–20,000 light curves per K2 Campaign. After identifying the transits of planets b and c, we re-ran TERRA on the photometry of K2-24 with the in-transit points removed and did not identify any additional transit candidates.

We obtained near-infrared adaptive optics images of K2-24 using NIRC2 on the 10 m Keck II Telescope on the night of 2015 April 01 UT. We used the 1024 × 1024 NIRC2 array and the natural guide star system; the target star was bright enough to be used as the guide star. The data were acquired using the narrow-band Br-γ filter using the narrow camera field of view with a pixel scale of 9.942 mas pixel⁻¹. The Br-γ filter has a narrower passband (2.13–2.18 μm), but a similar central wavelength (2.15 μm) compared the Ks filter (1.95–2.34 μm; 2.15 μm) and allows for longer integration times before saturation. A 3-point dither pattern was utilized to avoid the noisier lower left quadrant of the NIRC2 array. The 3-point dither pattern was observed three times with 1 co-add and a 5.5 second integration time for a total on-source exposure time of $3\times 3\times 5.5\;{\rm{s}}$ = 49.5 s.

The target star was measured with a resolution of 0.055 arcsec (FWHM). No other stars were detected within the 10 arcsec field of view of the camera. In the Br-γ filter, the data are sensitive to stars that have K-band contrast of ΔK = 4.2 at a separation of 0.1 arcsec and ΔK = 7.9 at 0.5 arcsec from the central star. We estimate the sensitivities by injecting fake sources with a signal-to-noise ratio of 5 into the final combined images at distances of N × FWHM from the central source, where N is an integer. Our combined NIRC2 image and contrast curve are shown in Figure 4.

We observed K2-24 with the High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) on the 10 m Keck Telescope I. Between 2015 June 24 and October 3, we obtained 32 spectra through an iodine cell mounted directly in front of the spectrometer slit. The iodine cell imprints a dense forest of absorption lines which serve as a wavelength reference. We also obtained a "template" spectrum without iodine. We used an exposure meter to achieve a constant signal to noise ratio of 110 per HIRES pixel on blaze near 550 nm. Exposure times were in the range 6–12 minutes. RVs were determined using standard procedures of the California Planet Search (CPS; Howard et al. 2010) including forward modeling of the stellar and iodine spectra convolved with the instrumental response (Marcy & Butler 1992; Valenti et al. 1995). The RVs are tabulated in Table 3. We also list the measurement uncertainty of each RV point, which range from 1.5 to 2.0 m s⁻¹ from the uncertainty on the mean RV of the ∼700 spectral chunks used in the RV pipeline.

We measured the strength of the Ca ii H & K lines and found that K2-24 is an inactive star. We see no emission reversal in the cores of these lines. Table 1 lists the median values of ${S}_{{\rm{HK}}}$ using the method of Isaacson & Fischer (2010) and $\mathrm{log}{R}_{{\rm{HK}}}^{\prime }$ computed using $B-V=0.673$ estimated from ${T}_{\mathrm{eff}}$ according to the relation from Valenti & Fischer (2005).

We searched for companions with separations smaller than $\approx 0.1$ arcsec, where our sensitivity to sources from AO imaging declines (see Figure 4). Adopting the methodology in Kolbl et al. (2015), we searched for spectroscopic binaries in our HIRES spectrum. We detect no secondary set of lines from a star having ${\rm{\Delta }}V\lt 5$ mag shifted by more than 15 km s⁻¹ relative to the lines of the primary star. Shifts of Δv ≳ 15 km s⁻¹ correspond to orbital separations of $\lesssim 4\;{\rm{AU}}$.

We analyzed our iodine-free template spectrum from HIRES using the SpecMatch spectrum synthesis code (Petigura 2015). SpecMatch is a general tool for extracting stellar ${T}_{\mathrm{eff}}$, $\mathrm{log}g$, [Fe/H], and $v\mathrm{sin}i$ by fitting high-resolution spectra. SpecMatch generates synthetic spectra at arbitrary ${T}_{\mathrm{eff}}$, $\mathrm{log}g$, [Fe/H], and $v\mathrm{sin}i$ by interpolating between LTE models of Coelho et al. (2005) and applying broadening kernels that account for line broadening due to stellar rotation and macroturbulence and the instrumental profile of the spectrometer. ${T}_{\mathrm{eff}}$, $\mathrm{log}g$, [Fe/H], and $v\mathrm{sin}i$ are adjusted in order to yield the best-matching spectrum in ${\chi }^{2}$ sense. We determined that K2-24 is a metal-rich G3 star having ${T}_{\mathrm{eff}}$ =  5743 ± 60  K, $\mathrm{log}g$ =  4.29 ± 0.07  dex, [Fe/H] =  0.42 ± 0.04  dex, and $v\mathrm{sin}i$ < 2 m s⁻¹. Our uncertainties in ${T}_{\mathrm{eff}}$ and [Fe/H] are based on comparisons with touchstone stars in the literature with stellar parameters from asteroseismology (Huber et al. 2013), LTE-modeling (Valenti & Fischer 2005; Torres et al. 2012), and Rossiter-McLaughlin measurements (Albrecht et al. 2012).

We converted spectroscopic parameters into physical stellar properties using the isochrones python package (Morton 2015), which provides a convenient interface to the Dartmouth suite of stellar isochrones (Dotter et al. 2008). K2-24 is slightly larger and more massive than the Sun: ${M}_{\star }$ =  1.12 ± 0.05  ${M}_{\odot }$ and ${R}_{\star }$ =  1.21 ± 0.11  ${R}_{\odot }$. We list the spectroscopic and derived physical properties of K2-24 in Table 1.

We analyzed K2 transit light curve using the same approach described by Crossfield et al. (2015). In brief, we fit each planet's transit separately using a minimization and Markov Chain Monte Carlo (MCMC) analysis (Foreman-Mackey et al. 2013) using the batman code (Kreidberg 2015) to model the light curves that assumes a linear transit ephemeris for each planet.

When modeling the transit photometry, we adopt a quadratic limb-darkening law. We used the LDTk Limb Darkening Toolkit (Parviainen & Aigrain 2015) to derive limb-darkening coefficients of u₁ = 0.568 ± 0.003 u₂ = 0.098 ± 0.005. We doubled the uncertainties associated with the limb-darkening parameters and incorporated them as Gaussian priors in the MCMC light curve analysis. All of the MCMC parameters show unimodal distributions. The transit profiles alone do little to constrain orbital eccentricity and give upper limits on ${e}_{b}\lt 0.78$ and ${e}_{c}\lt 0.81$ at 95% confidence. Figure 3 shows the K2-24 photometry and best fit models, and Table 2 summarizes the final values and uncertainties.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The transit profile constrains the mean stellar density if one assumes a circular orbit. Since we fit each planet separately, we obtain two independent measurements for ${\rho }_{\star ,\mathrm{circ}}$, ${1.00}_{-0.33}^{+0.21}$ g cm⁻³ and ${1.47}_{-0.23}^{+0.31}$ g cm⁻³. In addition we also have a spectroscopic estimate of ${\rho }_{\star }$ =  0.89 ± 0.23 g cm⁻³. All three estimates of mean stellar density are consistent at the 2-σ level. In our analysis, we have modeled the light curve as a single unblended star. While our AO and spectroscopic observations rule out stars with ${\rm{\Delta }}{Kp}\lesssim 5$ inside ∼4 AU and outside ∼20 AU, we have not covered parameter space entirely. There is a small possibility that our transit profiles could be diluted by an additional star, affecting primarily the derived planet radii. However, we confirm these planets without the need for statistical validation with RVs as described in the following section.

We detected RV variability matching the orbital periods and phases of K2-24 b and c that were measured from the K2 light curve. Measuring the masses of these planets (as described below) confirms their existence and rules out false positive scenarios.

We modeled the stellar RV time series as the sum of two Keplerian orbits. We considered both eccentric and circular orbits. Circular orbits require three parameters per planet: orbital period P, time of transit T₀, and the Doppler semi-amplitude K. In addition, we allowed for an arbitrary RV offset, γ, and a linear acceleration, ${dv}/{dt}$. To assess the quality of a given model we evaluated the log-likelihood, $\mathrm{ln}{ \mathcal L }$, according to the prescription given in Howard et al. (2014). This likelihood definition incorporates RV "jitter" (${\sigma }_{\mathrm{jit}}$), an additional RV uncertainty due to astrophysical and instrumental sources. To guard against non-physical values of K and ${\sigma }_{\mathrm{jit}}$, we parametrized the model using $\mathrm{log}K$ and $\mathrm{log}{\sigma }_{\mathrm{jit}}$. We imposed no prior on $\mathrm{log}{\sigma }_{\mathrm{jit}}$. Because P and T₀ are measured with exquisite precision from the K2 photometry, we held these parameters fixed during our RV analysis. While the K2-24 bc pair's proximity to resonance will result in strong dynamical interactions, we expect TTVs on the order of ∼3–6 hr, not an appreciable fraction of an orbital period (see Section 4.5.1). Therefore, we do not consider departures from strict Keplerian motion when modeling the RVs.

The red curve in Figure 5 shows the maximum likelihood model which we found using the Limited-Memory Broyden–Fletcher–Goldfarb–Shanno optimization routine (Byrd et al. 1995) as implemented in the scipy Python package (Jones et al. 2001). The bottom panels show the maximum likelihood models for each planet individually. We explored the likelihood surface using MCMC as implemented in the emcee Python package (Foreman-Mackey et al. 2013). Table 2 summarizes the median posterior values and the 14% and 86% quantiles. We detected the reflex velocities due to both planets. Assuming circular orbits, K2-24 b and K2-24 c have masses of 21.0 ± 5.4 ${M}_{\oplus }$ and 27.0 ± 6.9 ${M}_{\oplus }$, respectively.

Zoom In Zoom Out Reset image size

Eccentric models included two additional parameters per planet: e, eccentricity and ${\omega }_{\star }$, the longitude of periastron of the star's orbit. Following Eastman et al. (2013), we re-parametrized e and ${\omega }_{\star }$ as $\sqrt{e}\mathrm{cos}{\omega }_{\star }$ and $\sqrt{e}\mathrm{sin}{\omega }_{\star }$, which mitigates the Lucy-Sweeney bias toward non-zero eccentricity (Lucy & Sweeney 1971). The maximum likelihood model is shown as a blue dashed curve in Figure 5, and the MCMC posteriors are summarized in Table 2. When we included eccentricity in the models, the planets have masses 23.2 ± 5.3 ${M}_{\oplus }$ and 31.0 ± 6.4 ${M}_{\oplus }$, respectively. Interestingly, our eccentric models predict e_b =  ${0.24}_{-0.11}^{+0.11}$ while e_c is consistent with zero (<0.39 at 95% confidence).

We assessed the relative merits of the eccentric and circular models using the Bayesian Information Criterion (BIC). The BIC is defined as $\mathrm{BIC}=-2\mathrm{ln}{{ \mathcal L }}_{\mathrm{max}}+k\mathrm{ln}N$ where ${{ \mathcal L }}_{\mathrm{max}}$ is the maximized likelihood, k is the number of free parameters, and N is the number of observations (Schwarz 1978; Liddle 2004). For our RV time series N =  32 . When comparing two models, the model with the lower BIC is preferred. The BIC penalizes models with low likelihood and high complexity. BIC(circular) − BIC(eccentric) = −4.9. Because the best fit circular model has lower BIC, we adopt its associated best fit parameters as our preferred system parameters. However, we discuss the dynamical implications of eccentric orbits in Section 4.3.

We also observe a linear trend in the RVs of −23.9 ± 9.7 m s⁻¹ yr⁻¹. This trend is marginally significant and could indicate an additional body in the system. Following Winn et al. (2010), we consider the range of possible ${M}_{P}\mathrm{sin}i$ and a that could produce the observed trend. A body on a circular orbit with ${M}_{P}\ll {M}_{\star }$ induces a reflex acceleration of ${dv}/{dt}={{GM}}_{P}/{a}^{2}$. The mass and separation of the planet (or star) is given by

$$\begin{eqnarray*}&&{M}_{P}\mathrm{sin}i\sim 42\;{M}_{\oplus }{\left(\displaystyle \frac{a}{1\mathrm{AU}}\right)}^{2}.\end{eqnarray*}$$

We advocate for continued RV monitoring of K2-24 to determine whether the observed RV trend is due to an additional long-period planet.

Our circular orbital solution favors an RV "jitter" (${\sigma }_{\mathrm{jit}}$) of 3.4 m s⁻¹. Here, we assess whether that jitter is consistent with the ensemble of stars monitored with precision RVs. Isaacson & Fischer (2010) derived an empirical relationship between $\mathrm{log}{R}_{\mathrm{HK}}^{\prime }$, B-V color, and ${\sigma }_{\mathrm{jit}}$. The relationship in Isaacson & Fischer (2010) predicts that at star with B − V = 0.65¹⁶ and $\mathrm{log}{R}_{\mathrm{HK}}^{\prime }$ = −5.26 will have ${\sigma }_{\mathrm{jit}}\approx 2.0$ m s⁻¹. One explanation for the higher-than-expected jitter is the presence of additional short-period planets having K ≈ 1–2 m s⁻¹. Again, we encourage additional RV monitoring of K2-24 to search for additional short-period planets not detected in the K2 photometry. The detection such a planet would not only add to the dynamical richness of the K2-24 system, but would also lead to better constraints on the orbital parameters of K2-24 b and c.

K2-24 b & c are among only a handful of transiting sub-Saturns with well-measured masses. With two sub-Saturns in the same system, we have a rare chance to compare the possible compositions of these planets to each other and to the general population of sub-Saturns. We examine possible compositions with the interior and thermal evolution models of Lopez & Fortney (2014) which track the cooling and contraction of planets with H/He envelopes and allow us to convert measured masses, radii, and incident fluxes of K2-24 b & c into estimates of H/He mass fraction.

We modeled planets with solar metallicity H/He envelopes atop a fully differentiated Earth composition core. According to these models, we find that K2-24 b is $24\pm 8$% H/He by mass, while K2-24 c is $48\pm 9$% H/He by mass. K2-24 b & c then have core masses of $17.6\pm 4.3$ ${M}_{\oplus }$ and $16.1\pm 4.2$ ${M}_{\oplus }$, respectively. The uncertainty on the envelope fraction includes the observational uncertainties on planet mass, radius, age, and incident flux along with theoretical uncertainties such as the iron fraction and heat capacity of the rocky core. Our uncertainty in envelope mass is dominated by planet radius errors; our uncertainty in planet core mass is dominated by uncertainties in planet mass.

These planets are sufficiently large that our conclusions are insensitive to variations in the core composition. While pure water cores are likely unphysical, we repeated the above calculations for planets with 98% water cores in order to set a lower bound on the H/He envelope fraction. Using these models of K2-24 b & c, we found envelope mass fractions of $14\pm 5\%$ and $36\pm 8\%$ and core masses of $18.0\pm 4.9$ ${M}_{\oplus }$ and $17.2\pm 5.0$ ${M}_{\oplus }$, respectively. The effect of changing the assumed core composition from Earth-like to pure water has a small effect on the derived core masses (within the statistical uncertainties). We adopt $17.6\pm 4.3$ ${M}_{\oplus }$ and $16.1\pm 4.2$ ${M}_{\oplus }$ as the core masses of K2-24 b & c, respectively.

The inferred core-envelope structures of K2-24 b & c pose some challenges to explaining their formation. How did K2-24 c end up with twice as much gas as K2-24 b despite forming in the same disk with a similar core mass? Another challenge is explaining how K2-24 c is composed of half H/He gas, but somehow avoided runaway accretion as predicted in standard models of core accretion (e.g., Pollack et al. 1996; Lee et al. 2014).

While the different densities of planet pairs like Kepler-36b & c can be understood in terms of differing XUV-driven mass loss histories (Lopez & Fortney 2013), mass loss likely played only a minor role for K2-24 b & c. The planets are only modestly irradiated, and their cores are relatively massive compared to typical hot sub-Neptunes. Using the coupled thermal evolution and photo-evaporation model of Lopez & Fortney (2013), we find that both planets would have only been ≈1% more massive at an age of 10 Myr.

The fact that the K2-24 bc pair are near the 2:1 mean-motion resonance suggests they formed at larger orbital separations and experienced convergent inward migration (e.g., Murray & Dermott 2000; see also Deck & Batygin 2015). Formation at $\gtrsim 1\;{\rm{AU}}$ as opposed to their current locations (∼0.2 AU), together with inward migration, could explain the large inferred envelope fractions (Lee & Chiang 2015a). Lee & Chiang (2015b) derived analytic scaling relations for atmospheric accretion. For planets at ≳1 AU with dust-free atmospheres, they found that the gas-to-core mass ratio scales as ${M}_{\mathrm{core}}^{1}{T}_{\mathrm{eq}}^{-1.5}$ (Equation (24) of their paper), with ${T}_{\mathrm{eq}}$ equal to the equilibrium surface temperature. The dependence on ${T}_{\mathrm{eq}}$ arises because colder planets have lower opacities and therefore cool and accrete faster. Since K2-24 c presumably formed exterior to K2-24 b and had a lower ${T}_{\mathrm{eq}}$, these dust-free accretion models may explain, in a natural way, why the outer planet in the K2-24 bc pair has a more massive envelope.

In the case of dusty atmospheres, Lee & Chiang (2015b) found that the accreted gas fraction is independent of local disk temperature. Thus if the early atmospheres of K2-24 b & c were dusty, different formation locations could not alone explain their different envelope masses. Dusty envelopes might still be accommodated if K2-24 b initially formed with a more massive envelope but lost a large fraction of it to a giant impact. Recent studies of giant impacts have found that a collision with an equal mass impactor can reduce a planet's gas fraction by a factor of ∼2 (Inamdar & Schlichting 2015; Liu et al. 2015), thus providing an alternative explanation for the difference between the inferred gas fractions of K2-24 b and c.

With this one system, it is not possible to distinguish between the formation hypotheses of dust-free gas accretion versus giant impacts. However, it is intriguing that the longer period planet has the much larger envelope fraction, as we might expect from dust-free gas accretion and convergent migration (Lee & Chiang 2015a). As more sub-Saturns are found and characterized, ensemble properties should shed light on their formation. If stochastic processes like giant impacts determine envelope masses, we expect no correlation between orbital distance and envelope fraction. However, if gas accretion is governed by local disk properties, we should see correlations with orbital distance—as is arguably already observed by the well-known increase in the occurrence rate of Jupiter-mass gas giants beyond ∼1 AU (Cumming et al. 2008).

Explaining how K2-24 c roughly doubled in mass while accreting gas, yet somehow avoided runaway accretion is difficult in the context of standard models of core accretion (e.g., Pollack et al. 1996; Lee et al. 2014). The Lee & Chiang (2015a) scenario of dust-free gas accretion in a disk coupled with inward migration raises a concern of fine tuning since it requires that K2-24 c reach the threshold of runaway (its gas-to-core ratio is modeled to be $48\pm 9$%) without actually running away. However, we note that sub-Saturn-sized planets are not common outcomes of planet formation (2.9% of GK stars host a sub-Saturn with P < 100 days; Petigura et al. 2013) and near-resonant sub-Saturns are rarer still.

Given the current dataset, we are hesitant to claim a non-zero eccentricity for planet b. However, our eccentric model does have some precedent among previously discovered systems. GJ876c and b have orbital periods of 30.08 days and 61.12 days respectively, and, like K2-24 b and c, lie just outside 2:1 mean-motion resonance (Delfosse et al. 1998; Marcy et al. 1998, 2001). N-body fits to the GJ876 RVs show that planet c is moderately eccentric ($e=0.25591\pm 0.00093$) while planet b has a nearly circular orbit ($e=0.0292\pm 0.0015;$  Rivera et al. 2010). The high precision eccentricity measurements in this case are from the large Doppler amplitudes (∼100 m s⁻¹) and the detection of resonant interactions.

We consider here the dynamical implications of eccentric orbits, assuming the system has a long-lived orbital architecture. The system dynamics are governed by the two planets' masses, eccentricities, and longitudes of pericenter, ϖ.¹⁷ Instead of performing a uniform exploration of this six-dimensional parameter space, we consider systems drawn from our MCMC exploration of eccentric RV solutions. First, roughly 25% of the models in our MCMC chain satisfy ${a}_{c}(1-{e}_{c})\lt {a}_{b}(1+{e}_{b})$. Given that the distribution of planet ϖ is nearly uniform, many of these solutions correspond to crossing orbits.

We also considered whether the systems are Hill stable, using the full criterion, which is based on conservation of the quantity ${L}^{2}E$, where L and E are the total orbital angular momentum and energy of the system (Marchal & Bozis 1982; Milani & Nobili 1983). Roughly half of the MCMC realizations fail the Hill criterion. These include all solutions with ${e}_{b}\gtrsim 0.3$ or ${e}_{c}\gtrsim 0.3$. While eccentricity of planet b (e_b =  ${0.24}_{-0.11}^{+0.11}$) is likely less than 0.3, the 1-σ confidence interval extends past 0.3. In this case, the system must be in some type of resonant phase protection which prevents close approaches in order to be long-lived (e.g., Gladman 1993; Barnes & Greenberg 2007). On the other hand, orbits which satisfy the Hill criterion, though protected from collisions, are not necessarily long-lived as weak encounters can still lead to large and erratic variations in the orbital elements. To test this, we selected 100 planet masses, orbital eccentricities and longitudes of pericenter randomly from the MCMC chain. We integrated these initial conditions using a Wisdom-Holman mapping with a symplectic corrector employed (Wisdom & Holman 1991; Wisdom et al. 1996). Our timestep was 0.25 days, and we ensured that the fractional energy conservation was high (typically $\sim {10}^{-10}$). The integrations lasted for 10⁶ years, or $\approx 20$ million orbits of the inner planet.

Although roughly half of the orbits failed the Hill criterion, only 9 showed instability during the integrations (deviations in semimajor axes larger than 5% of the initial values). To understand why, we selected the orbits which failed the Hill criterion yet remained long-lived, and looked at the orbital evolution of the eccentricities and the angle ${\rm{\Delta }}\varpi ={\varpi }_{b}-{\varpi }_{c}$ on timescales of 3 × 10⁴ orbits of the inner planet. Roughly 70% show apparently regular evolution, and the majority exhibit libration of ${\rm{\Delta }}\varpi$, about zero or π, indicating a resonant protection mechanism. Although there is no preferred value of ${\rm{\Delta }}\varpi$ based on the RV data, at high eccentricities the 2:1 resonance is wide, and so it is not surprising that many are in resonance, albeit with large libration amplitudes.

The majority of the remaining 30% of the orbits appear chaotic, with erratic variation of eccentricities and alternations between circulating and libration of ${\rm{\Delta }}\varpi$ and/or the (mean-motion) resonant angles $2{\lambda }_{c}-{\lambda }_{c}-{\varpi }_{b}$ and $2{\lambda }_{c}-{\lambda }_{c}-{\varpi }_{c}$. It is interesting that despite this chaotic behavior, the effective lifetimes of these orbits are relatively long. We expect longer integrations would reveal unstable behavior.

In conclusion, this limited look into the long-term stability of the orbital solutions to the RV data suggests that orbits with large eccentricities are plausible, even if the system fails the Hill criterion, if the system is in resonance.

Here, we put the K2-24 system in the context of the other sub-Saturns. Figure 6 shows the densities and radii of planets having ${R}_{P}$ = 4–8 ${R}_{\oplus }$ and density measured to better than 50%, i.e., $\sigma (\rho )/\rho \lt 50\%$. The symbol colors represent the planet equilibrium temperature assuming zero albedo, and the symbol shapes indicate whether TTVs or RVs were used to measure planet mass. The K2-24 planets are labeled in bold. The K2-24 planets are fairly typical compared other sub-Saturn planets.

Zoom In Zoom Out Reset image size

The relative sizes and densities of the K2-24 bc pair are reminiscent of the Kepler-18 cd pair. Kepler-18c has a mass of 17.3 ± 1.9 ${M}_{\oplus }$, radius of 5.49 ± 0.26 ${R}_{\oplus }$, and a density of 0.59 g cm⁻³, similar to K2-24 b. Kepler-18 days has a mass of 16.4 ± 1.4 ${M}_{\oplus }$, a radius of 6.98 ± 0.33 ${R}_{\oplus }$, and a density of 0.27 ± 0.03 g cm⁻³. While Kepler-18b is smaller and less massive than K2-24 c, it has a similar density. Kepler-18 cd also lie near the 2:1 mean-motion resonance (Cochran et al. 2011).

While there are still relatively few sub-Saturns with well-measured masses and radii, there are some trends worth noting. Densities measured from TTVs tend to be lower than RV-measured densities. This trend was noted by Weiss & Marcy (2014) for planets smaller than 4 ${R}_{\oplus }$. Here, we offer some observational and astrophysical explanations. The Doppler semi-amplitude, K, depends primarily on planet mass, and has a weaker dependence on orbital period, eccentricity, and stellar mass. The TTV technique is also sensitive to planet mass; that sensitivity is amplified by a system's proximity to resonance. See Steffen (2015) for a more complete comparison of the sensitivities associated with TTVs and RVs. Thus, it is perhaps not surprising that the lowest density (i.e., lowest mass) sub-Saturns have more TTV than RV measurements. Lee & Chiang (2015a) offer a parallel astrophysical explanation: TTV measurements are most easily made for systems in or near mean-motion resonances; such resonant systems formed by convergent inward migration; the planets comprising a resonant system therefore formed at larger orbital distances where disk gas was colder, less dense, and optically thinner; such gas cools more rapidly and is therefore accreted more readily onto rocky cores, forming especially low-density ("super-puffy") planets (see also Section 4.2).

4.5.1. TTVs

Due to the proximity of K2-24 b & c to the 2:1 mean-motion resonance, coherent gravitational interactions between the planets will result in large TTVs. Lithwick et al. (2012) developed an analytic theory to describe TTVs near first order resonances, i.e., j+1:j resonances where j is an integer. We use this theory to estimate, within an order of magnitude, the amplitude of TTVs in this system. Following Lithwick et al. (2012), the normalized distance to resonance is given by

$$\begin{eqnarray*}&&{\rm{\Delta }}\equiv \displaystyle \frac{{P}^{\prime }}{P}\displaystyle \frac{j-1}{j}-1,\end{eqnarray*}$$

where P is the period of the inner planet and ${P}^{\prime }$ is the period of the outer planet. For the K2-24 bc pair, j = 2 and Δ = 0.014. Near resonance, TTVs are oscillatory with a "super period" given by

$$\begin{eqnarray*}&&{P}^{j}=\displaystyle \frac{{P}^{\prime }}{j{\rm{\Delta }}}.\end{eqnarray*}$$

For the K2-24 bc pair, ${P}^{j}$ is 1490 days or $\approx 4$ years.

Another important quantity that influences TTVs is ${Z}_{\mathrm{free}}$, a linear combination of the free complex eccentricities of the two planets.¹⁸ When $| {Z}_{\mathrm{free}}| \ll | {\rm{\Delta }}|$ or $| {Z}_{\mathrm{free}}| \gg | {\rm{\Delta }}|$ the amplitude of the TTV signal, $| V|$, is given by Equations (14) and (15) in Lithwick et al. (2012):

$$\begin{eqnarray}&&| V| \sim P\displaystyle \frac{{\mu }^{\prime }}{| {\rm{\Delta }}| }\left(1+\displaystyle \frac{| {Z}_{\mathrm{free}}| }{| {\rm{\Delta }}| }\right)\end{eqnarray} \tag{ 1 }$$

and

$$\begin{eqnarray}&&| {V}^{\prime }| \sim {P}^{\prime }\displaystyle \frac{\mu }{| {\rm{\Delta }}| }\left(1+\displaystyle \frac{| {Z}_{\mathrm{free}}| }{| {\rm{\Delta }}| }\right),\end{eqnarray} \tag{ 2 }$$

where $\mu ={M}_{P}/{M}_{\star }$. We consider a useful limiting case where $| {Z}_{\mathrm{free}}| \ll | {\rm{\Delta }}|$. In this case, the respective TTV amplitudes of planets b and c are

$$\begin{eqnarray}&&| V| \sim P\displaystyle \frac{{\mu }^{\prime }}{| {\rm{\Delta }}| }\sim 2.6\quad \mathrm{hr}\end{eqnarray} \tag{ 3 }$$

and

$$\begin{eqnarray}&&| {V}^{\prime }| \sim {P}^{\prime }\displaystyle \frac{\mu }{| {\rm{\Delta }}| }\sim 4.0\quad \mathrm{hr}.\end{eqnarray} \tag{ 4 }$$

Given that K2-24 is bright (V = 11.3) and that the planets are large, TTVs of this magnitude are easily detectable from the ground. We emphasize that the TTVs can be significantly larger when $| {Z}_{\mathrm{free}}| \approx e\gt | {\rm{\Delta }}| \approx 0.01$, and so the nominal TTV amplitude estimated above provides a rough lower limit to the TTV amplitude. Strictly speaking, however, there are some "coincidental" orbital configurations with $| {Z}_{\mathrm{free}}| \sim {\rm{\Delta }}$ where the TTV amplitudes could be smaller than the above estimates. However, such configurations are rare (see Lithwick et al. 2012). Given that $| {Z}_{\mathrm{free}}|$ depends on e_b, e_c, ${\varpi }_{b}$, and ${\varpi }_{c}$, observing and modeling the TTVs over an appreciable fraction of the 4 yr super period will place important constraints on the orbits of K2-24 b and c.

4.5.2. Transmission Spectroscopy