Exoplanet characterisation in the longest known resonant chain: the K2-138 system seen by HARPS

The detection of low-mass transiting exoplanets in multiple systems brings new constraints to planetary formation and evolution processes and challenges the current planet formation theories. Nevertheless, only a mere fraction of the small planets detected by Kepler and K2 have precise mass measurements, which are mandatory to constrain their composition. We aim to characterise the planets that orbit the relatively bright star K2-138. This system is dynamically particular as it presents the longest chain known to date of planets close to the 3:2 resonance. We obtained 215 HARPS spectra from which we derived the radial-velocity variations of K2-138. Via a joint Bayesian analysis of both the K2 photometry and HARPS radial-velocities (RVs), we constrained the parameters of the six planets in orbit. The masses of the four inner planets, from b to e, are 3.1, 6.3, 7.9, and 13.0 $\mathrm{M}_{\oplus}$ with a precision of 34%, 20%, 18%, and 15%, respectively. The bulk densities are 4.9, 2.8, 3.2, and 1.8 g cm$^{-3}$, ranging from Earth to Neptune-like values. For planets f and g, we report upper limits. Finally, we predict transit timing variations of the order two to six minutes from the masses derived. Given its peculiar dynamics, K2-138 is an ideal target for transit timing variation (TTV) measurements from space with the upcoming CHaracterizing ExOPlanet Satellite (CHEOPS) to study this highly-packed system and compare TTV and RV masses.


Introduction
Precise knowledge of both the mass and radius of planets is necessary before initiating more advanced studies. This is Full Table A.3 is only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc. u-strasbg.fr/viz-bin/cat/J/A+A/631/A90 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 bimodality 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 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 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 project 1 . 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 meanmotion 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 threebody 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 1 www.exoplanetexplorers.org (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.

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, 2016and 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-implementation 2 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)  flux-position systematics. 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 available 4 EPIC Variability Extraction and Removal for Exoplanet Science Targets (EVEREST) pipeline (Luger et al. 2016(Luger et al. , 2018 on the transitfiltered 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.

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-K2 large programme  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, showing significant 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 S MW . The periodograms associated to the RVs and different indices are given in Fig. A.1.

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 T eff = 5110 +204 −53 K. They also used the extinction A G = 0.0610 +0.2050 −0.0125 and the reddening E(BP − RP) = 0.0308 +0.1006 −0.0056 to determine the luminosity L = 0.517 +0.07 −0.07 L and radius R s = 0.92 +0.02 −0.07 R . 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.

Spectral analysis of the host star
The spectral analysis of the host star was performed to derive its fundamental parameters and the chemical composition of 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. (2012Deleuil et al. ( , 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], T eff = 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 t = 2.3 +0.44 −0.36 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 2.8 +3.8 −1.7 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], T eff = 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.

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 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 S MW time series. We also found a decreasing linear drift in the FWHM, H α and S MW 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 S MW . The priors and posteriors of all the hyperparameters are reported in Table A  Then, we derived the stellar rotation period from the log R HK index and the B−V colour as described in Mamajek & Hillenbrand (2008). We computed the log R HK from the S MW 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 log R HK = −4.76 ± 0.04. We deduced from it the Rossby number R 0 = 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 P rot = R 0 × τ c , we derived P rot = 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 P rot = 21.3 ± 3.2 d, assuming i = π 2 . 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 P rot = 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.

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 nearinfrared 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: A 2 (amplitude), P rot (stellar rotation period), λ 1 (coherent timescale) and λ 2 (relative weight between the periodic and decay terms). The matrix ∆t RV have elements ∆t i j = t i − t j 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 × 10 5 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: The first term is the square exponential kernel with the hyperparameters: A 1 (amplitude) and l (coherent time scale). The matrix ∆t K2 have elements ∆t i j = t i − t j 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 × 10 5 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.

Validity of the detections
We were able to derive the masses of planets b, c, d,  radii are compatible with Christiansen et al. (2018). The bulk densities are 4.9 +2.0 −1.8 , 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 this paper 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 uppermass 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 5.6 +2.9 −1.5 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.

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, as demonstrated 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 multitransiting 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 A90, page 8 of 17 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 observations of planets e and f, an analysis including TTVs would not allow to measure the masses of planets f and g.

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;Cuk 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).

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 c = 6.3 +1.1 −1.2 M ⊕ , M d = 7.9 +1.4 −1.3 M ⊕ and M e = 13.0 ± 2.0 M ⊕ with a precision of 20, 18 and 15%, respectively. For planet b, we have a strong constraint of M b = 3.1 ± 1.1 M ⊕ (precision 34%). The masses and radii derived lead to bulk densities of 4.9 +2.0 −1.8 , 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 M f < 8.7 M ⊕ and M g < 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 Observatory under 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).  (j)   −0.060 orbital eccentricity, planet mass and planet bulk density, the 99% credible intervals are given on the second lines. We assumed R = 695 508 km, M = 1.98842 × 10 30 kg, R ⊕ = 6 378 137 m, M ⊕ = 5.9736 × 10 24 kg and 1 AU = 149 597 870.7 km. The temperatures were derived assuming a zero albedo. The day-side temperature was computed assuming tidally synchronised rotation.