A Compact Multi-planet System with a Significantly Misaligned Ultra Short Period Planet

Since the failure of the second reaction wheel, the Kepler spacecraft has been re-purposed to observe a set of fields along the ecliptic. Each K2 campaign lasts ∼80 days (Howell et al. 2014), achieving similar precision to the original Kepler mission (Vanderburg et al. 2016a). K2-266 was observed during K2 Campaign 14 from UT 2017 June 02 until UT 2017 August 19, obtaining 3504 observations on a 30 minute cadence (see Figure 1). Following the strategy described in Vanderburg & Johnson (2014) and Vanderburg et al. (2016a), the light curves were extracted from the Kepler pipeline calibrated target pixel files from the Mikulski Archive for Space Telescopes,¹⁹ corrected for the K2 spacecraft-motion-induced systematics, and searched for transiting planet candidates. From our search of K2-266, we identified three super-Earth/sub-Neptune-sized transiting exoplanet candidates with periods of 0.66, 14.7, and 19.5 days, with signal-to-noise (S/N) values of 13.0, 114.6, and 111.5. In addition, some of us (M.H.K., M.O., H.M.S., and I.T.) performed a visual inspection of the light curve using the LCTOOLS²⁰ software (Kipping et al. 2015). From this visual inspection, we identified two additional Earth-sized exoplanet candidates with respective periods of 6.1 and 7.8 days and S/N values of 8.3 and 10.6. An additional visual inspection of the K2 light curve led to the identification of a sixth planet candidate at 56.7 days with an S/N value 6.6. The phase-folded light curves for each planet candidate is shown in Figure 2. We note that the two transits of this candidate overlap with other candidates in the system. The K2 light curve was reprocessed where all six planets were simultaneously fit along with the stellar variability and known K2 systematics. The corresponding light curve was flattened by dividing out the best-fit stellar variability using a spline fit with breakpoints every 0.75 days. The final light curve for K2-266, shown in Figure 1, has a 30 minute cadence noise level of 70 ppm, and a 6 hour photometric precision of 19 ppm.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Using the Tillinghast Reflector Echelle Spectrograph (TRES; Fűrész 2008)²¹ on the 1.5 m Tillinghast Reflector at the Fred L. Whipple Observatory (FLWO) on Mt. Hopkins, AZ we obtained eight observations of K2-266 between UT 2017 November 23 and UT 2018 April 10. TRES has a resolving power of λ/Δλ = 44000 and an instrumental radial velocity (RV) stability of 10–15 m s⁻¹. The spectra were optimally extracted, wavelength calibrated, and cross-correlated to derive relative RVs following the techniques described in Buchhave et al. (2010). We cross-correlate each spectrum, order by order, against the strongest observed spectrum, and fit the peak of the cross-correlation function summed across all orders to derive the relative RVs. Uncertainties are determined from the scatter between orders for each spectrum. We use RV standard stars to track the instrumental zero point over time, and apply these zero point shifts (typically < 15 m s⁻¹) to the relative RVs and propagate uncertainties in the zero point shifts to the RVs. This is why the strongest spectrum, correlated against itself, does not have an RV of 0 m s⁻¹. The final relative RVs are given in Table 2. Using the RV standards to set the absolute zero point of the TRES system, we also determine the RV of K2-266 on the IAU standard system to be 10.848 ± 0.066 km s⁻¹, where the uncertainty is dominated by the uncertainty in the shift from relative to absolute RV.

We refined the characterization of K2-266 by acquiring near-infrared spectra using TripleSpec on the 200'' Palomar Hale telescope on 2017 December 1. TripleSpec has a fixed slit of 1'' × 30'' slit, enabling simultaneous observations across the J, H, and K bands (1.0–2.4 microns) at a spectral resolution of 2500–2700 (Herter et al. 2008). Following Muirhead et al. (2014), we obtained our observations using a four-position ABCD pattern to reduce the influence of bad pixels on our resulting spectra. As in Dressing et al. (2017), we reduced our data using a version of the publicly available Spextool pipeline (Cushing et al. 2004) that was modified for use with TripleSpec data (available upon request from M. Cushing). We removed telluric contamination by observing an A0V star at a similar airmass and processing both our observations for both the A0V star and K2-266 with the xtellcor telluric correction package (Vacca et al. 2003).

After reducing the spectra, we estimated stellar properties by applying empirical relations developed by Newton et al. (2014, 2015) and Mann et al. (2013b, 2013a). Specifically, we estimated the stellar effective temperature and radius by measuring the widths of Al and Mg features using the publicly available, IDL-based tellrv and nirew packages (Newton et al. 2014, 2015). We then employed the stellar effective temperature–mass relation developed by Mann et al. (2013a) to infer the stellar mass from the resulting stellar effective temperature estimate. We also estimated the stellar metallicity ([M/H] and [Fe/H]) using the relations developed by Mann et al. (2013b). For more details about our TripleSpec analysis methods, see Dressing et al. (2017) and C. D. Dressing et al. (2018, in preparation).

The resulting stellar properties were T_eff = 4192 ± 77 K, ${M}_{\star }={0.67}_{-0.07}^{+0.08}\,{M}_{\odot }$, and R_⋆ = 0.63 ± 0.03. These values for the stellar mass are consistent with those estimated from our EXOFASTv2 analysis (see Table 3). However, the radius is ∼3σ different from the EXOFASTv2 fit using the broadband photometry and Gaia DR2 parallax.

To check for nearby stars (either physically associated companions or coincidental alignments) that may influence our results, we examined archival observations from National Geographic Society Palomar Observatory Sky Survey (NGS POSS) from 1952. The proper motion of K2-266 is μ_α = 56.9 mas and μ_δ = −68.8 mas, and has moved ∼6'' in the 66 years since the original POSS observations were taken. The present-day position of K2-266 is located right at the edge of the saturated point spread function of K2-266 in the original POSS plates. While the present-day position of K2-266 is not completely resolved in the POSS image, if there was a sufficiently bright background star at the present-day position of K2-266, we would expect to see some elongation of the POSS point spread function at that position. We see no evidence for such an elongation in POSS plates with either a red-sensitive or blue-sensitive emulsion. We estimate that we can rule out background stars at the present-day position of K2-266 down to a magnitude of about 19 in blue, and a magnitude of about 18 in the red. Figure 3 shows our archival imaging overlaid with the K2 photometric aperture used to extract the light curves.

Zoom In Zoom Out Reset image size

We used modern imaging from the Pan-STARRS data release to search for faint companions at distances greater than a few arcseconds from K2-266 (Flewelling et al. 2016). In the Pan-STARRS images, we identified one star located inside our best photometric aperture about nine magnitudes fainter than K2-266. In principle, if this star were a fully eclipsing binary (with 100% deep eclipses²² ), it could contribute a transit-like signal to the light curve of K2-266 with a depth of at most about 250 ppm. This is shallower than the transits of the two sub-Neptunes, but could in principle contribute the transits of the other four candidates. We therefore extracted the K2 light curve from a smaller aperture (shown in Figure 3 as a navy blue outline overlaid on the Pan-STARRS image of K2-266), which excludes the companion star detected in Pan-STARRS imaging. In the noisier light curve extracted from the smaller aperture, we find that the transits of the USP and 7.8 day planets are convincingly detected, but the transits of the two weaker candidates (at 6.1 and 56.7 days) do not convincingly appear (due to the increased noise in the light curve). We therefore cannot rule out a blended background eclipsing binary origin for at least one of those candidates.

We obtained high-resolution images of K2-266 using the Near Infrared Camera 2 (NIRC2) on the W. M. Keck Observatory. Two observations of each target were taken on UT 2017 December 28, one in the Br-γ filter and the other in the J-band (see Figure 3). NIRC2 has 9.942 mas pix⁻¹ pixel scale and 1024 × 1024 pixel array. The lower left quadrant of the array suffers from higher noise levels. To exclude this part of the detector, a three-point dither pattern was used. The final image shown in Figure 3 is created by shifting and co-adding the observations, after flat-fielding and sky subtraction. We see no other star in the 10'' field of view for K2-266. Our sensitivity to nearby companions is determined by injecting a simulated source with an S/N of 5. The final 5σ sensitivity curves as a function of spatial separation and the corresponding images for K2-266 in both Br-γ and J filters are shown in Figure 3.

