Issue 
A&A
Volume 608, December 2017



Article Number  A35  
Number of page(s)  13  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201731558  
Published online  05 December 2017 
Characterization of the K218 multiplanetary system with HARPS^{⋆}
A habitable zone superEarth and discovery of a second, warm superEarth on a noncoplanar orbit
^{1} Dept. of Astronomy & Astrophysics, University of Toronto, 50 St. George Street, M5S 3H4, Toronto, ON, Canada
email: cloutier@astro.utoronto.ca
^{2} Centre for Planetary Sciences, Dept. of Physical & Environmental Sciences, University of Toronto Scarborough, 1265 Military Trail, M1C 1A4, Toronto, ON, Canada
^{3} Institut de Recherche sur les Exoplanètes, Département de physique, Université de Montréal, CP 6128 Succ. Centreville, H3C 3J7, Montréal, QC, Canada
^{4} Observatoire Astronomique de l’Université de Genève, 51 chemin des Maillettes, 1290 Versoix, Switzerland
^{5} Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
^{6} Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150762 Porto, Portugal
^{7} Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, 4169007 Porto, Portugal
Received: 12 July 2017
Accepted: 6 September 2017
Aims. The bright M2.5 dwarf K218 (M_{s} = 0.36 M_{⊙}, R_{s} = 0.41 R_{⊙}) at 34 pc is known to host a transiting superEarthsized planet orbiting within the star’s habitable zone; K218b. Given the superlative nature of this system for studying an exoplanetary atmosphere receiving similar levels of insolation as the Earth, we aim to characterize the planet’s mass which is required to interpret atmospheric properties and infer the planet’s bulk composition.
Methods. We have obtained precision radial velocity measurements with the HARPS spectrograph. We then coupled those measurements with the K2 photometry to jointly model the observed radial velocity variation with planetary signals and a correlated stellar activity model based on Gaussian process regression.
Results. We measured the mass of K218b to be 8.0 ± 1.9M_{⊕} with a bulk density of 3.3 ± 1.2 g/cm^{3} which may correspond to a predominantly rocky planet with a significant gaseous envelope or an ocean planet with a water mass fraction ≳50%. We also find strong evidence for a second, warm superEarth K218c (m_{p,c}sini_{c} = 7.5 ± 1.3 M_{⊕}) at approximately nine days with a semimajor axis ~ 2.4 times smaller than the transiting K218b. After reanalyzing the available light curves of K218 we conclude that K218c is not detected in transit and therefore likely has an orbit that is noncoplanar with the orbit of K218b although only a small mutual inclination is required for K218c to miss a transiting configuration;  Δi ~ 1−2°. A suite of dynamical integrations are performed to numerically confirm the system’s dynamical stability. By varying the simulated orbital eccentricities of the two planets, dynamical stability constraints are used as an additional prior on each planet’s eccentricity posterior from which we constrain e_{b} < 0.43 and e_{c} < 0.47 at the level of 99% confidence.
Conclusions. The discovery of the inner planet K218c further emphasizes the prevalence of multiplanet systems around M dwarfs. The characterization of the density of K218b reveals that the planet likely has a thick gaseous envelope which, along with its proximity to the solar system, makes the K218 planetary system an interesting target for the atmospheric study of an exoplanet receiving Earthlike insolation.
Key words: techniques: radial velocities / methods: statistical / planets and satellites: detection / planets and satellites: fundamental parameters / planets and satellites: individual: K218
Table A.2 is also available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/qcat?J/A+A/608/A35
© ESO, 2017
1. Introduction
Exoplanets orbiting within their host star’s habitable zone may have surface temperatures that allow for the presence of liquid water on their surfaces, depending on the properties of the planetary atmosphere (Kasting et al. 1993). The presence of liquid water is a condition likely required to sustain extraterrestrial life. This implies that habitable zone exoplanets receive comparable levels of stellar insolation to what the Earth receives from the Sun. Habitable zone exoplanets therefore represent superlative opportunities to search for life outside of the solar system via the characterization of their atmospheric structure and composition via transmission spectroscopy for transiting exoplanets.
M dwarf host stars are ideal targets to probe potentially habitable exoplanetary atmospheres (e.g. Kaltenegger et al. 2011; Rodler & LópezMorales 2014). Transmission spectroscopy observations of transiting habitable zone (HZ) exoplanets around M dwarfs are favourable compared to those around Sunlike stars given the increased depth of the transit for a planet of a given size (e.g. Stevenson et al. 2010; Kreidberg et al. 2014). In addition, the orbital periods corresponding to the HZ are less than those around Sunlike stars (weeks for those in the HZ of M dwarfs, compared to 12 months) thus increasing the number of accessible transit events within a given observational baseline. M dwarfs are also known to frequently host multiple small planets (typically 2.5 planets per star with 0.5 ≤ r_{p}/R_{⊕} ≤ 4 and within 200 days; Dressing & Charbonneau 2015; Gaidos et al. 2016) thus enabling direct comparative planetology to be conducted on known multiplanet systems.
Montet et al. (2015) reported the detection of the HZ superEarth K218b originally proposed in the K2 light curve analysis of ForemanMackey et al. (2015). In these studies, two transit events were observed in Campaign 1 data from the repurposed Kepler spacecraft mission K2 whose field coverage only lasted for ~ 80 days. The existence of the planet was confirmed and uncertainties regarding its ephemeris were significantly reduced in Benneke et al. (2017; hereafter B17) who used followup transit observations with the Spitzer Space Telescope to detect an additional transit event. The now confirmed superEarth K218b orbits an M2.5 dwarf with a period of ~ 32.9 days placing it directly within the star’s habitable zone (Kopparapu et al. 2013). The measured radius of 2.38 R_{⊕} is suggestive of an extended H/He envelope (Valencia et al. 2013; Rogers 2015; Fulton et al. 2017) that may contain additional molecular species such as water and/or methane that could be detectable with the James Webb Space Telescope (JWST; Beichman et al. 2014). Owing to the proximity of the system (~ 34 pc, V = 13.5, I = 11.7, K = 8.9; Cutri et al. 2003; Zacharias et al. 2013), K218 is truly an attractive target for characterizing the atmosphere of a HZ superEarth with unprecedented precision in the JWSTera.
In this study we report the first measurement of the planetary mass of K218b using precision radial velocity measurements taken with the HARPS spectrograph. In this data we also find strong evidence for an additional planet of similar minimum mass whose orbit is interior to K218b but is not found to transit. In Sect. 2 we summarize the HARPS spectroscopic and K2 photometric observations used in our analysis, in Sect. 3 we analyze the periodic signals in the spectroscopic data and in Sect. 4 we discuss our various radial velocity modelling procedures. In Sect. 5 we present the results of our radial velocity analysis including the detection of a second superEarth K218c in the system which we show is nontransiting and therefore is not perfectly coplanar with K218b in Sect. 6. Lastly we perform a dynamical analysis of the twoplanet system in Sect. 7 to dynamically constrain the orbital eccentricities of the planets before concluding with a discussion in Sect. 8.
2. Observations
2.1. HARPS spectra
From April 2015 (BJD = 2 457 117.5) to May 2017 (BJD = 2 457 875.5), we collected 75 spectra of K218 (EPIC 201912552) with the highresolution (R = 115 000) HARPS spectrograph (Mayor et al. 2003; Pepe et al. 2004). The majority of exposure times were fixed to 1800 s with the exception of the following six epochs whose exposure times were modified to the following reported values: 2400 s (BJD2 450 000 = 7199.503915, 7200.503114), 1200 s (BJD2 450 000 = 7204.491167), and 900 s (BJD2 450 000 = 7810.806284, 7814.760772, 7815.759421). The online HARPS pipeline returned the extracted, wavelengthcalibrated spectra (Lovis & Pepe 2007). Initial radial velocity estimates were computed from the crosscorrelation of each spectrum with a numerical mask (Baranne et al. 1996; Pepe et al. 2002). Using each spectrum’s initial estimate, all spectra were shifted to a common reference frame by their corresponding barycentric correction such that spectral features originating from the target star become aligned while telluric features are shifted by minus the epoch’s barycentric correction. The median combination of these shifted spectra was then used to construct a custom reference spectrum at high signaltonoise (S/N). A telluric template was then constructed from the median combination of all residual spectra after removal of the high S/N reference stellar spectrum. The process of computing the median reference stellar spectrum was then repeated using the individual spectra with tellurics masked by the median telluric spectrum. We then computed precision radial velocities by performing a χ^{2}minimization of each spectrum with the reference spectrum (AstudilloDefru et al. 2015). Radial velocity uncertainties were then estimated directly on the reference spectrum (Bouchy et al. 2001).
Fig. 1 K2 photometric light curve after the removal of known unphysical spurious signals. The two transits of K218b are highlighted by the ong blue ticks with the expected times of midtransit for K218c highlighted with short green ticks (see Sect. 5). The solid yellow curve is the mean of the predictive GP distribution and the surrounding shaded regions mark its 99% confidence intervals. The upper left subpanel is a magnified view of the highlighted region to aid in the visualization of the data and the GP fit. 
From the extracted spectra we also derived a number of activity indicators including the time series of the Hα index which is sensitive to chromospheric activity and is computed following the definition in Bonfils et al. (2007). For the M dwarf K218 (V = 13.5; Henden & Munari 2014) the Hα index is favoured over the Ca ii H+K Mt. Wilson S index (Wilson 1968; Baliunas et al. 1995) due to the low S/N obtained in the blue. From the S index we derived $\mathrm{log}{\mathit{R}}_{\mathrm{HK}}^{\mathrm{\prime}}\mathrm{=}\mathrm{}\mathrm{5.247}\mathrm{\pm}\mathrm{0.318}$ (AstudilloDefru et al. 2017). Additionally we derived the full width at half maximum (FWHM) and bisector inverse slope (BIS) shape parameters of the crosscorrelation function which are modified by dark and/or bright active regions traversing the visible stellar surface. In Sect. 3 we use these ancillary time series to learn about the star’s activity simultaneously with our radial velocity measurements. All spectroscopic time series are reported in Table A.2.
2.2. K2 photometry
K218 was observed in longcadence mode during Campaign 1 of the K2 mission as part of the “Targeting M dwarfs with K2” proposal (GO1053^{1}, PI: B. Montet). The baseline of the K2 light curve is just 80 days but provides nearly continuous coverage between June 1st, 2014 (BJD = 2 456 810.5) and August 20th, 2014 (BJD = 2 456 890.5).
We obtained the full detrended light curve from the MAST^{2} data retrieval service. As a result of the loss of two reaction wheels onboard the Kepler spacecraft, photometric observations from the K2 mission exhibit a reduced pointing precision, and hence photometric precision, compared to the original Kepler mission. Raw K2 light curves must be detrended with the variable pointing of the spacecraft throughout the observing sequence. We selected the EVEREST reduction of the K2 light curve which performs this detrending correction (Luger et al. 2016).
The majority of the residual photometric variability following detrending of the light curve can be attributed to the intrinsic photometric variability of the star and two observed transits of K218b from Montet et al. (2015). Removal of the transit events provides a dataset that can be used to investigate the correlated photometric activity resulting from active regions traversing the visible stellar surface thus giving rise to the star’s observed photometric variability. For reference, the detrended light curve is shown in Fig. 1 along with our Gaussian process fit to the light curve (see Sect. 4.1 for an explanation of the fit).
3. Periodogram analysis
Accurate modelling of the stellar radial velocity (RV) variations requires knowledge of the strong periodicities present in the data. These signals include contributions from both orbiting planets and from the rotation of active regions present on the stellar surface which give rise to a correlated stellar RV signal which is modulated by the stellar rotation period and/or its harmonics. In the top panel of Fig. 2 we plot the LombScargle periodogram (Scargle 1982) of the raw RVs to determine which periodicities are present at high significance, that is, those with a low false alarm probability (FAP). In all LSperiodograms we calculated FAPs via bootstrapping with replacement using 10^{4} iterations and individually normalize each periodogram’s power by its standard deviation.
Fig. 2 Left column (top to bottom): LombScargle periodograms of the raw radial velocities (RV), the K2 photometry (Phot), the HARPS window function (WF), S index, Hα, full width half maximum (FWHM), and bisector inverse slope (BIS) time series. The orbital periods of K218b, K218c, the stellar rotation period, and its first harmonic are highlighted with vertical dashed lines. Right column: the false alarm probability as a function of normalized periodogram power for each time series. The FAP curve for the photometry spans very low power and is barely discernible in its subpanel. Suffice it to say that any signal with normalized power >10^{2} has a FAP ≪ 0.1%. The FAP curve for the S index exhibits FAP ≳ 30% for all power visible on its ordinate and therefore is only visible in the upper right of its subpanel. 
Two important features are detected in the LSperiodogram of the RVs. The first is a forest of peaks ranging from ~ 25–45 days with distinct peaks centred on both the orbital period of K218b (P_{b} ~ 33 days B17; FAP = 2.9%) and the approximate stellar rotation period from the K2 photometry (see Sect. 4.1 for fitting of the stellar rotation period P_{rot}~ 38.6 days; FAP = 0.1%). The second important feature is a pair of closely spaced peaks at ~ 9 days (FAP < 0.01%) which because of their similar period and power likely result from a single source. This feature at approximately nine days constitutes the strongest periodic signal in the periodogram of the raw RVs and is not observed in the periodograms of any of the ancillary time series, nor in the periodogram of the window function, all of which are shown in the remaining panels of Fig. 2. The aforementioned time series include the K2 photometry (see Fig. 1), the window function or time sampling of the HARPS observations, and four spectroscopic activity indicators: the S index, Hα index, FWHM, and the BIS of the crosscorrelation function. Together the presence of the strong ~ 9 d. signal in radial velocity and its absence elsewhere provides strong initial evidence for a second planet in the K218 system at ~ 9 d.
4. Joint modelling of planets and correlated RV activity
4.1. Training the GP activity model on ancillary time series
The K2 photometry of K218 exhibits quasiperiodic photometric variability with a semiamplitude of ~ 0.008 mag and a rotation period of P_{rot}~ 38.6 days as seen in Fig. 1. This makes K218 a moderately active early M dwarf in terms of its photometric variability (Newton et al. 2016) the origin of which likely results from the rotation of active regions across the projected stellar disk at or close to P_{rot} owing to the characteristically low amplitudes of differential rotation in M dwarfs (Kitchatinov & Olemskoy 2011). The observed photometric variability – or photometric activity – has a correlated manifestation in the variation of the star’s apparent radial velocity and certain spectroscopic indicators because it is a single physical process that is responsible for the activity in each time series (Aigrain et al. 2012).
In order to obtain accurate and selfconsistent detections of the planetary signals in radial velocity we must model the RV activity signal of K218 simultaneously with our planet model. The photometric stellar rotation period of ~ 38 days (see second panel in Fig. 2) is marginally detected in the LSperiodogram of the RVs. However P_{rot} is clearly discernible by eye in the K2 light curve (Fig. 1) and has significant power in the LSperiodogram of the K2 light curve (second panel in Fig. 2) although the power is spread over a wide range of periodicities. Because of this we consider in our first model – called Model 1 – the K2 photometric light curve, less the observed transits of K218b, to train our RV activity model whose covariance properties are common with the observed photometric variability. However two important caveats had to be considered when adopting the K2 photometry as our training set. The first being that the baseline of the photometry spans just 80 days implying that any temporal variation whose characteristic timescale is greater than this baseline will remain unconstrained or at best weakly constrained. Secondly the K2 photometry were obtained nearly eight months prior to our HARPS observations such that any evolution in the covariance structure of the stellar activity between observing sequences from say magnetic activity cycles, would not be captured in the training set. For these reasons we also considered the BIS time series as an alternative training set in a second round of modelling called Model 2. Being contemporaneous with the RV measurements, training on the BIS time series mitigates the two aforementioned issues. In place of the BIS we also tested training on the S index, Hα, and FWHM time series but find results consistent with training on the BIS. Following Faria et al. (2016) we also considered a joint activity + planet model but neglected any training of the activity model’s covariance structure in a third model; Model 3. Finally for comparison purposes we also considered a fourth model – called Model 4 – that neglected any contribution from stellar activity.
To implement this joint modelling procedure we followed Cloutier et al. (2017) by using a Gaussian process (GP) regression model to model the covariance between adjacent observations in our training sets where applicable (i.e. in Models 1 and 2). GP regression is an attractive method for modelling the stochastic processes that gives rise to observable RV activity signals as it is nonparametric and therefore independent of an assumed functional form of the signal. The GP prior was represented by a multivariate Gaussian distribution of functions described by a covariance matrix ${\mathit{K}}_{\mathit{ij}}\mathrm{=}{\mathit{k}}_{\mathit{ij}}\mathrm{\left(}{\theta}\mathrm{\right)}\mathrm{+}{\mathit{\sigma}}_{\mathit{i}}^{\mathrm{2}}{\mathit{\delta}}_{\mathit{ij}}$ with a function k_{ij}(θ) = k(t_{i},t_{j},θ) that parameterizes the covariance between values of the observable y(t) at the epochs t_{i} and t_{j} in t. The observable y(t) has associated uncertainties σ(t) which were added along the diagonal of the covariance matrix K in quadrature. The set of GP hyperparameters θ are unique to the chosen covariance function k_{ij}(θ) and are solved for in the training step. After solving for the GP hyperparameters and thus obtaining a unique GP prior distribution, the GP prior conditioned on the data y(t) becomes the predictive distribution. The mean function of the GP predictive distribution can be evaluated at previously unseen epochs t^{∗} using ${\mu}\mathrm{\left(}{{t}}^{\mathrm{\ast}}\mathrm{\right)}\mathrm{=}\mathit{K}\mathrm{\left(}{{t}}^{\mathrm{\ast}}\mathit{,}{t}\mathrm{\right)}\mathrm{\xb7}\mathit{K}\mathrm{(}{t}\mathit{,}{t}{\mathrm{)}}^{1}\mathrm{\xb7}{y}\mathrm{(}{t}\mathrm{)}\mathit{,}$(1)which we took to be our GP activity model of the RVs by evaluating Eq. (1) at t.
Because the stellar activity, and in particular the longterm photometric variation, is modulated by the stellar rotation period, we included a periodic term in our assumed covariance function k_{ij}(θ) with period equal to P_{rot}. We also included a radial component due to the stochastic temporal evolution of starspot lifetimes, spatial distributions, and contrasts thus forcing the covariance to not be strictly periodic. Explicitly the adopted covariance structure is parameterized by a quasiperiodic covariance kernel of the form ${\mathit{k}}_{\mathit{i,j}}\mathrm{\left(}{\theta}\mathrm{\right)}\mathrm{=}{\mathit{a}}^{\mathrm{2}}\mathrm{exp}\left[\mathrm{}\frac{\mathrm{}{\mathit{t}}_{\mathit{i}}\mathrm{}{\mathit{t}}_{\mathit{j}}{\mathrm{}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\lambda}}^{\mathrm{2}}}\mathrm{}{\mathrm{\Gamma}}^{\mathrm{2}}{\mathrm{sin}}^{\mathrm{2}}\left(\frac{\mathit{\pi}\mathrm{}{\mathit{t}}_{\mathit{i}}\mathrm{}{\mathit{t}}_{\mathit{j}}\mathrm{}}{{\mathit{P}}_{\mathrm{GP}}}\right)\right]\mathit{,}$(2)which is parameterized by four hyperparameters θ = (a,λ,Γ,P_{GP}): a the amplitude of the correlations, λ the exponential timescale, Γ the coherence scale of the correlations, and P_{GP} the periodic timescale of the correlations which we interpret as P_{rot}. We also included an additional scalar jitter parameter s which is added in quadrature to the diagonal elements of the covariance matrix K such that θ becomes (a,λ,Γ,P_{GP},s).
Using the Markov chain MonteCarlo (MCMC) ensemble sampler emcee (ForemanMackey et al. 2013) we sample the marginalized posterior probability density functions (PDFs) of the five hyperparameters assuming uniform priors on the logarithm of each hyperparameter and maximizing the Gaussian logarithmic likelihood function $\mathrm{ln}\mathrm{\mathcal{L}}\mathrm{=}\mathrm{}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{{y}}^{\mathit{T}}{\mathit{K}}^{1}{y}\mathrm{+}\mathrm{ln}\mathrm{det}\mathit{K}\mathrm{+}\mathit{N}\mathrm{ln}{\mathrm{2}\mathit{\pi}}^{\mathrm{)}}\mathit{,}$(3)where y is the vector of N observations in the training set. In Model 1 y = the binned photometric data points^{3} whereas y = BIS time series in Model 2.
The MCMC was initialized with 200 walkers and hyperparameter values (a,λ,Γ,P_{GP}) = (max(y− ⟨ y ⟩), 10^{2} days,1,36 days). We sampled the logarithmic hyperparameters up to ≈10 autocorrelation times to ensure adequate convergence of the chains. We also monitored the acceptance fraction for each walker and insist that it lies within 20–50%. The sampling of each hyperparameter’s marginalized posterior PDF commences following a burnin phase of the same duration. The resulting marginalized and joint posterior PDFs are shown in Fig. 3 along with kernel density estimations of each marginalized 1D distribution. From the posterior PDF of P_{GP} we measured a stellar rotation period of P_{rot} = $\mathrm{38.}{\mathrm{6}}_{0.4}^{\mathrm{+}\mathrm{0.6}}$ days.
Fig. 3 Marginalized and joint posterior PDFs of the logarithmic GP hyperparameters used to model the K2 photometry shown in Fig. 1. Kernel density estimations of each model parameter’s posterior are overlaid on the histograms with solid green lines. 
4.2. Joint modelling of RVs
We proceeded with modelling the RVs jointly with Keplerian solutions for both K218b and c plus a trained quasiperiodic GP to model the correlated RV residuals attributed to stellar activity. The marginalized posterior PDFs of the GP hyperparameters λ,Γ, and P_{GP} from training were then used as informative priors in the joint RV analysis which treat the remaining GP hyperparameters a and s as free parameters. We sampled the informative priors using the kernel density estimations of each hyperparameter’s PDF obtained during training. This methodology allows the model to learn the covariance structure of the stellar activity through observations which are independent of planetary sources and then apply that knowledge to the joint modelling of the RVs thus distinguishing between stellar activity and planetinduced Doppler shifts.
The RV modelling is again performed using emcee. In Models 1, 2, and 3 our RV model consisted of 16 parameters including the five GP hyperparameters discussed in Sect. 4.1, the systemic velocity of K218 γ_{0}, the orbital periods of the two planets P, their times of inferior conjunction T_{0}, their RV semiamplitudes K, and the MCMC jump parameters $\mathit{h}\mathrm{=}\sqrt{\mathit{e}}\mathrm{cos}\mathit{\omega}$ and $\mathit{k}\mathrm{=}\sqrt{\mathit{e}}\mathrm{sin}\mathit{\omega}$ describing each planet’s orbital eccentricity e and argument of periastron ω. This parameterization was chosen to minimize the correlation between e and ω as well as reduce the tendency for the MCMC sampler to favour higheccentricity solutions (Ford 2006). In Model 4 we only considered 11 model parameters as no GP activity model was included. Table 1 summarizes the adopted priors on each RV model parameter in each of the models considered in this study. We adopted noninformative priors for all Keplerian parameters other than the orbital period and time of midtransit of K218b which were wellconstrained by the transit light curve modelling in B17.
Fig. 4 The marginalized and joint posterior PDFs of the model parameters from Model 1 of the observed RVs. Model 1 of the observed RVs models the two planets with Keplerian orbital solutions and the residual RV activity signal with a GP regression model trained on the K2 photometry in Fig. 1. Kernel density estimations of the trained posteriors are shown in the histograms of the logarithmic GP hyperparameters λ,Γ, and P_{GP} (Cols. 2–4). 
Summary of RV models and adopted priors.
5. Results
5.1. Results from RV data analysis
Here we compare results from the four considered RV models. Each model contains Keplerian solutions for each of the two planets. Additionally Model 1 models the RV residuals with a quasiperiodic GP regression model that is trained on the K2 photometry in which the stellar rotation period is clearly detected (see Sect. 4.1). In this model the stellar rotation period, and hence the GP periodic term, is sufficiently distinct from the orbital period of K218b such that the two signals are not confused in our joint modelling and the measured semiamplitude of K218b is not misestimated. Model 2 models the RV residuals with a quasiperiodic GP regression model that is trained on the BIS time series which is contemporaneous with the RVs. Model 3 models the RV residuals with an effectively unconstrained quasiperiodic GP and Model 4 neglects any modelling of the RV residuals therefore assuming that they are uncorrelated.
The maximum a posteriori (MAP) values of each model parameter along with the 16th and 84th percentiles of the marginalized posterior PDFs are reported in Table A.1. The marginalized and joint posterior PDFs of the model parameters used in Model 1 are shown in Fig. 4. In Sect. 5.2 we shall see that Model 1 is the best predictor of the observed RVs and therefore we only show the results from Model 1 in Fig. 4.
We emphasize that all Keplerian model parameters for the two planets are consistent at 1σ across all four models. Recall that the stellar rotation period is only wellconstrained in Model 1 via the K2 photometry and yet the measured semiamplitudes of K218b are consistent in each of the four models. This further demonstrates that there is minimal confusion between the RV signals at the stellar rotation period (38.6 days) and at the orbital period of K218b (~ 32.93 days) both of which are not distinctly detected in the periodogram of the raw RVs (see top panel of Fig. 2) but appear to be hidden within a forest of peaks spanning periodicities between ~ 25–45 days. The consistency of all measured Keplerian parameters in each model also suggests that the RV residuals, following the removal of the two MAP Keplerian solutions, are weakly correlated because nearly identical RV solutions are obtained with and without a GP treatment of the RV residuals following the removal of our planet models. That is that K218 appears to be a spectroscopically quiet star with the majority of its observed RV variation being attributable to planetary companions. Being spectroscopically quiet is promising for the prospect of transmission spectroscopy of K218b; an observation that is significantly complicated by the presence of stellar activity. The quiet nature of K218 is highlighted by its low measured value of $\mathrm{log}{\mathit{R}}_{\mathrm{HK}}^{\mathrm{\prime}}\mathrm{=}\mathrm{}\mathrm{5.247}$.
5.2. RV model comparison
A formal model comparison between the four considered models was performed using timeseries crossvalidation to compute the likelihood of each model given various training and testing subsets of the observed RVs (Arlot & Celisse 2010). We split the RVs into chronological training sets with sizes ranging from 20 measurements to the size of the full dataset less one (i.e. 74 measurements). The model parameters for each of the four considered models were optimized on the training set and the likelihood of the corresponding model is evaluated on the testing set. The testing set was simply the next observation chronologically following the final observation in the training set. The resulting median likelihood and median absolute deviation for each model is reported at the bottom of Table A.1 and was used to distinguish which of our four RV models performs optimally on the prediction of unseen RV measurements and thus best fits the data without overfitting. Through timeseries crossvalidation we find that Model 1 is the best predictor of the observed RVs. In the remainder of this study we consider the results from Model 1 to be the measured values of the planets K218b and c.
To confirm that we have detected a second planet K218c in our RV data, we performed a second round of timeseries crossvalidation calculations. In these calculations we compared three RV models each containing 0, 1, or 2 planets. We considered K218b to be the lone planet in the one planet model. In each model we also considered a GP activity model that was trained on the K2 photometry as was Model 1 above. Following the same methodology as previously discussed we find median logarithmic likelihoods of lnℒ_{0} = −2.693 ± 0.056, lnℒ_{1} = −2.642 ± 0.047, and lnℒ_{2} = −2.566 ± 0.026. From this we find that lnℒ_{2}−lnℒ_{1} = 0.076 ± 0.054 > 0 therefore arguing that the two planet model is the best predictor of unseen RV measurements and confirming that our two planet model containing both K218b and c is the RV model most favoured by the data.
The contributions to the observed RVs from stellar activity and each planet were depicted in Fig. 5. Together these physical models account for all significant periodicities in the observed RVs. In Fig. 5 we show the raw RVs as well as the RVs corrected for the individual RV component and compare each time series to its MAP model. In panel b we can see the relative importance of the GP activity model at modelling the RV residuals following the removal of our two planet model. Owing to the spectroscopically quiet nature of K218, the contribution to the observed RVs from activity is relatively small yet still holds a periodic manifestation at the stellar rotation period, albeit with a small amplitude. The residual rms following the removal of all modelled contributions is 2.89 m s^{1}. This value is less than the median photon noise limit of the measured RVs of 3.56 m s^{1} suggesting that we have modelled all significant RV contributions. For comparison, the residual rms achieved in Model 4, which neglects any red noise modelling following the removal of the two Keplerian solutions, is 3.16 m s^{1}. This value is also less than the median photon noise limit suggesting that the GP regression modelling alone in Models 1, 2, and 3 does not result in overfitting of the data.
Fig. 5 Panel a: Raw RVs less the systemic velocity of K218. Panel b: RV contribution from stellar activity; raw RVs corrected by the MAP Keplerian orbital solutions for each detected planet. Panel c: RV contribution from K218b; raw RVs corrected for activity and K218c. Panel d: RV contribution from K218c; raw RVs corrected for activity and K218b. Panel e: the RV residuals. The solid curves in panels b, c, and d depict the mean GP activity model, and the MAP Keplerian models for K218b and c) respectively. The surrounding shaded region in panel b is the 68% confidence interval on the mean GP model. All RV units are in m s^{1}. 
The Model 1 MAP Keplerian orbital solutions for K218b and c, but with eccentricities fixed to zero, are shown in Fig. 6. Here we report circular orbital solutions given that with these data we can only place upper limits on each planet’s eccentricity rather than detect it directly. The RV data are phasefolded to each planet’s MAP orbital period and time of inferior conjunction and are corrected for stellar activity using the mean GP activity model trained on the K2 photometry and assuming a mean model equal to the superposition of the two MAP Keplerian solutions.
Fig. 6 Phasefolded RVs for each planet in the K218 planetary system (top: K218c and bottom: K218b). The RVs have been corrected for stellar activity with a quasiperiodic GP model trained on the K2 photometry. The solid curves depict the maximum aposteriori Keplerian orbital solutions for each planet with fixed circular orbits. 
6. Searching for transits of K218c
From our RV analysis in Sect. 5 we derived the approximate linear ephemeris of K218c. We can therefore predict the passage of K218c at inferior conjunctions within the mostly continuous K2 photometric monitoring shown in Fig. 1. In Fig. 1 we indicate the nine such passages of K218c. Given the comparable minimum masses of K218b and c, it is reasonable to expect that the two planets also have comparable radii (recall r_{p,b} ~ 2.38 R_{⊕}). Furthermore (Ciardi et al. 2013) argued that Kepler multiplanet systems with planet radii ≲3 R_{⊕} do not exhibit a size – semimajor axis correlation such that the inner K218c is not expected to have undergone significant atmospheric escape compared to K218b. However the two 10σ transits of K218b are clearly discernible byeye in the K2 photometry whereas the predicted transits of K218c are not. This suggests that either K218c is much smaller than K218b such that its resulting transit depth is below the threshold for detection, or that the orbit K218c is mutually inclined with that of K218b such that it misses a transit configuration.
Here we attempt to confirm that K218c is indeed not transiting in the K2 light curves. To do so we perform an MCMC sampling of the K218c transit model parameters (P_{c}, T_{0,c}, r_{p,c}/R_{s}, a_{c}/R_{s}, and impact parameter b_{c}) using the K2 photometry and following the removal of the known transits of K218b. A quadratic limbdarkening law is assumed with fixed parameters in the Kepler bandpass: a = 0.3695 and b = 0.3570. These values are interpolated from the tables of Claret & Bloemen (2011) based on the known K218 surface gravity and effective temperature. In each MCMC step we compute the corresponding transit model using the batman implementation (Kreidberg 2015) of the Mandel & Agol (2002) transit model. We assumed a circular orbit of K218c and adopt the same MCMC methodology utilized on the RV data in Sect. 4.2. The orbital period and time of inferior conjunction (i.e. time of midtransit) were sampled from their joint RV posterior which maintains their apparent correlation (see Fig. 4). Priors on the scaled planet radius, and impact parameter are assumed uniform. In this way the scaled planetary radius is uncorrelated with its measured minimum mass and the impact parameter is constrained to be  b_{c}  < 1 as is required for a transit to occur.
Based on our MCMC analysis we find that the values of r_{p,c}/R_{s} are consistent with zero, that is, no transit is detected in the K2 data. Assuming the most likely value of R_{s} = 0.411 R_{⊙}, we calculate a planet radius upper limit of r_{p,c} < 0.52 R_{⊕} at 99% confidence assuming that K218c is indeed transiting. However, if K218c were this size and transiting, albeit unknowingly due to its small size, the planet would have a bulk density of ≳295 g cm^{3} or ≳54 ρ_{⊕}; an unphysically large value given the compressibility of pure iron. Thus we conclude that K218c is not transiting in the K2 data and is therefore noncoplanar with K218b despite having a smaller orbital separation.
To visualize the data a selection of light curve models are compared to the phasefolded K2 photometry in Fig. 7. Models shown include a suite of r_{p,c} values including its upper limit derived from MCMC (0.52 R_{⊕}), 0.75 R_{⊕}, 1 R_{⊕}, and the radius of K218b (2.38 R_{⊕}) which is detected at the 10σ level in the K2 photometry (Montet et al. 2015).
Fig. 7 K2 light curve of K218 phasefolded to the maximum aposteriori orbital period and time of midtransit of K218c from Model 1. Four transit light curve models are overplotted for various illustrative values of the assumed size of K218c: the 99% upper limit on r_{p,c} (0.52 R_{⊕}), 0.75 R_{⊕}, 1 R_{⊕}, and the size of K218b (2.38 R_{⊕}). No transit of K218c is detected in the data. 
7. Dynamical stability and eccentricity restrictions
The nondetection of K218c in transit (see Sect. 6) suggests that its orbital plane is not perfectly coplanar with the outer transiting K218b whose semimajor axis is ~ 2.4 times greater than K218c’s. The orbital inclination of K218b is $\mathrm{89.578}{\mathrm{5}}_{0.0088}^{\mathrm{+}\mathrm{0.0079}}$ deg with a corresponding impact parameter of $\mathrm{0.60}{\mathrm{1}}_{0.011}^{\mathrm{+}\mathrm{0.013}}$ (B17). In order for the orbit of K218c to not pass in front of its host star its orbital inclination must be tilted either ≳1.41° or ≲−2.25° from the orbit of K218b depending on which hemisphere of the stellar disk its transit chord will traverse. Such a mutual inclination is consistent with the peak in the distribution of Kepler multiplanet mutual inclinations (Figueira et al. 2012; Fabrycky et al. 2014). If indeed the planetary angular momentum vectors are within only a few degrees and therefore nearly aligned then we can analytically evaluate their Hill stability given estimates of their orbital eccentricities and assuming an inclination correction factor that is close to unity (Gladman 1993). If we assume the simplest case of initially circular orbits then the system is strongly Hill stable given that the two planets are currently separated by ~ 23 mutual Hill radii.
Fig. 8 Smoothed 2D map depicting the fraction of stable systems as a function of each planet’s eccentricity based on a suite of dynamical integrations. The dashed black curve depicts the analytic condition for Hill stability from Gladman (1993) assuming the MAP masses and semimajor axes from Model 1 in Table A.1. 1D histograms depict the number of stable and unstable (yellow and red respectively) systems in eccentricity bins of width 0.05 and marginalized over all other dynamical parameters. The annotated numbers report each bin’s stability fraction in percentages. 
Accurate orbital eccentricities of small planets with precision radial velocities are notoriously difficult to measure. For example, the change in RV semiamplitude of a circular K218b compared to an eccentricity of 0.1 is ≲2 cm s^{1} (~0.5% of K_{b}) or 15 cm s^{1} (~5% of K_{b}) for an eccentricity of 0.3. The aforementioned values are both at least an order of magnitude smaller than the characteristic RV uncertainty of the HARPS measurements presented in this work. Given that the system is Hill stable at small eccentricities we can use dynamical simulations to constrain the orbital eccentricities of the planets insisting that the system remain stable throughout its simulated evolution.
To constrain the planet eccentricities we performed a suite of 10^{4} dynamical integrations wherein we sample linearly each planet’s e ∈ [ 0,1). In each simulation the orbital inclination of K218b is drawn from while the system’s mutual inclination is drawn from such that the planet inclinations remain uncorrelated with the orbital eccentricities thus permitting an unbiased assessment of the system’s stability across the Keplerian parameter space. We insist that K218c be nontransiting at the start of each simulation by setting Δi_{min,c} to be the minimum mutual inclination required for  b_{c}  > 1 and rejecting draws for which this is not true. The dispersion in sampled mutual inclinations is tuned such that the mode of the resulting distribution lies within ~ 1−2° (Fabrycky et al. 2014). Each planet’s initial semimajor axis, true anomaly, and absolute mass is drawn from a Gaussian distribution with a mean value equal to the parameter’s MAP value from Model 1 in Table A.1 and with standard deviations equal to its average measurement uncertainty. The stellar mass is drawn from . The ascending node longitudes and arguments of periapsis are both drawn from . The system is then integrated forward in time from the epoch of the first K2 photometric observation (BJD = 2 456 810.26222) for 10^{6} yr using the WisdomHolman symplectic integrator WHFast (Rein & Tamayo 2015) implemented in the opensource REBOUND Nbody package (Rein & Liu 2012). These integrations are not intended to provide a comprehensive overview of the system’s dynamical stability but rather are useful to show that the system can remain stable up to at least 1 Myr and provide constraints on the planet eccentricities.
We classify stable integrations as those in which the minimum distance between the planets never becomes less than their mutual Hill radius. The fraction of stable systems as a function of each planet’s eccentricity and marginalized over all other dynamical parameters is shown in Fig. 8. Strong correlations between the fraction of stable systems and dynamical parameters other than planet eccentricities was not apparent so we focus here on the effect of eccentricities only. At small eccentricities there is a large stable region wherein the fraction of systems that remain stable is ≳80% and the system is known to be Hill stable based on the analytic criterion. As we increase either planet’s eccentricity the fraction of stable systems decreases. This is also illustrated by further marginalizing over planet eccentricities and considering the onedimensional representations of each system’s stability fraction in the histograms shown in Fig. 8.
The RV analysis discussed Sect. 5 and our dynamical simulations provide two independent methods for constraining the orbital eccentricities of the K218 planets. We can therefore combine these independent results by using the dynamical stability fractions shown in Fig. 8 as an additional prior on the ith planet’s derived eccentricity posterior: ${\mathit{e}}_{\mathit{i}}\mathrm{=}{\mathit{h}}_{\mathit{i}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathit{i}}^{\mathrm{2}}$. To do this we resample each planet’s RV eccentricity posterior and accept draws with a probability equal to the stability fraction at that drawn eccentricity value ±0.025. This choice of bin width was varied between 0.01 and 0.1 and was found not to have a significant effect on the results. In this way numerous random samples from each planet’s RV eccentricity posterior are rejected due to the low corresponding stability fraction. This is especially true for large eccentricities wherein the system no longer satisfies the Hill stability criterion. From the modified eccentricity posteriors we can calculate the 99th percentiles and find that e_{b} < 0.43 and e_{c} < 0.47 at that confidence level. These are the eccentricity values reported in Table A.1 and represent a more stringent evaluation of each planet’s eccentricity than considering the RV data alone.
Fig. 9 K218b along with a sample of other small exoplanets in the planetary mass and radius space. Overlaid curves are twocomponent interior structure models of fullydifferentiated solid planets with mass fractions annotated for each curve. 
8. Discussion
With a set of 75 precision radial velocity measurements taken with the HARPS spectrograph we have obtained a robust mass measurement of the transiting HZ superEarth K218b and detected a second superEarth K218c. The orbit of the newly discovered K218c lies interior to that of K218b and yet the planet is nontransiting. This implies that the orbital planes of the planets are mutually inclined. In order for K218c to not be seen intransit the planetary system requires a mutual inclination of just ≳1.4° which is consistent with the observed distribution of mutually inclined multiplanet systems (Figueira et al. 2012; Fabrycky et al. 2014). Dynamical simulations of the system revealed that the oscillation timescale of the planets’ orbital inclinations is suggesting that it may take many years before K218c reaches a transiting configuration. Although exactly how long depends sensitivity on its current inclination which remains unknown. The discovery of RV planets in transiting M dwarf planetary systems further emphasizes the prevalence of multiple Earth to superEarthsized planets around nearby M dwarfs and that these additional planets can be uncovered with moderate RV followup (Cloutier et al. 2017). Multiplanet systems such as K218 provide unique opportunities to study planet formation processes around M dwarfs via direct comparative planetology.
The presence of a second planet in the K218 transiting system will result in mutual planetary interactions thus making the orbit of K218b nonKeplerian and possibly resulting in an observable transit timing variation (TTV). Assuming a mutual inclination of K218c that just misses a transiting configuration, we estimate the expected TTVs of K218b using the TTVFaster package (Deck et al. 2014; Agol & Deck 2016). We adopted the maximum aposteriori masses and orbital periods of the two planets from Model 1 and uniformly sample their eccentricities up the 99% upper limits reported in Table A.1. The remaining orbital parameters of the planets that are unconstrained by the RV data are sampled uniformly between 0 and 2π. We find a maximum TTV for K218b of ~ 40 s which is slightly less than, but of the same order as the uncertainty in its measured time of midtransit (~ 50 s; B17). Thus with photometric monitoring of at least comparable quality to the K2 photometry shown in Fig. 1, detecting TTVs in the K218 multiplanet system is unlikely to provide any significant new insight into the nature of the system. Indeed no significant TTVs were observed in the Spitzer light curves from B17.
With our measured mass of K218b the planet joins a select group of HZ planets with constraints on both its mass and radius. This represents a significant step towards searching for potentially habitable planets around stars earlier than ~M 4 (Dittmann et al. 2017). With its maximum a posteriori mass of 7.96 ± 1.91 M_{⊕} the bulk density of K218b (ρ_{p,b} = 3.3 ± 1.2 g cm^{3}) lies between that of an Earthlike rocky planet and a low density Neptunelike planet. The planet is therefore likely too large to be a terrestrial Earthlike planet (Valencia et al. 2007; Fulton et al. 2017). Including K218b on the exoplanet massradius diagram in Fig. 9 we find that the internal structure of K218b is consistent with a range of twocomponent solidplanet models (Zeng & Sasselov 2013) owing to the uncertainty in its measured radius and mass which is at the level of ~ 24%. In particular, the 1σ lower mass limit of K218b permits it to be a complete “water world” even though its upper mass limit is consistent with a largely rocky interior surrounded by a significant mass fraction of water ice. In this parameter space the physical parameters of K218b are most similar to the superEarths HD 97658b (Van Grootel et al. 2014) and Kepler 10c (Dumusque et al. 2014; Rajpaul et al. 2017) despite receiving ~ 65 and ~ 24 times less insolation than those two planets respectively. Furthermore K218b is of a similar mass to the habitable zone planet LHS 1140b (Dittmann et al. 2017) and receives a comparable level of insolation despite being ~ 1.6 times larger than LHS 1140b. Analyzing the massradius relationship of these small planets over a range of equilibrium temperatures is a critical step towards understanding which of these systems have retained significant atmospheric content thus making them more suitable to extraterrestrial life.
Distinguishing between K218b as a pure waterworld or a scaledup Earth with a significant gaseous envelope will likely require transmission spectroscopy followup observations either with highresolution spectrographs from the ground or from space with JWST. With J = 9.8, H = 9.1, and K = 8.9 (Cutri et al. 2003), we stress that K218 is currently the second brightest M dwarf with a transiting habitable zone planet behind the recently discovered LHS 1140b. In the coming years the sample of habitable zone M dwarf planets is expected to increase dramatically following the launch of TESS (Ricker et al. 2014), although the majority of TESS planets will be more distant than LHS 1140 (Sullivan et al. 2015).
Considering the prospect of observational followup of K218b in transmission spectroscopy, if we consider an atmosphere that is cloudfree and dominated by hydrogen, then spectral features from wellmixed nearIR absorbing species such as water would have amplitudes of $\mathrm{\Delta}\mathit{F}\mathit{/}\mathit{F}\mathrm{~}\mathrm{10}\mathit{H}{\mathit{r}}_{\mathrm{p}}\mathit{/}{\mathit{R}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{~}\mathrm{230}$ ppm where H = k_{B}T_{eq}/μg is the atmospheric scale height, T_{eq} is the planet’s equilibrium temperature set by the stellar insolation, μ is the mean molecular weight of the atmosphere, and g is the surface gravity (MillerRicci et al. 2009; Kaltenegger & Traub 2009). If instead the atmosphere is dominated by heavier elements, similar Earth’s (e.g. N_{2} + O_{2},μ = 29), then the transmission signal will be significantly smaller (~ 10 ppm) though potentially still detectable with JWST with several visits. Because of the brightness of its host star and the low bulk density of K218b, the system offers a unique opportunity to study superEarth atmospheres receiving Earthlike insolation in the JWSTera.
Binning the K2 photometry in one day bins results in N = 78 compared to the 3439 unbinned photometric observations thus drastically increasing the computational efficiency of the evaluating Eq. (3) in each step of the Markov chains.
Acknowledgments
R.C. is partially supported in this work by the National Science and Engineering Research Council of Canada and the Institute for Research on Exoplanets. R.C. would like to thank Dan Tamayo for his useful discussions regarding the dynamical simulations presented in this work and their interpretation. X.B., X.D., and T.F. acknowledge the support of CNRS/PNP (Programme national de planétologie). D.E. acknowledges support from the National Centre for Competence in Research “PlanetS” of the Swiss National Science Foundation (SNSF) and from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 724427). X.B., J.M.A. and A.W. acknowledge funding from the European Research Council under the ERC Grant Agreement No. 337591ExTrA. N.C.S. acknowledges support from Fundação para a Ciência e a Tecnologia (FCT) through national funds and from FEDER through COMPETE2020 by the following grants: UID/FIS/04434/2013 & POCI–01–0145FEDER–007672 and PTDC/FISAST/1526/2014 & POCI–01–0145FEDER–016886. NCS also acknowledges the support from FCT through Investigador FCT contract IF/00169/2012/CP0150/CT0002. Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS526555. Support for MAST for nonHST data is provided by the NASA Office of Space Science via grant NNX09AF08G and by other grants and contracts.
References
 Agol, E., & Deck, K. 2016, ApJ, 818, 177 [NASA ADS] [CrossRef] [Google Scholar]
 Aigrain, S., Pont, F., & Zucker, S. 2012, MNRAS, 419, 3147 [NASA ADS] [CrossRef] [Google Scholar]
 Arlot, S., & Celisse, A. 2010, Statist. Surv., 4, 40 [Google Scholar]
 AstudilloDefru, N., Bonfils, X., Delfosse, X., et al. 2015, A&A, 575, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AstudilloDefru, N., Delfosse, X., Bonfils, X., et al. 2017, A&A, 600, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baliunas, S. L., Donahue, R. A., Soon, W. H., et al. 1995, ApJ, 438, 269 [NASA ADS] [CrossRef] [Google Scholar]
 Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beichman, C., Benneke, B., Knutson, H., et al. 2014, PASP, 126, 1134 [NASA ADS] [CrossRef] [Google Scholar]
 Benneke, B., Werner, M., Petigura, E., et al. 2017, ApJ, 834, 187 [NASA ADS] [CrossRef] [Google Scholar]
 Bonfils, X., Mayor, M., Delfosse, X., et al. 2007, A&A, 474, 293 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bouchy, F., Pepe, F., & Queloz, D. 2001, A&A, 374, 733 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ciardi, D. R., Fabrycky, D. C., Ford, E. B., et al. 2013, ApJ, 763, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Claret, A., & Bloemen, S. 2011, A&A, 529, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cloutier, R., Doyon, R., Menou, K., et al. 2017, AJ, 153, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, 2MASS All Sky Catalog of point sources [Google Scholar]
 Deck, K. M., Agol, E., Holman, M. J., & Nesvorný, D. 2014, ApJ, 787, 132 [CrossRef] [Google Scholar]
 Dittmann, J. A., Irwin, J. M., Charbonneau, D., et al. 2017, Nature, 544, 333 [NASA ADS] [CrossRef] [Google Scholar]
 Dressing, C. D., & Charbonneau, D. 2015, ApJ, 807, 45 [NASA ADS] [CrossRef] [Google Scholar]
 Dumusque, X., Bonomo, A. S., Haywood, R. D., et al. 2014, ApJ, 789, 154 [NASA ADS] [CrossRef] [Google Scholar]
 Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ, 790, 146 [NASA ADS] [CrossRef] [Google Scholar]
 Faria, J. P., Haywood, R. D., Brewer, B. J., et al. 2016, A&A, 588, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Figueira, P., Marmier, M., Boué, G., et al. 2012, A&A, 541, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ford, E. B. 2006, ApJ, 642, 505 [NASA ADS] [CrossRef] [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [CrossRef] [Google Scholar]
 ForemanMackey, D., Montet, B. T., Hogg, D. W., et al. 2015, ApJ, 806, 215 [NASA ADS] [CrossRef] [Google Scholar]
 Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Gaidos, E., Mann, A. W., Kraus, A. L., & Ireland, M. 2016, MNRAS, 457, 2877 [NASA ADS] [CrossRef] [Google Scholar]
 Gladman, B. 1993, Icarus, 106, 247 [NASA ADS] [CrossRef] [Google Scholar]
 Gregory, P. C. 2005, ApJ, 631, 1198 [NASA ADS] [CrossRef] [Google Scholar]
 Henden, A., & Munari, U. 2014, Contributions of the Astronomical Observatory Skalnate Pleso, 43, 518 [Google Scholar]
 Kaltenegger, L., & Traub, W. A. 2009, ApJ, 698, 519 [NASA ADS] [CrossRef] [Google Scholar]
 Kaltenegger, L., Segura, A., & Mohanty, S. 2011, ApJ, 733, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Kitchatinov, L. L., & Olemskoy, S. V. 2011, MNRAS, 411, 1059 [NASA ADS] [CrossRef] [Google Scholar]
 Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, ApJ, 765, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Kreidberg, L. 2015, PASP, 127, 1161 [NASA ADS] [CrossRef] [Google Scholar]
 Kreidberg, L., Bean, J. L., Désert, J.M., et al. 2014, Nature, 505, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Lovis, C., & Pepe, F. 2007, A&A, 468, 1115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Luger, R., Agol, E., Kruse, E., et al. 2016, AJ, 152, 100 [Google Scholar]
 Mandel, K., & Agol, E. 2002, ApJ, 580, L171 [NASA ADS] [CrossRef] [Google Scholar]
 Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20 [NASA ADS] [Google Scholar]
 MillerRicci, E., Seager, S., & Sasselov, D. 2009, ApJ, 690, 1056 [NASA ADS] [CrossRef] [Google Scholar]
 Montet, B. T., Morton, T. D., ForemanMackey, D., et al. 2015, ApJ, 809, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Newton, E. R., Irwin, J., Charbonneau, D., et al. 2016, ApJ, 821, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Pepe, F., Mayor, M., Galland, F., et al. 2002, A&A, 388, 632 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pepe, F., Mayor, M., Queloz, D., et al. 2004, A&A, 423, 385 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rajpaul, V., Buchhave, L. A., & Aigrain, S. 2017, MNRAS, 471, L125 [NASA ADS] [CrossRef] [Google Scholar]
 Rein, H., & Liu, S.F. 2012, A&A, 537, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rein, H., & Tamayo, D. 2015, MNRAS, 452, 376 [NASA ADS] [CrossRef] [Google Scholar]
 Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, in Space Telescopes and Instrumentation 2014: Optical, Infrared, and Millimeter Wave, Proc. SPIE, 9143, 914320 [Google Scholar]
 Rodler, F., & LópezMorales, M. 2014, ApJ, 781, 54 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, L. A. 2015, ApJ, 801, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Scargle, J. D. 1982, ApJ, 263, 835 [NASA ADS] [CrossRef] [Google Scholar]
 Stevenson, K. B., Harrington, J., Nymeyer, S., et al. 2010, Nature, 464, 1161 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Sullivan, P. W., Winn, J. N., BertaThompson, Z. K., et al. 2015, ApJ, 809, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Valencia, D., Sasselov, D. D., & O’Connell, R. J. 2007, ApJ, 665, 1413 [NASA ADS] [CrossRef] [Google Scholar]
 Valencia, D., Guillot, T., Parmentier, V., & Freedman, R. S. 2013, ApJ, 775, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Van Grootel, V., Gillon, M., Valencia, D., et al. 2014, ApJ, 786, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Wilson, O. C. 1968, ApJ, 153, 221 [NASA ADS] [CrossRef] [Google Scholar]
 Zacharias, N., Finch, C. T., Girard, T. M., et al. 2013, AJ, 145, 44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zeng, L., & Sasselov, D. 2013, PASP, 125, 227 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Additional tables
Model parameters.
HARPS time series.
All Tables
All Figures
Fig. 1 K2 photometric light curve after the removal of known unphysical spurious signals. The two transits of K218b are highlighted by the ong blue ticks with the expected times of midtransit for K218c highlighted with short green ticks (see Sect. 5). The solid yellow curve is the mean of the predictive GP distribution and the surrounding shaded regions mark its 99% confidence intervals. The upper left subpanel is a magnified view of the highlighted region to aid in the visualization of the data and the GP fit. 

In the text 
Fig. 2 Left column (top to bottom): LombScargle periodograms of the raw radial velocities (RV), the K2 photometry (Phot), the HARPS window function (WF), S index, Hα, full width half maximum (FWHM), and bisector inverse slope (BIS) time series. The orbital periods of K218b, K218c, the stellar rotation period, and its first harmonic are highlighted with vertical dashed lines. Right column: the false alarm probability as a function of normalized periodogram power for each time series. The FAP curve for the photometry spans very low power and is barely discernible in its subpanel. Suffice it to say that any signal with normalized power >10^{2} has a FAP ≪ 0.1%. The FAP curve for the S index exhibits FAP ≳ 30% for all power visible on its ordinate and therefore is only visible in the upper right of its subpanel. 

In the text 
Fig. 3 Marginalized and joint posterior PDFs of the logarithmic GP hyperparameters used to model the K2 photometry shown in Fig. 1. Kernel density estimations of each model parameter’s posterior are overlaid on the histograms with solid green lines. 

In the text 
Fig. 4 The marginalized and joint posterior PDFs of the model parameters from Model 1 of the observed RVs. Model 1 of the observed RVs models the two planets with Keplerian orbital solutions and the residual RV activity signal with a GP regression model trained on the K2 photometry in Fig. 1. Kernel density estimations of the trained posteriors are shown in the histograms of the logarithmic GP hyperparameters λ,Γ, and P_{GP} (Cols. 2–4). 

In the text 
Fig. 5 Panel a: Raw RVs less the systemic velocity of K218. Panel b: RV contribution from stellar activity; raw RVs corrected by the MAP Keplerian orbital solutions for each detected planet. Panel c: RV contribution from K218b; raw RVs corrected for activity and K218c. Panel d: RV contribution from K218c; raw RVs corrected for activity and K218b. Panel e: the RV residuals. The solid curves in panels b, c, and d depict the mean GP activity model, and the MAP Keplerian models for K218b and c) respectively. The surrounding shaded region in panel b is the 68% confidence interval on the mean GP model. All RV units are in m s^{1}. 

