Dynamics and Formation of the Near-resonant K2-24 System: Insights from Transit-timing Variations and Radial Velocities

The vast majority of our current understanding of the masses and orbits of extrasolar planets is based on two techniques: radial velocities (RVs) and transit-timing variations (TTVs). Typically, RVs constrain ${M}_{p}\sin i$, the planet mass modulo (an unknown inclination angle). For high signal-to-noise data sets, deviations from sinusoidal RV curves can reveal orbital eccentricities, and for a few exceptional systems, non-Keplerian orbital dynamics have been observed (see, e.g., GJ876; Rivera et al. 2010; Nelson et al. 2016; Millholland et al. 2018). For transiting systems, the $\sin i$ ambiguity is negligible and RVs constrain planet mass and bulk composition directly. Such measurements have been made for planets as small as Earth (see, e.g., Kepler-78b; Howard et al. 2013; Pepe et al. 2013). Accordingly, RV mass measurements of transiting planets have helped reveal important trends in planetary bulk compositions, such as the onset of low-density envelopes above ${R}_{p}\approx 1.5{R}_{\oplus }$ (Marcy et al. 2014; Weiss & Marcy 2014; Rogers 2015).

While the early theoretical work on TTVs was developed a decade ago (Agol et al. 2005; Holman & Murray 2005), TTVs were not observed until NASA's Kepler mission provided high-precision, long baseline photometry (Holman et al. 2010). The TTV technique has achieved some remarkable results such as precise mass measurements of small planets in the Kepler-36 system (Carter et al. 2012), the discovery of a Laplace-like resonance in the Kepler-223 system (Mills et al. 2016), and mass measurements of non-transiting planets in the Kepler-88 system (Nesvorný et al. 2013).

While the RV and TTV techniques have been applied to many individual systems, only a handful of systems have benefited from joint analyses. Systems with TTVs have almost exclusively been discovered during the prime Kepler mission (Borucki et al. 2010; 2009–2013), which surveyed only 1/400 of the sky. While ≈40% of Kepler planets are in multi-planet systems (Rowe et al. 2014), planets typically need to be near mean-motion resonance to produce detectable TTVs. Holczer et al. (2016) reported TTVs for ≈260 Kepler planets, but most are too faint for precision RV measurements with current-generation instruments, which typically require host stars with $V\lesssim 13$ mag. As a result, fewer than 10 systems have mass constraints from both the TTV and the RV techniques (Mills & Mazeh 2017).

K2-24 has two known transiting planets, which were observed by Kepler during K2 operations (Howell et al. 2014). Petigura et al. (2016, P16 hereafter) reported mass measurements based on Keck/High Resolution Echelle Spectrometer (HIRES) RVs spanning one observing season. While P16 predicted TTV amplitudes of several hours based on their proximity to the 2:1 mean-motion resonance, the 80-day K2 baseline was too short to observe deviations from linear ephemerides.

Here, we present an extended RV time series and additional transit-timing measurements from Spitzer (Section 2). Our extended RV data set enables tighter constraints on the planet masses and reveals a third candidate planet in the system (Section 3). In Section 4, we perform a joint TTV/RV analysis, which provides improved constraints on planet masses, eccentricities, and core/envelope fractions (Section 5). In Section 6, we interpret the observed eccentricities in the context of system dynamics and formation scenarios, and we conclude in Section 7.

2.1. K2

K2-24 was observed during campaign 2 of the K2 mission from 2014 August 23 to 2014 October 13. To extract transit times, we used the photometry published in P16 and fit individual transits. We multiplied our transit model by a third-order polynomial to account for the long timescale variability seen in the photometry. For each transit, we first adopted the best-fit transit parameters from P16, which assumed linear ephemerides. We then fit the transit allowing the time of conjunction ${T}_{c}$ and the polynomial coefficients to vary. Figure 1 shows the K2 photometry along with the best-fit transit models.

Zoom In Zoom Out Reset image size

Care is required when assigning reasonable uncertainties to the measured transit times. K2 photometry contains correlated, non-Gaussian systematics that are mostly, but not entirely, removed during detrending.⁹ The derived transit times depend most sensitively on photometry collected during ingress or egress, which span one or two 30-minute long cadence measurements. Therefore, outliers have a significant effect on the derived transit times if they occur during ingress or egress. As an example, Benneke et al. (2017) found that a single outlier that occurred during one of the transits of K2-18b resulted in a $\approx 7\sigma$ error in the ephemeris reported in Montet et al. (2015).

We estimated the K2 transit-timing errors errors via bootstrap resampling. For each transit, we created 1000 realizations by randomly shuffling the residuals to the best-fit light curve and adding the shuffled residuals to the best-fit model. We then fit these bootstrap realizations using the methods described above and derived ${T}_{c}$ for each sample. We adopted the standard deviation of the resampled ${T}_{c}$ as the uncertainty on ${T}_{c}$. The bootstrapped uncertainties were roughly twice as large as the formal uncertainties, which assumed white and Gaussian distributed noise. Our measured transit times are listed in Table 1.

Table 1.  Transit Times

Instrument	Planet	i	Tc	$\sigma ({T}_{c})$
			days	days