We attempted to statistically validate each of the six candidates in K2-266, a process in which the probability of planethood is estimated. If the probability is above some threshold value, the candidate is upgraded to a validated planet. Our method of validation followed the approach taken by Mayo et al. (2018). In detail, we made use of vespa (Morton 2015), a Python package based on the work of Morton (2012). The vespa package calculates the false positive probability (FPP) of an exoplanet candidate by first simulating a population of synthetic stellar systems, each of which creates a transit signal due to a planet or eclipsing binary scenario. Then, vespa calculates the FPP by determining which synthetic systems are consistent with the input observations and calculating the fraction of those systems that correspond to an eclipsing binary scenario.

This determination is made based on inputs such as the sky position of the target, the transit signal, various stellar parameters, and contrast curves from any available high-resolution imaging. (A contrast curve relates the angular separation between the target star and an undetected companion to the maximum brightness for the putative companion.) In the case of K2-266, we provided as input to vespa the R.A. and decl., the phase-folded light curve of the candidate in question (with transits from other candidates removed), J, H, and K_S bandpass stellar magnitudes from 2MASS (Cutri et al. 2003; Skrutskie et al. 2006) and the Kepler magnitude, stellar parameters (${T}_{\mathrm{eff}}$, $\mathrm{log}g$, and $[\mathrm{Fe}/{\rm{H}}]$) calculated in Section 3, and contrast curves from two AO images.

After we subjected each of our six candidates to validation, we made two additional adjustments to their FPP estimates. First, there are eight spectra and corresponding RV measurements collected with TRES from 2017 November 23 to 2018 April 10. The RV measurements derived from the TRES spectra did not display any large variations indicative of a simple eclipsing binary, so we were able to eliminate that scenario. (Note that this is different from a background eclipsing binary or hierarchical eclipsing binary scenario, which we also consider.) By eliminating the possibility of a simple eclipsing binary, the probability of the planet scenario (and each false positive scenario) was increased such that the total probability remained at unity.

Second, according to Lissauer et al. (2012), the likelihood of one or more false positives decreases significantly when there is more than one candidate in a system. In the case of a system with more than two candidates, they estimate that a multiplicity boost factor of 50 is appropriate. As a result, we decreased the FPP for each candidate by a factor of 50.

After calculating FPP values for our six candidates, reducing the eclipsing binary scenario to zero probability, and including a multiplicity boost of 50, we found respective final FPP values of 3.02e − 05, 7.34e − 06, 9.40e − 06, 6.80e − 11, 1.16e − 12, and 4.90e − 06 for candidates K2-266.01, .02, .03, .04, .05, and .06. These values would each be low enough to easily validate all six candidates (e.g., Mayo et al. (2018) used an FPP threshold value of 1e − 4). However, given that we cannot rule out the possibility that the faint background star we identified in Section 2.4 is an eclipsing binary, we were only able to conclusively validate candidates K2-266.01, .03, .04, and .05, naming them K2-266 b, c, d, and e, respectively. We also refrain from validating candidates .02 and .06 because they have the lowest signal-to-noise ratios that do not pass our threshold (8.3 and 6.6, respectively). Validating such low-S/N candidates is challenging because it is difficult to prove that the weakest signals detected in Kepler or K2 data are astrophysical, and not the result of residual instrumental systematics or artifacts (Mullally et al. 2018).

The transit timing measurements for planets d and e listed in Table 3 can be used to derive dynamical constraints on their masses and orbits. In this section, we invert the planet pairs' TTVs to infer their masses, and in combination with the planet radii derived from out light-curve fitting, their densities. We model the planets' TTVs using the TTVFast code developed by Deck et al. (2014) and use the emcee package's (Foreman-Mackey et al. 2013) ensemble sampler, based on the algorithm of Goodman & Weare (2010), to sample the posterior distribution of the planetary masses and orbital elements. We model only the dynamical interactions of planets d and e, and ignore any perturbations from the other (potential) members of the system.²³ We assume planets d and e orbit in the same plane because small mutual inclinations have negligible influence on TTVs. We approximate the mid-transit uncertainties to be Gaussian-distributed about the median transit time determined by EXOFASTv2, with variances set to the larger of the two asymmetric error bars in Table 5. Our likelihood function is then computed based on the standard chi-squared statistic as $\mathrm{ln}{ \mathcal L }=-\tfrac{1}{2}{\chi }^{2}$. We impose a Gaussian prior with 0 mean and a variance of σ = 0.05 on the eccentricity vector components, ${e}_{i}\cos {\omega }_{i}$ and ${e}_{i}\sin {\omega }_{i}$, typical for eccentricities of multi-planet, sub-Neptune systems (Hadden & Lithwick 2014, 2017; Van Eylen & Albrecht 2015; Xie et al. 2016) We found that the inferred planet masses are insensitive to the assumed eccentricity priors after running additional MCMC fits with σ = 0.025 andσ = 0.1. We initialize our MCMC with 200 walkers and run for 120,000 iterations, saving all walker positions every 1000 iterations.

Figure 7 shows the observed timing variations of planets d and e, along with N-body solutions drawn from our MCMC posterior. From our TTV dynamical fit, we determine planet-star mass ratios of ${m}_{d}/{M}_{* }={3.9}_{-1.7}^{+2.5}\times {10}^{-5}$ for planet d and ${m}_{e}/{M}_{* }\,={6.3}_{-2.2}^{+2.8}\times {10}^{-5}$ for planet e. The TTVs yield no strong constraint on planet eccentricities and the posterior distributions largely mirror our assumed priors. We convert the dynamical constraints on planet-star mass ratios to constraints on the planetary masses and densities by combining the posterior samples from our TTV fit with posterior samples of our fit to stellar mass and planet radii computed with EXOFASTv2. The resulting posterior distributions of the planets' masses and densities are plotted in Figure 8. The inferred median and 1σ planet mass values are ${m}_{d}\,={8.4}_{-3.6}^{+5.4}\,{M}_{\oplus }$ and ${m}_{e}={13.6}_{-4.7}^{+6.1}\,{M}_{\oplus }$, and the densities are ${\rho }_{d}\,={2.7}_{-1.2}^{+1.8}\,{\rm{g}}/{\mathrm{cm}}^{3}$ and ${\rho }_{e}={5.6}_{-2.0}^{+2.6}\,{\rm{g}}/{\mathrm{cm}}^{3}$.

Our N-body dynamical model contains 10 free parameters, which are fit to only nine data points. This means that, at face value, our model is underconstrained and we are at risk of overfitting. To understand the origin of our dynamical mass constraints and ensure that they are not merely the result of overfitting, we analyze the TTVs using the analytic model of Hadden & Lithwick (2016). This analytic treatment reduces the dimensionality of the TTV model so that it is no longer underconstrained. Note that we adopt the masses derived from the more complete N-body model as our best-fit values; we present the analytic model simply as a consistency check to ensure that the N-body results are not overfitting the data because of poor MCMC convergence.

We use the formulas of Hadden & Lithwick (2016) to construct an analytic model for the TTVs of planets d and e as a function of planet periods P_i, initial transit times T_i, planet–star mass ratios μ_i, and the complex "combined eccentricity"

$$\begin{eqnarray}&&{ \mathcal Z }\approx \displaystyle \frac{1}{\sqrt{2}}({e}_{e}{e}^{i{\omega }_{e}}-{e}_{d}{e}^{i{\omega }_{d}})\,.\end{eqnarray} \tag{ 1 }$$

The analytic model reduces the total number of model parameters to eight by combining the planets' eccentricities and longitudes of perihelia into the single complex quantity, ${ \mathcal Z }$. The nth transit of planet d and e are modeled as

$$\begin{eqnarray}&&{t}_{i}(n)={T}_{i}+{{nP}}_{i}+\delta {t}_{{ \mathcal C },i}(n)+\delta {t}_{{ \mathcal F },i}(n)\end{eqnarray} \tag{ 2 }$$

where $\delta {t}_{{ \mathcal C }}\propto \mu ^{\prime}$, $\delta {t}_{{ \mathcal F }}\propto \mu ^{\prime} { \mathcal Z }$ with μ' being the perturbing planet's planet-star mass ratio. Expressions for $\delta {t}_{{ \mathcal C }}$ and $\delta {t}_{{ \mathcal F }}$ are given in Hadden & Lithwick (2016). We use the Levenberg-Marquardt minimization algorithm to fit our analytic model to the observed transit times and estimate uncertainties from the local curvature of the χ² surface (e.g., Press et al. 1992). The analytic fit, plotted in Figure 7, yields masses m_d/M_* = 3.0 ± 0.8 × 10⁻⁵ for planet d and m_e/M_* = 4.6 ± 1.0 × 10⁻⁵ for planet e, which are consistent with the N-body MCMC constraints.

