K2-138 g: Spitzer Spots a Sixth Planet for the Citizen Science System

We obtained a 0.38 to 0.7 μm spectrum of K2-138 using the Goodman spectrograph (Clemens et al. 2004) on the Southern Astrophysical Research Telescope (Program ID 2019A-0364; PI: K. Hardegree-Ullman), and a 0.7–2.4 μm spectrum using the SpeX spectrograph (Rayner et al. 2003) on the NASA Infrared Telescope Facility (Program ID 2017A-106; PI: K. Hardegree-Ullman). We followed the procedures outlined in Section 2.1–2.3 of Hardegree-Ullman et al. (2019) for observing the target and reducing the data. We compared the combined optical and infrared spectra between 0.38 and 1 μm to optical Sloan Digital Sky Survey (SDSS) spectral templates from Kesseli et al. (2017) following the procedures outlined in Section 4.1 in Hardegree-Ullman et al. (2020), which yielded a spectral type of G8 V, consistent with the spectral type found by Lopez et al. (2019). Figure 1 shows our 0.38 to 1 μm spectrum compared to G7 V, G8 V, and G9 V template spectra.

Zoom In Zoom Out Reset image size

In Table 1 we list the stellar parameters for K2-138 compiled from the Ecliptic Plane Input Catalog (EPIC; Huber et al. 2016), RAdial Velocity Experiment fifth data release (RAVE DR5; Kunder et al. 2017), Christiansen et al. (2018), Gaia DR2 (Gaia Collaboration et al. 2018; Bailer-Jones et al. 2018), the Transiting Exoplanet Survey Satellite (TESS) Input Catalog Candidate Target List (TIC CTL; Stassun et al. 2019), Lopez et al. (2019), and Hardegree-Ullman et al. (2020). The measured parameters are all consistent to within 1σ, except for one measurement of $\mathrm{log}(g)$ which is within 2σ. It is reassuring that different data sets and pipelines yield similar results, but when it comes to calculating planet parameters, small differences in stellar parameters can have a large impact. For large surveys and population studies of exoplanets, it is crucial to have a uniformly derived set of stellar parameters (e.g., Fulton et al. 2017; Berger et al. 2020; Hardegree-Ullman et al. 2020). For individual systems, however, it is typical to choose a single set of stellar parameters, which may be susceptible to systematic bias. Rather than cherry picking measurements from different references, we combined all the available measurements for T_eff, $\mathrm{log}(g)$, [Fe/H], and M_⋆. Instead of using a weighted mean, which would produce uncharacteristically small uncertainties, ¹⁹ we instead employed the following Monte Carlo method. For each measurement with symmetric uncertainties, we randomly drew 10⁴ values from a Gaussian distribution, and for asymmetric uncertainties we drew 10⁴ values from a split normal distribution. The posterior distributions were concatenated and we took the median, 16th, and 84th percentiles of the resultant distribution as our measurement and errors.

Table 1. K2-138 Stellar Parameters

Spectral Type	Teff	$\mathrm{log}(g)$	[Fe/H]	R⋆	M⋆	Distance	Reference
	K	$\mathrm{log}(\mathrm{cm}\,{{\rm{s}}}^{-2}$)	dex	R⊙	M⊙	pc