K2	b	0	2072.7954	0.0011
K2	c	0	2082.6248	0.0006
K2	b	1	2093.6806	0.0013
K2	b	2	2114.5654	0.0009
K2	c	1	2124.9879	0.0006
K2	b	3	2135.4505	0.0012
Spitzer	b	20	2490.6161	0.0011
Spitzer	c	10	2506.0002	0.0014
Spitzer	c	15	2717.5074	0.0015
Spitzer	b	31	2720.5049	0.0016

Note. Following a convention from the Kepler mission, times are given in ${\mathrm{BJD}}_{\mathrm{TBD}}\mbox{--}2454833$.

Download table as:  ASCIITypeset image

2.2. Spitzer

P16 used analytic approximations developed by Lithwick et al. (2012) to predict the expected TTVs of K2-24b and c. These approximations predicted anti-correlated sinusoidal TTVs having a "super-period" of roughly 4 years. Given the proximity of K2-24b and c to the 2:1 mean-motion resonance, P16 predicted large TTV amplitudes of several hours. However, the limited 80 day K2 baseline sampled only 5% of the TTV super-period, too small a fraction for TTVs to accumulate to detectable levels.

To cover a significant fraction of the expected TTV super-period, we used Spitzer to observe two additional transits of K2-24b on 2015 October 27 and 2016 June 13 and two additional transits of K2-24c on 2015 November 12 and 2016 June 10.¹⁰ The combined K2/Spitzer data set includes transit observations at three well-separated epochs, which is sufficient to constrain the mean transit period as well as the amplitude and phase of the approximately sinusoidal TTV signal.

When planning our 2015 Spitzer observations, we centered our observing sequence using the best-fit transit times of K2-24b and c based on the K2 data alone. To account for the substantial uncertainty due to TTVs, we observed K2-24b and c for 14 hr each. As shown in Figure 2, we observed a complete transit of K2-24b and a partial transit of K2-24c. We centered our 2016 Spitzer observations on the best-fit linear ephemeris that incorporated the K2 and 2015 Spitzer observations, and we observed K2-24b and c for 12 and 16 hr, respectively. Again, we observed a complete transit of K2-24b and a partial transit of K2-24c. In hindsight, after collecting the 2015 Spitzer transits, we should have performed a preliminary TTV model using plausible masses and eccentricities in order to better center our 2016 Spitzer observations.

Zoom In Zoom Out Reset image size

Following common practice, we included a 30-minute pre-observation sequence to mitigate the initial instrument drift in the science observations resulting from telescope temperature changes after slewing from the preceding target (Grillmair et al. 2012). To enhance the accuracy in positioning K2-24 on the IRAC detector, observations were taken in peak-up mode using the Pointing Calibration and Reference Sensor (PCRS) as a positional reference. We chose Spitzer/IRAC Channel 2 (4.5 μm) over Channel 1 (3.6 μm) because the instrumental systematics due to intra-pixel sensitivity variations are smaller (Ingalls et al. 2012). Our exposure times were set to 2 s to optimize the integration efficiency while remaining in the linear regime of the IRAC detector.

Following Benneke et al. (2017), we extracted multiple photometric light curves for each Spitzer data set using a wide range of fixed and variable aperture sizes. The purpose of extracting and comparing multiple photometric light curves is to choose the aperture that provides the lowest residual scatter and red-noise. We normalized the light curve by the median value and binned the data to a 60 s cadence. We found that this moderate binning did not affect the information content of the photometry, but provided more signal per data point allowing an improved correction of the systematics.

Raw aperture photometry from Spitzer contains large systematics due to the motion of the target star across the IRAC detector with percent-level intra-pixel sensitivity variations. To extract reliable transit times, we adopted the standard practice of modeling the Spitzer systematics and transit profile simultaneously. We used the pixel-level decorrelation (PLD) algorithm, first proposed by Deming et al. (2015), with modifications described in Benneke et al. (2017).

In our model, the following transit parameters were allowed to vary: transit midpoint ${T}_{c}$, planet-to-star radius ratio ${R}_{p}$/${R}_{\star }$, and impact parameter b. In addition, we parameterized the systematics in the Spitzer model using nine PLD coefficients, a white noise component, and two coefficients describing a polynomial trend of flux with time. Ideally, we would have allowed the transit duration ${T}_{14}$ to vary in our fits. However, because our Spitzer transit observations of K2-24c missed ingress, they could not meaningfully constrain ${T}_{14}$. For both K2-24b and c, we fixed ${T}_{14}$ to the value measured by P16 from K2 photometry. We explored the likelihood surface using Markov Chain Monte Carlo (MCMC). The maximum likelihood fits to the Spitzer photometry are shown in Figure 2, and the associated transit times are listed in Table 1.

2.3. Keck/HIRES Spectroscopy

We obtained  63  spectra of K2-24 using HIRES (Vogt et al. 1994) on the 10 m Keck I telescope between  2015 June 24  and  2017 October 03 . We collected 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 used an exposure meter to achieve a consistent signal-to-noise level of 110 per reduced pixel on blaze near 550 nm. We also obtained a "template" spectrum without iodine. The first 32 of these spectroscopic observations are described in P16.