The origins of the mass constraint can be qualitatively understood from the analytic model as follows: at conjunction, planets impart impulsive kicks to one another that change their instantaneous orbital periods. This effect is captured by the so-called "chopping" terms, $\delta {t}_{{ \mathcal C }}$, in Equation (2) (see also Nesvorný & Vokrouhlický 2014; Deck & Agol 2015). Indeed, we obtain nearly identical mass constraints from an analytic fit that drops the $\delta {t}_{{ \mathcal F }}$ terms from Equation (2) (and thereby further reduces the number of free parameters to six as the model no longer depends on the complex number ${ \mathcal Z }$). Because these $\delta {t}_{{ \mathcal F }}$ terms vary over a timescale much longer than the baseline of our observations, they are well-approximated by a linear trend and essentially degenerate with small changes to the T_i + nP_i terms in Equation (2). Thus, we have identified the the origin of our mass constraints with the measurement of the chopping signals, $\delta {t}_{{ \mathcal C },i}$, in the TTVs of planets d and e. Over the course of our observing baseline, planets d and e experience a single conjunction, at the time marked by a dashed line in Figure 7. The power of the TTV signal for constraining the planets' masses comes mainly from the impulsive changes in the planets' osculating orbital periods experienced at this conjunction causing the planets to arrive early (in the case of e) or late (in the case of d) at their next transits.

Next, we consider the dynamical stability of the system, along with the probability that all of the putative planets can be seen in transit. Although most Kepler multi-planet systems tend to have low mutual inclinations, this system is unique to date because there is a significant mutual inclination between the innermost planet and the other five. In the context of the known set of multi-planet systems, this appears significant. Ballard & Johnson (2016) found that the Kepler planet-hosting systems around cool stars appears to be drawn from two populations: one set of multi-transit systems and a second set of single-transit planets, which may also have statistically higher obliquities (Morton & Winn 2014) (the concept of these two populations is commonly called the "Kepler dichotomy"). One solution to the Kepler dichotomy is that the two populations are actually all multi-planet systems, and that the ones that appear to be singly-transiting are really systems with larger mutual inclinations, which can see seen as single-transit systems from a particular line of sight. Although most Kepler multi-planet systems (those with four or more planets) are fairly tightly confined to a roughly coplanar region, there is some precedent for multi-planet systems: Mills & Fabrycky (2017) found a two-planet system with a 24° mutual inclination. In cases like this, the question of how many planets in a multi-transiting system might be seen in transit at any one time becomes relevant (Brakensiek & Ragozzine 2016), as large mutual inclinations can lead to only a subset of the planets being seen from a given line of sight. K2-266 b is currently observed to have a grazing transit and a high mutual inclination with the remainder of the planets, the five of which reside in a roughly coplanar configuration. With an aim toward assessing where this system fits into the Kepler dichotomy, in this section we conduct an analysis of the transit likelihood for various numbers of planets in this system.

To test the dynamical and transit stability of these planets, we ran 250 N-body simulations of the system, drawing the initial orbital elements from the posteriors generated from the EXOFAST transit fit (more specifically, we draw a single link from the MCMC posterior at random for each of the 250 simulations, and then use all orbital elements from that link). We assigned the longitude of ascending node to be 2π, as it cannot be measured from the transit fits. The planetary masses are drawn from the posteriors provided by the EXOFAST fit, as are stellar mass and radius. For all calculated values of inclination, we broke the above/below solar midplane degeneracy by randomly assigning the value to be either greater or smaller than 90°. We also assume that the stellar obliquity is aligned with the plane containing the outer five planets (but there are no system-specific observations to support this assumption; instead, we make this assumption as a computational necessity, although we expect the stellar obliquity to be more aligned in multi-planet systems; see Morton & Winn (2014)). These 250 simulations were carried out using the hybrid Wisdom–Holman and Bulirsch–Stoer (B-S) integrator Mercury6 (Chambers 1999) for integration times of 10⁵ years, and with an initial time step of 8.5 minutes. Energy was conserved to better than one part in 10⁸ over the course of the simulations. When physical collisions occur, particles are removed from the simulation. The integration time was chosen to include many secular timescales of the system (Figure 9 shows that many periods of secular oscillations are included in the 10⁵ years time span).

In roughly 66% of our suite of 250 simulations, at least two planets in the system attain orbits that cross. In 23% of the simulations, the system experiences a true dynamical instability, in which a planet is ejected from the system or physically collides with another body. In the cases in which orbits cross, planets 0.02 and c are the culprits of the instability roughly 80% of the time. On the 10⁵ year integrations considered in this work, a sizable fraction (roughly a third) remain dynamically stable. As such, we cannot use dynamical arguments to argue against the existence of planet candidate 0.02, whose close orbit with planet c might otherwise be suspect.

As neither of the planet candidates can clearly be ruled out based on dynamical arguments, we next consider the dynamical evolution that leads to only a subset of the six known planets being seen in transit. We currently observe the system to be a six-planet system, but the innermost planet K2-266 b has a high measured impact parameter, meaning that is barely transiting. The simulations show significant inclination evolution over time for both K2-266 b and the other planets in the system. In Figure 9, we plot a representative case from our set of 250 simulations, where the semimajor axes and eccentricities of all planets remain confined relatively close to their currently measured values, but the orbital inclination of all six planets evolves.

One notable feature of the numerical simulations is the inclination evolution of all six planets. Secular evolution of the system causes the planetary orbits to evolve with time, resulting in configurations in which not all planets can be simultaneously seen in transit. Inspired by the present-day (apparently serendipitous) geometry, we extracted from the simulations the transit probability over time for varying numbers of planets in this system. The result of this analysis is presented in Figure 10 for two lines of sight. The first case considers the fixed line of sight corresponding to our current location (that of the solar system). The second case uses an optimized variable line of sight, which is recomputed at each time step of each integration to determine the largest multiplicity that can be observed from any location in the galaxy.

This analysis shows that observing six transiting planets in the system is rare, given the known components of the system. No matter which line of sight is considered, the system will appear to contain the six "known" planets a minority of the time. More commonly, the system will be seen as a five-planet system from the most favorable line of sight, and as a one- or two-planet system from our current line of sight. Most of the time, the inclinations of the outer five planets evolve and cause them to reside in non-transiting configurations. On the other hand, the probabilities are not vanishingly small. We expect to be able to reobserve a six-planet transiting system about 2.2% of the time in the future from our current line of sight, given the currently measured orbits of these planets. It is also important to note that we cannot be sure that the observed six planets are the only planets in this system: additional, non-transiting planets would alter the dynamics described here.

The sub-Neptunes K2-266 d and e, with periods of 19.482 days and 14.697 days, have a period ratio of 1.326, which is 0.59% away from the 4:3 mean motion commensurability. These planets reside in nearly the same orbital plane (with 8973 and 8974 inclinations). As such, the orbital periods of these planets suggest that they may reside in orbital resonance. However, true resonance is characterized by a librating resonance angle, and it is not clear from the orbital elements alone whether the resonance angle will librate or circulate. To determine the true resonance behavior of these two planets, we used the subset of the 250 simulations that did not experience orbit-crossing during the 10⁵ year integration time, and post-processed the results to search for resonances.

We have performed a resonance-finding algorithm to identify the time intervals in the simulations where the planets were in true resonance. We found regimes with nearly constant period ratio (constant P_e/P_d), generated arrays of resonance angles for all p:q resonances up to 29th order (while p ≤ 30), and automatically generated plots using the simulated orbital elements of planets K2-266 d and e for each resonance angle for each resonance order for all of the 250 simulations. Using the resulting resonance angles, we searched for librating behavior by breaking the time series into 5000 year intervals and searching for gaps in resonance angle space: Note that a circulating resonance angle will populate the entire 360° range of possible angles, whereas a librating angle will have gaps.

We find that, in our simulations, planets K2-266 d and e exhibit orbital resonances approximately 8.1% of the time over the simulations under consideration. The resonance angles populated in these cases have the forms

$$\begin{eqnarray}&&{\psi }_{1}=4{\lambda }_{o}-3{\lambda }_{i}-{\varpi }_{o},\end{eqnarray} \tag{ 3 }$$

and

$$\begin{eqnarray}&&{\psi }_{2}=4{\lambda }_{o}-3{\lambda }_{i}-{\varpi }_{i},\end{eqnarray} \tag{ 4 }$$

where λ is mean longitude and ϖ is longitude of pericenter. The subscripts denote the inner (i) and outer (o) planets. The four types of resonance behavior exhibited by this system, in order of occurrence rate, include: non-resonance, continuous resonance for the entire simulation lifetime, an initial condition close to resonance that loses the resonance as the system evolves, and very rarely, the attainment of a resonance after an initial period of non-resonance (see Ketchum et al. (2013) for a more detailed discussion of this process). We find that, for trials that start out in a resonant configuration, the typical libration width of the resonance generally is consistently around 190°. Although this width may evolve slightly as the simulation progresses, the resonances are not typically much deeper than this initial value.