⋯	5189 ± 156	4.557 ± 0.030	0.01 ± 0.12	0.807 ± 0.046	0.870 ± 0.052	182.50 ± 17.15	1
⋯	5228 ± 76	4.34 ± 0.12	0.26 ± 0.20	⋯	⋯	242.905 ± 94.869	2
K1 V ± 1	5378 ± 60	4.59 ± 0.07	0.16 ± 0.04	0.86 ± 0.08	0.93 ± 0.06	183 ± 17	3
⋯	${5110}_{-53}^{+209}$	⋯	⋯	${0.917}_{-0.070}^{+0.019}$	⋯	${202.585}_{-1.989}^{+2.028}$	4
⋯	5281 ± 129	4.53 ± 0.08	⋯	0.858 ± 0.049	0.91 ± 0.11	202.585 ± 2.009	5
G8	5350 ± 80	4.52 ± 0.15	0.14 ± 0.10	${0.86}_{-0.02}^{+0.03}$	0.93 ± 0.02	201.66 ± 6.38	6
G7	5303 ± 138	4.547 ± 0.150	0.035 ± 0.235	${0.841}_{-0.050}^{+0.056}$	${0.911}_{-0.280}^{+0.412}$	${202.585}_{-1.989}^{+2.028}$	7
G8 V ± 1	${5283}_{-163}^{+122}$	${4.538}_{-0.159}^{+0.090}$	${0.134}_{-0.176}^{+0.128}$	${0.839}_{-0.051}^{+0.054}$	${0.916}_{-0.096}^{+0.078}$	${202.585}_{-1.989}^{+2.028}$	8

Note. (1) EPIC (Huber et al. 2016), (2) RAVE DR5 (Kunder et al. 2017), (3) Christiansen et al. (2018), (4) Gaia DR2 (Gaia Collaboration et al. 2018; Bailer-Jones et al. 2018), (5) TIC CTL v8.01 (Stassun et al. 2019), (6) Lopez et al. (2019), (7) Hardegree-Ullman et al. (2020), (8) this work.

Download table as:  ASCIITypeset image

We computed a bolometric magnitude (M_bol) using the Gaia distance of ${202.585}_{-1.989}^{+2.028}$ pc reported by Bailer-Jones et al. (2018), accounting for interstellar reddening with the dustmaps package (Green et al. 2018), and applying a bolometric correction found using isoclassify (Huber et al. 2017). We computed the bolometric luminosity using ${L}_{\mathrm{bol}}={L}_{0}\,\times {10}^{-0.4{M}_{\mathrm{bol}}}$, where L₀ = 3.0128 × 10²⁸ W is the zero-point radiative luminosity (Mamajek et al. 2015). From the Stefan–Boltzmann law we derived a new stellar radius of ${0.839}_{-0.055}^{+0.060}\,{R}_{\odot }$ with our Monte Carlo averaged effective temperature and bolometric luminosity. Stellar parameters for K2-138 are listed in Table 1.

EPIC 245950175 (K2-138) was observed with K2 in 30 minute long cadence mode during Campaign 12 between 2016 December 15 and 2017 March 4, with a five day gap in the data about two-thirds of the way through the campaign due to a spacecraft safe-mode event. K2-138 was included in Campaign 12 Guest Observing Programs 12049, 12071, 12083, and 12122 (PI: E. Quintana, D. Charbonneau, A. Jensen, and A. Howard).

For our analysis, we used the instrumental systematics-corrected light curve (Figure 2) produced by the k2phot ²⁰ pipeline (Petigura et al. 2015; Aigrain et al. 2016). We compared the k2phot light curve to those produced by EVEREST (Luger et al. 2016) and K2SFF (Vanderburg & Johnson 2014), and found that the k2phot light curves had the lowest overall rms scatter and the fewest outliers. We first masked out data that was flagged in the k2phot pipeline as a thruster fire event or an outlier in background flux. Periodic transit signals were initially found by flattening the light curve with a Savitsky–Golay filter over a window of 101 points (∼50 hr), then running a box least-squares periodogram, iteratively masking out the higher signal-to-noise transits until there were no more convincing planet signals in the data. This search gave estimates of planet periods, transit times, and transit depths. Next, we made use of the exoplanet ²¹ toolkit to model stellar variability using a Gaussian process with a simple harmonic oscillator kernel, while simultaneously fitting planet transits as described in Foreman-Mackey et al. (2017). In order to simultaneously fit the K2 and Spitzer data (Section 2.4), we used the flattened the light curve by subtracting the Gaussian process stellar model without fitting out the planet transits, which is shown in the middle panel of Figure 2.