RVs were determined using standard procedures of the California Planet Search (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 measurement uncertainty of each RV point is derived from the uncertainty on the mean RV of the ∼700 spectral chunks used in the RV pipeline and ranges from  1.5 to  2.1  m s⁻¹. Table 2 lists the RVs and uncertainties.

Table 2.  Radial Velocities

Time	RV	σ(RV)	${S}_{\mathrm{HK}}$
days	m s−1	m s−1
2364.819580	0.85	1.68	0.132
2364.825101	1.72	1.52	0.130
2364.830703	9.99	1.59	0.132
2366.827579	−3.90	1.62	0.128
2367.852646	5.50	1.65	0.130
2373.888150	−3.77	1.78	0.094
2374.852412	−5.65	1.97	0.113
2376.863820	−6.09	1.79	0.131
2377.866073	−2.40	1.76	0.131
2378.834011	−1.33	1.60	0.131

Note. Radial velocities and uncertainties for K2-24. Times are given in ${\mathrm{BJD}}_{\mathrm{TBD}}\mbox{--}2454833$. We also provide the Mount Wilson ${S}_{\mathrm{HK}}$ activity index (Vaughan et al. 1978), which is measured to 1% precision.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Here we present our Keplerian analysis of the K2-24 RVs. The RVs exhibited ≈10 m s⁻¹ peak-to-trough variability that was not associated with the known ephemerides of K2-24b or c, which motivated searches for additional non-transiting planets. Figure 3 shows a Keplerian search using a modified version of the Two-dimensional Keplerian Lomb–Scargle (2DKLS) periodogram (O'Toole et al. 2009; Howard & Fulton 2016). When we measured the change in ${\chi }^{2}$ (periodogram power) between a three-planet fit and a two-planet fit, we found a peak at P = 420 days, with an empirical false alarm probability (eFAP) of 0.8%. While the eFAP was formally below the standard criterion of eFAP <1% for Doppler confirmation, a complete confirmation of this candidate would have required additional vetting such as an assessment of RV/activity correlations, which is beyond the scope of this work. We included this candidate our subsequent orbit fitting because it improved the quality of the RV fits to K2-24b and c.

Zoom In Zoom Out Reset image size

We analyzed the RV time series using the publicly available RV modeling package RadVel (Fulton et al. 2018). RadVel facilitates maximum a posteriori (MAP) model fitting and parameter estimation via MCMC. A Keplerian RV signal may be described by the orbital period P, time of inferior conjunction ${T}_{c}$, eccentricity e, longitude of periastron ω and Doppler semi-amplitude K, i.e., $\{P,{T}_{c},e,\omega ,K\}$. In our fitting and MCMC analysis, we adopted the following parameterization: $\{\mathrm{ln}P,{T}_{c},\sqrt{{e}_{}}\cos {\omega }_{},\sqrt{{e}_{}}\sin {\omega }_{},{K}_{}\}$. Our parameterization of e and ω enforces a uniform prior on eccentricity and prevents a Lucy–Sweeney bias toward nonzero eccentricities (Eastman et al. 2013; Fulton et al. 2018). Our preferred model consists of three Keplerians with eccentricities fixed to zero. We fixed the P and T_c of K2-24b and c to the P16 values. To aid convergence, we imposed a loose Gaussian prior on $\mathrm{ln}{P}_{d}$ of ${ \mathcal N }(\mathrm{ln}(440),1)$. Figure 4 shows the MAP model.

Zoom In Zoom Out Reset image size

Models with more free parameters will naturally lead to higher likelihoods at the expense of additional model complexity. To compare the quality of models of different complexity we used the Bayesian Information Criterion (BIC; Schwarz 1978). Models with smaller BIC are preferred. For the circular, three-planet model, BIC = 366.0. Models where candidate d is allowed to have a nonzero eccentricity were not favored by the BIC = 381.2. Models with only two planets on circular orbits were also disfavored by the BIC = 378.6.

To derive uncertainties on the model parameters, we used RadVel to sample the posterior probability via MCMC. RadVel automatically checks for convergence using the Gelman-Rubin statistic (Gelman & Rubin 1992). For K2-24b and c, our RV-only analysis yields masses of ${16.8}_{-3.1}^{+3.2}$ ${M}_{\oplus }$ and ${19.0}_{-3.8}^{+3.9}$ ${M}_{\oplus }$, respectively. We compare these masses to those determined by the joint TTV/RV analysis in Section 5. If candidate d is a bona fide planet, it has a mass of ${54}_{-14}^{+14}$ ${M}_{\oplus }$ and orbits at a distance of ${1.15}_{-0.05}^{+0.06}$ au. However, we do not treat candidate d in our subsequent analysis or discussion, because we have not performed a thorough confirmation and because it is decoupled dynamically from the inner two planets.

Even though the model with all three eccentricities set to zero was preferred in a BIC sense, we performed an analogous MCMC exploration with eccentric orbits to asses the extent to which the RVs alone constrain eccentricities. The RV data set only ruled out high eccentricity orbits, with upper limits of ${e}_{1}$ <  0.39and ${e}_{2}$ < 0.34 at 90% confidence.

As expected, the Spitzer observations revealed TTVs of several hours. In this section, we present an analysis of the transit times from K2 and Spitzer, folding in the constraints from RVs described in the previous section.

Lithwick et al. (2012, L12 hereafter) developed an analytical model for the TTVs that occur when two planets are near first-order mean-motion resonance (i.e., P₂:P₁ ≈ j:$j-1$, where j = 2, 3, ...). For a complete exposition of this formalism, see L12. Here, we provide a brief summary, in order to illustrate the type of constraints that the TTVs provide.

For planets near, but not in, first-order mean-motion resonance, L12 showed that their transit times, ${T}_{c,i}$, are described by a sinusoidal perturbation about a mean period, P:

$$\begin{eqnarray}&&{T}_{c,i}={T}_{c,0}+{Pi}+\mathrm{Re}(V)\sin {\lambda }_{j}+\mathrm{Im}(V)\cos {\lambda }_{j}.\end{eqnarray} \tag{ 1 }$$

Here, i is an integer index that labels the transit epoch, ${T}_{c,0}$ is the time of the first transit (i = 0), and V is the complex TTV amplitude. The longitude of conjunctions ${\lambda }_{j}$, is an angle that advances linearly with time and is given by

$$\begin{eqnarray}&&{\lambda }_{j}=j{\lambda }_{2}-(j-1){\lambda }_{1},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\lambda }_{1}=\displaystyle \frac{2\pi }{{P}_{1}}(t-{T}_{c,1}),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\lambda }_{2}=\displaystyle \frac{2\pi }{{P}_{2}}(t-{T}_{c,2}).\end{eqnarray} \tag{ 4 }$$