Dynamical systems are often chaotic and we would like to quantify the chaotic behavior of K2-266. The system, as observed, has six planets in a compact configuration with a relatively large mutual inclination between the innermost planet and the others. Our numerical simulations, described above, indicate that while the outer planets (K2-266 d, e, and .06) are generally dynamically stable, the middle planets (K2-266.02 and d) can experience scattering or other non-periodic time evolution, potentially leading to orbit crossing. Non-periodic behavior of this nature can be indicative of chaos.

The evolution of a chaotic planetary system is extremely sensitive to its initial conditions. Chaos is often parameterized by the Lyapunov exponents of the system, which determine the rate of exponential divergence of orbits with similar initial conditions. In contrast, general observational uncertainties in the orbital elements can also lead to non-chaotic divergence if initial orbital elements are drawn from different locations of the posteriors (in cases where ensembles of simulations are used to sample the uncertainties). Either sufficiently large observational errors or the sufficiently rapid onset of chaos will thus make both numerical integrations and analytic explorations less certain.

To test the chaotic behavior of the K2-266 system, we ran 400 integrations of this system, with the orbital parameters and masses drawn from the posteriors generated by the transit fit. Each simulation was carried out using the Rebound N-Body integration package (Rein & Liu 2012), where we used a total integration time of 1000 years, the IAS15 integrator (Rein & Tamayo 2015), and an initial time step of eight minutes. For each of these integrations, we evaluate the chaotic nature of the initial conditions by employing the Mean Exponential Growth factor of Nearby Orbits (MEGNO) indicator (Cincotta et al. 2003), implemented in the Rebound N-body code. For chaotic trajectories, the MEGNO indicator, Y, grows linearly in time at a rate of 1/t_Ly where t_Ly is the Lyapunov time, while for regular trajectories it asymptotically approaches Y = 2. We compute MEGNO values for the 400 draws from the posteriors of this system at times between 1 and 1000 years. These realizations provide a good sample of the parameter space spanned by the observational posteriors. Of the 400 realizations, only 4.5% can be categorized as unambiguously regular at the end of the integration (where the criterion for regularity is taken to be MEGNO < 4). Moreover, we find no strong correlations between planetary parameters and MEGNO values using our current simulation set. We attempted to trace chaotic behavior using the period ratio of the resonant planets d and e (as done in Figure 3 of Deck et al. 2012), the ratio perihelion/aphelion of .02 and c, and by using the orbital elements of planet b, but no trends emerged. This finding is likely due to the high multiplicity and tightly packed nature of the system: there are multiple equally-important sources of dynamical chaos. In Figure 11, we plot the median MEGNO indicator value for these 400 simulations considered in this section at periodic intervals in the 1000 year integrations. A MEGNO indicator less than four denotes regular orbits, while a MEGNO of four or greater denotes measured chaos. The median MEGNO indicator reaches four (denoting the measurable onset of chaos) at roughly 100 years.

From this MEGNO analysis, we know that the majority of orbits allowed by the transit posteriors are chaotic. It is important to note that, while some of these chaotic posterior points are likely destined to experience instabilities based on our 10⁵ year numerical simulations, a chaotic system does not necessarily mean a dynamically unstable system, or even a particularly active system (the planetary orbits in our solar system are known to be chaotic, as is the Kepler-36 planetary system; see Deck et al. (2012)). Chaos implies that similar initial conditions will diverge over some timescale, so precise future predictions of planetary orbits can no longer be made. Specifically, for two given sets of similar initial conditions, integrations of both cases could result in systems that are dynamically stable and continuously transit, but the values of the phase space variables (including planet locations) can diverge over time if the system is chaotic. One implication of this analysis is that a large amount of uncertainty in forward integrations comes from chaos, rather than only from the uncertainty of the transit posteriors.

To check that EPIC 248435395 is actually an associated companion to K2-266, we directly compare the Gaia Data Release 2 proper motions and distances for both systems (Gaia Collaboration et al. 2016, 2018). K2-266 has a Gaia parallax of 12.87 ± 0.06 mas, corresponding to a distance of 77.8 ± 0.6 pc, while EPIC 248435395 has a parallax of 12.85 ± 0.06 mas and a distance of 77.7 ± 0.6 pc. All systematic uncertainties on the Gaia DR2 parallax should be <0.1 mas. Therefore, the two stars are at the same distance. The Gaia DR2 proper motions for K2-266 are μ_α = 56.9 mas yr⁻¹ and μ_δ = −68.8 mas yr⁻¹. These are very similar to the proper motions for EPIC 248435395, which are μ_α = 53.2 mas yr⁻¹ and μ_δ = −72.7 mas yr⁻¹ (see Table 1). The Gaia parallaxes and proper motions are consistent with the two stars being a widely separated binary.

Using this distance, the ∼42'' separation would correspond to an orbital semimajor axis of ∼3200 au. Using our determined mass of EPIC 248435395 of 0.561 $\,{M}_{\odot }$ and K2-266 of 0.649 $\,{M}_{\odot }$ this would correspond to an orbital period 0.16 Myr. This period would result in a ∼4 mas motion from the binary orbit, within the expected micro-arcsecond expected astrometric precision (Gaia Collaboration et al. 2018). However, EPIC 248435395 will move ∼045 over the nominal five-year Gaia mission, and it may be difficult to differentiate the contribution from the binary orbit from the star's proper motion. Moreover, as we show in the following section, EPIC 248435395 is likely a binary itself, which could further confuse the astrometric solution.

The ∼42'' separation means that the two systems were well-resolved by K2. We were therefore able to produce a separate light curve for the companion star. As we did for K2-266 and following the strategy described in Vanderburg & Johnson (2014) and Vanderburg et al. (2016a), we searched EPIC 248435395 for possible transit signals. From this search, we identify a transit signal at 8.5 days with an S/N of 7.0 around EPIC 248435395, below the typical Kepler S/N threshold of 7.1 and after K2 of 9. The final light curve for EPIC 248435395 has a 30 minute noise level of 151 ppm and a 6 hour photometric precision of 42 ppm.

We also inspected the NGS POSS and Pan-Stars images of EPIC 248435395 in a manner similar to that described in Section 2.4 (see Figure 3). Because it shares a common proper motion with K2-266, it has also moved about 6 arcseconds since it was imaged by POSS in 1952. Because EPIC 248435395 is fainter than K2-266, its saturated point spread function does not extend all the way to the star's current-day position, so we are able to rule out background objects down to the POSS limiting magnitude of about 20 in the red, and 19 in the blue. We identified no nearby companions in Pan-STARRS imaging. We also obtained high-resolution images in the Br-γ filter and the J-band of EPIC 248435395 using NIRC2 on UT 2017 December 28 (see Section 2.5 and Figure 3). We see no evidence for any additional companions from our "patient" and AO imaging.

We also obtained two spectra of EPIC 248435395 using TRES, on UT 2018 February 7 and UT 2018 April 3. They were analyzed according to the same procedures described in Section 2.2. The absolute RV of EPIC 248435395 is 12.84 ± 0.31 km s⁻¹, about 2 km s⁻¹ different from K2-266. Given their projected separation, this is slightly larger than one would expect if the two stars are bound. However, we note that this is unlikely to be the true systemic velocity of EPIC 248435395. The large error on the velocity is not because of poor RV precision, but because our two RVs exhibit a 450 m s⁻¹ variation; EPIC 248435395 is likely to be a binary itself, and we do not know the full amplitude of variation or its orbital phase. We conclude that the RVs of the components of the wide pair are not inconsistent with expectations for a bound system; K2-266 is likely to be the primary star in a hierarchical triple system. The TRES RVs for EPIC 248435395 are shown in Table 2.

Table 2.  Relative Radial Velocities for K2-266 and EPIC 248435395

${\mathrm{BJD}}_{\mathrm{TDB}}$	RV (m s−1)	σRV (m s−1)	Target
2458081.028322	−30.1	17.7	K2-266
2458090.980928	−2.5	32.0	K2-266
2458106.965261	−3.2	43.6	K2-266
2458107.923621	−55.5	28.2	K2-266
2458211.644609	−26.3	21.2	K2-266
2458212.644636	−36.5	25.6	K2-266
2458213.663584	−23.9	33.3	K2-266
2458218.808648	−30.0	32.5	K2-266
2458156.911319	425.9	46.9	EPIC 248435395
2458211.675712	−21.6	47.9	EPIC 248435395

