EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 631 (November 2019) A&A, 631 (2019) A90 Full HTML
Open Access
Issue
A&A
Volume 631, November 2019
Article Number A90
Number of page(s) 17
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201936267
Published online 28 October 2019

© T. A. Lopez et al. 2019

Licence Creative Commons
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 Introduction

Precise knowledge of both the mass and radius of planets is necessary before initiating more advanced studies. This is particularly true for internal structure modelling which, beyond the degeneracies inherent to this type of investigation, requires a precise determination of the planetary bulk density (Seager et al. 2007; Brugger et al. 2017). This is also true for atmospheric studies which rely on mass measurements to derive the scale height, which is an essential parameter to interpret observations in transmission spectroscopy (e.g. Seager et al. 2009). Few low-mass planets have both a precise radius and mass. One of the reasons is that the Kepler mission focussed on faint stars (V > 13), prohibiting subsequent Doppler monitoring from the ground for most of them.

Understanding the compositions of small planets allows us to put strong constraints on their formation path, and even helps to understand the formation of the solar system (e.g. Raymond et al. 2018). For example, Kepler provides evidence for a bi-modality in the radius distribution of small planets, with a gap around 1.6 R⊕ (Rogers 2015; Fulton et al. 2017; Fulton & Petigura 2018). Moreover, there are several examples of planets with a similar radius but radically different compositions: Kepler-11f and Kepler-10c have a radius of 2.61 ± 0.25 R⊕ and 2.35 ± 0.05 R⊕ respectively,whereas the first is gaseous and the second is rocky (Lissauer et al. 2011; Dumusque et al. 2014). This has also been observed for planets within the same system (Bonomo et al. 2019) and planets with sizes in the range of 1.5 R⊕–4 R⊕ can have compositions ranging from gaseous (dominated by H–He) to rocky (dominated by iron and silicates), and even to icy (water-worlds with a large amount of volatiles) (Léger et al. 2004; Rogers et al. 2011; Rogers 2015). The transition between rocky and gaseous planets is not thoroughly understood and more well-characterised low-mass planets are required. The K2 mission provides us with a large number of small planet candidates around relatively bright stars, allowing precise spectroscopic follow-up from the ground to determine their masses (e.g. Osborn et al. 2017; Barros et al. 2017; Santerne et al. 2018; Lam et al. 2018).

Among these, some are multiple systems of small transiting planets. They are really valuable as they highlight differences in planetary formation and evolution within the same environment, allowing us to perform comparative studies (Bonomo et al. 2019). These systems are often compact and coplanar (Lissauer et al. 2011; Fabrycky et al. 2014; Winn & Fabrycky 2015) and bring further constraints to the evolution and migration models. Indeed, they are believed to have migrated to their current position rather than formed in-situ (e.g. Ormel et al. 2017). In addition, some of these tightly-packed systems present resonances, which are a consequence of inward migration when drifting planets pile-up in or near mean-motion resonances, thus creating chains of resonances (Raymond et al. 2008; Fabrycky et al. 2014; Ormel et al. 2017).

This is the case of the K2-138 planetary system. K2-138 is a moderately bright star (V = 12.2) observed in the campaign 12 of the K2 mission. Four planets were initially identified by citizen scientists involved in the Exoplanet Explorers project1. A detailed analysis of the light curve revealed the presence of potentially six planets. The first five form a compact inner system with radii ranging from 1.6 R⊕ to 3.3 R⊕ and periods between 2.35 and 12.76 days (Christiansen et al. 2018). For the outer planet, only two transits are available in the K2 data, indicating a potential period of 42 days. Additional Spitzer observations confirm its planetary nature (Hardegree-Ullman et al., in prep.). The five inner planets form a chain of near first-order 3:2 mean-motion resonance, and all the planets lie just outside these resonances. This is the longest known chain to date (Fabrycky et al. 2014). The five inner planets also are in a chain of three-body Laplace resonances (Christiansen et al. 2018), which is similar to the TRAPPIST-1 system (Luger et al. 2017).

The brightness of K2-138 makes it an interesting target for ground-based spectroscopic follow-up. The precise mass measurement of these small planets is within reach of an instrument like the High Accuracy Radial velocity Planet Searcher (HARPS) on a 3.6 m telescope. In this paper, we present radial velocity (RV) observations acquired in 2017 and 2018, allowing us to constrain precisely the masses of three of the six planets in the system, and to put strong constraints on the others. We used these measurements to predict transit timing variations (TTVs) for the purpose of observing them with the CHaracterizing ExOPlanet Satellite (CHEOPS). The observations and data reduction are described in Sect. 2. The spectral analysis, activity characterisation, and the joint analysis of the radial velocities and photometry are reported in Sect 3. Finally, we discuss the validity of the results and the magnitude of transit timing variations in Sect. 4.

2 Observations and data reduction

2.1 K2 Photometry