The time it takes ${\lambda }_{j}$ to advance by 2π is known as the super-period P_j, which is given by

$$\begin{eqnarray}&&{P}_{j}\equiv \displaystyle \frac{{P}_{2}}{j| {\rm{\Delta }}| },\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}\equiv \displaystyle \frac{{P}_{2}}{{P}_{1}}\displaystyle \frac{j-1}{j}-1.\end{eqnarray} \tag{ 6 }$$

For the K2-24bc pair, ${\rm{\Delta }}=0.013$ and P_j = 1595 days. The complex TTV amplitudes are given by

$$\begin{eqnarray}&&{V}_{1}={P}_{1}\displaystyle \frac{{\mu }_{2}}{\pi {j}^{2/3}{(j-1)}^{1/3}{\rm{\Delta }}}\left(-f-\displaystyle \frac{3}{2}\displaystyle \frac{{Z}_{\mathrm{free}}^{* }}{{\rm{\Delta }}}\right)\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{V}_{2}={P}_{2}\displaystyle \frac{{\mu }_{1}}{\pi j{\rm{\Delta }}}\left(-g+\displaystyle \frac{3}{2}\displaystyle \frac{{Z}_{\mathrm{free}}^{* }}{{\rm{\Delta }}}\right),\end{eqnarray} \tag{ 8 }$$

respectively, where μ is the planet-star mass ratio, and f and g are order unity scalar coefficients that depend j and Δ and are given in L12. For the K2-24bc, f = −1.16 and g = 0.38. ${Z}_{\mathrm{free}}^{* }$ is the complex conjugate of the following linear combination of the planets complex eccentricities:

$$\begin{eqnarray}&&{Z}_{\mathrm{free}}={{fz}}_{\mathrm{free},1}+{{gz}}_{\mathrm{free},2},\end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray}&&z=e\cos \varpi +{ie}\sin \varpi .\end{eqnarray} \tag{ 10 }$$

Our full TTV model contains the following free parameters: $\{{P}_{1},{T}_{c,1},{\mu }_{1},{P}_{2},{T}_{c,2}$, ${\mu }_{2},\mathrm{Re}({Z}_{\mathrm{free}}),\mathrm{Im}({Z}_{\mathrm{free}})\}$.

We incorporated Gaussian priors of ${\mu }_{1}=48\pm 9$ ppm and ${\mu }_{2}=53\pm 11$ ppm based on our RV analysis in Section 3. We confirmed that Gaussian priors were appropriate by checking that the RV-only constraints on μ₁ and μ₂ are well-described by normal distributions, with negligible covariance (Pearson r = 0.09).

We explored the range of plausible planet masses and orbits given the measured transit times using the Affine-Invariant MCMC sampler of Goodman & Weare (2010). We found that employing parallel tempering dramatically reduced the number of iterations needed for convergence (Earl & Deem 2005). We let 16 walkers evolve for 50000 iterations at five different temperatures, discarding the first 10000 iterations as burn in. We verified that the chains were well-mixed by computing the autocorrelation length scale τ for each chain at each temperature and confirming that τ is much smaller than the number of iterations.

In Figure 5, we display the measured and modeled transit times with respect to an adopted reference linear ephemeris. The models sampled from the posterior are a good fit to the observed transit times and gradually diverge from one another after the last Spitzer measurement. To facilitate future observations of K2-24b and c, we include the predicted transit times and uncertainties through 2025 in the Appendix.

Figure 6 shows the two-parameter joint posterior distributions. Note the strong covariance between μ₁ and μ₂. As expected, the TTVs enabled a tight constraint on the planet mass ratio of ${M}_{p,2}/{M}_{p,1}$ =  ${0.81}_{-0.02}^{+0.03}$. As a point of comparison, the RV-only fits constrained the mass ratio to ${M}_{p,2}/{M}_{p,1}\,={1.10}_{-0.26}^{+0.34}$, which is consistent at the 1σ level.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Note also the strong covariance between μ and ${Z}_{\mathrm{free}}$. The priors on μ₁ and μ₂ help to break the μ–${Z}_{\mathrm{free}}$ degeneracy, and we detect significant nonzero real imaginary components of ${Z}_{\mathrm{free}}$. While ${Z}_{\mathrm{free}}$ only constrains linear combinations of the eccentricities, we could infer that (1) at least one of the planets has a nonzero eccentricity and (2) the eccentricities are likely $| {Z}_{\mathrm{free}}| \sim 0.08$. Recall that the RV analysis in Section 3 only provided upper limits of ${e}_{1}\lt 0.39$ and ${e}_{2}\lt 0.34$. Because the TTVs constrain only linear combinations of the e₁ and e₂, we cannot rule out high eccentricity solutions. However, as we discuss in Section 5, these solutions are unlikely given the low eccentricities typically observed in compact Kepler multi-planet systems.