Zoom In Zoom Out Reset image size

We observed K2-138 with the Infrared Array Camera (IRAC) on the Spitzer Space Telescope (DDT 13253; PI: J. Christiansen) at the predicted transit time of the putative sixth planet. We used Channel 2 (4.5 μm) since it is less affected by intrapixel sensitivity variations than Channel 1 (3.6 μm). The observation began with a 30 minute pre-observation stare, which was discarded in the analysis but included in the Astronomical Observation Request to allow the telescope and instruments to settle after slewing (Grillmair et al. 2012). To minimize the pixel-phase effect and achieve a pointing accuracy to within ∼0.1 pixel, we conducted pre-observations in peak-up mode using the Pointing Calibration and Reference Sensor (Ingalls et al. 2012).

Observations of K2-138 were conducted between 2018 March 15 and 16 for a total duration of 11 hr centered near the predicted time of mid-transit from the K2 ephemeris of the sixth planetary signal found by Christiansen et al. (2018). Individual frame exposure times were set to 2 s to stay in the linear regime of the detector for this bright target. The subarray mode was used to minimize readout times and data volume. In total, 19,840 individual frames were taken.

We performed centroiding and aperture photometry using photutils (Bradley et al. 2019), fitting a 2D Gaussian to each image. To select the optimal aperture radius, we computed photometry from the centroid positions using fixed radii between 1.5 and 3.0 pixels in 0.1 pixel increments. Background levels were found by the method described by Knutson et al. (2011). This process entails masking out the regions within a radius of 12 pixels from the centroid along with the central two rows and columns, then finding the median background value of the pixels after clipping 3σ outliers.

We modeled systematics in the Spitzer light curves using pixel-level decorrelation (PLD; Deming et al. 2015). PLD has become a premier technique for correcting Spitzer systematics in planet transit analyses (e.g., Beichman et al. 2016; Benneke et al. 2017; Dressing et al. 2018; Feinstein et al. 2019; Livingston et al. 2019; Berardo et al. 2019), and was developed to account for intrapixel sensitivity variations that produce intensity fluctuations in the photometry. In our analysis, we used PLD to model the Spitzer systematics simultaneously with the exoplanet system parameters. The full model is described by

$$\begin{eqnarray}&&S(t)=\displaystyle \frac{{\displaystyle \sum }_{i=1}^{n}{w}_{i}{D}_{i}(t)}{{\displaystyle \sum }_{i=1}^{n}{D}_{i}(t)}+m\cdot t+{M}_{\mathrm{tr}}(\theta ,t),\end{eqnarray} \tag{ 1 }$$

where w_i are individual time-independent pixel weights in the n selected pixels in the region centered on the star, D_i (t) is the observed flux (or counts) in the individual pixels of the selected region for each time step t, m is the slope of a linear temporal ramp, and ${M}_{\mathrm{tr}}(\theta ,t)$ is the transit model with model parameters θ. The first part of this equation normalizes the individual pixel intensities so their sum at each time step is unity. In our analysis we used a 3 × 3 pixel region centered around the brightest pixel.

Using emcee (Foreman-Mackey et al. 2013), we simultaneously model the K2 and Spitzer data, computing the posterior probability distributions for six transiting planets and the Spitzer systematics. To model the transits we used batman (Kreidberg 2015), which solves the analytic equations for an exoplanet transit as derived in Mandel & Agol (2002). We computed posterior probability distributions for the mid-transit times T₀, orbital periods P, the ratios of planet-to-star radii R_p /R_⋆, the scaled semimajor axes a/R_⋆, impact parameters b, two sets of quadratic limb-darkening coefficients q₁ and q₂ (one set for K2 and another for Spitzer), and the nine pixel weights w_i and linear slope m from Equation (1). We performed an autocorrelation analysis ²² to ensure chain convergence. Due to our large set of parameters, we used 500 walkers and 250,000 steps.