K2-138 was observed for 79 days during the campaign 12 of the K2 mission (Howell et al. 2014) between December 15, 2016 and March 04, 2017 in long cadence mode. On February 01, 2017 the spacecraft entered safe mode resulting in a 5.3 days gap. Five distinct transiting signals were identified and analysed as described in Christiansen et al. (2018) leading to the discovery of five sub-Neptune planets with sizes ranging from 1.6 R⊕ to 3.3 R⊕ in a chain of near 3:2 resonances. Additionally, a potential sixth planet was identified with a 42 days period. Hardegree-Ullman et al. (in prep.) were able to confirm it using observations from the Spitzer Space Telescope. False-positive scenarios are discussed and ruled out in Christiansen et al. (2018). A sketch of the system is shown in Fig. 1 along with the planet distances to the resonances.

Our group independently discovered K2-138 using a new transit search routine that was optimised for multi-planetary systems. The code is based on a Python re-implementation2 of the Box-fitting least squares (BLS; Kovács et al. 2002) algorithm. To search for multi-planetary systems, the code repeats the BLS transit search several times per light curve after removing the model of transits found in the previous BLS search. For each campaign we usually apply the BLS three times to all the light curves reduced with the POLAR pipeline. The candidates that pass the vetting procedure were re-analysed with our transit searching code with a higher number of BLS calls to allow the discovery of more planets in the system. For K2-138 the first search revealed two candidates with periods 5.4 days and 8.3 days. Other transits were also clearly seen in the light curve by eye. For the second search we performed the BLS six times and found two clearer candidates with periods 3.6 days and 2.35 days and a possible candidate with a period of 12.8 days. As a very interesting planetary transiting system, we started the RV follow-up with HARPS early on after its identification by our group.

We used the Planet candidates from OptimaL Aperture Reduction (POLAR) pipeline (Barros et al. 2016) to reduce the pixel data and extract the light curve. The POLAR pipeline is summarised as follows: We downloaded the calibrated pixel data (pixel files) from the Mikulski Archive for Space Telescopes (MAST)3 and performed a photometric extraction using the adapted CoRoT imagette pipeline from Barros et al. (2014), which uses an optimal aperture following the point spread function. The centroid positions were then computed using the Modified Moment Method (Stone 1989) to correct the flux-positionsystematics. This correction was performed using a self-flat-fielding method similar to Vanderburg & Johnson (2014). The first transit (planet d) was removed due to the significantly higher noise level at the very beginning of the campaign. We also used the open source freely available4 EPIC Variability Extraction and Removal for Exoplanet Science Targets (EVEREST) pipeline (Luger et al. 2016, 2018) on the transit-filtered light curve. EVEREST rely on an extension of the pixel level decorrelation used on Spitzer data (Deming et al. 2015). We used both the POLAR and EVEREST light curves for a complete combined analysis with the radial velocities. We report the results from the EVEREST reduction as it has slightly better K2 residuals and allows for a better convergence of the analysis. Nevertheless, they are consistent with the POLAR ones. Figure 2 shows the extracted light curve from EVEREST. We also corrected the K2 light curve for spot and faculae modulation: the light curve was first detrended using a second order polynomial regression to remove low-frequency variations, outliers were removed with a sigma-clipping and the activity was filtered as described in Sect. 3.3. Finally, we performed the transit fitting jointly with the radial velocity analysis as described in details in Sect. 3.3.

thumbnail Fig. 1

K2-138 planetary system. All distances are drawn to scale. Planetary radii are enlarged by a factor 50 for readability. Orbits, assumed circular, are shown in black whereas the location of the 3:2 resonances are shown in dashed green.
Open with DEXTER

2.2 HARPS radial-velocity follow-up

We collected 215 spectra of K2-138 over 79 nights between September 25, 2017 and September 04, 2018 with the HARPS spectrograph (R ~115 000) mounted on the 3.6 m Telescope at ESO La Silla Observatory (Mayor et al. 2003). These observations were conducted as part of the ESO-K2large programme (ID 198.C-0169) with the aim of detecting the masses of the five transiting planets known at this time. We used an exposure time of 1800s, giving a mean signal-to-noise ratio (S/N) of 28.5 per pixel at 550 nm. The spectra were reduced using the HARPS pipeline. Radial velocities (RVs) were derived by cross-correlating the spectra with a numerical template of a G2V star (Baranne et al. 1996; Pepe et al. 2002) and the photon noise was estimated as described in Bouchy et al. (2001). The average photon noise is 3.92 m s−1. Finally, the full width half maximum (FWHM) and the bisector inverse slope (BIS) of the cross-correlation functions were measured and their uncertainties computed as described in Santerne et al. (2015). We rejected 21 measurements that were affected by Moon contamination, showingsignificant anomalies in their RV and FWHM. The remaining 194 measurements were used for the analysis described in Sect. 3.3 and are reported in Table A.3. We also measured indices sensitive to stellar activity: Hα, Na I D and Ca II from which we derived SMW. The periodograms associated to the RVs and different indices are given in Fig. A.1.