Download table as:  ASCIITypeset image

Following a procedure similar to the one used for K2-266 (see Section 3), we perform a fit of the exoplanet candidate around EPIC 248435395 using EXOFASTv2 (Eastman 2017). Within this analysis, we simultaneously fit the flattened K2 light curve (see Figure 12) and SED (see Figure 4) for EPIC 248435395. To constrain the mass of EPIC 248435395 within the fit, we use a Gaussian prior of $0.584\pm 0.015\,{M}_{\odot }$ from Mann et al. (2015), but with the uncertainties inflated to 5%. We constrain the radius within the fit using the broadband photometry shown in Table 1. We set a starting point on the ${T}_{\mathrm{eff}}$ and $[\mathrm{Fe}/{\rm{H}}]$ of the host star to be 3699 K and −0.006 dex from the EPIC catalog (Huber et al. 2016). We also enforce the same upper limit on the V-band extinction from the Schlegel et al. (1998) dust maps of 0.0548 as we did for the fit of K2-266. We use the Gaia DR2 parallax (12.85 mas) with a conservative uncertainty of 0.1 mas as prior in the fit. The final determined system parameters are shown in Tables 3 and 7. We note that this is not a confirmed or validated planet.

Table 3.  Median Values and 68% Confidence Intervals for the Stellar Parameters of the K2-266 and EPIC 248435395 from EXOFASTv2

Stellar	Parameters	K2-266	EPIC 248435395
M*	Mass (M⊙)	0.686 ± 0.033	0.581 ± 0.029
R*	Radius (R⊙)	${0.703}_{-0.022}^{+0.024}$	${0.649}_{-0.035}^{+0.040}$
L*	Luminosity (L⊙)	0.1502 ± 0.0057	${0.0617}_{-0.0028}^{+0.0027}$
ρ*	Density (cgs)	${2.79}_{-0.30}^{+0.29}$	${2.99}_{-0.51}^{+0.56}$
$\mathrm{log}g$	Surface gravity (cgs)	${4.581}_{-0.037}^{+0.032}$	${4.577}_{-0.056}^{+0.052}$
Teff	Effective temperature (K)	${4285}_{-57}^{+49}$	3570 ± 110
[Fe/H]	Metallicity	$-{0.12}_{-0.42}^{+0.40}$	$-{0.19}_{-0.81}^{+0.43}$
Av	V-band extinction	${0.029}_{-0.019}^{+0.018}$	${0.028}_{-0.019}^{+0.018}$
σSED	SED photometry error scaling	${5.0}_{-1.1}^{+1.6}$	${5.5}_{-1.2}^{+1.8}$
πa	Parallax (mas)	12.960 ± 0.100	12.928 ± 0.098
d	Distance (pc)	${77.16}_{-0.59}^{+0.60}$	${77.35}_{-0.58}^{+0.59}$

Note.

^aThe MIST Isochrones were not used in the EXOFASTv2 fit for K2-266.

Download table as:  ASCIITypeset image

Table 4.  Median Values and 68% Confidence Intervals for Planetary Parameters of K2-266 from EXOFASTv2