The resultant K2 light curve fit using the median values of the posteriors is shown in Figure 2, with the phase-folded light curves shown in Figure 3. Figure 4 shows a clear transit event in the Spitzer data for K2-138 g, confirming the existence of a sixth planet in the K2-138 system. We note that Christiansen et al. (2018) obtained high-resolution AO imaging of the K2-138 system, ruling out nearby stellar companions that could contaminate or mimic a planet signal. Further, since K2-138 is a multiplanet system, it is more likely that additional transit-like signals come from another planet (validation by multiplicity, e.g., Lissauer et al. 2014; Sinukoff et al. 2016). Table 2 lists all the derived planet parameters for the K2 and Spitzer data. We compare the K2 and Spitzer light curves for K2-138 g in Figure 5. The transit durations for the light curves are nearly identical, but the transit depth posterior distributions show a slightly larger radius in the Spitzer data, although the difference is <1σ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2. K2-138 Planet Parameters

Planet	Period	T0	T14	Rp	Rp	a	a	b	i	F	Teq a	TSM
	d	BJD-2457700	hr	R⋆	R⊕	R⋆	au		°	F⊕	K
b	${2.35321}_{-0.00036}^{+0.00037}$	${40.3713}_{-0.0059}^{+0.0061}$	${1.947}_{-0.168}^{+0.150}$	${0.01628}_{-0.00111}^{+0.00116}$	${1.49}_{-0.14}^{+0.15}$	${8.68}_{-0.58}^{+0.60}$	${0.0338}_{-0.0031}^{+0.0034}$	${0.35}_{-0.23}^{+0.21}$	${87.72}_{-1.54}^{+1.52}$	${424.51}_{-98.15}^{+126.30}$	${1157}_{-74}^{+78}$	${2.62}_{-0.98}^{+1.88}$
c	${3.56015}_{-0.00023}^{+0.00022}$	${40.3210}_{-0.0028}^{+0.0029}$	${2.339}_{-0.087}^{+0.089}$	${0.02463}_{-0.00078}^{+0.00091}$	${2.26}_{-0.17}^{+0.18}$	${11.39}_{-0.70}^{+0.53}$	${0.0442}_{-0.0039}^{+0.0039}$	${0.28}_{-0.18}^{+0.17}$	${88.59}_{-0.97}^{+0.92}$	${248.95}_{-55.02}^{+70.82}$	${1012}_{-61}^{+65}$	${25.91}_{-7.14}^{+10.33}$
d	${5.40484}_{-0.00051}^{+0.00049}$	${43.1569}_{-0.0035}^{+0.0037}$	${2.764}_{-0.086}^{+0.097}$	${0.02678}_{-0.00081}^{+0.00093}$	${2.46}_{-0.18}^{+0.19}$	${14.81}_{-0.86}^{+0.64}$	${0.0575}_{-0.0049}^{+0.0050}$	${0.23}_{-0.16}^{+0.17}$	${89.11}_{-0.75}^{+0.62}$	${147.12}_{-32.19}^{+41.30}$	${888}_{-53}^{+57}$	${23.25}_{-6.44}^{+8.74}$
e	${8.26147}_{-0.00053}^{+0.00052}$	${40.6451}_{-0.0021}^{+0.0022}$	${2.967}_{-0.070}^{+0.076}$	${0.03385}_{-0.00094}^{+0.00110}$	${3.11}_{-0.22}^{+0.24}$	${20.40}_{-1.57}^{+1.20}$	${0.0792}_{-0.0077}^{+0.0076}$	${0.38}_{-0.21}^{+0.15}$	${88.92}_{-0.55}^{+0.61}$	${77.79}_{-17.73}^{+23.51}$	${757}_{-47}^{+52}$	${23.27}_{-6.28}^{+8.27}$
f	${12.75758}_{-0.00121}^{+0.00120}$	${38.7018}_{-0.0035}^{+0.0044}$	${3.243}_{-0.111}^{+0.130}$	${0.02982}_{-0.00122}^{+0.00130}$	${2.73}_{-0.21}^{+0.23}$	${27.32}_{-2.11}^{+1.95}$	${0.1064}_{-0.0105}^{+0.0108}$	${0.48}_{-0.19}^{+0.12}$	${88.99}_{-0.35}^{+0.43}$	${43.02}_{-10.03}^{+13.24}$	${653}_{-42}^{+45}$	${23.45}_{-12.83}^{+32.85}$
gK2	${41.96645}_{-0.00665}^{+0.00603}$	${73.8541}_{-0.0044}^{+0.0044}$	${4.035}_{-0.172}^{+0.156}$	${0.03472}_{-0.00177}^{+0.00171}$	${3.18}_{-0.26}^{+0.28}$	${58.59}_{-4.64}^{+4.84}$	${0.2287}_{-0.0231}^{+0.0249}$	${0.73}_{-0.07}^{+0.05}$	${89.29}_{-0.11}^{+0.11}$	${9.29}_{-2.23}^{+2.92}$	${445}_{-30}^{+31}$	${18.08}_{-10.83}^{+29.41}$
gSpitzer b	⋯	${493.5155}_{-0.0065}^{+0.0103}$	${4.054}_{-0.175}^{+0.159}$	${0.03759}_{-0.00230}^{+0.00214}$	${3.44}_{-0.31}^{+0.32}$	⋯	⋯	⋯	⋯	⋯	⋯	⋯