2.3 Gaia astrometry and contamination

We used the Gaia Data Release 2 (Gaia Collaboration 2018) parallax, taking into account the correction from Schönrich et al. (2019), to derive the stellar distance. K2-138 has a corrected parallax of 4.962 ± 0.065 mas, leading to a distance of 201.66 ± 6.38 pc. Table 1 gives the astrometric properties of K2-138. Andrae et al. (2018) derived the effective temperature from the three-band photometry and the parallax, giving K. They also used the extinction and the reddening to determine the luminosity and radius . The stellar parameters from Gaia DR2 are consistent with the distance estimate, effective temperature and stellar radius which are derived in the joint Bayesian analysis in Sect. 3.3 and via the spectral analysis. In complement to the planetary transit validation in Christiansen et al. (2018), we also checked the presence of contaminants resolved by Gaia. There is only one source already reported in Christiansen et al. (2018), within 14′′ with a G mag of 20.1, giving Δm = 8.1. The contamination is therefore negligible, at most at the level of 10−3, assuming the contaminant falls within the K2 aperture.

3 Data analysis and results

3.1 Spectral analysis of the host star

The spectral analysis of the host star was performed to derive its fundamental parameters and the chemical compositionof the photosphere. We used the HARPS spectra corrected from systemic velocity and planetary reflex-motion. We removed the spectra with a S/N lower than ten in order 47 (550 nm) and the ones contaminated by the Moon (S/N above 1.0 in fibre B). We then co-added the spectra in a single 1D spectrum. Then, the spectral analysis was performed using two sets of tools.

Firstly, we proceeded as described in Deleuil et al. (2012, 2014). The single 1D spectrum was carefully normalised. The spectral analysis was performed using the Versatile Wavelength Analysis (VWA) package (Bruntt et al. 2010). The spectral parameters were determined by fitting the Iron lines until the derived Fe I and Fe II abundances minimised the correlation with both the equivalent width and the excitation potentials. The surface gravity was determined from the pressure-sensitive lines Mg I b, Na I D and the calcium lines at 612.2 nm, 616.2 nm and 643.9 nm. We found log g = 4.52 ± 0.15 [cgs], Teff = 5350 ± 80 K, υmicro = 0.90 ± 0.10 km s−1, υmacro = 1.9 ± 0.1 km s−1, υ sin i⋆ = 2.5 ± 1.0 km s−1, [M/H] = 0.15 ± 0.04, and [Fe/H] = 0.14 ± 0.10. The abundances of elements with isolated lines were derived. They are shown in Table A.1. For the Lithium, we do not find any significant Li I line. We computed the stellar age Gyr using the calibration with the chromospheric emission (Donahue 1993; Wright et al. 2004). The age derived from the joint analysis (Sect 3.3) is Gyr. Both are consistent. The final stellar parameters are given in Table 1. They are all fully consistent with those obtained within the joint analysis with PASTIS, as described in Sect 3.3. They are also consistent with the results in Christiansen et al. (2018).

As an independent check of the stellar parameters, we employed the same tools as used for stars in the SWEET-Cat catalogue (Santos et al. 2013). The method is in essence the same as in the first part (equivalent widths of Fe I and Fe II lines, and iron excitation and ionisation equilibrium). Local thermodynamical equilibrium is assumed in the framework of the 2014 version of the code MOOG (Sneden 1974), together with the grid of Kurucz ATLAS9 plane-parallel model atmospheres (Kurucz 1993) and the equivalent width measurements of the iron lines from the code ARES (Sousa et al. 2015). The derived values are log g = 4.327 ± 0.069 [cgs], Teff = 5277 ± 37 K, υmicro = 0.838 ± 0.060 km s−1 and [Fe/H] = 0.069 ± 0.024. The log g value, known to be systematically biased, was corrected using a calibration based on asteroseismic targets (Mortier et al. 2014). The corrected value is log g = 4.51 ± 0.07 [cgs]. All these values agree with the previous ones, within uncertainties.

thumbnail Fig. 2

K2-138 light curve and associated models. Top: extracted K2 light curve, using the EVEREST pipeline, of K2-138 with in red the Gaussian process regression trained on the out-of-transit light curve. The positions of transits are marked with coloured lines at the bottom of the top figure. Bottom: light curve after subtraction of the Gaussian process fit. The best six-planet photometric model of the transits is superimposed in red.
Open with DEXTER
thumbnail Fig. 3

HARPS radial velocity time series obtained in two seasons (in red). The six-planets Keplerian model and the Gaussian process regression are shown in black with a 68.3% confidence interval in grey.
Open with DEXTER

3.2 Activity and stellar rotation