In the text 
Fig. 6 Phasefolded RVs for each planet in the K218 planetary system (top: K218c and bottom: K218b). The RVs have been corrected for stellar activity with a quasiperiodic GP model trained on the K2 photometry. The solid curves depict the maximum aposteriori Keplerian orbital solutions for each planet with fixed circular orbits. 

In the text 
Fig. 7 K2 light curve of K218 phasefolded to the maximum aposteriori orbital period and time of midtransit of K218c from Model 1. Four transit light curve models are overplotted for various illustrative values of the assumed size of K218c: the 99% upper limit on r_{p,c} (0.52 R_{⊕}), 0.75 R_{⊕}, 1 R_{⊕}, and the size of K218b (2.38 R_{⊕}). No transit of K218c is detected in the data. 

In the text 
Fig. 8 Smoothed 2D map depicting the fraction of stable systems as a function of each planet’s eccentricity based on a suite of dynamical integrations. The dashed black curve depicts the analytic condition for Hill stability from Gladman (1993) assuming the MAP masses and semimajor axes from Model 1 in Table A.1. 1D histograms depict the number of stable and unstable (yellow and red respectively) systems in eccentricity bins of width 0.05 and marginalized over all other dynamical parameters. The annotated numbers report each bin’s stability fraction in percentages. 

In the text 
Fig. 9 K218b along with a sample of other small exoplanets in the planetary mass and radius space. Overlaid curves are twocomponent interior structure models of fullydifferentiated solid planets with mass fractions annotated for each curve. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.