Notes.

^a Equilibrium temperatures were computed assuming a Bond albedo of 0.3. ^b Period, a, b, and i were computed jointly with the K2 data.

Download table as:  ASCIITypeset image

The ratio of orbital periods between successive K2-138 planets are: c:b = 1.513, d:c = 1.518, e:d = 1.529, f:e = 1.544, and g:f = 3.290. In order to determine how close to 3:2 resonance these planets are, we estimated the mean-motion resonance widths using the program from Volk & Malhotra (2020), ²³ which is based on the analytical derivations of resonance widths in the single-planet limit from Murray & Dermott (2000). For this calculation, we used the masses of planets K2-138 b, c, d, and e from Lopez et al. (2019), and estimated masses of K2-138 f and g (${6.72}_{-3.86}^{+8.04}$ and ${8.94}_{-5.91}^{+12.89}{M}_{\oplus }$) from mass–radius relationships (Ning et al. 2018). The results of this calculation, out to fourth-order mean-motion resonances, are shown in Figure 6. Within the upper and lower planet mass limits, K2-138 b, c, d, and e are near (within a few half-widths) their mutual 3:2 resonances at low eccentricity, but the outer pair of planets are not near any low-order resonances.

Zoom In Zoom Out Reset image size

The sizeable gap between K2-138 f and g leads to speculation that there could be additional nontransiting planets in the system. Indeed, Gilbert & Fabrycky (2020) suggest ∼20% of high-multiplicity planet systems host additional planets in the gaps between detected planets. Each consecutive planet pair of K2-138 has period ratios that slip further away from 3:2, and assuming the orbital period ratios continued at 1.544 (f:e), planets could be expected with orbital periods near 19.70, 30.42, and 46.98 days. However, if the K2-138 planets were all in perfect 3:2 resonance with planet b, they would orbit at 17.87, 26.80, and 40.21 days. Without additional data, we are unable to conclude whether or not K2-138 g would be near a 3:2 resonance with a planet in the gap.