The stellar rotation period was first derived using the K2 light curve. We computed the autocorrelation function (ACF; Fig. 4) as described in McQuillan et al. (2013). It shows a dominant periodicity at 24.44 ± 2.35 d. Then, we trained a Gaussian process (GP) with a quasi-periodic kernel (Eq. (1) in Sect. 3.3) on the six hours-binned, transit-filtered light curve. We used a uniform prior on the GP period, between 15 and 40 days. All the priors are listed in Table A.4. We found a period of d.

We also used the spectroscopic data to constrain the stellar rotation. The BIS periodogram (Fig. A.1) shows a peak with a false alarm probability (FAP) below 10% at around 12.5 d, which is close to half the rotation period. The Hα index and FWHM periodogram (Fig. A.1) show a significant peak at 25 days of period, below 0.1% FAP. In contrast, there is no clear signal in the SMW time series. We also found a decreasing linear drift in the FWHM, Hα and SMW times series (Fig. A.1). This may be indicative of a magnetic cycle (Lovis et al. 2011a). We then trained a GP on the FWHM, BIS, Hα, Na I D and SMW. The priors and posteriors of all the hyperparameters are reported in Table A.4. We obtained loose constraints of respectively d, 22 ± 13 d, d, 22 ± 12 d, and 23 ± 12 d.

Then, we derived the stellar rotation period from the index and the B−V colour as described in Mamajek & Hillenbrand (2008). We computed the from the SMW measurements and the APASS B−V colour (Henden et al. 2011), following the calibrations from Noyes et al. (1984) and the correction from Lovis et al. (2011b) assuming an Iron abundance of 0.14 ± 0.10 and a B−V = 0.839 ± 0.062 corrected from the extinction. This led to an averaged value of . We deduced from it the Rossby number R0 = 1.51 ± 0.11 (Mamajek & Hillenbrand 2008). From the B−V colour, we derived the convective turnover time τc = 20.6 ± 1.7 d (Noyes et al. 1984), and using the relation Prot = R0 × τc, we derived Prot = 31 ± 4 d. This value is compatible, within errors, with the ones derived above.

Finally, we derived the projected rotational velocity of K2-138 using calibrations from the FWHM. We obtained υ sin i⋆ = 2.2 ± 1.2 km s−1 which we combined with the stellar radius from Gaia DR2 (Sect. 2.3) to derive an upper limit on the stellar rotation period Prot = 21.3 ± 3.2 d, assuming .

In the end, as most of the indicators tend to show a stellar rotation period at 25 days, we modelled the activity-induced RVs using a Gaussian process regression as described in the next section with Prot = 24.7 ± 2.2 d as prior of the stellar rotation period. This prior is drawn from the symmetrised value obtained with the GP learning on the light curve alone, taking a conservative width.

Table 1

Stellar properties of K2-138.

thumbnail Fig. 4

Autocorrelation function of K2-138 light curve. The dashed line in black shows the dominant period, at 24.4 days, corresponding to the stellar rotation. It was computed on the smoothed ACF, following recommendations from McQuillan et al. (2013).
Open with DEXTER

3.3 PASTIS analysis

Similarly to previous studies (Osborn et al. 2017; Barros et al. 2017; Santerne et al. 2018; Lam et al. 2018), we used the Bayesian software PASTIS (Díaz et al. 2014; Santerne et al. 2015) to jointly analyse the HARPS radial velocities, the K2 light curve and the spectral energy distribution (SED). For the SED, we used the magnitudes in optical from the APASS survey and in near-infrared from the 2MASS and AllWISE surveys (Henden et al. 2015; Munari et al. 2014; Cutri et al. 2014). The list of magnitudes is given in Table 1. The RVs were modelled with Keplerian orbits and the stellar activity with a GP. We did not filtered the activity using de-correlation of the spectroscopic proxies as we did not find any significant correlations between them and the RVs (the highest correlation with RVs is 0.14 for the FWHM). We used a quasi-periodic kernel with the following covariance matrix: (1)

The first term is the quasi-periodic kernel with the hyperparameters: A2 (amplitude), Prot (stellar rotation period), λ1 (coherent timescale) and λ2 (relative weight between the periodic and decay terms). The matrix ΔtRV have elements Δtij = ti − tj from the RV time series. The second term is the identity matrix I multiplied by the data uncertainty vector σ plus an extra source of uncorrelated noise σj (jitter).

The transits were modelled with the jktebop package (Southworth 2008) using an oversampling factor of 30 to account for the long integration time of the data (Kipping 2010). The SED was modelled using the BT-Settl library of stellar atmosphere models (Allard et al. 2012).

A Markov chain Monte Carlo (MCMC) method is implemented in PASTIS. It was used to derive the system parameters and their uncertainties. The spectroscopic parameters were converted into physical stellar parameters using the Dartmouth evolution tracks (Dotter et al. 2008) at each step of the chains. Similarly, the limb darkening coefficients, assuming a quadratic law, were computed using the stellar parameters and tables from Claret & Bloemen (2011).