Planetary Parameters:		b	K2-266.02	c	d	e	K2-266.06
P	Period (days)	0.658524 ± 0.000017	${6.1002}_{-0.0017}^{+0.0015}$	${7.8140}_{-0.0016}^{+0.0019}$	${14.69700}_{-0.00035}^{+0.00034}$	19.4820 ± 0.0012	${56.682}_{-0.018}^{+0.019}$
RP	Radius ($\,{R}_{\oplus }$)	${3.3}_{-1.3}^{+1.8}$	${0.646}_{-0.091}^{+0.099}$	${0.705}_{-0.085}^{+0.096}$	${2.93}_{-0.12}^{+0.14}$	${2.73}_{-0.11}^{+0.14}$	${0.90}_{-0.12}^{+0.14}$
TC	Time of conjunction (${\mathrm{BJD}}_{\mathrm{TDB}}$)	${2457945.7235}_{-0.0030}^{+0.0032}$	${2457913.413}_{-0.011}^{+0.013}$	${2457907.5812}_{-0.012}^{+0.0099}$	2457944.8393 ± 0.0012	2457938.5410 ± 0.0013	${2457913.436}_{-0.013}^{+0.014}$
T0	Optimal conjunction time (${\mathrm{BJD}}_{\mathrm{TDB}}$)	${2457949.6747}_{-0.0030}^{+0.0032}$	${2457943.9143}_{-0.0064}^{+0.0066}$	2457946.6510 ± 0.0064	2457944.8393 ± 0.0012	2457938.5410 ± 0.0013	${2457913.436}_{-0.013}^{+0.014}$
a	Semimajor axis (au)	${0.01306}_{-0.00021}^{+0.00020}$	${0.05761}_{-0.00093}^{+0.00090}$	0.0679 ± 0.0011	${0.1035}_{-0.0017}^{+0.0016}$	${0.1249}_{-0.0020}^{+0.0019}$	${0.2546}_{-0.0041}^{+0.0040}$
i	Inclination (degrees)	${75.32}_{-0.70}^{+0.62}$	${87.84}_{-0.46}^{+0.84}$	${88.28}_{-0.41}^{+0.81}$	${89.46}_{-0.25}^{+0.32}$	${89.45}_{-0.18}^{+0.25}$	${89.40}_{-0.14}^{+0.26}$
e	Eccentricity	⋯	${0.051}_{-0.036}^{+0.051}$	${0.042}_{-0.030}^{+0.043}$	${0.047}_{-0.032}^{+0.043}$	${0.043}_{-0.030}^{+0.036}$	${0.31}_{-0.17}^{+0.11}$
${\omega }_{* }$	Argument of periastron (degrees)	⋯	${88}_{-62}^{+60}$	87 ± 61	87 ± 62	${89}_{-58}^{+57}$	${83}_{-59}^{+57}$
Teq	Equilibrium temperature (K)	1515 ± 18	${721.7}_{-8.4}^{+8.7}$	${664.5}_{-7.7}^{+8.0}$	${538.3}_{-6.3}^{+6.5}$	${490.1}_{-5.7}^{+5.9}$	${343.3}_{-4.0}^{+4.1}$
MP	Mass ($\,{M}_{\oplus }$)	${11.3}_{-6.5}^{+11}$	${0.209}_{-0.089}^{+0.15}$	${0.29}_{-0.11}^{+0.17}$	${9.4}_{-2.0}^{+2.9}$	${8.3}_{-1.8}^{+2.7}$	${0.70}_{-0.30}^{+0.87}$
K	RV semi-amplitude (m s−1)	${10.3}_{-5.9}^{+10.}$	${0.094}_{-0.040}^{+0.067}$	${0.119}_{-0.046}^{+0.073}$	${3.17}_{-0.69}^{+0.99}$	${2.53}_{-0.56}^{+0.84}$	${0.158}_{-0.068}^{+0.20}$
logK	Log of RV semi-amplitude	${1.01}_{-0.37}^{+0.31}$	$-{1.03}_{-0.24}^{+0.23}$	-0.92 ± 0.21	${0.50}_{-0.11}^{+0.12}$	${0.40}_{-0.11}^{+0.12}$	$-{0.80}_{-0.24}^{+0.35}$
${R}_{P}/{R}_{* }$	Radius of planet in stellar radii	${0.043}_{-0.017}^{+0.023}$	${0.0084}_{-0.0011}^{+0.0012}$	${0.0092}_{-0.0011}^{+0.0012}$	${0.03827}_{-0.00057}^{+0.00071}$	${0.03564}_{-0.00061}^{+0.00067}$	${0.0117}_{-0.0015}^{+0.0018}$
$a/{R}_{* }$	Semimajor axis in stellar radii	${4.00}_{-0.15}^{+0.14}$	${17.63}_{-0.66}^{+0.60}$	${20.80}_{-0.78}^{+0.71}$	${31.7}_{-1.2}^{+1.1}$	${38.2}_{-1.4}^{+1.3}$	${77.9}_{-2.9}^{+2.7}$
δ	Transit depth (fraction)	${0.0018}_{-0.0012}^{+0.0025}$	${0.000071}_{-0.000018}^{+0.000022}$	${0.000085}_{-0.000019}^{+0.000023}$	${0.001465}_{-0.000043}^{+0.000055}$	${0.001270}_{-0.000043}^{+0.000048}$	${0.000136}_{-0.000033}^{+0.000046}$
Depth	Flux decrement at mid-transit	${0.000596}_{-0.000091}^{+0.000068}$	${0.000071}_{-0.000018}^{+0.000022}$	${0.000085}_{-0.000019}^{+0.000023}$	${0.001465}_{-0.000043}^{+0.000055}$	${0.001270}_{-0.000043}^{+0.000048}$	${0.000136}_{-0.000033}^{+0.000046}$
τ	Ingress/egress transit duration (days)	${0.00686}_{-0.00045}^{+0.00051}$	${0.00116}_{-0.00032}^{+0.00042}$	${0.00134}_{-0.00033}^{+0.00045}$	${0.00569}_{-0.00045}^{+0.00067}$	${0.00601}_{-0.00063}^{+0.00073}$	${0.00281}_{-0.00096}^{+0.0019}$
T14	Total transit duration (days)	${0.01389}_{-0.00085}^{+0.0012}$	${0.083}_{-0.012}^{+0.014}$	${0.094}_{-0.014}^{+0.016}$	${0.1420}_{-0.0015}^{+0.0014}$	${0.1527}_{-0.0015}^{+0.0014}$	${0.143}_{-0.023}^{+0.024}$
TFWHM	FWHM transit duration (days)	${0.00695}_{-0.00042}^{+0.00062}$	${0.082}_{-0.013}^{+0.015}$	${0.092}_{-0.014}^{+0.016}$	0.1362 ± 0.0015	${0.1466}_{-0.0016}^{+0.0015}$	0.140 ± 0.024
b	Transit impact parameter	${1.011}_{-0.024}^{+0.027}$	${0.64}_{-0.25}^{+0.13}$	${0.60}_{-0.28}^{+0.14}$	${0.29}_{-0.17}^{+0.12}$	${0.36}_{-0.16}^{+0.10}$	${0.64}_{-0.33}^{+0.17}$
bS	Eclipse impact parameter	⋯	${0.68}_{-0.26}^{+0.13}$	${0.64}_{-0.30}^{+0.14}$	${0.31}_{-0.18}^{+0.12}$	${0.379}_{-0.16}^{+0.100}$	${0.86}_{-0.33}^{+0.30}$
${\tau }_{S}$	Ingress/egress eclipse duration (days)	⋯	${0.00131}_{-0.00038}^{+0.00055}$	${0.00148}_{-0.00039}^{+0.00057}$	${0.00613}_{-0.00048}^{+0.00064}$	${0.00641}_{-0.00052}^{+0.00069}$	${0.0037}_{-0.0037}^{+0.0035}$
${T}_{S,14}$	Total eclipse duration (days)	⋯	${0.084}_{-0.016}^{+0.019}$	${0.096}_{-0.017}^{+0.020}$	${0.1485}_{-0.0051}^{+0.0100}$	${0.1593}_{-0.0052}^{+0.010}$	${0.13}_{-0.13}^{+0.10}$
${T}_{S,\mathrm{FWHM}}$	FWHM eclipse duration (days)	⋯	${0.083}_{-0.017}^{+0.020}$	${0.094}_{-0.018}^{+0.021}$	${0.1424}_{-0.0052}^{+0.0099}$	${0.1528}_{-0.0053}^{+0.011}$	${0.13}_{-0.13}^{+0.11}$
${\delta }_{S,3.6\mu {\rm{m}}}$	Blackbody eclipse depth at 3.6 μm (ppm)	${210}_{-140}^{+290}$	${0.41}_{-0.11}^{+0.14}$	${0.299}_{-0.069}^{+0.088}$	${1.24}_{-0.12}^{+0.15}$	${0.511}_{-0.057}^{+0.070}$	${0.00160}_{-0.00044}^{+0.00062}$
${\delta }_{S,4.5\mu {\rm{m}}}$	Blackbody eclipse depth at 4.5 μm (ppm)	${280}_{-180}^{+380}$	${0.94}_{-0.24}^{+0.31}$	${0.76}_{-0.17}^{+0.22}$	${4.25}_{-0.35}^{+0.44}$	${2.05}_{-0.19}^{+0.24}$	${0.0135}_{-0.0036}^{+0.0050}$
${\rho }_{P}$	Density (cgs)	${1.77}_{-0.88}^{+1.9}$	${4.27}_{-0.66}^{+0.79}$	${4.51}_{-0.67}^{+0.81}$	${2.03}_{-0.43}^{+0.64}$	${2.21}_{-0.47}^{+0.70}$	${5.27}_{-0.87}^{+1.5}$
loggP	Surface gravity	${3.00}_{-0.13}^{+0.14}$	2.69 ± 0.12	2.75 ± 0.11	${3.026}_{-0.100}^{+0.12}$	${3.03}_{-0.10}^{+0.12}$	${2.93}_{-0.12}^{+0.18}$
Θ	Safronov Number	${0.0046}_{-0.0015}^{+0.0018}$	${0.00191}_{-0.00065}^{+0.00094}$	${0.00285}_{-0.00086}^{+0.0012}$	${0.0340}_{-0.0072}^{+0.010}$	${0.0387}_{-0.0082}^{+0.012}$	${0.0203}_{-0.0070}^{+0.017}$
$\langle F\rangle$	Incident flux (109 erg s−1 cm−2)	${1.197}_{-0.055}^{+0.059}$	${0.0612}_{-0.0028}^{+0.0031}$	${0.0441}_{-0.0020}^{+0.0022}$	${0.01898}_{-0.00086}^{+0.00093}$	${0.01304}_{-0.00060}^{+0.00064}$	${0.00286}_{-0.00024}^{+0.00025}$
TP	Time of periastron (${\mathrm{BJD}}_{\mathrm{TDB}}$)	${2457945.7235}_{-0.0030}^{+0.0032}$	${2457913.43}_{-0.95}^{+0.93}$	2457907.6 ± 1.2	2457944.8 ± 2.3	2457938.6 ± 2.9	${2457913.3}_{-5.2}^{+5.0}$
TS	Time of eclipse (${\mathrm{BJD}}_{\mathrm{TDB}}$)	${2457946.0528}_{-0.0030}^{+0.0032}$	2457916.46 ± 0.17	${2457903.67}_{-0.18}^{+0.19}$	${2457937.49}_{-0.40}^{+0.41}$	2457948.28 ± 0.44	${2457885.2}_{-8.7}^{+8.6}$
TA	Time of ascending node (${\mathrm{BJD}}_{\mathrm{TDB}}$)	${2457945.5589}_{-0.0030}^{+0.0032}$	${2457911.942}_{-0.076}^{+0.13}$	${2457905.682}_{-0.080}^{+0.14}$	${2457941.28}_{-0.17}^{+0.28}$	${2457933.82}_{-0.19}^{+0.32}$	${2457903.3}_{-4.2}^{+3.4}$
TD	Time of descending node (${\mathrm{BJD}}_{\mathrm{TDB}}$)	${2457945.8881}_{-0.0030}^{+0.0032}$	${2457914.884}_{-0.13}^{+0.077}$	${2457909.478}_{-0.14}^{+0.078}$	${2457948.40}_{-0.28}^{+0.17}$	${2457943.26}_{-0.33}^{+0.19}$	${2457923.6}_{-3.5}^{+4.3}$
$e\cos {\omega }_{* }$		⋯	$-{0.000}_{-0.045}^{+0.044}$	0.000 ± 0.037	-0.000 ± 0.043	$-{0.000}_{-0.035}^{+0.036}$	0.00 ± 0.24
$e\sin {\omega }_{* }$		⋯	${0.026}_{-0.022}^{+0.046}$	${0.022}_{-0.018}^{+0.038}$	${0.025}_{-0.020}^{+0.035}$	${0.025}_{-0.020}^{+0.034}$	${0.19}_{-0.15}^{+0.16}$
V/VC		⋯	${0.973}_{-0.044}^{+0.022}$	${0.978}_{-0.037}^{+0.018}$	${0.975}_{-0.034}^{+0.020}$	${0.975}_{-0.033}^{+0.020}$	${0.80}_{-0.12}^{+0.15}$
${M}_{P}\sin i$	Minimum mass ($\,{M}_{\oplus }$)	${10.9}_{-6.3}^{+11}$	${0.208}_{-0.089}^{+0.15}$	${0.29}_{-0.11}^{+0.17}$	${9.4}_{-2.0}^{+2.9}$	${8.3}_{-1.8}^{+2.7}$	${0.70}_{-0.30}^{+0.87}$
${M}_{P}/{M}_{* }$	Mass ratio	${0.000049}_{-0.000028}^{+0.000051}$	${0.00000091}_{-0.00000039}^{+0.00000065}$	${0.00000126}_{-0.00000049}^{+0.00000077}$	${0.0000413}_{-0.0000091}^{+0.000013}$	${0.0000362}_{-0.0000081}^{+0.000012}$	${0.0000031}_{-0.0000013}^{+0.0000039}$
$d/{R}_{* }$	Separation at mid-transit	${4.00}_{-0.15}^{+0.14}$	${16.98}_{-0.92}^{+0.81}$	${20.17}_{-1.0}^{+0.87}$	${30.7}_{-1.7}^{+1.4}$	${37.0}_{-1.7}^{+1.6}$	${59}_{-11}^{+13}$
PT	A priori non-grazing transit prob	${0.2389}_{-0.0092}^{+0.010}$	${0.0584}_{-0.0027}^{+0.0033}$	${0.0491}_{-0.0020}^{+0.0026}$	${0.0313}_{-0.0014}^{+0.0018}$	${0.0261}_{-0.0011}^{+0.0012}$	${0.0166}_{-0.0030}^{+0.0039}$
${P}_{T,G}$	A priori transit prob	${0.2616}_{-0.0098}^{+0.011}$	${0.0594}_{-0.0027}^{+0.0034}$	${0.0500}_{-0.0021}^{+0.0026}$	${0.0338}_{-0.0015}^{+0.0020}$	${0.0280}_{-0.0012}^{+0.0013}$	${0.0170}_{-0.0030}^{+0.0039}$
PS	A priori non-grazing eclipse prob	⋯	0.0546 ± 0.0025	0.0465 ± 0.0020	${0.02941}_{-0.00093}^{+0.0013}$	${0.0245}_{-0.0010}^{+0.0012}$	${0.0113}_{-0.0016}^{+0.0017}$
${P}_{S,G}$	A priori eclipse prob	⋯	0.0555 ± 0.0026	${0.0473}_{-0.0020}^{+0.0021}$	${0.0317}_{-0.0010}^{+0.0015}$	${0.0263}_{-0.0012}^{+0.0013}$	${0.0116}_{-0.0017}^{+0.0018}$

