Confirming the 3:2 Resonance Chain of K2-138

Since we run our simulations to explore a large range of mass and orbital properties, we analyze our results to see which planetary parameters lead to the libration of the resonance angles. We first compare the distributions of planetary mass, eccentricity, and orbital period for simulations where each angle is librating to distributions of the same parameter from simulations where each angle is not librating using a two-sample Kolmogorov–Smirnov test. In this test, the null hypothesis is that the two samples (parameter from librating simulations and parameter from circulating simulations) are drawn from the same population, and the resulting p-value is the probability that the null hypothesis is correct. Therefore, a small p-value ($p<0.05$) allows us to reject this hypothesis and suggests that the two distributions are statistically distinct.

We find that the libration of the three-body angle between the inner three planets depends on the eccentricities and orbital periods of the three planets but only depends on the mass of K2-138b. Similarly, the libration of the three-body angle between the outer three planets depends on the eccentricities of K2-138d and f, the orbital periods of the three planets, and the mass of K2-138f. The three-body angle between the middle planets and the four two-body angles also depend on some mixture of masses, eccentricities, and orbital periods. We show Kernel Density Estimations of the distributions of some of the more interesting dependencies in Figure 2 and report the p-values resulting from our K-S tests involving mass and orbital period in Table 3.

For each angle and planet property pair with a small K-S p-value, we report the median and 68.3% confidence interval of the parameter in the simulations where the angle is librating in Table 3. To further quantify which planet properties are statistically distinct between our simulations with librating resonance angles and the initial distributions, we perform a T-test on each parameter with a small K-S p-value. The T-test null hypothesis states that there is no statistical difference between two groups, and a small p-value indicates that an observed difference is not due to chance. We find that the masses of K2-138b, K2-138c, and K2-138f are statistically distinct between our simulations with librating resonance angles and the estimates from Lopez et al. (2019). Our results suggest that these three planets are slightly more massive than the radial velocity data can constrain, with masses of $3.46^{+1.07}_{-1.01}$, 6.53${}^{+1.17}_{-1.13}$, and 2.65${}^{+1.90}_{-1.62}~{}M_{\oplus}$, respectively. If we constrain the mass of K2-138f using the two-body angle between it and K2-138e, we would recover an even larger mass of 3.01${}^{+2.58}_{-1.98}~{}M_{\oplus}$, but the libration of this angle also depends on the orbital period of K2-138e.

The periods of all five planets shift slightly within the beginning of our simulations, even when they do not lock into resonance, and this shift results in a distribution of final orbital periods that is statistically distinct from those measured by Lopez et al. (2019) for all planets except K2-138f. Because of this, we compare the orbital periods of the planets between simulations with librating resonance angles and simulations without libration via the T-test. We recover large p-values for all orbital periods except K2-138e, suggesting that our small K-S p-values are due to chance. For the two-body angle between K2-138e and K2-138f to librate, K2-138e requires a slightly larger period of $8.262\pm 0.002$ days than that measured by Lopez et al. (2019) ($8.2615\pm 0.0002$ days).

The differences in the masses of K2-138b, K2-138c, and K2-138f and in the orbital period of K2-138e could explain why we do not recover a five-planet resonance chain in all of our simulations. We discuss potential reasons why K2-138f is dynamically decoupled in the majority of our simulations below in Section 4.2.

Ideally, we would constrain the masses and orbital parameters of the planets by fitting the TTVs of the planets or by photodynamically fitting the system. The TTVs, although estimated to be on the order of 2–5min (Christiansen et al., 2018; Lopez et al., 2019), have so-far been illusive, indicating that the cadence or precision of the data is insufficient to tease out the perturbations. We summarize our photodynamic fitting efforts in the Appendix.

We hypothesize three reasons why K2-138f is not part of the resonance chain in the majority of our simulations:

1. 1.

   The planet is in resonance, but its orbital parameters, mass, and/or the masses or orbits of the other planets fall in the part of parameter space that is narrower than the data currently constrain (H1)

2. 2.

   The planet was once in resonance but this resonance has since been broken or disrupted (H2)

3. 3.

   The planet is not and was never in resonance (H3)

We explore each of these hypotheses below.