The complete list of priors of fitted parameters is given in Table A.4. For the stellar temperature, surface gravity and Iron abundance, we used normal priors centred on the values from the spectral analysis. For the period and transit epoch, we used normal priors centred on the values from Christiansen et al. (2018). For the systemic distance to Earth, we used a normal prior centred on the Gaia Data Release 2 (Gaia Collaboration 2018) distance value. For the orbital inclination, we used a sine distribution and for the orbital eccentricity, we used a truncated normal distribution with width σ = 0.083 as reported by Van Eylen et al. (2019) for small multi-transiting exoplanets (R < 6 R⊕). For the other parameters, we used uniform uninformative priors with ranges shown in Table A.4.

The parameter space was explored in several steps and the priors reported in Table A.4 are applicable to the last step of the analysis. First, we normalised the light curve for the joint analysis. For this, we ran 20 MCMCs with 3 × 105 iterations on the photometric data alone, modelling the six transiting planets starting from the values reported in Christiansen et al. (2018), and the stellar modulation with a GP, using the following square exponential kernel: (2)

The first term is the square exponential kernel with the hyperparameters: A1 (amplitude) and l (coherent time scale). The matrix ΔtK2 have elements Δtij = ti − tj from the photometric time series. The second term is similar to the one in Eq. (1). We used a square exponential kernel instead of a quasi-periodic kernel as it requires fewer parameters and the goal at this stage was to fit the light curve and not to derive the stellar rotation period as already done in Sect. 3.2 on the binned light curve. We checked the convergence with a Kolmogorov-Smirnov test, removed the burn-in phase and merged the remaining chains. Then, we subtracted the resulting model of the six planets, removing the transits. Finally, the prediction at observed times on the residuals (activity modulation and residual systematics) of the GP was subtracted from the complete light curve, including the transits. The result is shown in Fig. 2, bottom part. We used the resulting normalised light curve in the final analysis with radial velocities. For this, we ran 96 MCMCs with 6 × 105 iterations. We there also checked the convergence with a Kolmogorov-Smirnov test, removed the burn-in phase and merged the remaining chains. We used the results as starting points for a new iteration of the analysis, to increase the number of chains with the same posterior probability distributions. The results of the final run are shown in Table A.4, including stellar parameters derived in the joint analysis. The RVs and transits in phase at the periods of each planets are shown respectively in Figs. 6 and 7. The radial velocities, together with the six-planet model and the activity GP regression, are shown in Fig. 3. We also ran a full analysis using the PARSEC evolution tracks (Bressan et al. 2012) to check whether the parameters were consistent, which is the case, as it is also reported in Table A.4.

thumbnail Fig. 5

Lomb-Scargle periodogram of HARPS radial velocity. The peak position marked by the purple line corresponds to a period of 25 days (stellar rotation). The orbital periods of the six planets are shown with coloured vertical lines. False alarm probability levels are shown with horizontal lines.
Open with DEXTER

4 Discussion

4.1 Validity of the detections

We were able to derive the masses of planets b, c, d, and e: Mb= 3.1 ± 1.1 M⊕, M⊕, M⊕ and Me = 13.0 ± 2.0 M⊕ with a precision of 34, 20, 18 and 15%, respectively. The radii are compatible with Christiansen et al. (2018). The bulk densities are , 2.8 ± 0.7, 3.2 ± 0.7, and 1.8 ± 0.4 g cm−3, respectively,ranging from Earth to Neptune-like values, as shown in the mass-radius plane in Fig. 8. This brings constraints on their internal structures. Considering their densities and that their mass are ≲ 10 M⊕, these planets likely have a rocky core and a substantial atmospheric layer, composed of volatiles. Such a layer is not taken into account in current models of super-Earth interiors (Brugger et al. 2017). The modelling of the planets is therefore beyond the scope of thispaper and will be the subject of a forthcoming study. The detection of multiple relatively low-density planets around a metallic star is also consistent with previous studies showing an increase in frequency and upper-mass boundary of Neptune-like planets with the metallicity of the host star (Courcol et al. 2016). For planets f and g, we have upper limits at 99% of 8.7 M⊕ and 25.5 M⊕ on the masses, leading to upper limits of 2.1 and 5.1 g cm−3 on the densities.

The activity level in the RVs, as fitted by the GP is m s−1. Unfortunately, this cannot be compared with the photometric activity level since there is a one-year gap between the photometric and spectroscopic observations and the activity level could have changed. We can exclude a stellar origin for the detected signals as the stellar rotation period does not correspond to the periods of the signals at the exception of planet f which we do not detect as it is likely absorbed by the GP regression, being close to half the rotation period (e.g. Damasso et al. 2018). As such, the upper-limit provided for planet f is likely to be underestimated, and the mass may be higher. Assuming a density for planet f in the range [1.79, 4.85] g cm−3, which is the density range of the other planets detected in the system, we can estimate masses between 8.0 M⊕ and 21.6 M⊕, for the derived radius. This range of values is above the upper-limit we obtained for planet f, in line with a likely absorption of the signal by the GP. In this regard, it could have been beneficial to have simultaneous spectropolarimetric observations to optimise the activity filtering as in Hébrard et al. (2016). It may have prevented the non-detection of planet f. Besides, we are not able to search for the presence of planets in the gaps of the broken resonance chain for the same reason.