Download table as:  ASCIITypeset images: 1 2

Table 5.  Median Values and 68% Confidence Intervals for the Additional Parameters of K2-266 from EXOFASTv2

Wavelength Parameters:		Kepler
u1	Linear limb-darkening coeff	${0.630}_{-0.10}^{+0.064}$
u2	Quadratic limb-darkening coeff	${0.107}_{-0.054}^{+0.082}$
Telescope Parameters:		TRES
γrel	Relative RV offset (m s−1)	$-{29.9}_{-8.6}^{+7.4}$
σJ	RVjitter (m s−1)	${0.00}_{-0.00}^{+19}$
${\sigma }_{J}^{2}$	RV jitter variance	$-{130}_{-140}^{+480}$
Transit Parameters:

Planet	Transit Date	Added Variance	Transit Mid Time	Baseline flux
		σ2 × 10−10	${\mathrm{BJD}}_{\mathrm{TDB}}$	F0
b, c, d, g	Full K2 LC	${7.6}_{-2.7}^{+2.9}$	N/A	1.0000035 ± 0.0000030
e	K2 UT 2017-06-10	${19}_{-18}^{+26}$	2457915.44761 ± 0.00106	0.999981 ± 0.000016
f	K2 UT 2017-06-14	${10}_{-15}^{+21}$	2457919.05628 ± 0.00088	1.000012 ± 0.000014
e	K2 UT 2017-06-25	$-{3}_{-11}^{+15}$	2457930.13813 ± 0.00101	${1.000002}_{-0.000012}^{+0.000013}$
f	K2 UT 2017-07-04	${25}_{-18}^{+26}$	2457938.54211 ± 0.00108	${0.999986}_{-0.000017}^{+0.000016}$
e	K2 UT 2017-07-10	${16}_{-17}^{+25}$	2457944.83597 ± 0.00137	${0.999971}_{-0.000015}^{+0.000016}$
f	K2 UT 2017-07-23	${54}_{-27}^{+41}$	2457958.03343 ± 0.00112	${0.999986}_{-0.000021}^{+0.000020}$
e	K2 UT 2017-07-25	${13}_{-14}^{+21}$	2457959.52590 ± 0.00090	${1.000000}_{-0.000015}^{+0.000014}$
e	K2 UT 2017-08-08	${30}_{-19}^{+26}$	2457974.23919 ± 0.00114	1.000009 ± 0.000017
f	K2 UT 2017-08-12	${5}_{-15}^{+23}$	2457977.50614 ± 0.00120	1.000007 ± 0.000015

Download table as:  ASCIITypeset image

The origin of Neptune-sized planets is not clear, yet they appear to be one of the most common types of planets. The planets range in size from 2–6 $\,{R}_{\oplus }$ and have been discovered orbiting >25% of all stars (Howard et al. 2012; Fressin et al. 2013; Buchhave et al. 2014; Fulton et al. 2017). Recent statistical studies of the observed amplitude of transmission spectral features of warm Neptunes show a correlation with equilibrium temperature or its bulk H/He mass fraction (Crossfield & Kreidberg 2017). However, there are only a small number of Neptune-sized planets that are amenable to transmission spectroscopy with current facilities like the Hubble Space Telescope (HST).

To understand if the planets orbiting K2-266 would be viable targets for transmission spectroscopic measurements, we follow the technique described in Vanderburg et al. (2016a) to calculate their expected atmospheric scale height and S/N per transit. Using data from NASA's Exoplanet archive (Akeson et al. 2013), we also calculate the atmospheric scale height and S/N for all planets with R_p < 3 $\,{R}_{\oplus }$ (see Table 6), updating this table from Rodriguez et al. (2018). The calculations are done in the H-band to understand their accessibility using the Wide Field Camera 3 instrument on HST, as well as the future feasibility to observe them with the suite of instruments on the upcoming James Webb Space Telescope (JWST). Purely based on the inferred sizes of the planets, it is expected that K2-266 b, d, and e might all have thick gaseous atmospheres (Weiss & Marcy 2014), but the uncertainty in the size of planet b (due to the grazing transit configuration) and its proximity to its host star makes this unclear (see Figure 6 for the probability distribution of planet b's radius from our global fit). While our transit fit indicates that the most probable size of planet b is about three times the size of the Earth, virtually all known USP planets known are smaller than 2 $\,{R}_{\oplus }$ (Sanchis-Ojeda et al. 2014; Winn et al. 2017)

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 6.  The Best Confirmed or Validated Planets for Transmission Spectroscopy with R_P < 3 $\,{R}_{\oplus }$

Planet	Period (days)	RP($\,{R}_{\oplus }$)	S/Na	Reference	Discovery
GJ 1214 b	1.58	2.85	1.00	Charbonneau et al. (2009)	MEarth
55 Cnc e	0.74	1.91	0.41	Dawson & Fabrycky (2010)	RVs
HD 97658 b	9.49	2.35	0.27	Dragomir et al. (2013)	RVs
TRAPPIST-1 f	9.21	1.04	0.24	Gillon et al. (2017)	Spitzer
K2-136 c	17.31	2.91	0.19	Ciardi et al. (2018), Livingston et al. (2018), Mann et al. (2018)	K2
GJ 9827 b	1.21	1.75	0.18	Niraula et al. (2017), Rodriguez et al. (2018)	K2
K2-167 b	9.98	2.82	0.16	Mayo et al. (2018)	K2
K2-266 e	14.70	2.93	0.15	This work	K2
GJ 9827 d	6.20	2.10	0.15	Niraula et al. (2017), Rodriguez et al. (2018)	K2
HIP 41378 b	15.57	2.90	0.14	Vanderburg et al. (2016a)	K2
HD 3167 b	0.96	1.70	0.14	Vanderburg et al. (2016b), Christiansen et al. (2017)	K2
K2-233 d	24.37	2.65	0.13	David et al. (2018)	K2
K2-266 f	19.48	2.73	0.12	This work	K2
K2-28 b	2.26	2.32	0.12	Hirano et al. (2016)	K2
K2-199 c	7.37	2.84	0.12	Mayo et al. (2018)	K2
K2-155 c	13.85	2.60	0.11	Diez Alonso et al. (2018), Hirano et al. (2018)	K2
Kepler-410 A b	17.83	2.84	0.10	Van Eylen et al. (2014)	Kepler
HD 106315 b	9.55	2.40	0.10	Rodriguez et al. (2017), Crossfield et al. (2017)	K2

Note.