In the previous section, we presented a joint TTV/RV analysis of the K2-24 system. Here, we provide an updated assessment of planet properties based on our combined TTV/RV analysis in Section 4 and compare them to those presented in P16, which only included RVs. Orbital eccentricities are substantially improved over P16, and we also find modest gains in planet mass precision and constraints on the planets' core/envelope structures.

5.1. Planet Mass

5.2. Core/Envelope Structure

Petigura et al. (2017) examined the distribution of core masses ${M}_{\mathrm{core}}$ and envelope masses ${M}_{\mathrm{env}}$ in a sample of 20 sub-Saturns (${R}_{p}$ = 4–8 ${R}_{\oplus }$), which included K2-24b and c. Planets in this size range are well-approximated by a two-component model consisting of a high-density core and a thick envelope of near solar composition H/He. Lopez & Fortney (2014) constructed a grid of model planets having different ${M}_{\mathrm{core}}$ and ${M}_{\mathrm{env}}$ and computed their radii given different levels of stellar irradiation. For each planet in the sample, Petigura et al. (2017) used the Lopez & Fortney (2014) grid to derive the range of ${M}_{\mathrm{core}}$ and ${M}_{\mathrm{env}}$ consistent with the observed planet mass and radii.

For K2-24b and c, Petigura et al. (2017) derived envelope fractions of ${f}_{\mathrm{env},{\rm{b}}}$ = ${28}_{-6}^{+7}$% and ${f}_{\mathrm{env},{\rm{c}}}$ = ${57}_{-10}^{+9}$%. We repeated this analysis using the updated planet masses and radii and found ${f}_{\mathrm{env},{\rm{b}}}$ = ${26}_{-3}^{+3}$% and ${f}_{\mathrm{env},{\rm{c}}}$ = ${52}_{-3}^{+5}$%. Our new values are consistent with Petigura et al. (2017), but with smaller formal uncertainties. This stems mainly from the improved stellar radius (see Table 3) and from the fact that, in the sub-Saturn size range, radius alone is a good proxy for envelope fraction (Lopez & Fortney 2014).

Table 3.  K2-24 System Parameters

Parameter	Value	Notes
Stellar Parameters
${T}_{\mathrm{eff}}$ (K)	5625 ± 60	(A)
$\mathrm{log}g$ (dex)	4.29 ± 0.05	(A)
[Fe/H] (dex)	+0.34 ± 0.04	(A)
K (mag)	9.18 ± 0.02	(B)
${\pi }_{\star }$ (mas)	5.84 ± 0.05	(C)
${M}_{\star }$ (${M}_{\odot }$)	1.07 ± 0.06	(D)
${R}_{\star }$ (${R}_{\odot }$)	1.16 ± 0.04	(D)
Model Parameters
P1 (days)	${20.88977}_{-0.00035}^{+0.00034}$	(E)
${T}_{c,1}$ (BJD−2454833)	${2072.8855}_{-0.0053}^{+0.0055}$	(E)
μ1 (ppm)	${53.3}_{-5.2}^{+5.2}$	(E)
P2 (days)	${42.3391}_{-0.0012}^{+0.0012}$	(E)
${T}_{c,2}$ (BJD−2454833)	${2082.4485}_{-0.0079}^{+0.0078}$	(E)
μ2 (ppm)	${43.4}_{-4.7}^{+4.8}$	(E)
$\mathrm{Re}({Z}_{\mathrm{free}})$	${0.038}_{-0.003}^{+0.004}$	(E)
$\mathrm{Im}({Z}_{\mathrm{free}})$	${0.070}_{-0.007}^{+0.008}$	(E)
Derived Parameters
$| {Z}_{\mathrm{free}}|$	${0.080}_{-0.007}^{+0.009}$	(F)
${M}_{p,2}/{M}_{p,1}$	${0.81}_{-0.02}^{+0.03}$	(F)
${M}_{p,1}$ (${M}_{\oplus }$)	${19.0}_{-2.1}^{+2.2}$	(F)
${M}_{p,2}$ (${M}_{\oplus }$)	${15.4}_{-1.8}^{+1.9}$	(F)
${R}_{p,1}$ (${R}_{\oplus }$)	${5.4}_{-0.2}^{+0.2}$	(F)
${R}_{p,2}$ (${R}_{\oplus }$)	${7.5}_{-0.3}^{+0.3}$	(F)
${\rho }_{1}$ (g cm−3)	${0.64}_{-0.10}^{+0.12}$	(F)
${\rho }_{2}$ (g cm−3)	${0.20}_{-0.03}^{+0.04}$	(F)
e1	${0.06}_{-0.01}^{+0.01}$	(G)
e2	<0.07 (90% conf.)	(G)