While K2-138’s resonance chain is interesting, we now question how the resonances formed. Although long-scale migration — forming the planets more widely spaced and further from their star than observed — often results in resonance pairs (e.g., Snellgrove et al., 2001; Papaloizou & Terquem, 2005; Rein & Tamayo, 2015) and resonance chains (Cossou et al., 2013) and is often used as the explanation of such chains (e.g., Kepler-223, Mills & Fabrycky, 2017), resonance chains can potentially form in situ²²2Here, we use the term “in situ” to distinguish the evolution after the giant impact phase from long-scale migration (e.g., Lee & Peale, 2002) and to align with other works exploring the in situ formation of super-Earths (e.g., Dawson et al., 2015). The local growth of these planets does not require that the planetesimals that form them be local, as they could have also accumulated by radial drift. , with only small changes to their semi-major axes. MacDonald & Dawson (2018) explored three different pathways for forming resonance chains, including long-scale migration and two pathways consistent with in situ formation: short-scale migration (where the planets form outside of resonance and small shifts to their orbital periods lock them into resonance) and eccentricity damping. In addition, Morrison et al. (2020) find that close-in super-Earths and mini-Neptunes, such as those in K2-138, can lock into resonance chains due to dissipation from a depleted gas disk and maintain resonance once the gas disk is fully dissipated. Since all currently confirmed resonance chains are consistent with both in situ formation and long-scale migration (MacDonald & Dawson, 2018), we test all three different pathways for forming the resonance chain of K2-138.

Following the methods of MacDonald & Dawson (2018), we simulate the formation of this resonance chain via long-scale migration, short-scale migration, and eccentricity damping only. Here, long-scale migration (e.g., Mills & Fabrycky, 2017) assumes that the planets form far from their star and migrate inwards until they reach their currently observed semi-major axes; short-scale migration (e.g., MacDonald et al., 2016) assumes that the planets do not undergo significant changes to their semi-major axes during or after the giant impact phase; and eccentricity damping only (e.g., Dong & Dawson, 2016) assumes constant angular momentum. For these three pathways, we apply an inward migration force and/or eccentricity damping forces on timescales $\tau_{a}\sim 10^{4}$–$10^{6}$ and $\tau_{e}\sim 10^{3}$–$10^{5}$ years, respectively, following the prescription in Papaloizou & Larwood (2000). We draw these timescales from independent log-uniform distributions (MacDonald et al., 2021). We apply damping forces to only the outer planet³³3We do not know the migration rates for the planets since they depend on the conditions of the disk, so we implicitly assume that the migration of the inner planets is on a much longer timescale than that of the outer planet. for long-scale and short-scale migration, and damp the eccentricities of all planets in the eccentricity damping only simulations. We start all simulations with all planets wide of their observed commensurability and with no librating resonance angles. For each formation pathway, we run 500 simulations, damping the semi-major axes and eccentricities where applicable using the modify_orbits_forces routine in the REBOUNDx library. We use the stellar and planetary properties as defined in Table 1.

We find that we are able to form this resonance chain via all three pathways described above, as each suite of simulations contains numerous sets of initial conditions that lead to fully librating resonance chains. We summarize the resulting centers and amplitudes for the three-body resonance angles in Table 5. Since both short-scale migration and eccentricity damping are consistent with in situ formation, we find that the resonance chain of K2-138 could have formed in situ. We show an example simulation from the short-scale migration suite in Figure 3. We caution against using the values in Table 5 to draw any stronger conclusions, as this study is only able to say what is possible and not what is more likely.

The compositions of planets can also be used to constrain the formation and evolution of a system. As a first-order approximation of the compositions of K2-138’s planets, we first explore the planet’s bulk densities, comparing them to Earth-like and less dense compositions. We draw 1000 pairs of mass and radius estimates for the confirmed planets from normal distributions of the parameters in Lopez et al. (2019), while obeying the 99% credibility intervals of density. In Figure 4, we plot 200 of these samples with mass-radius curves for compositions of pure water, Earth-like, and 1% H/He envelopes. We see that K2-138b is consistent with a terrestrial composition while K2-138c, K2-138d, and K2-138e are less dense and require large volatile layers. K2-138f also has a low density, below 2.068 $\mathrm{g\,cm^{-3}}$, and most likely requires the largest atmosphere envelope.