thumbnail Fig. 6

HARPS radial velocities in phase at each of the six orbital periods (from top to bottom: planet b to g).
Open with DEXTER
thumbnail Fig. 7

Phase-folded transit light curves of planets (from top to bottom) b, c, d, e, f, and g. The best transit model is shown in red. The time scale is not the same between the first three and last three plots. The coverage of the last transit is not as good as in Christiansen et al. (2018) due to the differences in the pipeline used to correct for K2 systematics.
Open with DEXTER
thumbnail Fig. 8

Mass-radius diagram of small planets with masses up to 22 M⊕ and radius up to 5 R⊕. From top to bottom, the lines denotes different compositions for solid planets in between pure water and pure iron (Brugger et al. 2017). We superimposed the known planets (from the NASA Exoplanet Archive, updated on 2019 May 05) in this mass-radius range where the grey-scale depends on the precision on the mass and radius. We mark the positions of the planets K2-138b, c, d, e, f and g. The uncertainties are the 68.3% confidence interval. For planets f and g, which arenot detected, we overplot the 99% confidence interval on the mass in dashed black.
Open with DEXTER

4.2 Transit timing variations

K2-138 is a remarkable dynamical system with five planets in a chain close to the 3:2 resonance. In addition, these planets form a chain of three-body Laplace resonances (Christiansen et al. 2018). This makes K2-138 an ideal target to study transit timing variations. We estimated the transit timing variations using the TTVFaster code (Agol & Deck 2016). We first assumed zero eccentricities for the six planets and took the median masses from Table A.2. We found amplitudes of the order of 2.0, 4.1, 7.3,4.5, 6.4, and 0.02 min for planets from b to g, respectively. Then, assuming median eccentricities from Table A.2, we obtained amplitudes of the order 34, 42, 66, 37, and 41 min for planets b to f (Fig. 9). These amplitudes are excluded by K2 observations which do not show any significant variations at the level of 8–10 min, asdemonstrated by Christiansen et al. (2018). Therefore, the planets are likely close to circular orbits, which is compatible with the posterior distributions we obtained. This is also in line with results showing that tightly packed multi-transiting planets have low eccentricities (Van Eylen et al. 2019). Finally, using the upper limits derived for the masses of planets f and g, we computed the corresponding TTVs of planets e and f. Planet g does not impact TTVs of planet f in a significant (detectable) way as it is far from resonance, being located after a double gap in the chain of resonance. Therefore, it is not possible to use TTVs of planet f to constrain the mass of planet g even if there are undetected and non-transiting planets filling the gap in the chain. Still assuming circular orbits, TTVs of planet e reach an amplitude of the order eight minutes only for masses of planet f above 12 M⊕, which is above the upper limit we derived. However, as we cannot exclude an absorption of the signal of planet f by the GP, it is more conservative to include masses for planet f up to around 12 M⊕. Using the upper limit mass for planet f, the predicted TTVs of planet e are of the order of 6.4 min, and using the median mass they are of 4.5 min. Both are well within reach of space mission like CHEOPS which has a cadence of up to 60 s. Observations of transits of planet e using CHEOPS would be beneficial to constrain further the mass of planet f through a photodynamical analysis (Barros et al. 2015). Additionally, more transit observations, with a higher cadence, of planets b, c and d would allow to constrain the masses of planets c, d and e from TTVs which would make K2-138 an interesting benchmark system for comparing RV and TTV masses. This would allow to better calibrate the two mass measurement techniques. Nevertheless, at this stage, and without more precise photometric observationsof planets e and f, an analysis including TTVs would not allow to measure the masses of planets f and g.

thumbnail Fig. 9

Theoretical transit timing variations predicted using TTVFaster, assuming zero eccentricities for all the planets. Top: planet b (in blue), planet c (in red), planet d (in green). Bottom: planet e (in orange), planet f (in purple), planet g (in brown).
Open with DEXTER

4.3 A favourable system to search for co-orbitals

Multi-planetary systems with a relatively large number of planetary components have been studied in theoretical works as potential environments for the formation of co-orbital planet pairs (e.g. Cresswell & Nelson 2006). In these configurations, two planet-like objects share the same orbital path captured on the gravitational well of each other in different possible 1:1 mean motion resonance (MMR) configurations (e.g. Laughlin & Chambers 2002). Although theoretical results allow these configurations to exist in nature (e.g. Laughlin & Chambers 2002; Ćuk et al. 2012; Leleu et al. 2019a) and different works have searched for them (e.g. Madhusudhan & Winn 2009; Janson 2013; Lillo-Box et al. 2018a,b), none has yet been detected but some candidates have been announced (e.g. Leleu et al. 2019b).