Multiplanet systems have been found to be highly coplanar (e.g., Fabrycky et al. 2014; Zhu et al. 2018; Gilbert & Fabrycky 2020), however, the more distant a planet orbits, the closer to 90° inclination it must be to be in a transiting geometry. Assuming orbital periods of 19.70 and 30.42 days, planets around K2-138 would need to be at inclinations above 889 and 892, respectively, for us to observe them in transit. Even within the solar system, the planets are nearly coplanar, yet they still have mutual inclinations between 033 and 63 (Winn & Fabrycky 2015).

We further explore the possibility of planets within the gap between planets f and g using DYNAMITE, ²⁴ which uses population statistics to predict previously undetected planets (Dietrich & Apai 2020). This model takes inputs of stellar parameters (radius, mass, and temperature) and known planet parameters (inclination, radius, and period), and yields probability distributions where the population models predict a planet or planets might exist. We considered four different scenarios as inputs to DYNAMITE, which are shown in Figure 7: (a) all currently known/detected K2-138 planets, (b) the removal of K2-138 c and e, (c) all K2-138 planets with a planet injected at 19.70 days, and (d) all K2-138 planets with a planet injected at 30.42 days. In scenario (a), DYNAMITE predicted a planet or planets to be within the gap between K2-138 f and g, and beyond K2-138 g. The model accurately predicted the locations of K2-138 c and e in scenario (b). From our planet injection tests in scenarios (c) and (d), the models predicted a planet near 30 and 20 days, respectively.

Zoom In Zoom Out Reset image size

Due to its distance from first-order resonance, TTV measurements for K2-138 g would be difficult. Using TTVFaster (Agol & Deck 2016) we computed TTV amplitudes of ${2.23}_{-0.43}^{+0.38}$, ${4.70}_{-1.11}^{+1.13}$, ${8.57}_{-1.46}^{+1.40}$, ${6.92}_{-2.18}^{+3.84}$, ${7.17}_{-1.13}^{+1.23}$, and ${0.09}_{-0.05}^{+0.12}$ minutes for K2-138 b, c, d, e, f, and g, respectively. Our inputs to TTVFaster were the masses of planets K2-138 b, c, d, and e from Lopez et al. (2019) and the aforementioned estimated masses of K2-138 f and g. We also assumed zero eccentricity. Our average six minute (1σ) K2 timing precision was insufficient to measure TTVs for this system, however, higher cadence (one minute) observations with CHEOPS should improve the timing precision enough to allow detection of TTVs of the inner five planets.

In measuring the masses of the inner four K2-138 planets, Lopez et al. (2019) did not identify additional planets, though additional signals might have been absorbed by their Gaussian process to fit out stellar activity at the 5.6 m s⁻¹ level. The mass measurement of K2-138 f was hindered by its orbital period of 12.8 days, near half of the 24.7 day stellar rotation period. The stellar rotation period might also hinder detection of a planet in a orbit near the next 3:2 resonance beyond planet f around 20 days. Future planet searches and mass measurements for this system would likely benefit from simultaneous photometric and RV observations, as Kosiarek & Crossfield (2020) suggest this could enhance the precision of RV measurements. Lopez et al. (2019) were unable to reliably measure a mass of K2-138 g either, but assuming a mass of ${8.94}_{-5.91}^{+12.89}\,{M}_{\oplus }$, we predict an RV semi-amplitude of ${1.79}_{-1.18}^{+2.56}\,{\rm{m}}\,{{\rm{s}}}^{-1}$. If there were planets between f and g of similar masses to the other planets in the system, we would expect them to have RV semi-amplitudes between 1.5 and 2.5 m s⁻¹, which would make them similarly difficult to detect due to stellar activity levels.