^aThe predicted signal-to-noise ratios relative to GJ 1214 b. All values used in determining the S/N were obtained from the NASA Exoplanet Archive (Akeson et al. 2013). If a system did not have a reported mass in the NASA Exoplanet Archive or it was not a 2σ result, we used the Weiss & Marcy (2014) Mass–Radius relationship to estimate the planet's mass. ^bOur calculation for the S/N of 55 Cnc e assumes a H/He envelope because it falls just above the pure rock line determined by Zeng et al. (2016). However, 55 Cnc e is in a USP orbit, making it unlikely that it would hold onto a thick H/He envelope. We do not include K2-266 b due to its grazing configuration.

Download table as:  ASCIITypeset image

Table 7.  Median Values and 68% Confidence Intervals for Planetary Parameters of EPIC 248435395 from EXOFASTv2

Planetary Parameters:		EPIC 248435395.01
P	Period (days)	${8.53467}_{-0.00063}^{+0.00073}$
RP	Radius ($\,{R}_{\oplus }$)	${10.6}_{-8.5}^{+2.7}$
TC	Time of conjunction (${\mathrm{BJD}}_{\mathrm{TDB}}$)	${2457907.9884}_{-0.0039}^{+0.0034}$
T0	Optimal conjunction time (${\mathrm{BJD}}_{\mathrm{TDB}}$)	${2457925.0578}_{-0.0028}^{+0.0025}$
a	Semimajor axis (au)	0.0684 ± 0.0012
i	Inclination (degrees)	${87.22}_{-0.23}^{+0.30}$
Teq	Equilibrium temperature (K)	${530.3}_{-7.7}^{+7.6}$
${R}_{P}/{R}_{* }$	Radius of planet in stellar radii	${0.147}_{-0.12}^{+0.042}$
$a/{R}_{* }$	Semimajor axis in stellar radii	${22.7}_{-1.4}^{+1.3}$
δ	Transit depth (fraction)	${0.022}_{-0.021}^{+0.014}$
Depth	Flux decrement at mid-transit	${0.00096}_{-0.00026}^{+0.00020}$
τ	Ingress/egress transit duration (days)	0.0154 ± 0.0016
T14	Total transit duration (days)	${0.0317}_{-0.0026}^{+0.0061}$
TFWHM	FWHM transit duration (days)	${0.0158}_{-0.0013}^{+0.0032}$
b	Transit impact parameter	${1.119}_{-0.14}^{+0.043}$
${\delta }_{S,3.6\mu {\rm{m}}}$	Blackbody eclipse depth at 3.6 μm (ppm)	${21}_{-20}^{+15}$
${\delta }_{S,4.5\mu {\rm{m}}}$	Blackbody eclipse depth at 4.5 μm (ppm)	${74}_{-71}^{+49}$
$\langle F\rangle$	Incident flux (109 erg s−1 cm−2)	0.0179 ± 0.0010
TP	Time of periastron (${\mathrm{BJD}}_{\mathrm{TDB}}$)	${2457907.9884}_{-0.0039}^{+0.0034}$
TS	Time of eclipse (${\mathrm{BJD}}_{\mathrm{TDB}}$)	${2457912.2558}_{-0.0036}^{+0.0031}$
TA	Time of ascending node (${\mathrm{BJD}}_{\mathrm{TDB}}$)	${2457905.8548}_{-0.0040}^{+0.0035}$
TD	Time of descending node (${\mathrm{BJD}}_{\mathrm{TDB}}$)	${2457910.1221}_{-0.0037}^{+0.0032}$
$d/{R}_{* }$	Separation at mid-transit	${22.7}_{-1.4}^{+1.3}$
PT	A priori non-grazing transit prob	${0.0387}_{-0.0039}^{+0.0043}$
${P}_{T,G}$	A priori transit prob	${0.0499}_{-0.0045}^{+0.0039}$
Wavelength Parameters:		Kepler
u1	Linear limb-darkening coeff	${0.35}_{-0.14}^{+0.16}$
u2	Quadratic limb-darkening coeff	${0.37}_{-0.14}^{+0.11}$
Transit Parameters:		dat UT p248-43-53 (Kepler)
${\sigma }^{2}$	Added Variance	${0.00000000090}_{-0.00000000054}^{+0.00000000057}$
F0	Baseline flux	1.0000000 ± 0.0000026

Download table as:  ASCIITypeset image

K2-266 b is a particularly interesting target for atmospheric follow-up because of its status as a USP exoplanet. Lopez (2017) suggests that, in order for USP planets to have radii larger than ∼2 R_⊕, they should have formed with high-metallicity, water-rich envelopes, and are likely to remain water-rich today. In addition, the theoretical models of Owen & Wu (2013) and Owen & Wu (2017) suggest that planets with the radial size and orbital period of K2-266 b reside near the boundary in parameter space where photoevaporation becomes an important source of mass loss. In addition, if K2-266 b has a typical magnetic field strength, its close proximity allows for interactions between the magnetospheres of the star and the planet (Adams 2011). Combinations of mass measurements of K2-266 b, refined radius measurements, and atmospheric constraints on water abundance could be used together to paint a complete picture of where in the disk this planet originated and when it reached its current-day location.

K2-266 d and e are two of the best sub-Neptunes for atmospheric characterization, and their longer periods (as compared to the others in Table 6) provide a valuable opportunity for a comparative atmospheric study between hot and warm sub-Neptunes. Interestingly, the Near-IR and IR brightness of K2-266 and EPIC 248435395 should allow for high-S/N observations using short exposure time for all four instruments on JWST: Near Infrared Camera (NIRCam), Near Infrared Imager and Slitless Spectrograph (NIRISS), Near-Infrared Spectrograph (NIRSpec), and the Mid-Infrared Instrument (MIRI) (Beichman et al. 2014; Kalirai 2018).

The dynamics of the K2-266 system is characterized by several remarkable features: The innermost planet (K2-266 b) is highly inclined relative to the rest of the planets, the orbits of planet candidate K2-266.02 and validated planet K2-266 c are in close proximity, the two sub-Neptunes (K2-266 d, e) are either in or extremely close to a mean-motion resonance, and the outer candidate K2-266.06 has a moderately eccentric orbit. Taken together, these factors place the planetary system orbiting K2-266 in a particularly unique realm.

Among the exoplanetary systems discovered thus far, only a small number have been determined to be in true resonance (Rivera et al. 2010; Lissauer et al. 2011; Carter et al. 2012; Barclay et al. 2013; Mills et al. 2016; Luger et al. 2017; Millholland et al. 2018; Shallue & Vanderburg 2018). From our numerical simulations initialized with orbital elements from the transit fit, in the stable simulations in which no planetary orbits cross, we find that planets K2-266 d and K2-266 e are in true resonance for 8.1% of the time, as characterized by a librating resonance angle. This significant fraction makes K2-266 another member of the short list of stars hosting systems containing potentially resonant exoplanets. Additional transits, which will improve orbital period precision, would enable future refinement of this resonance fraction.

From the time-evolution of the MEGNO indicator, we find that the average draw from the posterior becomes noticeably chaotic after roughly 100 years. The high mutual inclination between validated planet K2-266 b and the rest of the planets is also intriguing. As shown in Figure 10, because of the high present-day inclination of K2-266 b, this system is rarely (perhaps 5% of the time) in a configuration where all six planets can be seen simultaneously from our current line of sight. A smaller number of planets are expected to be seen in transit most of the time. Similarly, it is possible that the system hosts more than six planets, but we are seeing only six of them in transit at the current epoch. Tighter limits on the planetary posteriors will allow for a more precise determination of the future transit probability for each (known) planet, and will place constraints on any additional bodies in the system. In future work, a numerical survey of the parameter space subtended by the measured posteriors might allow for additional constraints on planetary parameters based on dynamical stability limits.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The formation of misaligned orbits, such as that of K2-266 b, remains an open problem. Petrovich et al. (2018) proposes that most USP planets form through non-linear secular interactions ("secular chaos"). In this scenario, the proto-USP planet starts with an orbital period of 5–10 days, is excited to high eccentricity, and is subsequently tidally captured onto a short-period orbit. Note that the process that leads to high eccentricity (e.g., planet–planet scattering) could also produce inclined orbits. As a result, one potential signature of this process could be a high mutual inclination for the USP planet relative to the other planets, as seen in this system. However, most companions to a USP planet generated in this manner would generally have orbital periods of 10 days or larger, whereas the K2-266 system has two planets with shorter periods of only six and seven days. For this secular chaos mechanism to form this system, the six- and seven-day planets would need to migrate inward after the eccentricity excitement and subsequent circularization of the USP planet; yet, as K2-266 b may have an unevaporated atmosphere based on current radius estimates, a dynamical history allowing it to form further from the star (subject to less photoevaporation) seems favorable.