The outcome of studies that analyse the dynamical interactions between planetary embryos during the formation and first stages of the evolution of multi-planetary systems suggests that in a relatively large percentage of the systems pairs of embryos end up being captured in these 1:1 resonances. In particular, Cresswell & Nelson (2008) found that in 30% of their simulations the final system contained one pair of co-orbital planets, while Leleu et al. (2019a) also found this same result for 13% of their simulated systems (including disk dissipation). The systems found by Leleu et al. (2019a) ending up in co-orbital configurations were in a vast majority trapped in 3:2 and 4:3 MMR with a third planet and demonstrated that resonant chains can stabilise co-orbital configurations that would be unstable otherwise.

Given its architecture with multiple planets in MMR, K2-138 represents an excellent system to look for these co-orbital configurations in further studies. A dedicated analysis of the K2 light curve and radial velocity in this regard is, however, out of the scope of this paper and will be studied in future works by following procedures similar to those described in Lillo-Box et al. (2018a) and Janson (2013).

5 Conclusion

In this paper, we report the characterisation of the planets around K2-138. Their parameters were derived via a Bayesian combined analysis of both K2 photometry and HARPS radial-velocities. We computed the following masses for planets c, d and e: M⊕, M⊕ and Me = 13.0 ± 2.0 M⊕ with a precision of 20, 18 and 15%, respectively. For planet b, we have a strong constraint of Mb = 3.1 ± 1.1 M⊕ (precision 34%). The masses and radii derivedlead to bulk densities of , 2.8 ± 0.7, 3.2 ± 0.7, and 1.8 ± 0.4 g cm−3 ranging from Earth to Neptune-like densities for planet b to e. For the two outer planets, we were not able to detect their masses, mostly due to the level of stellar activity. The upper masses derived are Mf < 8.7M⊕ and Mg < 25.5 M⊕ at 99% confidence interval, though the one for planet f should be taken with caution as its signal may have been partly absorbed by the activity filtering. The timing precision obtained with K2 does not allow to use TTVs to bring further constraints on the masses of planets f and g. Additional observations of transits in this system, with a higher cadence than K2, would allow both to constrain the mass of planet f and to use K2-138 as a benchmark system to compare RV and TTV masses. Finally, thanks to its dynamical peculiarities, K2-138 is also a good candidate to search for co-orbital bodies.

Acknowledgements

We thank Jessie L. Christiansen and Kevin Hardegree-Ullman for the useful discussions on TTVs and on the Spitzer transit detection of planet g. We are grateful to the pool of HARPS observers who conducted part of the visitor-mode observations at La Silla Observatory: Julia Seidel, Amaury Triaud, David Martin, Nicola Astudillo, Jorge Martins, Vedad Hodzic. Based on observations made with ESO Telescopes at the La Silla Paranal Observatoryunder programme ID 198.C-0.168. This publication makes use of The Data & Analysis Center for Exoplanets (DACE), which is a facility based at the University of Geneva (CH) dedicated to extrasolar planets data visualisation, exchange and analysis. DACE is a platform of the Swiss National Centre of Competence in Research (NCCR) PlanetS, federating the Swiss expertise in Exoplanet research. The DACE platform is available at https://dace.unige.ch. This research was made possible through the use of the AAVSO Photometric All-Sky Survey (APASS), funded by the Robert Martin Ayers Sciences Fund. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This paper includes data collected by the K2 mission. Funding for the K2 mission is provided by the NASA Science Mission directorate. This research has made use of the Exoplanet Follow-up Observation Program website, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This research has made use of NASA’s Astrophysics Data System Bibliographic Services. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France. The original description of the VizieR service was published in A&AS, 143, 23. D.J.A. acknowledges support from the STFC via an Ernest Rutherford Fellowship (ST/R00384X/1). S.C.C.B. acknowledges support from Fundação para a Ciência e a Tecnologia (FCT) through Investigador FCT contract IF/01312/2014/CP1215/CT0004. O.D.S.D. is supported in the form of work contract (DL 57/2016/CP1364/CT0004) funded by national funds through FCT. This work was supported by FCT – Fundação para a Ciência e a Tecnologia through national funds and by FEDER – Fundo Europeu de Desenvolvimento Regional through COMPETE2020 – Programa Operacional Competitividade e Internacionalização by these grants: UID/FIS/04434/2019; PTDC/FIS-AST/28953/2017 & POCI-01-0145-FEDER-028953 and PTDC/FIS-AST/32113/2017 & POCI-01-0145-FEDER-032113.

Appendix A Supplementary tables and figures

Table A.1

Chemical abundances of host star, relative to Sun, for main elements with number of lines used for each element.