Note. (A) Brewer et al. (2016). (B) 2MASS (Skrutskie et al. 2006). (C) Gaia DR2 (Gaia Collaboration et al. 2018). (D) Derived from A, B, and C using the methodology described in Fulton & Petigura (2018). (E) See Section 4. (F) Derived from the posterior samples of (D) and (E). (G) Same as F, but with the eccentricity prior described in Section 5.

Download table as:  ASCIITypeset image

One challenge in explaining the formation of K2-24c is to determine how the planet acquired such a large envelope, while avoiding runaway accretion. As a point of reference, in the canonical core accretion models of Pollack et al. (1996), Saturn forms first as a ≈12 ${M}_{\oplus }$ core that accretes H/He from the protoplanetary disk. At the crossover mass (i.e., when ${M}_{\mathrm{env}}\approx {M}_{\mathrm{core}}$ or when ${f}_{\mathrm{env}}\approx 50 \%$), runaway accretion begins and Saturn quickly grows to its final mass.

One way to resolve the ${f}_{\mathrm{env}}\approx 50 \%$ problem is to imagine that the disk dissipated right as K2-24c approached the runaway phase. While impossible to rule out, this scenario requires special timing of planet formation and is thus a priori unlikely. More likely, the inferred structure of K2-24c points to an incomplete understanding of core-nucleated accretion and motivates further theoretical explanations of planet conglomeration in the sub-Saturn mass regime.

5.3. Eccentricity

By combining TTVs and RVs, we achieved significantly tighter constraints on eccentricity than those from either technique alone. The full RV data set only provided weak upper limits on the planet eccentricities of ${e}_{1}$ <  0.39  and ${e}_{2}\lt 0.34$. The TTVs, in contrast, constrained ${\mu }_{1}{Z}_{\mathrm{free}}$ and ${\mu }_{2}{Z}_{\mathrm{free}}$ (Equations (7) and (8)). Because RVs constrain planet mass directly, they break some of the μ–${Z}_{\mathrm{free}}$ degeneracy inherent to a TTV-only analysis.

Our TTV/RV model provided the following constraints on $\mathrm{Re}({Z}_{\mathrm{free}})$ and $\mathrm{Im}({Z}_{\mathrm{free}})$:

$$\begin{eqnarray*}\mathrm{Re}({Z}_{\mathrm{free}}) & = & {{fe}}_{1}\cos {\varpi }_{1}+{{ge}}_{2}\cos {\varpi }_{2}={0.038}_{-0.003}^{+0.004}\\ \mathrm{Im}({Z}_{\mathrm{free}}) & = & {{fe}}_{1}\sin {\varpi }_{1}+{{ge}}_{2}\sin {\varpi }_{2}={0.070}_{-0.007}^{+0.008}.\end{eqnarray*}$$

These constraints amount to lines in the ${e}_{1}\cos {\varpi }_{1}$-${e}_{2}\cos {\varpi }_{2}$ and ${e}_{1}\sin {\varpi }_{1}$-${e}_{2}\sin {\varpi }_{2}$ planes with slopes determined by f and g. Because TTVs only constrain linear combinations of e₁ and e₂ there are still significant e₁–e₂ degeneracies, even after folding in the RV constraints. Figure 7 shows the large range of e₁ and e₂ consistent with our TTV/RV analysis. Note, however, that e₁ and e₂ cannot both be zero. Our analysis does not formally exclude high eccentricity solutions. These solutions, however, are disfavored for stability reasons and because TTV-active systems are observed to have eccentricities of a few percent.

Zoom In Zoom Out Reset image size

Various groups have characterized the distribution of eccentricities among large numbers of Kepler multi-planet systems, modeling eccentricities as a Rayleigh distribution parameterized by a mean eccentricity $\langle e\rangle$. Studies of TTV-active multi-planet systems have found $\langle e\rangle$ = 0.01–0.03 (Wu & Lithwick 2013; Hadden & Lithwick 2014). Analyses of transit durations in multi-planet systems where the host stars have well-measured densities have found $\langle e\rangle$ = 0.05–0.07 (Van Eylen & Albrecht 2015; Xie et al. 2016). That TTV-active systems exhibit lower $\langle e\rangle$ than the more general class of multi-planet systems suggests a distinct formation pathway.

Under the assumption that K2-24 is drawn from the population of TTV-active Kepler multi-planet systems, we applied a Rayleigh prior on eccentricity $\langle e\rangle$ = 0.03. Figure 7 shows the joint distribution of e₁ and e₂ including this prior. The eccentricity of K2-24c assumes the prior distribution. Solutions with nonzero e₁ are favored because ${e}_{1}\sim 0$ requires ${e}_{2}\sim 0.2$, which is strongly disfavored by our prior. For the remainder of the paper, we adopt e₁ = ${0.06}_{-0.01}^{+0.01}$ and ${e}_{2}\lt 0.07$ (90% conf.). We discuss the dynamical origins of these eccentricities in Section 6.

Here, we explore the dynamical origins of the K2-24 system architecture. In Section 6.1, we discuss how the system evolves on secular timescales. In Section 6.2, we consider several formation scenarios and assess whether they are consistent with the observed eccentricities.

6.1. Secular Evolution

While K2-24b and c are near the 2:1 mean-motion resonance, they cannot be locked in resonance. Resonant locking generally requires that $e\gtrsim {{\rm{\Delta }}}^{2}/\mu$, and for both planets ${{\rm{\Delta }}}^{2}/\mu \sim 3$. Therefore, the long-term dynamical evolution of K2-24b and c is dominated by secular interactions. The coplanar secular evolution of the planets' eccentricities may be visualized as trajectories in the e-${\rm{\Delta }}\varpi$ plane, where ${\rm{\Delta }}\varpi$ is the angle between the apses.¹¹