We note that the outer five planets of K2-138 are all sub-Neptunes similar in size, and planet b is likely a rocky super-Earth with a density of ${5.01}_{-2.00}^{+2.73}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. Common sizing of multiplanet systems has previously been found for Kepler systems (e.g., Millholland et al. 2017; Wang 2017; Weiss et al. 2018; Gilbert & Fabrycky 2020). From a planet formation standpoint, Adams et al. (2020) found that energy optimization occurs when planets are nearly equal in mass for low-mass (super-Earth/sub-Neptune) planet systems, which is consistent with what we see with K2-138. Though, we note that the outer planets of K2-138 have larger radii than the inner planets, a trend consistent with the findings of Ciardi et al. (2013), Millholland et al. (2017), Kipping (2018), and Weiss et al. (2018), possibly the result of enhanced photoevaporation closer to the star. We plot the planet radii with respect to incident stellar flux for the K2-138 planets, compared to the population of the K2 planets (shown as density contours) from Hardegree-Ullman et al. (2020) in Figure 8. K2-138 b has an incident flux (F_⊕) over 400 times higher than Earth and is the only planet in the system below the planet radius valley. The other planets in the system receive less than 250 F_⊕, apparently low enough to retain an atmosphere.

Zoom In Zoom Out Reset image size

From Figure 8, it appears that many of the K2-138 planets are inflated relative to their counterparts with similar incident stellar flux. If the system was relatively young, we would expect the planets to still be undergoing mass loss. Lopez et al. (2019) computed an age of ${2.3}_{-0.36}^{+0.44}\,\mathrm{Gyr}$ for K2-138 based on chromospheric emission, and ${2.8}_{-1.7}^{+3.8}\,\mathrm{Gyr}$ from their joint RV, light curve, and spectral energy distribution analysis. Similarly, we input photometry, stellar parameters, and a rotation period of 24.7 days into the isochrone fitting with gyrochronology package stardate ²⁵ and compute an age of 2.8 ± 0.3 Gyr, consistent with Lopez et al. (2019). However, a visual assessment of the raw flux in Figure 2 and a Lomb–Scargle periodogram yields a significant peak corresponding to a period of ∼12.5 days. This rotation period corresponds to a younger age closer to ∼0.9 Gyr. We note that even at this younger age, it is unlikely that the planets are still undergoing significant mass loss since this process occurs within the first few hundred megayears (Lopez et al. 2012).

Another possibility for these relatively large planets is tidally induced radius inflation (Millholland 2019). The Kepler mission unveiled a statistical overabundance of planet pairs just outside of first-order mean-motion resonances, specifically 2:1 and 3:2 (Lissauer et al. 2011; Millholland & Laughlin 2019). As noted in Section 3.1, most of the planet pairs of K2-138 fall just outside a 3:2 resonance. Tidal forces from the host star can push planets into near-resonant configurations, but host-star tides alone cannot explain how all the energy from this process is dissipated to keep planets in this configuration. Millholland & Laughlin (2019) showed that obliquity tides may be the source of energy dissipation that helps sculpt these near-resonant systems. Consequently, these tidal forces heat the planet interiors, leading to atmospheric inflation (Millholland 2019). We posit that K2-138 is a strong candidate for planet radius inflation due to obliquity tides.

To date, there have only been nine other exoplanet systems with six or more confirmed planets, ²⁶ including RV discovered systems HD 10180 (six planets), HD 219134 (six planets), and HD 34445 (six planets), and transiting systems Kepler-11 (six planets), Kepler-20 (six planets), Kepler-80 (six planets), Kepler-90/KOI-351 (eight planets), TRAPPIST-1 (seven planets), and TOI-178 (six planets). Perhaps most similar to K2-138, however, is the HD 158259 system, with five confirmed planets and a sixth candidate outer planet (Hara et al. 2020), all near 3:2 orbital mean-motion resonances. Four of the confirmed planets were detected in RV data with the SOPHIE spectrograph, and the innermost planet was found to be transiting in TESS data. The outermost candidate planet orbits every 17.4 days, close to the stellar rotation period, complicating confirmation of this planet. The five innermost planets of K2-138 and HD 158259 are each located at nearly identical distances to their host stars, as shown in Figure 9. We estimated HD 158259 planet radii for the five nontransiting planets using the planet masses from Hara et al. (2020) and the mass–radius relationships of Chen & Kipping (2017). These nontransiting planets are also all similar-sized sub-Neptunes, with estimated radii larger than 2 R_⊕—again consistent with the aforementioned common sizing of multiplanet systems. Each respective planet in HD 158259 is slightly smaller than its counterpart K2-138 planet, which could be the result of HD 158259 being a larger host star (1.08 ± 0.10 M_⊙) that is more efficient at stripping away planetary atmospheres by intense irradiation (Ehrenreich et al. 2015). We note, however, that there are significant uncertainties in planetary mass–radius relationships. Without transit data, it is difficult to test whether or not this system undergoes tidal radius inflation as mentioned in Section 3.2.