thumbnail Fig. A.1

From left to right and top to bottom: time series and associated Lomb-Scargle periodogram of (a and b) radial velocity (RV); (c and d) bisector span (BIS); (e and f) full width half maximum (FWHM); (g and h) Hα index; (i and j) S index (SMW). The peak position marked by the purple line corresponds to a period of 25 d (stellar rotation). The orbital periods of the six planets are marked in the RV periodogram. For the FWHM, we fitted a linear drift of −32.39 m s−1 yr−1. For the Hα index, the drift is − 7.28 ± 1.34. For the SMW index, the drift is − 57.83 ± 4.85 (reference epoch is RJD_TDB 55 500 for all of them).
Open with DEXTER
Table A.2

System parameters of K2-138 obtained from PASTIS.

Table A.3

Radial velocity data.

Table A.4

Parameters used in the analysis.

References

All Tables

Table 1

Stellar properties of K2-138.

In the text
Table A.1

Chemical abundances of host star, relative to Sun, for main elements with number of lines used for each element.

In the text
Table A.2

System parameters of K2-138 obtained from PASTIS.

In the text
Table A.3

Radial velocity data.

In the text
Table A.4

Parameters used in the analysis.

In the text

All Figures

thumbnail Fig. 1

K2-138 planetary system. All distances are drawn to scale. Planetary radii are enlarged by a factor 50 for readability. Orbits, assumed circular, are shown in black whereas the location of the 3:2 resonances are shown in dashed green.
Open with DEXTER
In the text
thumbnail Fig. 2

K2-138 light curve and associated models. Top: extracted K2 light curve, using the EVEREST pipeline, of K2-138 with in red the Gaussian process regression trained on the out-of-transit light curve. The positions of transits are marked with coloured lines at the bottom of the top figure. Bottom: light curve after subtraction of the Gaussian process fit. The best six-planet photometric model of the transits is superimposed in red.
Open with DEXTER
In the text
thumbnail Fig. 3

HARPS radial velocity time series obtained in two seasons (in red). The six-planets Keplerian model and the Gaussian process regression are shown in black with a 68.3% confidence interval in grey.
Open with DEXTER
In the text
thumbnail Fig. 4

Autocorrelation function of K2-138 light curve. The dashed line in black shows the dominant period, at 24.4 days, corresponding to the stellar rotation. It was computed on the smoothed ACF, following recommendations from McQuillan et al. (2013).
Open with DEXTER
In the text
thumbnail Fig. 5

Lomb-Scargle periodogram of HARPS radial velocity. The peak position marked by the purple line corresponds to a period of 25 days (stellar rotation). The orbital periods of the six planets are shown with coloured vertical lines. False alarm probability levels are shown with horizontal lines.
Open with DEXTER
In the text
thumbnail Fig. 6

HARPS radial velocities in phase at each of the six orbital periods (from top to bottom: planet b to g).
Open with DEXTER
In the text
thumbnail Fig. 7

Phase-folded transit light curves of planets (from top to bottom) b, c, d, e, f, and g. The best transit model is shown in red. The time scale is not the same between the first three and last three plots. The coverage of the last transit is not as good as in Christiansen et al. (2018) due to the differences in the pipeline used to correct for K2 systematics.
Open with DEXTER
In the text
thumbnail Fig. 8

Mass-radius diagram of small planets with masses up to 22 M⊕ and radius up to 5 R⊕. From top to bottom, the lines denotes different compositions for solid planets in between pure water and pure iron (Brugger et al. 2017). We superimposed the known planets (from the NASA Exoplanet Archive, updated on 2019 May 05) in this mass-radius range where the grey-scale depends on the precision on the mass and radius. We mark the positions of the planets K2-138b, c, d, e, f and g. The uncertainties are the 68.3% confidence interval. For planets f and g, which arenot detected, we overplot the 99% confidence interval on the mass in dashed black.
Open with DEXTER
In the text
thumbnail Fig. 9

Theoretical transit timing variations predicted using TTVFaster, assuming zero eccentricities for all the planets. Top: planet b (in blue), planet c (in red), planet d (in green). Bottom: planet e (in orange), planet f (in purple), planet g (in brown).
Open with DEXTER
In the text
thumbnail Fig. A.1

From left to right and top to bottom: time series and associated Lomb-Scargle periodogram of (a and b) radial velocity (RV); (c and d) bisector span (BIS); (e and f) full width half maximum (FWHM); (g and h) Hα index; (i and j) S index (SMW). The peak position marked by the purple line corresponds to a period of 25 d (stellar rotation). The orbital periods of the six planets are marked in the RV periodogram. For the FWHM, we fitted a linear drift of −32.39 m s−1 yr−1. For the Hα index, the drift is − 7.28 ± 1.34. For the SMW index, the drift is − 57.83 ± 4.85 (reference epoch is RJD_TDB 55 500 for all of them).
Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.