We simulated plausible long-term evolutions of K2-24b and c by taking 1000 draws from the posterior samples from Section 5 and integrating them for 10000 years with the Mercury N-body integrator (Chambers 1999). These integrations revealed several qualitative apsidal outcomes: circulation, libration about ${\rm{\Delta }}\varpi =0^\circ$ (aligned apses), and libration about ${\rm{\Delta }}\varpi =180^\circ$ (anti-aligned apses). Indeed, the observational data is not yet precise enough to conclusively determine which of these regimes the systems actually occupies. We show representative examples of circulation and libration in Figure 8. Inspection of these solutions shows that while at present time ${e}_{1}$ is likely larger than ${e}_{2}$, at other phases of the secular cycle ${e}_{2}$ may be larger than ${e}_{1}$.

Zoom In Zoom Out Reset image size

6.2. Origin of Eccentricities

Here, we consider several plausible mechanisms for exciting eccentricities, and assess whether they are consistent with the observed eccentricities of K2-24b and c.

6.2.1. Self-excitation

We first considered the possibility that the eccentricities are self-excited, since gravitational interactions between two planets on initially circular orbits will pump eccentricities up to a certain value. To simulate this, we performed an integration with Mercury using representative planet masses and setting the initial eccentricity to zero. As expected, the planets gained some eccentricity but never exceeded e = 0.005. Eccentricities smaller than 0.005 are excluded by the data (see Figure 7), implying that some other process is required to explain the observed eccentricities.

6.2.2. Divergent Migration Through Resonance

A well-known mechanism to excite eccentricity is divergent migration through mean-motion resonance. In this scenario planets begin interior to resonance with zero eccentricity. As shown in Batygin & Morbidelli (2013), migration through resonance corresponds with a separatrix crossing, after which the planets emerge with nonzero eccentricities and anti-aligned apses (${\rm{\Delta }}\varpi =180^\circ$). As shown in Batygin (2015), the exited relic eccentricities are set by the planet-star mass ratios μ and initial eccentricities, which are usually assumed to be small.

In models of early solar system evolution by Tsiganis et al. (2005), such a resonance crossing is used to trigger the onset of a transient dynamical instability. We note that divergent migration could be driven by gravitational scattering with a planetesimal disk (Minton & Levison 2014).

Figure 9 shows the time evolution of a simulation where K2-24b and c are adiabatically driven through resonance using fictitious forces. During the resonant crossing, eccentricities are quickly excited to ${e}_{1}$ = 0.03 and ${e}_{2}$ = 0.02. In this scenario, ${\rm{\Delta }}\varpi$ is driven to 180 deg, and the libration amplitude is very small. Given that this mechanism produces planets that are stationary in the e–${\rm{\Delta }}\varpi$ plane, we can directly compare the present day e to the predicted values from divergent migration.

Zoom In Zoom Out Reset image size

In Figure 7, we compare the predicted eccentricities to our present day constraints. Eccentricities of $({e}_{1},{e}_{2})=(0.03,0.02)$ are disfavored by the data, both with and without the Rayleigh prior on eccentricity. Moreover, the mechanism that drives divergent migration (e.g., planetesimal scattering) is also likely to damp eccentricities. Therefore, $({e}_{1},{e}_{2})=(0.03,0.02)$ corresponds to upper bounds on the eccentricities the planets could acquire through this mechanism. This tension disfavors divergent resonant crossing as the sole explanation for the planet eccentricities, but future measurements of $e$ and ϖ for both planets would shed additional light on this interpretation.

6.2.3. Disk-driven Stochastic Excitation

Another mechanism that excites eccentricities is stochastic interactions between young planets and a turbulent disk (Adams et al. 2008). Density fluctuations within a turbulent protoplanetary disk cause eccentricities to grow approximately like a random walk, with $\mathrm{rms}(e)\propto \sqrt{t}$. One mechanism to drive density fluctuations is the magnetorotational instability (MRI). In the limit of ideal MRI-driven turbulence, Okuzumi & Ormel (2013) showed that the growth of e can be constructed from analytical arguments:

$$\begin{eqnarray*}\mathrm{rms}(e) & \sim & 0.033{\left(\displaystyle \frac{\alpha }{0.01}\right)}^{1/2}\left(\displaystyle \frac{{\rm{\Sigma }}}{{10}^{3}\,{\rm{g}}\,{\mathrm{cm}}^{2}}\right){\left(\displaystyle \frac{a}{0.1\mathrm{au}}\right)}^{2}\\ & & \times \,{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{-1}{\left(\displaystyle \frac{n}{100{\mathrm{days}}^{-1}}\right)}^{1/2}{\left(\displaystyle \frac{t}{10\mathrm{Myr}}\right)}^{1/2}\end{eqnarray*}$$

where α is Shakura–Sunyaev viscosity parameter, σ is the surface density, and n is the mean-motion. This equation suggests that if planets are embedded in a gas disk for a significant fraction of a 10 Myr disk lifetime, as they must have been to capture their H/He envelopes, they can acquire the several percent eccentricities we observe today.

In order to illustrate this process, we performed a Mercury integration where we subjected the planets to appropriately scaled stochastic velocity kicks over a period of 2 × 10⁵ year. The simulation setup was identical to that of Batygin & Adams (2017). The resulting evolution is shown in Figure 9. Note that unlike in the case of divergent migration through resonance, the apsidal offset ${\rm{\Delta }}\varpi$ takes on a broad range of values, resulting in an observable distinction between the two dynamical excitation mechanisms.

6.2.4. Summary

We considered three mechanisms for exciting planet eccentricities: self-excitation, divergent migration, and stochastic pumping. We found that self-excitation cannot explain the present day eccentricities. Divergent migration produces eccentricities that are qualitatively similar to the values observed today, although the predicted eccentricities are formally inconsistent with our measured values. Stochastic pumping can account for the present day eccentricities.

We stress that this is not an exhaustive analysis of excitation mechanisms. Among the mechanisms considered, however, stochastic pumping remains the most plausible explanation, given the data. Divergent migration predicts specific values for ${e}_{1}$, ${e}_{2}$, and ${\rm{\Delta }}\varpi$ which can be corroborated with future observations. For example, measurements of secondary eclipse times place tight constraints on $e\cos \omega$. When combined with the constraints from this paper, such measurements would constrain e and ϖ separately.

We have presented a joint TTV/RV analysis of the K2-24 system based on RVs from Keck/HIRES and transit observations with K2 and Spitzer. Our analysis provides new constraints on planet masses and core/envelope structure. Importantly, we leveraged the synergies between TTV and RV measurements to provide tight constraints on planet eccentricities of ${e}_{1}\sim {e}_{2}\sim 0.08$. Assuming the planets are drawn from the ensemble of Kepler multi-planet systems, we found a small, but significantly nonzero eccentricity of ${0.06}_{-0.01}^{+0.01}$ for K2-24b, and we ruled out eccentricities larger than  0.07  for K2-24c. These eccentricities are relics of the planets' past formation histories, and we found that stochastic interactions with a gas disk are a viable explanation for the observed dynamical state.

Future advances in the exoplanet census and RV instruments will expand the number of systems amenable to similar studies. Next-generation RV facilities at large telescopes such as VLT/ESPRESSO (González Hernández et al. 2017), Keck/KPF (Gibson et al. 2016), and GMT/GCLEF (Szentgyorgyi et al. 2016) will enable RV measurements of a large sample of faint Kepler planet hosts, including many TTV-active systems. Also, ESA's PLATO mission (Rauer 2013) will conduct a transit survey over ≈2000 deg² for 2–3 years and add to the sample of planets with long baseline photometry.

Proceeding along an orthogonal direction, NASA's TESS mission (Ricker et al. 2014) will soon survey the entire sky, casting a wide net for planets around bright stars. These bright stars will be more amenable to RV follow-up than our current sample from Kepler and K2. One challenge is the limited baseline of TESS observations. During a nominal two-year mission, most of the sky would receive 27 days of TESS observations. While this will be sufficient to detect near-resonant systems, the baseline is too short to adequately sample TTV super-periods, which are typically measured in years. Extensions to TESS that would allow for subsequent transit measurements of known planets would therefore be exceedingly valuable.

We thank the anonymous reviewer for helpful suggestions that improved this manuscript. E.A.P. acknowledges support from Hubble Fellowship grant HST-HF2-51365.001-A 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. This work made use of NASA's Astrophysics Data System Bibliographic Services.

We thank Spitzer Science Center Director Tom Soifer for awarding discretionary time. This work is based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. Support for this work was provided by NASA through an award issued by JPL/Caltech.

Some of the data presented herein were 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 NASA. We thank Andrew Howard, Lea Hirsch, Lauren Weiss, Howard Isaacson, and Molly Kosiarek for their assistance with the Keck/HIRES observations. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.

Software: Astropy (Astropy Collaboration et al. 2013), batman (Kreidberg 2015), ChainConsumer (Hinton 2016), emcee (Foreman-Mackey et al. 2013), isoclassify (Huber et al. 2017), lmfit (Newville et al. 2014), matplotlib (Hunter 2007), Mercury (Chambers 1999), numpy/scipy (van der Walt et al. 2011), pandas (McKinney 2010), and RadVel (Fulton et al. 2018).

Facilities: Keck:I (HIRES) - , Spitzer - , Kepler - .

Table 4 lists the predicted transit times and uncertainties for K2-24b and c up to 2025.

Table 4.  Predicted Transit Times

Planet	i	UTC Date	Tc	$\sigma ({T}_{c})$
			days	days
b	0	2014 Sep 05	2072.7962	0.0007
c	0	2014 Sep 15	2082.6250	0.0006
b	1	2014 Sep 26	2093.6803	0.0006
b	2	2014 Oct 17	2114.5650	0.0005
c	1	2014 Oct 27	2124.9877	0.0006
b	195	2025 Oct 31	6146.5011	0.0556
c	96	2025 Oct 31	6146.7839	0.1006
b	196	2025 Nov 20	6167.3924	0.0567
b	197	2025 Dec 11	6188.2829	0.0578
c	97	2025 Dec 12	6189.1158	0.1048

Note. Predicted transit times for K2-24b and c, where i, is an index that labels individual transits. Times are given in ${\mathrm{BJD}}_{\mathrm{TBD}}\mbox{--}2454833$.

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image