Zoom In Zoom Out Reset image size

We qualitatively compared the orbital spacing (a/R_⋆) of these high-multiplicity systems with transiting planets (Figure 10). In addition to the K2-138 system, there is a sizeable gap between the outermost detected transiting planets of the Kepler-11, Kepler-20, and Kepler-80 systems. HD 219134 has two transiting planets and four nontransiting planets detected via RV measurements, again with a large gap between the two outermost planets. This large outermost planet gap is also present in the RV system HD 34445. Notably, a nontransiting planet was identified in the gap between outer planets Kepler-20 f and d with RV data (Buchhave et al. 2016). As noted in Section 3.1, planets orbiting further out must be closer to i = 90° to be in a transiting geometry, but RV and TTV data may uncover unseen planets. We encourage further investigations of this outer planet gap feature in high-multiplicity planet systems in order to disambiguate whether it is caused by observational biases or planet formation processes.

Zoom In Zoom Out Reset image size

We computed the transmission spectroscopy metric (TSM) for the K2-138 planets as defined by Equations (1), (2), and (3) of Kempton et al. (2018). The TSM is the expected signal-to-noise for a 10 hr observing program with the James Webb Space Telescope (JWST)/Near Infrared Imager and Slitless Spectrograph (NIRISS). For planets b, c, d, and e, we used the planet masses measured by Lopez et al. (2019), and for planets f and g, we used our estimated masses. The equilibrium temperature was calculated assuming zero albedo and full day–night heat redistribution. The resultant TSM values are listed in Table 2, and are ∼20 for the outer planets, falling to 2.62 for the innermost planet b. These values are well below the recommended threshold of TSM >90 for high-quality atmospheric characterization of sub-Neptune-sized planets. For now, the K2-138 planets are unlikely to be selected as high-priority targets for JWST observations.

The European Space Agency Atmospheric Remote-sensing Infrared Exoplanet Large-survey (ARIEL) space mission aims to gather transmission spectra of 1000 exoplanets during its four-year mission (expected to launch in 2028) in order to study their composition, formation, and evolution. Edwards et al. (2019) compiled a list of potential targets for ARIEL, taking into account currently known stellar and planet parameters. K2-138 falls very near the average star system considered for this target list. The inner five planets of K2-138 also fall within the range of planets considered for the potential target list, but very few planets with orbital periods beyond ∼20 days will likely be considered, all but ruling out observations of K2-138 g. However, since K2-138 contains five similarly sized sub-Neptunes with a ∼500 K range of equilibrium temperatures from warm to temperate, these planets might provide a unique test bed for comparative sub-Neptune atmosphere studies.

We have confirmed the existence of K2-138 g, solidifying K2-138 as the largest K2 multiplanet system. K2-138 g breaks the continuous near 3:2 mean-motion resonance of the inner five planets, but the sizeable gap between K2-138 f and g hints at the possibility there could be additional nontransiting planets in this system. We encourage future observations of this potential key benchmark system to (1) constrain TTVs of the inner planets, (2) enable more precise masses and potential discovery of additional planets with simultaneous photometric and RV measurements, and (3) facilitate comparative atmospheric studies of warm to temperate sub-Neptune planets.