We model the interiors of the four inner planets (b-e) using the planet structure code MAGRATHEA⁴⁴4We used a preliminary version; the final version will be available at https://github.com/Huang-CL/Magrathea. (Huang et al. submitted, 2022), which calculates the pressure, density, temperature, and radius of a spherically symmetric planet with defined mass in each differentiated layer in 1-D using the fourth order Runge-Kutta method. We assume a surface pressure of one bar and use MAGRATHEA’s default model. The default model uses equations of state for solid (hexagonal-close-packed⁵⁵5Under the pressure and temperature of Earth’s inner core, the iron most likely takes this structure (Vočadlo et al., 1999)) and liquid iron in the core, bridgmanite and post-perovskite silicate in the mantle, and water, ice-VII, ice-VII’, and ice-X in the hydrosphere (Oganov & Ono, 2004; Sakai et al., 2016; Dorogokupets et al., 2017; Smith et al., 2018; Grande et al., 2019). We model the atmosphere as an ideal gas with a mean molecular weight of 3 $\mathrm{g~{}cm^{-3}}$, similar to a hydrogen-helium mixture.

To mitigate the degeneracy between water mass fraction and atmospheric mass fraction with core mass fraction, we separate our analysis of each planet into two suites of 1000 models: one where we explore the water mass fraction and one where we explore the atmospheric mass fraction. For our water analysis, we limit the hydrosphere to liquid and ice phases, although the planets’ equilibrium temperatures would suggest a vapor layer. We assume an isentropic temperature profile with a surface temperature of 300 K (e.g., Hakim et al., 2018). For our atmosphere analysis, we assume an ideal gas using an isentropic temperature profile and set the surface temperatures equal to the equilibrium temperatures derived from Lopez et al. (2019). A real gas would require less atmosphere mass, while a less steep temperature profile would require more. It is important to note that models such as ours are inherently limited for large, hot planets with interior conditions beyond our experimentally-determined equations of state, but although limited, our models are efficient at exploring the range of possible interior solutions.

For each of the 1000 samples of mass and radius, we use a secant method and vary the mass percentage in each layer until the simulated radius matches the sample to 0.01%. We calculate the water mass fraction (WMF) uniformly across the range of core mass fractions while fixing the core-to-mantle mass ratio. For the atmosphere mass fraction (AMF), we fix the core-to-mantle mass ratio to 1:2 similar to Earth and calculate the AMF uniformly across WMF.

We show the resulting WMF in Figure 5. We find that K2-138b most likely has a WMF between 9.0 and 47%, depending on the core-to-mantle mass ratio, and an estimated water fraction of 24.3${}^{+39.0}_{-22.0}$% with an Earth-like core-to-mantle mass ratio. 32% of the samples result in hydrosphere-free solutions, meaning the planet could be composed of only core and mantle and no water. With Earth-like core-to-mantle mass ratios, only 18.3% of the K2-138c samples and 25.5% of the K2-138d samples have three-layer solutions with less than 90% WMF, suggesting appreciable atmospheres. For reference, models of Neptune suggest at least 80% mass in a water-dominated fluid layer (Scheibe et al., 2019).

In Figure 6, we show the AMF needed to match samples of the planets’ masses and radii across WMF. Across all possible water mass fractions, we predict AMF of the four planets of: 1.7${}^{+9.3}_{-0.9}\times 10^{-3}$%, 5.3${}^{+7.9}_{-3.8}\times 10^{-4}$%, 7.0${}^{+13}_{-5.6}\times 10^{-4}$%, 0.022${}^{+0.016}_{-0.011}$%. K2-138b’s mass and radius could be explained with a hot, inflated H/He atmosphere layer, but many solutions require atmospheric masses of less than 10^-6 M_⊕, 1.25% of Venus’s atmospheric mass. With WMF = 0, K2-138c and K2-138d require an AMF of over 0.001%. The pressure and temperature under the atmosphere of K2-138c with 50% WMF is around 10 bar and 4000 K which, in our model, creates a small layer of liquid water before transitioning to high-pressure ices in the model. This temperature and pressure suggest that the water would be gaseous, but understanding this boundary requires a model that couples the atmosphere and interior. K2-138e requires an AMF around 15-30 times more than K2-138c and d at zero WMF.

Aside from the similar inferred compositions of K2-138c and K2-138d, the possible compositions of the planets of K2-138 have little overlap. Their densities decrease and inferred volatile content increases with orbital period. Other systems with resonances, such as TRAPPIST-1, and other compact multi-planet systems have inter-planetary similarities in their sizes and masses and therefore in their inferred compositions (Weiss et al., 2018; Millholland & Winn, 2021; Agol et al., 2021). This similar sizing within systems, especially when paired with the regular orbital spacing these systems also exhibit, could be telling of the formation history and/or the subsequent dynamical evolution (e.g., Adams, 2019; MacDonald et al., 2020; Mishra et al., 2021), and K2-138’s lack of intra-system uniformity could be just as telling.