Issue 
A&A
Volume 653, September 2021



Article Number  A41  
Number of page(s)  37  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202140728  
Published online  06 September 2021 
Warm terrestrial planet with half the mass of Venus transiting a nearby star^{★,}^{★★}
^{1}
Instituto de Astrofísica e Ciências do Espaço, CAUP, Universidade do Porto,
Rua das Estrelas,
4150762,
Porto, Portugal
email: olivier.demangeon@astro.up.pt
^{2}
Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto,
Rua Campo Alegre,
4169007,
Porto, Portugal
^{3}
Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências da Universidade de Lisboa,
Campo Grande,
1749016 Lisboa, Portugal
^{4}
Departamento de Física da Faculdade de Ciências da Universidade de Lisboa,
Edifício C8,
1749016 Lisboa, Portugal
^{5}
Observatoire de Genève, Université de Genève,
Chemin Pegasi, 51,
1290
Sauverny, Switzerland
^{6}
Physics Institute, University of Bern,
Sidlerstrasse 5,
3012
Bern, Switzerland
^{7}
Instituto de Astrofísica de Canarias (IAC),
Calle Vía Láctea s/n,
38205
La Laguna,
Tenerife, Spain
^{8}
Departamento de Astrofísica, Universidad de La Laguna (ULL),
38206
La Laguna,
Tenerife, Spain
^{9}
Consejo Superior de Investigaciones Cientícas, Spain
^{10}
Centro de Astrobiología (CSICINTA),
Crta. Ajalvir km 4,
28850
Torrejón de Ardoz,
Madrid, Spain
^{11}
INAF – Osservatorio Astronomico di Trieste,
via G. B. Tiepolo 11,
34143
Trieste, Italy
^{12}
INAF – Osservatorio Astrofisico di Torino,
via Osservatorio 20,
10025
Pino Torinese, Italy
^{13}
Fundación G. Galilei – INAF (Telescopio Nazionale Galileo),
Rambla J. A. Fernández Pérez 7,
38712
Breña Baja,
La Palma, Spain
^{14}
INAF – Osservatorio Astronomico di Brera,
Via E. Bianchi 46,
23807
Merate, Italy
^{15}
INAF – Osservatorio Astronomico di Palermo,
Piazza del Parlamento 1,
90134
Palermo, Italy
^{16}
Institute for Fundamental Physics of the Universe,
Via Beirut 2,
34151
Miramare,
Trieste, Italy
^{17}
European Southern Observatory,
Alonso de Córdova 3107,
Vitacura,
Región Metropolitana, Chile
^{18}
European Southern Observatory,
KarlSchwarzschildStrasse 2,
85748,
Garching b. München, Germany
^{19}
Department of Physics, and Institute for Research on Exoplanets, Université de Montréal,
Montréal
H3T 1J4, Canada
Received:
4
March
2021
Accepted:
12
July
2021
In recent years, the advent of a new generation of radial velocity instruments has allowed us to detect planets with increasingly lower mass and to break the one Earthmass barrier. Here we report a new milestone in this context by announcing the detection of the lowestmass planet measured so far using radial velocities: L 9859 b, a rocky planet with half the mass of Venus. It is part of a system composed of three known transiting terrestrial planets (planets b–d). We announce the discovery of a fourth nontransiting planet with a minimum mass of 3.06_{−0.37}^{+0.33} M_{⊕} and an orbital period of 12.796_{−0.019}^{+0.020} days and report indications for the presence of a fifth nontransiting terrestrial planet. With a minimum mass of 2.46_{−0.82}^{+0.66} M_{⊕} and an orbital period 23.15_{−0.17}^{+0.60} days, this planet, if confirmed, would sit in the middle of the habitable zone of the L 9859 system. L 9859 is a bright M dwarf located 10.6ṗc away. Positioned at the border of the continuous viewing zone of the James Webb Space Telescope, this system is destined to become a corner stone for comparative exoplanetology of terrestrial planets. The three transiting planets have transmission spectrum metrics ranging from 49 to 255, which undoubtedly makes them prime targets for an atmospheric characterization with the James Webb Space Telescope, the Hubble Space Telescope, Ariel, or groundbased facilities such as NIRPS or ESPRESSO. With an equilibrium temperature ranging from 416 to 627 K, they offer a unique opportunity to study the diversity of warm terrestrial planets without the unknowns associated with different host stars. L 9859 b and c have densities of 3.6_{−1.5}^{+1.4} and 4.57_{−0.85}^{+0.77} g cm^{−3}, respectively, and have very similar bulk compositions with a small iron core that represents only 12 to 14% of the total mass, and a small amount of water. However, with a density of 2.95_{−0.51}^{+0.79} g cm^{−3} and despite a similar core mass fraction, up to 30% of the mass of L 9859 d might be water.
Key words: techniques: radial velocities / techniques: photometric / planets and satellites: detection / planets and satellites: terrestrial planets / planets and satellites: composition / stars: individual: L 9859
Full Table B.1 is only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/cat/J/A+A/653/A41
© ESO 2021
1 Introduction
In the pastyears, radial velocity (RV) instruments such as HARPS^{1} (Mayor et al. 2003), HARPSN^{1} (Cosentino et al. 2012), and more recently, CARMENES^{1} (Quirrenbach et al. 2014) and ESPRESSO^{1} (Pepe et al. 2021), have demonstrated that it is now possible to detect planets with masses similar to the mass of the Earth using RV s (e.g., AstudilloDefru et al. 2017a; Rice et al. 2019; Zechmeister et al. 2019; Suárez Mascareño et al. 2020). These results represent an important achievement in the quest for life outside the Solar System. However, it is important to keep pushing toward lower masses and longer periods to ensure that we remain capable to measure the mass of a transiting Earth analog in the habitable zone of a bright host star.
The detection of biosignatures on an exoplanet depends on our capability of studying its atmosphere; this currently relies on transit spectroscopy (e.g., Kaltenegger 2017). Spacebased transit surveys such as Kepler/K2 (Borucki et al. 2010; Howell et al. 2014) and TESS^{1} (Ricker et al. 2015) and even groundbased surveys such as TRAPPIST^{1} (Gillon et al. 2011) have revealed hundreds of transiting terrestrial planets (e.g., Batalha et al. 2013). However, the community still has to detect and study the atmosphere of one of them (Kreidberg et al. 2019). A large fraction of the known terrestrial planets are part of multiplanetary systems (Lissauer et al. 2011). Multiplanetary systems are laboratories for a variety of studies: Planetplanet interactions (e.g., Barros et al. 2015), planetary formation and migration (e.g., Rein 2012; Albrecht et al. 2013; Delisle 2017), and/or comparative planetology (e.g., Mandt et al. 2015; Millholland et al. 2017). The discovery and accurate characterization of a system with multiple transiting terrestrial planets amenable to transit spectroscopy would thus represent a crucial milestone.
The L 9859 system, also known as the TESS Object of Interest 175 (TOI175) system, is a multiplanetary system announced by Kostov et al. (2019, hereafter K19) as composed of three transiting exoplanets with radii ranging from 0.8 to 1.6 Earth radii (R_{⊕}). The host star is a bright (magK = 7.1, Cutri et al. 2003, magV = 11.7, Zacharias et al. 2012) nearby (10.6194 pc, Gaia Collaboration 2018; BailerJones et al. 2018) Mdwarf star (Gaidos et al. 2014). One interesting particularity of this system is its location, with a right ascension (RA) of 08:18:07.62 and declination (DEC) of −68:18:46.80, at the border of the continuous viewing zone (~200 days per year) of the James Webb Space Telescope (JWST, Gardner et al. 2006). This system is thus a prime target for a comparative study of rocky planet atmospheres within the same system (Greene et al. 2016; Morley et al. 2017).
The HARPS spectrograph (Mayor et al. 2003) was used to carry out an RV campaign to measure the masses of these three planets (Cloutier et al. 2019, hereafter C19). The masses of the two outer planets were constrained to 2.36 ± 0.36 and 2.24 ± 0.53 Earth masses (M_{⊕}), leading to bulk densities of 5.3 ± 1.2 and 3.2 ± 1.2 g cm^{−3} for planet c and d, respectively. C19 were unable to constrain the mass of the inner planet b and delivered an upper limit of 1.01 M_{⊕} (with a 95% confidence level). The PFS spectrograph (Crane et al. 2006, 2008, 2010) was also used to attempt measuring the mass of the three planets. With only 14 PFS measurements, Teske et al. (2020) derived masses of 1.32 ± 0.73, 1.24 ± 0.95, and 2.11 ± 0.72 M_{⊕} for planets b, c, and d, respectively. These mass estimates are less precise, but roughly compatible with those of C19. Due to the low number and the lower precision of the PFS data, we did not include these measurements in our analysis.
We report here the results of a followup RV campaign with ESPRESSO (Pepe et al. 2021) aimed at refining the mass of the planets in the L 9859 system. In Sect. 2, we present the RV and photometric data sets. In Sect. 3, we characterize the host star. We describe our analysis of the data sets in Sect. 4. Finally, in Sects. 5 and 6 we discuss the particularities and the importance of this system.
2 Datasets
Our analysis of the L 9859 system relies on RV and photometric time series. The RV s were obtained with the HARPS (Sect. 2.1.1) and ESPRESSO (Sect. 2.1.2) instruments. The light curve (LC) was acquired by TESS (Sect. 2.2) space telescope.
2.1 Highresolution spectroscopy
2.1.1 HARPS
C19 obtained 165 spectra with HARPS, which is installed at the 3.6 telescope of the ESO La Silla Observatory (programs 198.C0838, 1102.C0339, and 0102.C0525) between October 17, 2018 (barycentric Julian date, BJD = 2458408.5), and April 28, 2019 (BJD = 2458601.5). HARPS is a fiberfed crossdispersed echelle spectrograph operating in a temperature and pressureregulated vacuum chamber. It covers wavelengths from 380 to 690 nm with an average spectral resolution of R = 115 000. We obtained the RV s from C19 and refer for details of the observations and their processing to this publication. However, we caution that in order to reproduce the results presented by C19, in particular the RV time series and its generalized LombScargle periodogram (GLSP, Zechmeister & Kürster 2009), we had to exclude four measurements obtained at 2 458 503.795048, 2 458 509.552019, 2 458 511.568314, and 2 458 512.581045 BJD. We identified these measurements with a 4 σ iterative sigma clipping. These measurements were excluded from all the analyses in this paper. All measurements were obtained with fiber B pointed at the sky (no simultaneous observation of a calibration source). One hundred and forty measurements were obtained with an exposure time of 900 s, resulting in an average signaltonoise ratio (S/N) of 41 per resolution element at 650 nm. For the remaining 21 measurements, the exposure time varied from 500 to 1800 s, resulting in a median S/N of 49. The RV s were extracted from the spectra through template matching (AstudilloDefru et al. 2017b). Their median precision (1 σ uncertainty) is 2.08 m s^{−1}.
In addition to the RV measurements, C19 provided the measurement of several stellar activity indicators: the full width at half maximum (FWHM) of the crosscorrelation function (CCF), the bisector span of the CCF (BIS), the depth of the H_{α}, H_{β}, H_{γ} lines, the depth of the sodium doublet NaD, and the Sindex based on the depth of the Ca II H & K doublet. All these indicators are sensitive to chromospheric or photospheric activity.
2.1.2 ESPRESSO
We obtained 66 spectra with ESPRESSO, which is installed at the VLT telescopes of the ESO Paranal Observatory between November 14, 2018 (BJD = 2458436.5), and March 4, 2020 (BJD = 2458912.5), as part of the ESPRESSO Guaranteed Time Observation (programs 1102.C0744, 1102.C0958, and 1104.C0350). ESPRESSO (Pepe et al. 2021) is also a fibrefed highresolution echelle spectrograph operating in a temperature and pressureregulated vacuum chamber. It covers wavelengths from 380 to 788 nm with an average spectral resolution of R = 140 000 in its single UT highresolution mode (HR21, slowreadout mode) that was used for these observations. All measurements were obtained with the sky on fiber B. All measurements were obtained with a 900 s exposure time, resulting in an average S/N of 70 per resolution element at 650 nm. The RV s were extracted from the spectra using version 2.2.1 of the ESPRESSO pipeline DataReduction Software (DRS)^{2}. It computes the CCF of the skysubtracted spectra with a stellar line mask to estimate the RV (Baranne et al. 1996). In this case, the mask was optimized for stars of spectral type M2 V. The CCF was then fit with an inverted Gaussian model. The parameters of the profile are the continuum level; the center of the Gaussian profile, which provides the measurement of the RV ; and its FWHM. Finally, the amplitude provides a measure of the contrast of the CCF. The uncertainties on the measured RV s are computed using the algorithms described in Bouchy et al. (2001) and reflect the photonoiselimited precision. The uncertainties on the FWHM are estimated as the double of RV uncertainties. In addition to the RV, FWHM, and contrast measurements, we computed several activity indicators: the BIS (Queloz et al. 2001), the depth of the H_{α} line, the sodium doublet (NaD, Díaz et al. 2007), and the Sindex (Lovis et al. 2011; Noyes et al. 1984).
From the 66 measurements, we discarded three measurements, obtained at 2 458 645.496, 2 458 924.639, and 2 458 924.645 BJD_{TDB}, due to their high RV uncertainties (identified through an iterative 4 σ clipping). An inspection of the night reports indicates that these measurements were obtained under poor observing conditions: Strong wind, poor seeing, and cirri and bright moon for the first measurement. The last measurement was even interrupted by high winds. The median precision (1 σ uncertainty) obtained on the ESPRESSO RV s is 0.8 m s^{−1} (a factor 2.6better than the HARPS RV s). At about the middle of our RV campaign, in June 2019, the fiberlink of ESPRESSO was replaced. This resulted in an increased throughput, but required us to consider an RV offset between the data taken beforeand after this intervention (Pepe et al. 2021).
2.2 Highprecision photometry with TESS
L 9859 (TIC 307 210 830, TOI175) was observed by TESS in short cadence (2 min) in 9 sectors (2, 5, 8, 9, 10, 11, 12, 28, and 29) with cameras 4 and 3. These observations correspond to ~ 243 days of noncontinous observations taken between August 22, 2018 (BJD = 2458352.5), and September 22, 2020 (BJD = 2459114.5). We downloaded the LC s from the Mikulski Archive for Space Telescopes (MAST) using the python package astroquery. The LC data products provided by the TESS pipeline (Jenkins et al. 2016) provide two LC s, the simple aperture photometry SAP LC and the presearch dataconditioned simple aperture photometry PDCSAP LC (Smith et al. 2012; Stumpe et al. 2014). In contrast to the SAP LC, the PDCSAP LC is detrended using common basis vectors computed over all stars observed on the same CCD. For our analyses, we exclusively used the PDCSAP LC. From the LC, we removed the data points whose quality flags where showing the bits 1, 2, 3, 4, 5, 6, 8, 10, and 12 following the example provided by the TESS team. Following a procedure inspired by K19, we detrended the LC from the residual stellar activity signal and instrument noise using a Gaussian process (GP). We masked all the transits of the three planets using the ephemerides and transit durations provided by K19 and fit the resulting LC with a GP model using the celerite Python package (ForemanMackey et al. 2017; ForemanMackey 2018) and a mean shift between sectors. The functional form of the kernel we used was the one of a damped harmonic oscillator chosen for its flexibility and smooth variations, allowing us to model the unknown mixture of stellar activity and residual instrumental noise. Its equation is (1)
where Q, the quality factor, was fixed to , S_{0} is the amplitude, and ω_{0} is the angular frequency corresponding to the break point in the power spectral density of this kernel.
The fit was performed using an affineinvariant ensemble sampler for MCMC (Goodman & Weare 2010) implemented in the Python package emcee (ForemanMackey et al. 2013), which samples the posterior probability density function. We used a multidimensional Gaussian distribution for the likelihood. For the prior s, we used a uniform prior between −20 and 15 on ln S_{0}, and we obtained a posterior providing an estimate of ppm, using the median and the 68% confidence interval. For ω_{0}, we used a uniform prior between −20 and 15 on lnω_{0}, and we obtained an estimate of (in lnday^{−1}). We did not attribute prior s to the offset between sectors, and the retrieved values are compatible with the values provided in Table C.1. We used 32 walkers (for 11 free parameters) and performed a first run of 500 iterations as burnin. The initial positions for this first run were drawn from the prior for S_{0} and ω_{0} and set to 0 for the offset between sectors. After this first run, we reset the emcee sampler and performed a second run of 2000 iterations, which started from the last positions of the previous run. After this second run, we examined the histogram of the acceptance fraction of the chains to identify chains that had significantly lower acceptance fractions than the others. A lower acceptance fraction implies a stronger correlation between consecutive iterations, which will increase the sampling error of the posterior PDF inferred from the histograms of the chains. We also examined the histogram of the logarithm of the posterior probability of the chain (estimated by the average of this value computed over the last 1% of the iterations ofeach chains). The objective was to understand whether all the chains have converged toward regions of the parameter space that have similar posterior probability density values. In this case, both histograms are monomodal, indicating that all chains have similar acceptance fractions and sample regions of the parameter space with similar posterior probability density values. We confirmed that all the chains converged and converged to the same region of the parameter space using the Geweke criterion (Geweke 1992). All the chains indeed converged to the same region of the parameter space after the first 750 iterations of the second run. We further confirmed that the remaining parts of the chains converged and were long enough by computing the integrated autocorrelation time using the method implemented in emcee and verifying that it was ten times shorter than the remaining number of iterations.
We normalized the LC by dividing it by the best GP model, whose parameters values are the median values of the converged MCMC chains. Finally, we cut the LC to keep only data points within 1.5 transit durations on either side of each midtransit time. This reduced the number of data points and the computation time.
3 Characterization of the M dwarf L 9859 A
According to K19, L 9859 A is an M3V star. The derivation of accurate stellar properties through highresolution spectroscopy for M stars is complicated because blended lines prevail. We thus used several approaches to characterize L 9859 A in order to assess and discuss the homogeneity and the accuracy of the outcomes. This analysis is presented in detail in Appendix A, and we summarize the results in this section.
3.1 Stellar atmospheric parameters
To derive the stellar parameters, effective temperature (T_{eff}), surface gravity (log g), and metallicity ([Fe∕H]), we chose to fit the combined spectrum of L 9859 A constructed using 61 ESPRESSO spectra (S∕N = 1063 at 7580 Å) with the latest version of the spectral synthesis code STEPARSYN (Tabernero et al. 2018, 2021, see Appendix A.1.1 for more details). We adopted the estimates provided by STEPARSYN except for the uncertainty on T_{eff}, which we identified as underestimated (see Appendix A.1). We enlarged this uncertainty to encompass the best values provided by the other methods within 1 σ. The set of adopted estimates is provided in Table H.1.
3.2 Stellar modeling: Mass, radius, and age
Thanks to the high precision and accuracy of Gaia parallactic distances (10.6194 ± 0.0032 pc inferred from the GaiaDR2 parallax by BailerJones et al. 2018) and the wellsampled photometric spectral energy distribution (SED, see Appendix A.1.3), we can derive a reliable estimate of the absolute bolometric luminosity of L 9859: 0.01128 ± 0.00042 L_{⊙}. Added to our estimate of T_{eff} (see Sect. 3.1, Appendix A.1 and Table H.1), we infer the radius of L 9859 A to be using the StefanBoltzmann law. This agrees well (better than 1 σ) with the literature value derived by K19 from the massradius relations for M and K dwarfs of Boyajian et al. (2012). We derived the mass of L 9859 A using the Virtual Observatory SED Analyzer online tools^{3} (VOSA, Bayo et al. 2008, see Appendix A.2 for more details). VOSA derives the mass by comparing the measured T_{eff} and bolometric luminosity to BTSettl evolutionary tracks (Allard et al. 2012). Finally, we determined the age of L 9859 A using the photometry and distance provided by Gaia (see Appendix A.2 for more details). We compared the location of L 9859 A in the colormagnitude diagram (see Fig. 1) to mean sequences of stellar members of the β Pictoris moving group (~20 Myr, MiretRoig et al. 2020), the TucanaHorologium moving group (~ 45 Myr, Bell et al. 2015), the Pleiades open cluster (~120 Myr, Gossage et al. 2018), and the field (possible ages in the range 0.8–10 Gyr). This comparison allowed us to infer that L 9859 has an age consistent with that of the field. This age estimate is confirmed by our kinematics analysis, which indicates that L 9859 A is a thindisk star that does not belong to any known young moving group (see Appendix A.2 and Table A.4). The adopted radius, mass, and ages of L 9859 A are provided in Table H.1.
Fig. 1
Absolute magnitude (in the G Gaia bandpass) vs. color (magnitude difference between the Gaia bands G_{BP} and G_{RP}): L 9859 A (TOI175) is located in the Gaia colormagnitude diagram together with the mean sequences of young clusters and moving groups (Luhman 2018) and the main sequence of stars (Cifuentes et al. 2020). The error bars of L 9859 A are smaller than the symbol size. The gray area represents the 1 σ dispersion of field M dwarfs. 
3.3 Stellar Mg and Si abundances
Stellar abundances of Mg and Si are valuable constraints for modeling the interior of planets (see Sect. 5.3). However, deriving individual abundances of M dwarfs from visible spectra is a very difficult task (e.g., Maldonado et al. 2020). We estimated the abundances of Mg and Si following the procedure described in Adibekyan et al. (2017). From the APOGEE DR16 (Jönsson et al. 2020), we selected cool stars (T_{eff} < 5500 K, the choice of this temperature limit does not have a significant impact) with metallicities similar to that of L 9859 A within 0.05 dex. We considered only stars with the highest S∕N (> 500) spectra to guarantee the high quality of the extracted parameters and abundances of these stars. Because L 9859 A is a member of the Galactic thindisk population (see Table A.4), only stars belonging to the thindisk population were selected. The selection of the thindisk stars was based on the [Mg/Fe] abundance of the APOGEE stars (see, e.g., Adibekyan et al. 2012). With these constraints, we obtained a sample of about 1000 thindisk stars with properties similar to those of our target. The mean abundances of Mg and Si of these stellar analogs were adopted as proxy for the empirical abundances, and their standard deviation (startostar scatter) was adopted as the uncertainty (see Table H.1).
3.4 Stellar rotation and activity periods
As mentioned in Sects. 2.1.1 and 2.1.2, the HARPS and ESPRESSO instruments give access to the time series of several activity indicators. These activity indicators are sensitive to variations in the stellar chromosphere, but not to the presence of planets in the system. They are therefore ideal for identifying periodicities that arise from stellar chromospheric activity. To identify these periods, we computed the GLSP of all available activity indicators, see Fig. 2. This figure also includes the GLSP of the RV measurements.
The GLSP s of the ESPRESSO activity indicators suggest that the rotation period (P_{rot}) of L 9859 A is days, measured on the highest peak of the FWHM GLSP, in agreement with C19. The GLSP s of the FWHM, the contrast of the CCF, and the Sindex all show peaks at this period with a falsealarm probability (FAP) below 0.1%. The FAP levels were computed using the analytical relation described in Zechmeister & Kürster (2009) for the Zechmeister–Krster (ZK) normalization. Our GLSP sof the HARPS activity indicators are consistent with those presented by C19. The GLSP s of the BIS, the Sindex, and H_{α} show peaks with an FAP below 0.1%, but not at the same period. However, as noted by C19, the peak with the highest significance, which is found in the GLSP of H_{α}, is close to 80 days. This period is used by C19 as an estimate of the rotation period.
Photometric time series can also provide insight into the stellar rotation periods. The appearance and disappearance of dark and bright active regions due to stellar rotation, such as spots and plages, produce a modulation in the LC. To investigate the rotational modulation in the TESS LC, we first fit the PDCSAP TESS LC with a GP and an offset for each sector. Using the retrieved offsets between the sectors, we computed the GLSP of the TESS LC presented in Fig. 3 (see also Appendix C). The three highest peaks in this periodogram in order of decreasing amplitude are at 93, 115, and 79 days. The 79day periodicity is a confirmation of the 80day period identified in the GLSP s of the spectroscopic time series presented in Fig. 2. However, the 93 and 115day periodicities are absent in these periodograms.
Overall, the spectroscopic and photometric time series all exhibit power at a period of 80 days. This is thus our best guess for the rotation period of L 9859. However, the power spectrum of all these stellar activity indicators depicts a complex activity pattern thatdoes not seem to be fully described by only one periodicity and its harmonics.
Fig. 2
GLSP of the RV and activity indicators from ESPRESSO (a) and HARPS (b) data. The last row for both instruments presents the window function. The vertical dotted lines indicate from right to left the orbital period of the planets b, c, d, and e, the planetary candidate 05, and half and the full stellar rotation period (assumed here to be 80 days). The horizontal lines indicate the amplitude levels corresponding to 10% (dashed line), 1% (dotdashed line) and 0.1% (dotted line) of the FAP. The amplitudes of the GLSP s are expressed using the ZK normalization described in Zechmeister & Kürster (2009, Eq. (5)). The FAP levels are computed using the analytical relation also described in Zechmeister & Kürster (2009) for this normalization. We display the GLSP of the BIS for completeness and comparison with C19, but we caution that the reliability of BIS measurements from CCF s for M dwarfs is uncertain (Rainer et al. 2020). 
Fig. 3
GLSP of the TESS LC. The format of the this figure is identical to that in Fig. 2. In particular, the power of the GLSP is normalized using the ZK normalization. The highest peak in this periodogram is for a period of 93 days. 
4 Radial velocity and lightcurve modeling
4.1 Search for additional planets in the L 9859 system
K19 and C19 confirmed the presence of three transiting planets in the L 9859 system. Using the new sectors from TESS and the new RV data from ESPRESSO, we wish to improve the precision of the planetary parameters and search for additional planets.
The GLSP of the HARPS RV data (see Figs. 2b) shows six peaks above an FAP of 10% at about 3.7 (orbital period of planet c), 7.6 (orbital period of planet d), 13, 15, 23, and 40 (~ P_{rot}∕2) days. The GLSP of the ESPRESSO RV data (see Figs. 2a) shows two narrow peaks above an FAP of 10% at about 13 and 23 days. The fact that the two peaks identified in the ESPRESSO data are also present in the HARPS data and are not an obvious fraction of the stellar rotation period indicates that there might be two additional planets in the system.
Due to the high computational cost linked to the analysis of the nine TESS sectors, we divided our analysis into three steps. In the first step (Sect. 4.1.1), we analyze the TESS LC alone in order to refine the properties of the three known transiting planets and in particular their ephemerides. In the second step (Sect. 4.1.2), we use these ephemerides as prior for the analysis of the highresolution spectroscopy data. The main objective of this second step is to assess the presence of additional planets in system L 9859 (Sect. 4.1.3). Finally, in a third step (Sect. 4.2), we perform a final joint analysis of the RV s and the LC to obtain the final parameters of the system.
4.1.1 LC analysis
To model the planetary transits, we used a modified version^{4} of the Python package batman^{5} (Kreidberg 2015). The parameters used for each planets are the orbital period P, the time of inferior conjunction (t_{ic}), the products of the planetary eccentricity by the cosine and sine of the stellar argument of periastron (e cos ω and e sin ω), the ratio of the planetary radius to that of the star (R_{p}∕R_{*}), and the cosine of the planetary orbital inclination (cos i_{p}). The model also included the stellar density (ρ_{*}). For the limbdarkening law, we used the four coefficients of the nonlinear model (u_{1,TESS}, u_{2,TESS}, u_{3,TESS}, and u_{4,TESS}). To this set of parameters, we added one additive jitter term (σ_{TESS}) for the photometry in all TESS sectors to account for a possible underestimation of the error bars (Baluev 2009).
To infer the values of these parameters, we maximized the posterior probability density function (PDF) of the model as prescribed by the Bayesian inference framework (e.g., Gregory 2005). The likelihood functions we used were multidimensional Gaussians. To obtain robust error bars, we explored the parameter space with an affineinvariant ensemble sampler for MCMC implemented in the Python package emcee (ForemanMackey et al. 2013). We adapted the number of walkers to the number of free parameters in our model. As a compromise between speed and efficiency, we used ⌈n_{free} × 2.5 × 2⌉∕2 walkers, where n_{free} is the number of free parameters and ⌈ ⌉ is the ceiling function. This allowed us to have an even number of walkers that was at least twice (~ 2.5 times) the number of free parameters, as suggested by the authors of emcee. The initial values of each walker were obtained from the output of a maximization of the posterior PDF made with the Nelder–Mead simplex algorithm (Nelder & Mead 1965) implemented in the Python package scipy.optimize. The initial values for the Nelder–Mead simplex maximization were drawn from the prior s of the parameters. The objective of this premaximization was to start the emcee exploration closer to the best region of the parameter space and thus reduce its convergence period. Our experience is that this usually results in a reduction of the overall computational time because the Nelder–Mead simplex algorithm usually converges faster than emcee.
The prior PDF assumed for the parameters were noninformative and are given in Table H.1 (column prior), along with references justifying their use when needed (column Source prior). Along with the posterior PDF provided in the same table, this allowed a qualitative assessment of the impact of the prior on the posterior (inferred values). A detailed description of the reasons for the choice of each prior is given in Appendix D.
To choose the initial values for the analysis, those used to start the preminimization, we usually use values drawn from the prior s. However, here, we did not analyze the full TESS LC, only small portions of it around the location of the transits (see Sect. 2.2). Consequently, drawing initial values from noninformative prior s would very likely cause the simulated transits to fall outside of the selected portions of the LC and make the optimization impossible. To prevent this, we drew the initial values for P, t_{ic}, R_{p} ∕R_{*}, and cos i_{p} from the posterior PDF s obtained by K19.
We used 50 000 MCMC iterations and analyzed the chains using the same procedure as is described in Sect. 2.2. The posterior distributions of the parameters of the three transiting planets were then used as prior s for the analysis of the RV s.
4.1.2 RV analysis
Our model of the RV s is composed of three main components: The planetary model, the stellar activity model, and the instrumental model. In their analysis of the HARPS data, C19 demonstrated the importance of stellar activity mitigation for this system. They inferred an amplitude of ~ 7 m s^{−1} for the stellar activity signal compared to ≲2 m s^{−1} for the semiamplitude of the three planetary Keplerians. We thus paid particular care to the stellar activity mitigation and used two different approaches. The first approach is similar to the one used by C19. We fitted the RV data using Keplerians for the planetary signals and a GP with a quasiperiodic kernel for the stellar activity. The mathematical expression of the kernel of this GP is (2)
where A_{rv} is the amplitude of the covariance, τ_{decay} is the decay timescale, P_{rot} is the period of recurrence of the covariance, and γ is the periodic coherence scale (e.g., Grunblatt et al. 2015). We used the Python package george^{5} (Ambikasaran et al. 2015) for the implementation. For the interpretation of the results, it is valuable to understand the impact of these hyperparameters on the stellar activity model that this kernel produces (e.g., Angus et al. 2018; Haywood et al. 2014). A_{rv} scales with the amplitude of the stellar activity signal. P_{rot} indicates its main periodicity and is considered as a measure of the stellar rotation period (Angus et al. 2018). τ_{decay} and γ are two indicators of the coherence of the stellar activity signals. τ_{decay} governs the aperiodic coherence, the coherence between one period and the next periods. It is considered a measure of the timescale of growth and decay of the active regions (Haywood et al. 2014). If it is longer than P_{rot}, the stellar activity pattern will change slowly from one rotation period to the next. γ controls the periodic coherence, that is, the coherence of the signal within a stellar rotation period. It is considered an indicator of the number of active regions. Thelarger γ, the weaker the correlation between two points within a rotation period. γ governs the complexity of the harmonic content of the stellar activity signal (Angus et al. 2018).
For the second approach, we used the same model, but we jointly fit the RV s and the FWHM values that accompany each RV measurement. The FWHM was fit with a GP with a quasiperiodic kernel. This kernel is independentof the one used for the RV, but it uses the same hyperparameters, except for the amplitude (A_{FWHM}). This approach, inspired by Suárez Mascareño et al. (2020) and subsequently LilloBox et al. (2020), relies on the assumption that the variations in FWHM are solely due to stellar activity and that their periodicity and coherence are the same as the stellar activity component of the RV. Under these assumptions, the joint fit of the RV and FWHM data sets allows constraining the hyperparameters of the quasiperiodic kernel better. In contrast to a first fit of the FWHM s followed by a second fit of the RV s using the marginalized posterior of the first fit as prior for the second, this approach preserves the correlation between the hyperparameters.
For the planetary model, we used a constant systemic velocity (v_{0}) and one Keplerian function per planet in the system. The parameters of each Keplerian are the semiamplitude (K) of the RV signal, and similarly to Sect. 4.1.1, the orbital parameters P, t_{ic}, e cosω, and e sin ω. The Keplerians were implemented using the Python packages radvel^{5} (Fulton et al. 2018).
For the instrumental model, as mentioned in Sect. 2.1.2, we considered three instruments in our model because of the fiberlink change of ESPRESSO: HARPS, ESPRESSO before (pre), and ESPRESSO after the intervention (post). We used ESPRESSO_{pre} as RV reference, meaning that v_{0} was measured with the data coming from this instrument. We modeled the RV offsets with the other two instruments with two offset parameters (ΔRV_{HARPS∕pre} and ΔRV_{post∕pre}). The FWHM is also subject to offsets between instruments, and our model includes a constant level for each instrument (C_{pre}, C_{post}, and C_{HARPS}). Finally, for the RV and FWHM and for each instrument, we considered one additive jitter parameter to account for a potential underestimation of the measurement errors due to underestimated or even nonconsidered noise sources (Baluev 2009) (σ_{RV,pre}, σ_{RV,post}, σ_{RV,HARPS}, σ_{FWHM,pre}, σ_{FWHM,post}, and σ_{FWHM,HARPS}).
To infer the values of these parameters, we performed a preminimization followed by an MCMC exploration as described in Sect. 4.1.1. The only difference was that this time, the initial values were all drawn from the prior s. The prior PDF s assumed for the parameters are given in Table H.1, except for the prior of P and t_{ic} of the three transiting planets. For these, we used the posterior PDF s of our analysis of the TESS LC (provided in a footnote of Table H.1). A detailed description of the reasons for the choice of each prior is given in Appendix D.
4.1.3 Evidence for additional planets in the L 9859 system
We analyzed our RV data with six different models by varying the number of planets in the system from three to five and the stellar mitigation approach with or without the FWHM data (see Sect. 4.1.2). After each analysis, we inspected the output of the fit using plots such as those provided in Figs. 4 and 5. Figure 4 shows the RV time series including the data from both instruments, the best planetary plus activity model, and the residuals of this model. Figure 5 displays the GLSP of the combined RV data, the residuals, the planetary and stellar activity models sampled at the same times as the RV data, and the window function (WF).
Extensive outputs are shown and discussed in Appendix F. From the fit of the three planets model (see Figs. F.1 and F.2), the GLSP of the residuals displays a narrow peak at 13 days, which we consider to be a strong insight for the presence of a fourth planet in the L 9859 system at this period. For the analysis with four planets, we adopted a noninformative prior for the orbital period of the potential fourth planet (see Table H.1). However, to speed up convergence, we drew its initial values from a Gaussian distribution with a mean of 13 days and a standard deviation of 1 day. As shown in Fig. 5, the GLSP of the residuals of the fourplanet model shows two narrow peaks around 1.743 and 2.341 days. These two peaks are aliases of one another. Because transit signals in the TESS LC are absent at these periods, we did not explore the possibility of a planet at these periods. However, the peak at 23 days in the GLSP of the RV s appears to be absorbed by the stellar activity model despite the absence of signal at 23 days in the GLSP s of the FWHM and other activity indicators. We thus performed another analysis with five planets. We again set a noninformative prior for the orbital period of the potential fifth planet (see Table H.1), but we drew its initial values from a Gaussian distribution with a mean of 23 days and standard deviation of 1 day. The fit converged toward a significant detection of the semiamplitude of a fifth Keplerian signal.
Table 1 regroups the Bayesian information criterion (BIC) values computed for all the models we tested. However, the BIC is not necessarily adapted for our analysis because our models are nonlinear and our prior s uninformative but relatively complex (see Appendix D). Consequently, we also computed the Bayesian evidence () of our models using the Perrakis algorithm (Perrakis et al. 2014) using the Python implementation bayev^{6} (Díaz et al. 2014). We computed the logarithm of based on 5000 sets of parameters values and repeated the process 150 times. From these 150 computations, we extracted the median and the 68% confidence interval (using the 16th and 84th percentiles) and report these values in Table 1. The Bayesian evidence agrees with the BIC values. According to both criteria, the fourplanet model is favored and obtains the best values (minimum for the BIC and maximum for the Bayesian evidence). The only difference is in the absolute difference between the four and the fiveplanet models. The BIC values of the fiveplanet model is significantly higher (ΔBIC = 3 for the RV+FWHM analysis), while the Bayesian evidence of these two models is very similar ().
We thus conclude that our additional ESPRESSO RV campaign allows us to identify one additional planet in the L 9859 system: a fourth planet, hereafter planet L 9859 e, with an orbital period of 12.80 days. We also identify a planetary candidate, a potential fifth planet, hereafter planet 5, with an orbital period of 23.2 days. We show in Sect. 4.3 that these two additional planets do not transit.
Finally, retrieving the relevant information on L 9859 from the new Gaia Early Data Release 3 (EDR3), we note that an astrometric excess noise of 0.171 mas is reported, and the reduced unit weight error (RUWE) statistics has a value of 1.27. At G = 10.6 mag. The star is not so bright as to be strongly affected by unmodeled systematics due to the limited calibration. The Gaia EDR3 astrometry information (particularly RUWE) can thus be interpreted as providing weak evidence for the possible existence of an unresolved, massive outer companion (e.g., Belokurov et al. 2020; Penoyre et al. 2020). However, no longterm trend is observed in our RV analysis.
Fig. 4
Outcome of the fit of the four planets model: top left: RV time series along with the best model (solid green line), which includes the planetary signals and best prediction from the GP stellar activity model. The 1 σ uncertainties from the GP prediction are also displayed (shaded green area). For this plot, we subtracted the systemic velocity and the instruments offsets from the RV data (see values in Table H.1). Bottom left: Time series of the residuals of the best model. Right: Zoom on a small portion of the time series to better visualize the shorttimescale variations. 
Fig. 5
Outcome of the fit of the fourplanet model: GLSP s of the RV time series (top) and of the planetary (second) and stellar activity (third) models sampled at the same time as the RV data, GLSP of the time series of the residuals (fourth) andthe window function (bottom). The vertical lines on the GLSP s correspond to the orbital periods of planets b, c, d, and e, half and the full rotation period (estimated at 80 days) from right to left. 
Comparison of different models of the RV s of the L 9859 system.
4.2 Joint analysis of RV and photometry data
For the joint analysis of the RV and photometry data, we only fit the best model identified by the RVonly analysis due to the much highercomputational time associated with the data of the nine TESS sectors: The four planets plus stellar activity model on the RV and FWHM data sets.
4.2.1 Keplerian model
The model of the RV, FWHM, and LC data as well as the inference process is similar to the models used in Sects. 4.1.2 and 4.1.1. The prior PDF assumed for the parameters is given in Table H.1 and discussed in Appendix D. The initial values were drawn from the prior PDF s with a few exceptions. For P, t_{ic}, R_{p} ∕R_{*}, and cos i_{p} of the three transiting planets, we used the posterior PDF obtained by K19 to draw the initial values. For P of the two exterior planets, we used Gaussian prior s with a standard deviation of 1 day and a mean value of 13 and 23 days for planet e and planetary candidate 5, respectively.
From our MCMC exploration, we extracted the estimates of the model parameters using the median of the converged iterations as best model values and their 16th and 84th percentiles as the boundaries of the 68% confidence level intervals. We also derived estimates for secondary parameters. As opposed to the model parameters (also called main or jumping parameters) described in the previous sections, secondary parameters are not used in the parameterization chosen for our modeling and are not necessary to perform the MCMC exploration. However, they provide quantities that can be computed from values of the main parameters and are of interest for describing the system. The secondary parameters that we computed are ΔF∕F the transit depth, i the orbital inclination, e the eccentricity, ω the argument of periastron, a the orbital semimajor axis, M_{ref} the mean anomaly at a given reference time (set as BTJD = 1354, the time of the first TESS measurement), b the impact parameter, D14 the outer transit duration (duration between first and fourth contact), D23 the inner transit duration (duration between second and third contact), R_{p} the planetary radius, M_{p} the planetary mass, F_{i} the incident flux at the top of the planetary atmosphere, and T_{eq} the equilibrium temperature of the planet (assuming an albedo of 0). After the full MCMC analysis, we drew for each iteration a mass, a radius, and an effective temperature value for the star using Gaussian distributions whose mean and standard deviation were set according the results of our stellar analysis (see Sect. 3 and Table H.1). We then consistently computed the value of all the secondary parameters at each iteration of the emcee exploration, which provided us with chains for the secondary parameters. Finally, we estimated their best model values and 68% confidence intervals with the same method as for the main parameters.
4.2.2 Dynamical stability and parameters of the L 9859 system
In compact multiplanetary systems such as L 9859, the assumption of longterm stability of the system can bring strong constraints on the planetary masses and orbital properties. Both K19 and C19 performed Nbody dynamical simulations with the objective of constraining the orbital eccentricity of the planets in this system. Both studies provide compatible conclusions: the eccentricity of planets c and d should be 0.1 or lower. As only the three inner planets were known at the time, the discovery of a fourth planet in this system requires revisiting this question. To do this, we used the framework implemented in the spock Python package (Tamayo et al. 2021, 2020, 2016). spock has been developed specifically to assess the stability of compact multiplanetary systems. It performs a short, and thus relatively inexpensive, Nbody simulations (10^{4} orbits of theinner planet) using the Python package rebound (Rein & Liu 2012). This simulation is then used to compute metrics based on established stability indicators (see Tamayo et al. 2020, and references therein). These metrics are then provided to a machinelearning algorithm that estimates the probability that the simulated system is stable on the long term (typically 10^{9} orbits of the inner planet). According to spock, the probability that the system described by the best model parameters inferred from our joint analysis of the RV, FWHM, and photometry data is stable is 0. This means that the simulated system becomes unstable during the short Nbody simulation (within 10^{4} orbits of theinner planet). This stresses the importance of considering the dynamical stability for this system.
Following the procedure described in Tamayo et al. (2021), we used spock to compute the probability of stability of the 10^{5} versions of the L 9859 systems described by the last 10^{5} converged MCMC iterations of our joint analysis. For these computations, we used the WHFast symplectic integrator (Rein & Tamayo 2015) of rebound. We set a maximum distance of 0.4 AU (~ 6 times the semimajor axis of planet e), meaning that all simulations that led to one of the planets traveling 0.4 AU away from the barycenter of the system were stopped and their probability of stability was set to 0. For each MCMC iteration we considered, we provided for the Nbody simulation the mass of the star, the masses of the planets, and their orbital elements orbital period, semimajor axis, inclination, eccentricity, argument of periastron passage, mean anomaly at the beginning of the simulation (set as 1354 BTJD, the time of the first TESS measurement where BTJD = BJD−2 457 000), and the longitude of ascending node. All these quantities, except for the longitude of the ascending node, are either main or secondary parameters of the model (see Sect. 4.2). Their values were thus taken directly from the MCMC chains or from their associated secondary parameters chains. For the longitudes of the ascending node, we drew values from a uniform distribution between 0 and 2 π.
With the probability of longterm stability estimated for the last 10^{5} iterations of our MCMC analysis of the joint fit of the data, we selected the iterations for which the probability of stability is higher than 40% (as in Tamayo et al. 2020). This left us with only 1588 iterations. From these iterations and using their probability of stability as weight, we computed the weighted median and the weighted 16th and 84th percentile that we used as the best model values and the boundaries of the 68% confidence interval, respectively, as suggested by Tamayo et al. (2021). These estimates now describe a system with a high probability of longterm stability, and they are reported in Table H.1. The phasefolded data (RV and photometry) and the best model are displayed in Figs. 6 and 7. The main impact of the longterm dynamical stability condition is on the eccentricity of planet c, which decreases from to . The eccentricities of the other planets remain unchanged or decrease slightly, but they are well within 1 σ of the previous estimates. The other parameters of the system are all compatible with their previous estimates at better than 1 σ. With these updated estimates, the eccentricities of the three transiting planets satisfy the constraints derived by both K19 and C19 from their respective Nbody simulations.
Finally, in order to assess whether planets c and d are in mean motion resonance, we performed an additional Nbody simulation for each iteration of the system with a probability of longterm stability higher than 40%. For each iteration, we started the simulation using the parameter values found in the MCMC chains or the associated secondary parameters chains, as before. We used rebound and the WHFast symplectic integrator with a time step of 10^{−4} year∕2π (which corresponds to ~ 10^{3} time steps per orbits of planet c). We integrated each simulation for the duration of our observations, 560 days between the beginning of the TESS observations and the last ESPRESSO point. For each time step, we calculated the 2:1 resonant angles (θ) of planet c and d, whose equation is (e.g., Quillen & French 2014, Eq. (1))
where λ is the mean longitude. As explained in Delisle (2017), if planets c and d are in mean motion resonance, their resonant angles should librate around a constant value. Following a procedure already used by Hara et al. (2020), we computed the derivative of the resonant angles using the finite difference approximation and averaged their value over the duration of the simulation. The normalized histogram of the 1588 values of the average derivatives of the resonant angles obtained is not compatible with zero and indicates that planets c and d are not in mean motion resonance.
Fig. 6
Phasefolded HARPS and ESPRESSO RV s, best model (top) and residuals (bottom) for the four planets. The HARPS data, presented in Sect. 2.1.1, are displayed with empty blue circles, and the ESPRESSO data, presented in Sect. 2.1.2, are displayed with orange circles. The filled orange circles are for the data taken before the fiber change of ESPRESSO. The empty orange circles are for data taken after the change. For clarity, the error bars of the HARPS and ESPRESSO data points are not displayed. For this plot, the stellar activity model has been subtracted from each data point. The points with error bars in red correspond to averages of the data within evenly spaced bins in orbital phase, whose size is 0.07 orbital period. The best model is shown with a green line. Before the subtraction of the stellar activity model the RMS of the RV data is 3.5, 3.4, and 3.2 m s^{−1} for HARPS, ESPRESSO_{pre}, and ESPRESSO_{post}, respectively. After the subtraction of the stellar activity model, it is 2.9, 2.5, and 2.3 m s^{−1} for HARPS, ESPRESSO_{pre}, and ESPRESSO_{post}, respectively. Finally, after subtraction of the planetary model, the RMS of the residuals is 1.8, 1.2, and 0.7 m s^{−1} for HARPS, ESPRESSO_{pre}, and ESPRESSO_{post}, respectively. 
4.3 Three transiting planets
Our RV analysis (see Sect. 4.1.3) concluded with the existence of a fourth planet and a planetary candidate. They have not been reported before. Assuming that all planets in the system are coplanar, we can infer an orbital inclination of degrees and predict the impact parameter of planet e () and candidate 5 (). From these impact parameter distributions, we estimate a probability of 4.8% and 0.11% that planet e and planetary candidate 5 transit their host star, respectively. Using the nine TESS sectors and the best ephemerides inferred from our analysis, we do not detect any sign of transit from either planet e or planetary candidate 5 (see Fig. E.1 and Appendix E for more details of the analysis we performed).
5 Discussion
5.1 Stellar activity modeling and mitigation
Stellar activity mitigation is a current focus of the exoplanet community due to its impact on the detection and characterization of lowmass planets, both in RV (e.g., Dumusque et al. 2017) and transit photometry (e.g., Barros et al. 2020). For this analysis, we used a GP with a quasiperiodic kernel to account for the important stellar activity imprint on the RV data that were already identified by C19. We analyzed the data with two slightly different approaches (see Sect. 4.1.2): one used a GP on the RV data alone, and the other used the time series of a stellar activity indicator (here the FWHM) that was fit simultaneously with the RV. The motivation for the second approach is to place stronger constraints on the hyperparameters of the GP. In the case of L 9859, we have shown in Sect. 4.1.3 that the two approaches provide similar answers for the preferred model. A comparison of the posterior PDF of all common parameters to the two approaches shows that they also provide compatible estimates (within 1 σ).
Fig. 7
Phasefolded TESS LC, best model (top) and residuals (bottom) for the three transiting planets. The data presented in Sect. 2.2 are displayed in black. For clarity, the error bars are not displayed. The points with error bars in red correspond to averages of the data within evenly spaced bins in orbital phase whose size corresponds to 5 min. The best model is shown with a black line. The standard deviation of the raw and binned residuals is indicated above each residual plot. 
5.2 Fourplanet system hosting the smallest planet measured through RV
The additional six sectors we analyzed compared to K19 improved the characterization of the three transiting planets presented by K19 and C19 (see Table H.1). The ephemerides of the three planets are improved by factors ~ 2 and ~ 10 for the time of transit and the orbital period, respectively. The relative precisions of the radius ratios (R_{p} ∕R_{*}) are also improved by a factor ~2 for the two inner planets and by a factor ~4 for planet d.
We also improved the mass determinations for these three planets. We derived the mass of planet b with 40% relative precision (C19 only provided an upper limit). With an RV semiamplitude of and a mass of M_{⊕} (half the mass of Venus), L 9859 b currently is the lowestmass exoplanet measured through RV^{7}. It represents a new milestone that illustrates the capability of ESPRESSO to yield the masses of planets with RV signatures of about 10 cm s^{−1} in multiplanetary systems even when there is stellar activity. The relative precision of the RV semiamplitude of the other two previously known planets is also improved by a factor ~ 1.5 for planet c and by a factor of ~2 for planet d. We obtain a relative mass precision of 11% and 14% for planets c and d, respectively, which is the best precision for the mass measurement of superEarths around M dwarfs that can currently be achieved (Suárez Mascareño et al. 2020; LilloBox et al. 2020).
For the three transiting planets, we achieve bulk densities with relative precision of 46, 21, and 24% for planets b, c, and d, respectively.Considering the size and mass of these planets and the difficulties associated with a precise characterization of the mass and radius of M dwarfs, these density measurements are references for the field. Figure 8 shows these three planets in the massradius diagram and in the context of the known exoplanet population. These three planets are located below the radius gap (Fulton et al. 2017; Fulton & Petigura 2018; Cloutier & Menou 2020) and appear to be mostly rocky (see Sect. 5.3).
We also expand the view of this system with the discovery of a fourth planet and a planetary candidate. These planets do not transit, but with minimum masses of and M_{⊕}, they are probably both rocky planets or water worlds (also called ocean worlds, e.g., Adams et al. 2008). With an equilibrium temperature of K, the planetary candidate 5, if confirmed, would orbit in the habitable zone of its parent star.
Fig. 8
Massradius diagram of the smallplanet population. Each point represents a confirmed exoplanet for which the mass and radius are measured with a relative precision better than 50%. These data have been extracted from exoplanet.eu (Schneider et al. 2011). The shape of the points indicates the technique used to measure the mass of the planet: circles for RV, and squares for transit timing variations. The color of the point reflects the intensity of the incident flux. The level of transparency of the error bars indicates the relative precision of the planetary bulk density. The better the precision, the more opaque the error bars. The three transiting planets in the L 9859 system are labeled and appear circled in black. The labeled blue stars indicate the Solar System planets. The colored dashed lines are the massradius models from Zeng et al. (2016). The gray region indicates the maximum collision stripping of the mantle. The shaded horizontal blue line represent the radius gap (Fulton et al. 2017). L 9859 b is in a sparsely populated region of the parameter space and currently the lowestmass planet with a mass measured through RV. Lower planetary masses have all been measured through transittiming variation, e.g., for Trappist1 h (Gillon et al. 2017) to the left of L 9859 b. This plot has been produced using the code available at https://github.com/odemangeon/massradius_diagram. 
5.3 Internal composition of the three transiting superEarths
We performed a Bayesian analysis to determine the posterior distribution of the internal structure parameters of the planets. The method follows Dorn et al. (2015, 2017), and has been used in Mortier et al. (2020), Leleu et al. (2021), and Delrez et al. (2021). The model consists of two parts. The first is the forward model, which provides the planetary radius as a function of the internal structure parameters (iron molar fraction in the core, Si and Mg molar fraction in the mantle, mass fraction of all layers, age of the planet, and irradiation from the star), and the second part is the Bayesian analysis, which provides the posterior distribution of the internal structure parameters based on the observed radii, masses, and stellar parameters (in particular, its composition). The details of the analysis performed along with additional outputs are provided in Appendix G.
Figure 9 provides the ternary diagrams representing the posterior distributions of the composition of the three transiting planets in system L 9859. Furthermore, Figs. G.1 to G.3 provide the detailed posterior distributions of the most important parameters (mass fractions and composition of the mantle) of each planet. The three planets are characterized by small iron cores (12–14% in mass), which reflects the small iron abundance (compared to Si and Mg) in the star. According to the Bayesian analysis, the two innermost planets are likely to have a small mass fraction of water (the mode of the distribution is at 0) and a low gas mass, if they have any gas at all. Interestingly, the internal structure parameters of L 9859 d according to the Bayesian analysis are substantially different: the mode of the watermass fraction distribution is at ~ 0.3, whereas the mode of the gas mass peaks at ~10^{−6} M_{⊕}. Because the Bayesian analysis provides the joint distribution of all planetary parameters, we can easily compute the probability that the mass fraction of gas and water is higher in L 9859 d than in L 9859 b and L 9859 c. Based on our model, the values are 79.3% and 72.0% for gas and water, respectively, for planet d versus planet b. These values are 79.6% and 79.1% for gas and water, respectively, for planet d versus planet c. Planet d therefore appears to be likely richer in gas and water. On the other hand, planets b and c are very similar in composition. We emphasize finally that these numbers result from the Bayesian analysis, and they therefore depend on the assumed prior s that we took to be as uninformative as possible.
Our modeling favors a dry and hydrogen and heliumfree model for planet b and c. The posterior distributions of their gas and water content peak at 0, but the 3 σ confidence interval still allows for a water mass fraction of up to ~25% (see Figs. G.1 and G.2). In order to understand how promising planets b, c, and even d are for atmospheric characterization, we need to understand whether these warm planets (T_{eq} between ~ 400 and ~ 600K) could retain a waterdominated atmosphere. Providing a robust answer to this question requires modeling the complex phase diagram of water (e.g., French et al. 2009; Mousis et al. 2020; Turbet et al. 2020), the radiative transfer in a waterdominated atmosphere irradiated by an M star including potential runaway greenhouse effects (e.g., Arnscheidt et al. 2019), and the hydrodynamic escape of water potentially assisted by ultraviolet photolysis (e.g., Bourrier et al. 2017). This analysis is beyond the scope of this paper. However, we can study the example of the TRAPPIST1 system (Gillon et al. 2017; Luger et al. 2017) for comparison. Turbet et al. (2020) stressed the impact of irradiation on a waterdominated atmosphere. If the received irradiation is higher than the runaway greenhouse irradiation threshold (e.g., Kasting et al. 1993), which should be the case for TRAPPIST1 b to d (Wolf 2017), water should be in a steamed phase instead of a condensed phase, as classically assumed. In this case, the estimated water content of the planets decreases by several orders of magnitude. The authors further concluded that planets with masses lower than 0.5 M_{⊕} that are more irradiated than the runaway greenhouse irradiation threshold are probably unable to retain more than a few percent of water by mass due to an efficient hydrodynamic escape. For the TRAPPIST1 system, Bourrier et al. (2017) followed a theoretical study from Bolmont et al. (2017), however, in the attempt to assess the water loss experienced by the planets during their lifetime. The authors concluded that planets g and planets closer in could have lost up to 20 Earth oceans through hydrodynamic escape. However, they noted that depending on the exact efficiency of the photolysis, even TRAPPIST1 b and c could still harbor significant amounts of water.
L 9859 b is similar in mass and radius to TRAPPIST1 d. However, it is significantly more irradiated (T_{eq} = 288 ± 5.6 K for TRAPPIST1 d Gillon et al. 2017). L 9859 b might thus undergo or have undergone efficient hydrodynamic escape. L 9859 c and d are more massive than any of the TRAPPIST1 planets, but they are also more irradiated (T_{eq} = 400.1 ± 7.7 K for TRAPPIST1 b Gillon et al. 2017), and the comparison is thus less pertinent. They are likely to have undergone runaway greenhouse effect, but their higher masses might inhibit atmospheric escape. A more detailed study and observational evidence are thus required to reliably assess the nature and content of the atmosphere of the transiting planets in the L 9859 system.
Fig. 9
Ternary diagrams showing the internal composition (mass fractions of the gas (H and He), the volatile (water) and the refractory elements) for the three transiting planets in system L 9859. 
6 Conclusion: L 9859, a benchmark system for superEarth comparative exoplanetology around an M dwarf
Multiplanetary systems are ideal laboratories for exoplanetology because they offer the unique possibility of comparing exoplanets that formed in the same protoplanetary disk and are illuminated by the same star. According to the exoplanet archive^{8} (Akeson et al. 2013), we currently know 739 multiplanetary systems. A large fraction of them (~ 60%) were discovered by the Kepler survey (Borucki et al. 2010; Lissauer et al. 2011). From a detailed characterization and analysis of the properties of the Kepler multiplanetary systems, Weiss et al. (2018, hereafter W18) extracted the socalled “peas in a pod” configuration. The authors observed that consecutive planets in the same system tend to have similar sizes. The planets also appear to be preferentially regularly spaced. The authors also noted that the smaller the planets, the tighter their orbital configuration. For Table 2, we computed the metrics identified by W18 for the L 9859 system along with the distributions of these metrics derived by the authors from their sample. From Table 2, we conclude that system L 9859 closely follows the peas in a pod configuration. Most systems in the W18 sample have FGK host stars. For example, none of the 51 host stars that host four or more planets have a mass lower than 0.6 M_{⊙}. The fact that the L 9859 system also follows the peas in a pod configuration thus further strengthens the universality of this configuration and the constraints that it brings on planet formation theories. The only trend observed by W18 that the L 9859 system does not display is the positive correlation between the equilibrium temperature difference of consecutive planets and their radius ratio. Furthermore, assuming a v sin i_{*} of 1 km s^{−1}, the semiamplitude of the expected Rossiter–McLaughling effect (e.g., Queloz et al. 2000) is 40 cm s^{−1}, 1 m s^{−1}, and 37 m s^{−1} for planets b, c,and d, respectively. These amplitudes might at least for planets c and d be within the reach of highresolution spectrographs such as ESPRESSO. This would give us access to the spinorbital angle in this system, which would further constrain itsarchitecture and the possible mechanisms of its formation and migration.
The fact that L 9859 A is an M dwarf sets this system apart among multiplanetary systems. According to the exoplanet archive and the recent literature, only seven multiplanetary systems arecurrently confirmed (including L 9859) around M dwarfs for which the planetary masses and radius of at least two planets have been measured. The other six are TRAPPIST1 (Gillon et al. 2017), LTT3780 (Cloutier et al. 2020), TOI1266 (Demory et al. 2020), LHS1140 (LilloBox et al. 2020), K2146 (Hamann et al. 2019), and Kepler138 (JontofHutter et al. 2015). With a V magnitude of 11.7 and a distance of 10.6 pc, L 9859 is the brightest and closest of these systems.
Finally, according to the transmission spectrum metric (TSM, Kempton et al. 2018), with values of 49, 37, and 255 for planets b, c, and d, respectively, the three transiting planets in system L 9859 are comfortably above the thresholds proposed by Kempton et al. (2018) for a superEarth atmospheric characterization with the JWST. This threshold is 12 for planets with radii smaller than 1.5 R_{⊕}, like planets b and c, and 92 for planets with radii between 1.5 and 2.75 R_{⊕}, like planet d. Figure 10 shows the TSM values for the wellcharacterized smallplanet population. L 9859 b and c are the two planets with the highest TSM value below 1.5 R_{⊕}, and L 9859 d has the second highest above this value. These three planets thus belong to the most favorable warm to temperate (T_{eq} < 650 K) superEarths (R_{p} < 1.5 R_{⊕}) for an atmospheric characterization. Furthermore, L 9859 is located at the border of the continuous viewing zone (~ 200 days per year) of the JWST, making it a golden system for atmospheric characterization and comparative planetology. Even if the TSM is specifically tailored to the JWST, these planets are also suitable for transmission spectroscopy with other facilities such as ESPRESSO, the HST (Sirianni et al. 2005), NIRPS (Bouchy et al. 2017), or Ariel (Tinetti et al. 2016).
Peas in a pod statistics in system L 9859.
Fig. 10
Transmission spectrum metric vs. planetary radius diagram of the smallplanet population. Each point represents a confirmed exoplanets with mass and radius measured with a relative precision better than 50%. These data have been extracted from the exoplanet archive. The shape of the points indicates the technique used to measure the mass of the planet: circles for RV, and squares for transit timing variations. The color of the point reflects the equilibrium temperature of the planet. The level of transparency of the error bars indicates the relative precision of the planetary bulk density. The better the precision, the more opaque the error bars. The three transiting planets in system L 9859 are labeled and appear circled in black. We also display the names of the other planets with the highest transmission spectrum metrics. This plot has been produced using the code available at https://github.com/odemangeon/massradius_diagram. 
Acknowledgements
The authors acknowledge the ESPRESSO project team for its effort and dedication in building the ESPRESSO instrument. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This work has made use of data from the European Space Agency (ESA) mission (Gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC). 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 NASA/IPAC Infrared Science Archive, which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology. This publication makes use of VOSA, developed under the Spanish Virtual Observatory project supported by the Spanish MINECO through grant AyA201784089. VOSA has been partially updated by using funding from the European Union’s Horizon 2020 Research and Innovation Programme, under Grant Agreement no 776 403 (EXOPLANETSA). This work was supported by FCT  Fundação para a Ciência  through national funds and by FEDER through COMPETE2020  Programa Operacional Competitividade e Internacionalização by these grants: UID/FIS/04434/2019; UIDB/04434/2020; UIDP/04434/2020; PTDC/FISAST/32113/2017 and POCI010145FEDER032113; PTDC/FISAST/28953/2017 and POCI010145FEDER028953; PTDC/FISAST/28987/2017 and POCI010145FEDER028987; PTDC/FISAST/30389/2017 and POCI010145FEDER030389. O.D.S.D. is supported in the form of work contract (DL 57/2016/CP1364/CT0004) funded by FCT. J.I.G.H. acknowledges financial support from Spanish Ministry of Science and Innovation (MICINN) under the 2013 Ramón y Cajal programme RYC201314 875. J.I.G.H., A.S.M., R.R., and C.A.P. acknowledge financial support from the Spanish MICINN AYA201786 389P. A.S.M. acknowledges financial support from the Spanish Ministry of Science and Innovation (MICINN) under the 2019 Juan de la Cierva Programme. This work has been carried out with the support of the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). The authors acknowledge the financial support of the SNSF and in particular YA and JH acknowledge the SNSF for supporting research through the grant 200020_19203. N.J.N. acknowledges support form the following projects: CERN/FISPAR/0037/2019, PTDC/FISOUT/29048/2017. R. A. is a Trottier Postdoctoral Fellow and acknowledges support from the Trottier Family Foundation. This work was supported in part through a grant from FRQNT. A.A.K. and S.G.S. acknowledge the support from FCT in the form of the exploratory projects with references IF/00028/2014/CP1215/CT0002, IF/00849/2015/CP1273/CT0003. A.S. acknowledges support from the Italian Space Agency (ASI) under contract 201824HH.0. The financial contribution from the agreement ASIINAF n.201816HH.0 is gratefully acknowledged. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project FOUR ACES; grant agreement no. 724427). Used Simbad, Vizier, exoplanet.eu. Most of the analyses presented in this paper were performed using the Python language (version 3.6) available at python and several scientific packages: Numpy (van der Walt et al. 2011), Scipy (Virtanen et al. 2020), Pandas (McKinney 2010), Ipython (Pérez & Granger 2007), Astropy (Astropy Collaboration 2013, 2018) and Matplotlib (Hunter 2007).
Appendix A Characterization of the M dwarf L 9859 A
A.1 Atmospheric parameters of L 9859 A: Detailed description of the different methods
In addition to the derivation made by K19, we applied three different methods to derive the T_{eff}, log g and [Fe∕H] of L 9859 A.
A.1.1 Spectral synthesis with STEPARSYN
We employed the BTSettl model grid (Allard et al. 2012), the radiative transfer code turbospectrum (Plez 2012), and a VALD3based line list (Ryabchikova et al. 2015). The stellar atmospheric parameters of our selected set of synthetic spectra span between 2600 and 4500 K in T_{eff}, 4.0 to 6.0 dex in log g, and 1 to +0.5 dex in [Fe/H]. In addition, we took the instrumental broadening into account by means of a Gaussian kernel (R = 140 000). We used the latest version of the steparsyn code (Tabernero et al. 2018, 2021) to infer the stellar parameters. We fit the combined spectrum of L 9859 A that was constructed using 61 ESPRESSO spectra (S∕N = 1063 at 7580 Å). We selected the TiO band system at 7050 Å together with some Fe I and Ti I lines (see Marfil et al. 2020) to fit the observations. The latest version of STEPARSYN relies on emcee (ForemanMackey et al. 2013), a Markov chain Monte Carlo (MCMC) method used to fully sample the underlying distribution of the stellar parameters of L 9859. In addition to the T_{eff}, [Fe∕H], and log g values shown in Table A.2, the method also provides an estimate for the quadratic sum of the macroturbulence (ζ) and the stellar equatorial spin velocity projected on the plane of the sky (v sini): km s^{−1}.
A.1.2 Machinelearning regression with odusseas
The ODUSSEAS software (AntoniadisKarnavas et al. 2020) receives a 1D spectrum and its resolution as input. The pseudoequivalent widths are measured and used as input for a supervised machinelearning algorithm (ridge regression model) that is used to derive the spectroscopic parameters T_{eff} and [Fe∕H]. The implementation relies on the machinelearning Python package scikit learn. The training and testing sets were taken from a reference sample of 65 HARPS spectra with associated T_{eff} and [Fe∕H] derived by Casagrande et al. (2008) and Neves et al. (2012). When the spectra provided as input did not have the same resolution as the HARPS spectra from the reference sample, the spectra with the highest resolution were degraded (by convolution) to the lowest of the two resolutions. The estimates of T_{eff} and [Fe∕H] result from the average of 100 determinations obtained by randomly shuffling and splitting the training and testing groups. The reported uncertainties are the wide uncertainties of the machinelearning models at this resolution, after taking the intrinsic uncertainties of the reference sample parameters during the machinelearning process into consideration. The estimates provided by this method are also reported in Table A.2.
A.1.3 Spectral energy distribution fitting with vosa
The VOSA (Bayo et al. 2008) online tool estimates the T_{eff}, [Fe∕H], log g, extinction (A_{V}), and alpha enhancement by fitting the photometric SED with theoretical models. It also computes the total flux (F_{tot}) by integrating over the best template and then uses the distance to infer the luminosity (L). VOSA offers a wide variety of stellar models. We chose the BTSettl model (Allard et al. 2012) for its treatment of dust and clouds, which is important for lowmass stars. Because of the small distance of 10.6194 pc (inferred from GAIA parallaxes, BailerJones et al. 2018), we fixed the extinction to 0. The photometric measurements we used for the photometric SED are listed in Table A.1. The T_{eff}, [Fe∕H], and log g provided by this analysis are provided in Table A.2. Additionally, the fitting procedure inferred an alphaelement enhancement ([α/Fe]) of dex and a luminosity of L = 0.01128 ± 0.00042 L_{⊙}.
Broadband photometry of L 9859
A.1.4 K19 approach
K19 estimated T_{eff} and log g from two mostly independent derivations. T_{eff} was derived using the StefanBoltzman law. The required bolometric luminosity was estimated from V and Kband photometry using empirical bolometric correction relations (Pecaut & Mamajek 2013; Mann et al. 2015, erratum). For the radius, they used 0.312 ± 0.014 R_{⊙} derived from the massluminosity relation of Benedict et al. (2016) and the massradius relation of Boyajian et al. (2012). [Fe∕H] was derived from SED fitting (Stassun et al. 2017; Stassun & Torres 2016). This procedure also yielded an estimate of T_{eff} that was compatible within 1 σ with the previous one, but was not preferred by the authors.
A.1.5 Choice of the adopted set of atmospheric parameters
Table A.2 compiles the four estimates of the spectroscopic parameters of L 9859 A obtained with the four approaches presented above. It makes sense to separate them into two groups: the VOSA and K19 estimates, which rely on the photometric SED, on one side and the spectral synthesis and machinelearning estimates, which rely on the highresolution ESPRESSO spectra, on the other. For T_{eff}, the SED based estimates are similar in terms of best values and uncertainties. They are both compatible within 1 σ with the two ESPRESSO based estimates. However, the latter are 2.5 times more precise. The ESPRESSO based estimates provide similar uncertainties, but are only compatible at 2.25 σ. We do not currently know of any study that demonstrates the higher accuracy of one of the two ESPRESSO based approaches for M stars. Consequently, we did not exclude any of these estimates as an obvious outlier. However, judging from the data, the spectral synthesis and machinelearning uncertainties appear to be underestimated. For [Fe∕H], the four estimates are compatible within 1.6 σ. As expected, the spectral synthesis and machinelearning methods provide more precise estimates with uncertainties up to five times better. Finally, the two log g estimates provided bythe spectral synthesis and VOSA approaches are compatible within 1 σ. The spectral synthesis method provides a more accurate estimate because it uses data with high spectral resolution.
Different approaches to the spectroscopic parameters of L 9859
In this paper, which focuses on the characterization of the planets in system L 9859, we need to conclude with one final set of T_{eff}, [Fe∕H], and log g estimates. In order to keep a physically selfconsistent set of estimates, we decided to use the best values inferred by one method as the final best values for the three spectroscopic parameters. The use of highresolution spectroscopy data, which offers the possibility of directly characterizing the chromospheric lines, is clearly an asset for inferring [Fe∕H] and log g compared to the use of the photometric SED. The larger wavelength coverage toward the infrared offered by the SED can provide important constraints for the inference of T_{eff}. However, the T_{eff} estimates provided by the SED based methods include the estimates of the methods based on high spectral resolution within 1 σ. We thus decided to use one of the two methods based on high spectral resolution to obtain our set of best values. We chose the spectral synthesis method because of the lack of a benchmark analysis demonstrating the accuracy of the relatively recent machinelearning approach and because it does not provide an estimate for log g. The only exception was the uncertainties on the T_{eff}, which we identified as underestimated. We chose to enlarge this uncertainty to encompass the best values provided by the other three methods within 1 σ, leading to the adopted values and uncertainties provided in Table H.1.
A.2 Stellar modeling: Mass, radius, and age
The derivation of the radius of L 9859 A is already presented in detail in Sect. 3.2, but for the derivation of its mass, we again used several methods. The first method relies on our estimate of log g (see Sect. 3.1), which combined with our radius estimates provides a mass of . The second method relies on the stellar density retrieved by K19 from the fit of the transits of the three transiting planets. Combined with our radius estimates, it provides a mass of . Our third method relies on the massluminosity relation in K band of Mann et al. (2019). From the absolute K magnitude of 6.970 ± 0.019 mag, obtained from the observed magnitude provided by the 2MASS catalog (Cutri et al. 2003) and the distance provided by the Gaia collaboration (BailerJones et al. 2018), we obtain a mass of 0.290 ± 0.020 M_{⊙}. The fourth approach is based on the recently published studies of M dwarfs by Cifuentes et al. (2020, see in particular Table 6). The authors performed a comprehensive analysis of 1843 nearby bright lowmass star using SED photometry. They derived bolometric luminosities, effective temperature, radius, and mass for this sample. The masses are based on Schweitzer et al. (2019). They thus provide an equivalence between absolute bolometric luminosity, effective temperature, radius, and mass. Our bolometric luminosity estimate would indicate a radius of 0.343 ± 0.082 R_{⊙} and a mass of 0.338 ± 0.087 M_{⊙}. Our effective temperature estimate would indicate a radius of 0.433 ± 0.086 R_{⊙} and a mass of 0.432 ± 0.090 M_{⊙}. For our fifth approach, we used the VOSA online tools that are described in Appendix A.1.3. VOSA derives the stellar mass of 0.273 ± 0.030 M_{⊙} by comparing the measured T_{eff} and bolometric luminosity to evolutionary tracks (BTSettl model Allard et al. (2012) for consistency with our analysis of the photometric SED). Finally, K19 provided an estimate of the mass of 0.313 ± 0.014 M_{⊕} using the massluminosity relation for M dwarfs of Benedict et al. (2016). They derived the luminosity from Kband observations.
Table A.3 gathers all these estimates of the radius and mass of L 9859 A. From these, we also computed the resulting stellar densities through Monte Carlo simulations. We drew 100,000 samples of stellar mass and radius from normal distributions with mean and standard deviation as provided by the estimates from the corresponding row of Table A.3. When the error bars were asymmetric, we used the average of the upper and lower uncertainties as standard deviation. From these 100,000 samples, we computed 100,000 stellar density values. We then computed the estimate of the stellar density using the 50th, 16th and 84th percentiles. The relative precision on the stellar density provides us with a lower limit on the relative precision that we can achieve for the planetary density (see Table H.1). The absolute value of the stellar density will also impact the measured planetary densities and thus is of particular interest for modeling their interior (see Sect. 5.3). All stellar density estimates agree within 1 σ. However, the dispersion of best values shows that the one inferred from Cifuentes et al. (2020) using the T_{eff} is clearly off. The associated mass and radius are also significantly above all others. This might be due to the scale of the Cifuentes et al. (2020) study. The table from which we derive our estimates is a summary of the properties of around 2000 stars, which might be relevant for a large sample, but might fail to accurately represent a specific case such as that of L 9859 A. We thus discarded this estimate. We also note that the Cifuentes et al. (2020) estimates that are based on the bolometric luminosity (instead of the T_{eff}) agree well with the others.
Mass, radius, and density of L 9859 derived with different approaches
The remaining radius estimates agree within 1 σ, but their uncertainties vary by a factor of up to ~ 6 between the K19 estimate and the one from Cifuentes et al. (2020) based on the bolometric luminosity. As already mentioned in the previousparagraph, because of the scale of the Cifuentes et al. (2020) study, their uncertainties are probably overestimated. The uncertainties of the other two estimates differ by less than a factor two. We adopted the values derived from the StefanBoltzmann law because they are based on first principles.
The mass estimates also agree within 1 σ, but their uncertainties vary by a factor of up to ~ 10. Compared with the dispersion of the best values, the K19 uncertainty appears to be underestimated. The log g and stellar density based values and the Cifuentes et al. (2020) uncertainties, in contrast, appear to be overestimated. In between the two remaining estimates, VOSA and Mann et al. (2019), we adopted the estimate derived with VOSA. The VOSA tools provided T_{eff}, [Fe∕H], and log g values in good agreement with the value we adopted. We also used VOSA to derive the bolometric luminosity used to derive L 9859 A radius. The VOSA mass estimate thus provides a physically consistent set of stellar parameters. The final set of adopted values and uncertainties is provided in Table H.1.
To determine the age of L 9859 A, we used the accurate photometry and distance provided by Gaia. We constructed the colormagnitude diagram shown in Fig. 1, where we also depict the wellknown empirically determined mean sequences of stellar members of the β Pictoris moving group (~20 Myr, MiretRoig et al. 2020), the TucanaHorologium moving group (~45 Myr, Bell et al. 2015), the Pleiades open cluster (~120 Myr, Gossage et al. 2018), and the field (possible ages in the range 0.8–10 Gyr). These sequences were taken from Luhman (2018) and Cifuentes et al. (2020) and were derived by employing Gaia data; therefore the direct comparison with L 9859 A is feasible without any systematic effect. From its location in the Gaia colormagnitude diagram, we infer that L 9859 likely has an age that is consistent with that of the field (our target lies below the mean field sequence of M dwarfs). We did not correct L 9859 A data for interstellar extinction because based on its optical and infrared photometry (Table A.1) and optical HARPS and ESPRESSO spectroscopy, there is no evidence of strong or anomalous absorption. The field age is consistent with the measured mass and radius of the star, and the actual position of L 9859 A below the bottom borderline of the 1 σ dispersion of the field sequence also agrees with a slightly subsolar metallicity. Finally, the kinematics of L 9859 A can also provide indications about its age. Using the RV systemic velocity, the Gaia parallax, the RA/DEC coordinates and proper motions, we derived the UVW velocities of L 9859 A (see Table A.4). L 9859 A appears to belong to the thin disk and does not belong to any know young moving group. Therefore it is kinematically older than the oldest moving group currently known, that is, its age is older than 800 Myr.
Kinematics of L 9859 A
Appendix B Measuring radial velocities and activity indicators
Table B.1 provides the measurements of the RV s and activity indicators from the ESPRESSO spectrograph used in this paper. For the RV s and activity indicators measurements from the HARPS spectrograph, we refer to C19.
ESPRESSO RV, FWHM, BIS, contrast, S_{index}, H_{α}, NaD, and BERV measurements for L 9859
Appendix C Rotational modulation in photometric time series
In order to address the presence of stellar activity induced modulation in the TESS data, we first attempted to fit the LC with GP and mean offsets for each sector. The GP was implemented with the celerite Python package, as in Sect. 2.2, but this time, the functional form of the kernel was designed to model quasiperiodic signal. Its equation is (C.1)
and it is taken from ForemanMackey et al. (2017, eq. 56). P_{rot} is an estimator of the stellar rotation period, L is the correlation timescale, B is a positive amplitude term, and C is a positivefactor. We performed the fit by maximizing the log likelihood with emcee. We used 32 walkers. For each walker, we first maximized the log likelihood using the LBFGSB algorithm (Morales & Nocedal 2011; Zhu et al. 1997; Byrd et al. 1995) implemented in the scipy.optimize Python package. Then we performed a first exploration of 5000 iterations followed by a second exploration of 10 000 iterations starting at the last position of the first iteration.
The posterior PDF of the main hyperparameters (B, L and P_{rot}) is presented in Fig. C.1. The rotation period is poorly constrained ( days). It is also worth noting that the retrieved amplitude and timescales are low and short: ppm for the amplitude, and days for the timescales. In particular, appears too low to be physical because the timescale is expected to be on the same order of magnitude as or higher than the rotation period (Angus et al. 2018; Haywood et al. 2014). The timescale and the amplitude are strongly correlated. The short timescale and low amplitude can thus be tentatively explained by this degeneracy, which would result in a strong underestimation of both quantities.
Because the rotation period is only poorly determined, we used the GLSP as a more modelindependent approach to determining the rotational modulation in the TESS LC. TESS is designed for highprecision relative photometry (as opposed to highprecision absolute photometry). The photometry can thus show offsets between each sector that would strongly affect the GLSP of the LC. We thus used the offsets derived from our GP fit (see Table C.1) to realign the different sectors before we computed the GLSP. As the retrieved offsets are several orders of magnitude higher than the amplitude of the GP signal, we can assume that they are independent of the exact model used to describe the stellar activity (hyperparameters and choice of kernel). The result of the GLSP analysis is presented in Fig. 3 and discussed in Sect. 3.4.
Photometric offset derived for each TESS sector
Fig. C.1
posterior distributions of the main hyperparameters of the rotational kernel (Equation ((C.1))) 
Appendix D Choice of priors
The prior PDF used for the analyses described in Sects. 4.1.1, 4.1.2 and 4.2 are provided in Table H.1 (column prior). In this appendix, we explain the reasons for the choice of each prior.
D.1 Priors used for the TESS LC analysis (Sect. 4.1.1)
For the instrumental prior, the TESS additive jitter term (σ_{TESS}), we adopted a uniform distribution between zero and five times the median value of the reported error bars. The orbital parameters e cosω and e sin ω were assigned a joint prior. A joint prior consists of a transformation between two sets of parameters to define the prior on the new set of parameters instead. In this case, ecosω and esinω were converted into e and ω. For the prior PDF of e, as recommended by Kipping (2013), we used a beta distribution with the following values for the two shape parameters: a = 0.867 and b = 3.03. For the prior PDF of ω, we used a uniform distribution between − π and π. The remaining planetary parameters, P, t_{ic}, R_{p} ∕R_{*}, and cosi_{p}, were also assigned a joint prior. This joint prior, which we call transiting prior, also includes the stellar density ρ_{*}. Its main objective is to exclude regions of the parameter space where the three transiting planets are not transiting. It performs two changes of coordinates. It first computes the impact parameter (b) from P, ρ_{*}, and cos i_{p} (assuming a circular orbit), effectively converting the parameter cosi_{p} into b. Then it computes the orbital phase (ϕ) from P and t_{ic}. For this conversion, we need to define a reference time that corresponds to ϕ = 0. We chose this reference time to be the floored value of the first ESPRESSO observation, t_{ref} = 1436 BTJD. Then t_{ic} = t_{ref} + Pϕ. We thus transformed the set of parameters ρ_{*}, P, t_{ic}, R_{p} ∕R_{*}, and cos i_{p} into the new set of parameters ρ_{*}, P, ϕ, R_{p} ∕R_{*}, and b. To ρ_{*}, we assigned as prior the posterior of the K19 analysis. To P, we assigned a Jeffreys distribution between 0.1 day and the time span of the RV observations (~ 520 days). To avoid degenerate values of t_{ic} separated by a multiple of the period, we chose as prior a uniform distribution between zero and one for ϕ. For R_{p} ∕R_{*}, we assigned a uniform distribution between 10^{−3} and 1. For the prior of b, we used a uniform distribution between 0 an 2 in order to allow grazing transiting, but we imposed the condition that b < 1 + R_{p}∕R_{*} to ensure that the configuration is transiting.
Finally, for the prior on the limbdarkening coefficients, we used Gaussian PDF s whose first two moments were defined using the Python package ldtk^{11} (Parviainen & Aigrain 2015). Using a library of synthetic stellar spectra, it computes the limbdarkening profile of a star that is observed in a given spectral bandpass (specified by its transmission curve), and defined by its T_{eff}, log g, and [Fe∕H]. Provided the values and error bars for these stellar parameters (see Sect. 3.1) and the spectral bandpass of TESS, ldtk uses an MCMC algorithm to infer the mean and standard deviation of the Gaussian PDF s for the coefficients of a given limbdarkening law (nonlinear in our case). ldtk relies on the library of synthetic stellar spectra generated by Husser et al. (2013). It covers the wavelength range from 500 Å to 5.5 μm and the stellar parameter space delimited by 2 300 K ≤ T_{eff} ≤ 12 000 K, 0.0 ≤ log g ≤ +6.0, − 4.0 ≤ [Fe∕H] ≤ +1.0, and − 0.2 ≤ [α∕Fe] ≤ +1.2. This parameter space is well within the requirements of our study (see Table A.2).
D.2 Priors used for the RV analysis (Sect. 4.1.2)
Regarding the instrumental prior s, the prior PDF of the offsets between the RV instruments (ΔRV_{HARPS∕pre}, ΔRV_{post∕pre}) are Gaussian distributions with means equal to the difference of the median values of the data sets and variances equal to the sum of their variances. The prior PDF s of the constant levels of the FWHM (C_{pre}, C_{post}, and C_{HARPS}) are Gaussian distributions with means equal to the median values of each data set and variances equal to their variances. The prior PDF of the additive jitter parameters (σ_{RV,pre}, σ_{RV,post}, σ_{RV,HARPS}, σ_{FWHM,pre}, σ_{FWHM,post}, and σ_{FWHM,HARPS}) are uniform distributions between zero and five times the median values of the reported error bars for each data set. Regarding the star related prior s, the prior PDF of the systemic velocity (v_{0}) is a Gaussian with the mean equal to the median value of the RV data taken by ESPRESSO before the fiber change and a variance equal to its variance. The other parameters are the hyperparameters of the quasiperiodic kernels. The prior PDF s of the two amplitudes (A_{RV}, A_{FWHM}) are uniform between zero and the maximum of the peaktopeak values of the joint data sets taken by the three instruments. For the period of recurrence (P_{rot}), the prior PDF chosen is a Jeffreys distribution between 5 days and the time span of our observations (~ 520 days). Considering the age and the spectral type of L 9859, 5 days appears to be a good lower limit for the rotation period. This prior comfortably encompasses the estimate of ~ 80 days made by C19 based on the periodogram of the H_{α} measurements. For the decay timescale (τ_{decay}), we chose a Jeffreys distribution between 2.5 days and five times the time span of observations. This upper limit was set to prevent the GP from producing stellar activity models that would be completely coherent over the time span of our observations. In other words, we imposed that the stellar activity signal is quasiperiodic and not periodic. The objective was to avoid that the GP reproduces planetary signals. Furthermore, we imposed that the decay timescale was superior to half of the period of recurrence. This condition, suggested by Angus et al. (2018) and Haywood et al. (2014), prevents the GP from producing stellar activity signals that are too incoherent and thus close to white noise. In these cases, the GP signal and the additive jitter terms start to become degenerate. The prior PDF of the periodic coherence scale (γ) was uniform between 0.05 and 5. The typical value for γ in the literature is thought to be 0.5 (Dubber et al. 2019). This prior is designed to explore one order of magnitude below and above this typical values.
Regarding the planetary prior s, the prior PDF of K is uniform between 0 and the maximum of the peaktopeak values of the RV data sets taken by the three instruments. For the ephemeris parameters (P and t_{ic}) and for the three known transiting planets, we used as prior s the posterior s of our analysis of the TESS LC (see notes ‡ at the end of Table H.1).
For the nontransiting planets that we identified in the GLSP, we used a joint prior. This joint prior converts P and t_{ic} into P and ϕ, similarly to what was done within the transiting joint prior in Appendix D.1. The reference time used, which corresponds to ϕ = 0, is the same (t_{ref} = 1436 BTJD). We chose a uniform distribution between zero and one for ϕ. For P, we used a Jeffreys distribution between 0.1 day and the time span of the RV observations (~ 520 days). Finally, the last two parameters are ecosω and esinω. We used the same joint prior as in Appendix D.1, which results in a beta distribution with shape parameters a = 0.867 and b = 3.03 for the prior PDF of e (Kipping 2013), and a uniform distribution between − π and π for ω.
D.3 Priors used for the joint analysis of the RV and LC data (Sect. 4.2)
The prior s used for this analysis are the same as were used for the analysis of the TESS LC (see Appendix D.1). For the parameters that are not present in this analysis, we used the same prior s as we used for our analysis of the RV data (see Appendix D.2). All prior s are mentioned in Table H.1.
Appendix E Searching for the transits of planet e and planetary candidate 5
We searched the TESS data for previously unreported planetary transit signal including planet e and planetary candidate 5. We used a procedure similar to Barros et al. (2016). For this analysis, we did not use the LC detrended with a GP described in Sect. 2.2 because the flexibility of the GP might alter the transit signals. Instead we detrended each sector separately by dividing the LC by a spline interpolation of third degree. We used a knot every 0.5 days. Combined with an iterative 3 σ clipping to identify outliers, this allowed us to better preserve unidentified transits signals in the detrended LC (Barros et al. 2016). Then we removed the transits of the three known transiting planets by cutting out data within a window of two transit durations centered on the predicted transit time. The additional 0.5 transit duration before and after transit allows accounting for errors in the ephemerides or unknown transittiming variations. We then performed a box leastsquares (BLS) search (Kovács et al. 2002) to find periodicities between 0.5 and 40 days. The resulting periodogram is shown in Fig. E.2, with the highest peak corresponding to 1.049 days, which is probably due to aliases linked to the rotation of Earth. Phasefolding the LC at this period does not show a typical transit signature. No other significant peaks are seen in the BLS periodogram, including at the periods of the candidate planets detected in RV (see Sect. 4.1.3). We also performed a transit search using the TLS software (Hippke & Heller 2019) and obtained the same conclusion.
To confirm the absence of transit signal for planet e and planetary candidate 5, we phasefolded the TESS LC using the ephemeris of Table H.1. We do not observe any transit signal in either case. We show that if the planetary radii are similar to those of the other transiting planets, the transit signal would have been clear in the TESS LC.
Fig. E.1
Phasefolded TESS LC assuming the best model ephemerides of planet e (a) and planetary candidate 5 (b). The black points are the TESS data points at the original cadence. The red line is the data binned in phase using bins of 15 min. The dashed pink and brown lines are the expected transit signal assuming that the planets have the same radius as planet d (see Table H.1). 
Fig. E.2
Periodogram provided by the BLS search in the TESS data. The dashed vertical pink and brown links indicate the orbital period of planet e and planetary candidate, respectively. There is no significant power at these periods. 
Appendix F Evidence for additional planets in the L 9859 system
As mentioned in Sect. 4.1.3, in order to assess the presence of additional planets in the L 9859 system, we first performed the two analyses that included only the three previously known planets. Figures F.1 and F.2 follow the same format as Figs. 4 and 5. Figure F.1 shows the RV time series including the data from both instruments, the best three planets plus activity model, and the residuals of this fit. Figure F.2 displays the GLSP of the combined RV data and the residuals, the GLSP s of the planetary and stellar activity model sampled at the same times as the RV time series, and the WF. The GLSP of the combined RV s in Fig. F.2 shows two narrow peaks with an FAP below 0.1 % at the two periods, 13 and 23 days, which were previously identified as potential additional planetary signals. The GLSP of the residuals displays a narrow peak at 13 days.
The analyses with four planets converges toward a significant detection of the semiamplitude of a fourth Keplerian signal. Similarly to Fig. F.1, Fig. 4 shows the time series and the best model, and Fig. 5 shows the GLSP s. We also performed an iterative GLSP analysis in Fig. F.5. This allows showing the peak on the GLSP that corresponds to planet b, which is invisible in other figures. The GLSP of the residuals after the subtraction of the model for planet b (also shown in Fig. 5) shows two peaks around 1.743 and 2.341 days and no peak around 23 days. The analysis of the TESS LC did not show transit signals at 1.743 or 2.341 days, so that we did not pursue a planetary origin for these peaks. However, the GLSP of the activity model does show a peak around 23 days. This indicates that the signal at 23 days might be generated by stellar activity. Based on our stellar activity model, which analyzes the FWHM data simultaneously with the RV data, we can also analyze the behavior of this activity indicator. Figures F.3 and F.4 show similar information to that in Figs. 4 and 5, but for the FWHM data. There is no significant power around 23 days in the GLSP of the combined FWHM data or in those of the stellar activity model and the residuals. Similarly, the GLSP s of all the other activity indicators (see Figs. Fig. 2 and Fig. 3) do not display significant power around 23 days. The analysis of the activity indicator does not confirm the stellar activity origin of the 23day signal.
Consequently, we performed other analyses with five planets. The fits converge toward a significant detection of the semiamplitude of a fifth Keplerian signal. Figures F.6 and F.7 show the time series, the best model, and the GLSP s. The GLSP of the stellar activity model still displays power around 23 days, but it is less significant and has a much more flattened profile compared with the fourplanet analyses (Fig. 5).
Fig. F.1
Outcome of the fit of the threeplanet model: The format of this figure is identical to that in Fig. 4, but isdescribed again here for convenience. (Top left) RV time series along with the best model (solid green line)that includes the planetary signals and best prediction from the GP stellar activity model. The 1 σ uncertainties from the GP prediction are also displayed (shaded green area). For this plot, we subtracted the systemic velocity and the instruments offsets from the RV data (see values in Table H.1). (Bottom left) Time series of the residuals of the best model. (Right) Zoom on a small portion of the time series to better visualize the shorttimescale variations. 
Fig. F.2
Outcome of the fit of the threeplanet model: The format of this figure is identical to that in Fig. 5, but is described again here for convenience. GLSP s of the RV time series (top) and of the planetary (second) and stellar activity (third) models sampled at the same times as the RV data. GLSP of the time series of the residuals (fourth) and the window function (bottom). The vertical lines on the GLSP s correspond to the orbital periods of planets b, c, and d, half and thefull rotation period (estimated at 80 days) from right to left. 
Fig. F.3
Outcome of the fit of the fourplanet model regarding the FWHM: The structure of this figure is similar to that of Fig. 4 or Fig. F.1, but the FWHM data and model are displayed instead of the RV data and model. 
Fig. F.4
Outcome of the fit of the fourplanet model regarding the FWHM: The structure of this figure is similar to that of Fig. 5 or Fig. F.2, but the FWHM data and model are displayed instead of the RV data and model. 
Fig. F.5
Iterative GLSP for the four planets model: GLSP of the RV data (top) and the window function (bottom). The GLSP s of the data shown in the previous row minus the model for planets c, e, and d, the stellar activity model, and the model for planet b are displayed in the second, third, fourth, fifth, and sixth row, respectively. 
Fig. F.6
Outcome of the fit of the fiveplanet model: The format of this figure is identical to that in Figs. 4 and F.1. 
Fig. F.7
Outcome of the fit of the fiveplanet model: The format of this figure is identical to that in Figs. 5 and F.2. 
Appendix G Internal composition of three transiting superearths
As explained in Sect. 5.3, our framework for modeling the interior of the three transiting planets is composed of a forward model and a Bayesian retrieval. In the forward model, each planet is made of four layers: an iron or sulfur inner core, a mantle, a water layer, and a gas layer. We used for the core the equation of state (EOS) of Hakim et al. (2018), for the silicate mantle, the EOS of Sotin et al. (2007), and the water EOS was taken from Haldemann et al. (2020). These three layers constitute the solid part of the planets. The thickness of the gas layer (assumed to be made of pure H and He) was computed as a function of the stellar age, mass, and radius of the solid part and irradiation from the star using the formulas of Lopez & Fortney (2014).
In the Bayesian analysis part of model, we proceeded in two steps. We first generated 150000 synthetic stars, their mass, radius, effective temperature, age, and composition ([Si/H], [Fe/H], and [Mg/H]), as well as the associated error bars, which were taken at random following the stellar parameters quoted above. For each of these stars, we generated 1000 planetary systems for which we varied the internal structure parameters of all planets, and we assumed that the bulk Fe/Si/Mg molar ratios are equal to the stellar ratios. We then computed the transit depth and RV semiamplitude for each of the planets and retained models that fit the observed data within the error bars. With this procedure, we included the fact that all synthetic planets orbit a star with exactly the same parameters. Planetary masses and radii are correlated by the fact that the fitted quantities are the transit depth and RV semiamplitude, which depend on the stellar radius and mass. In order to take this correlation into account, it is therefore important to fit the planetary system at once, and not each planet independently.
The prior s used in the Bayesian analysis are the following: The mass fraction of the gas envelope is uniform in log, the mass fraction (relative to the solid planet, i.e., excluding the mass of gas) of the inner core, mantle, and water layer are uniform on the simplex (the surface on which they add up to one). Finally, we constrain the mass fraction of water to be 50 % at most (Thiabaud et al. 2014; Marboeuf et al. 2014). The molar fraction of iron in the inner core is uniform between 0.5 and 1, and the molar fraction of Si, Mg, and Fe in the mantle is uniform on the simplex (they add up to one).
The posterior distributions of the most important parameters (mass fractions and composition of the mantle) of each planet in L 9859 are shown in Figs. G.1 to G.3.
Fig. G.1
Corner plot showing the main internal structure parameters of L 9859 b. We show the mass fraction of the inner core, the mass fraction of water, the Si and Mg mole fraction in the mantle, the Fe mole fraction in the inner core, and the mass of gas (log scale). The values at the top of each column are the mean and 5% and 95% quantiles. 
Appendix H Additional table
Parameter estimates of the planetary system L 9859.
References
 Adams, E. R., Seager, S., & ElkinsTanton, L. 2008, ApJ, 673, 1160 [Google Scholar]
 Adibekyan, V. Z., Sousa, S. G., Santos, N. C., et al. 2012, A&A, 545, A32 [Google Scholar]
 Adibekyan, V., Figueira, P., & Santos, N. C. 2017, in EWASS Special Session, 4 (2017): Star–Planet Interactions (EWASSSS42017) [Google Scholar]
 Akeson, R. L., Chen, X., Ciardi, D., et al. 2013, PASP, 125, 989 [Google Scholar]
 Albrecht, S., Winn, J. N., Marcy, G. W., et al. 2013, ApJ, 771, 11 [Google Scholar]
 Allard, F., Homeier, D., & Freytag, B. 2012, Philos. Trans. Roy. Soc. London A, 370, 2765 [Google Scholar]
 Ambikasaran, S., ForemanMackey, D., Greengard, L., Hogg, D. W., & O’Neil, M. 2015, IEEE Trans. Pattern Anal. Mach. Intell., 38 [Google Scholar]
 Angus, R., Morton, T., Aigrain, S., ForemanMackey, D., & Rajpaul, V. 2018, MNRAS, 474, 2094 [Google Scholar]
 AntoniadisKarnavas, A., Sousa, S. G., DelgadoMena, E., et al. 2020, A&A, 636, A9 [Google Scholar]
 Arnscheidt, C. W., Wordsworth, R. D., & Ding, F. 2019, ApJ, 881, 60 [Google Scholar]
 Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33 [Google Scholar]
 Astropy Collaboration (PriceWhelan, A. M., et al.) 2018, AJ, 156, 123 [Google Scholar]
 AstudilloDefru, N., Díaz, R. F., Bonfils, X., et al. 2017a, A&A, 605, L11 [EDP Sciences] [Google Scholar]
 AstudilloDefru, N., Forveille, T., Bonfils, X., et al. 2017b, A&A, 602, A88 [Google Scholar]
 BailerJones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58 [Google Scholar]
 Baluev, R. V. 2009, MNRAS, 393, 969 [Google Scholar]
 Baranne, A., Queloz, D., Mayor, M., et al. 1996, A&AS, 119, 373 [Google Scholar]
 Barros, S. C. C., Almenara, J. M., Demangeon, O., et al. 2015, MNRAS, 454, 4267 [Google Scholar]
 Barros, S. C. C., Demangeon, O., & Deleuil, M. 2016, A&A, 594, A100 [Google Scholar]
 Barros, S. C. C., Demangeon, O., Díaz, R. F., et al. 2020, A&A, 634, A75 [EDP Sciences] [Google Scholar]
 Batalha, N. M., Rowe, J. F., Bryson, S. T., et al. 2013, ApJS, 204, 24 [Google Scholar]
 Bayo, A., Rodrigo, C., Barrado Y Navascués, D., et al. 2008, A&A, 492, 277 [Google Scholar]
 Bell, C. P. M., Mamajek, E. E., & Naylor, T. 2015, MNRAS, 454, 593 [Google Scholar]
 Belokurov, V., Penoyre, Z., Oh, S., et al. 2020, MNRAS, 496, 1922 [Google Scholar]
 Benedict, G. F., Henry, T. J., Franz, O. G., et al. 2016, AJ, 152, 141 [Google Scholar]
 Bolmont, E., Selsis, F., Owen, J. E., et al. 2017, MNRAS, 464, 3728 [Google Scholar]
 Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977 [Google Scholar]
 Bouchy, F., Pepe, F., & Queloz, D. 2001, A&A, 374, 733 [Google Scholar]
 Bouchy, F., Doyon, R., Artigau, É., et al. 2017, The Messenger, 169, 21 [NASA ADS] [Google Scholar]
 Bourrier, V., de Wit, J., Bolmont, E., et al. 2017, AJ, 154, 121 [Google Scholar]
 Boyajian, T. S., von Braun, K., van Belle, G., et al. 2012, ApJ, 757, 112 [Google Scholar]
 Byrd, R. H., Lu, P., Nocedal, J., & Zhu, C. 1995, SIAM J. Sci. Comput., 16, 1190 [Google Scholar]
 Casagrande, L., Flynn, C., & Bessell, M. 2008, MNRAS, 389, 585 [Google Scholar]
 Cifuentes, C., Caballero, J. A., CortésContreras, M., et al. 2020, A&A, 642, A115 [Google Scholar]
 Cloutier, R., & Menou, K. 2020, AJ, 159, 211 [Google Scholar]
 Cloutier, R., AstudilloDefru, N., Bonfils, X., et al. 2019, A&A, 629, A111 [Google Scholar]
 Cloutier, R., Eastman, J. D., Rodriguez, J. E., et al. 2020, AJ, 160, 3 [Google Scholar]
 Cosentino, R., Lovis, C., Pepe, F., et al. 2012, Proceedings of the SPIE: Groundbased and Airborne Instrumentation for Astronomy IV, 8446, 84461V [Google Scholar]
 Crane, J. D., Shectman, S. A., & Butler, R. P. 2006, Proceedings of the SPIE: Groundbased and Airborne Instrumentation for Astronomy, 6269, 626931 [Google Scholar]
 Crane, J. D., Shectman, S. A., Butler, R. P., Thompson, I. B., & Burley, G. S. 2008, Proceedings of the SPIE: Groundbased and Airborne Instrumentation for Astronomy II, 7014, 701479 [Google Scholar]
 Crane, J. D., Shectman, S. A., Butler, R. P., et al. 2010, Proceedings of the SPIE: Groundbased and Airborne Instrumentation for Astronomy III, 7735, 773553 [Google Scholar]
 Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog: II/246 [Google Scholar]
 Delisle, J.B. 2017, A&A, 605, A96 [Google Scholar]
 Delrez, L., Ehrenreich, D., Alibert, Y., et al. 2021, Nat. Astron., 5, 775 [Google Scholar]
 Demory, B.O., Pozuelos, F. J., Gómez Maqueo Chew, Y., et al. 2020, A&A, 642, A49 [Google Scholar]
 Díaz, R. F., Cincunegui, C., & Mauas, P. J. D. 2007, MNRAS, 378, 1007 [Google Scholar]
 Díaz, R. F., Almenara, J. M., Santerne, A., et al. 2014, MNRAS, 441, 983 [Google Scholar]
 Dorn, C., Khan, A., Heng, K., et al. 2015, A&A, 577, A83 [Google Scholar]
 Dorn, C., Venturini, J., Khan, A., et al. 2017, A&A, 597, A37 [Google Scholar]
 Dubber, S. C., Mortier, A., Rice, K., et al. 2019, MNRAS, 490, 5103 [Google Scholar]
 Dumusque, X., Borsa, F., Damasso, M., et al. 2017, A&A, 598, A133 [Google Scholar]
 ForemanMackey, D. 2018, Res. Notes Am. Astron. Soc., 2, 31 [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
 ForemanMackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220 [Google Scholar]
 French, M., Mattsson, T. R., Nettelmann, N., & Redmer, R. 2009, Phys. Rev. B, 79, 054107 [Google Scholar]
 Fulton, B. J., & Petigura, E. A. 2018, AJ, 156, 264 [Google Scholar]
 Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109 [Google Scholar]
 Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. 2018, PASP, 130, 044504 [Google Scholar]
 Gaia Collaboration 2018, VizieR Online Data Catalog: I/345 [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [Google Scholar]
 Gaidos, E., Mann, A. W., Lépine, S., et al. 2014, MNRAS, 443, 2561 [Google Scholar]
 Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Sci. Rev., 123, 485 [Google Scholar]
 Geweke, J. 1992, Bayesian Stat., 4, 641 [Google Scholar]
 Gillon, M., Jehin, E., Magain, P., et al. 2011, EPJ Web of Conferences: Detection and Dynamics of Transiting Exoplanets, St. Michel l’Observatoire, France, eds. F. Bouchy, R. Díaz & C. Moutou, 11, 06002 [Google Scholar]
 Gillon, M., Triaud, A. H. M. J., Demory, B.O., et al. 2017, Nature, 542, 456 [Google Scholar]
 Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65 [Google Scholar]
 Gossage, S., Conroy, C., Dotter, A., et al. 2018, ApJ, 863, 67 [Google Scholar]
 Greene, T. P., Line, M. R., Montero, C., et al. 2016, ApJ, 817, 17 [Google Scholar]
 Gregory, P. C. 2005, Bayesian Logical Data Analysis for the Physical Sciences: A Comparative Approach with ‘Mathematica’ Support, ed. P. C. Gregory (Cambridge University Press) [Google Scholar]
 Grunblatt, S. K., Howard, A. W., & Haywood, R. D. 2015, ApJ, 808, 127 [Google Scholar]
 Hakim, K., Rivoldini, A., Van Hoolst, T., et al. 2018, Icarus, 313, 61 [Google Scholar]
 Haldemann, J., Alibert, Y., Mordasini, C., & Benz, W. 2020, A&A, 643, A105 [Google Scholar]
 Hamann, A., Montet, B. T., Fabrycky, D. C., Agol, E., & Kruse, E. 2019, AJ, 158, 133 [Google Scholar]
 Hara, N. C., Bouchy, F., Stalport, M., et al. 2020, A&A, 636, L6 [Google Scholar]
 Haywood, R. D., Collier Cameron, A., Queloz, D., et al. 2014, MNRAS, 443, 2517 [Google Scholar]
 Henden, A. A., Templeton, M., Terrell, D., et al. 2016, VizieR Online Data Catalog: II/336 [Google Scholar]
 Hippke, M., & Heller, R. 2019, A&A, 623, A39 [Google Scholar]
 Howell, S. B., Sobeck, C., Haas, M., et al. 2014, PASP, 126, 398 [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Husser, T.O., Wendevon Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6 [Google Scholar]
 Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, Proceedings of SPIE: Software and Cyberinfrastructure for Astronomy IV, 9913, 99133E [Google Scholar]
 Jönsson, H., Holtzman, J. A., Allende Prieto, C., et al. 2020, AJ, 160, 120 [Google Scholar]
 JontofHutter,D., Rowe, J. F., Lissauer, J. J., Fabrycky, D. C., & Ford, E. B. 2015, Nature, 522, 321 [Google Scholar]
 Kaltenegger, L. 2017, A&ARv, 55, 433 [Google Scholar]
 Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108 [Google Scholar]
 Kempton, E. M.R., Bean, J. L., Louie, D. R., et al. 2018, PASP, 130, 114401 [Google Scholar]
 Kipping, D. M. 2013, MNRAS, 434, L51 [Google Scholar]
 Kostov, V. B., Schlieder, J. E., Barclay, T., et al. 2019, AJ, 158, 32 [Google Scholar]
 Kovács, G., Zucker, S., & Mazeh, T. 2002, A&A, 391, 369 [Google Scholar]
 Kreidberg, L. 2015, PASP, 127, 1161 [Google Scholar]
 Kreidberg, L., Koll, D. D. B., Morley, C., et al. 2019, Nature, 573, 87 [Google Scholar]
 Leleu, A., Alibert, Y., Hara, N. C., et al. 2021, A&A, 649, A26 [Google Scholar]
 LilloBox, J., Figueira, P., Leleu, A., et al. 2020, A&A, 642, A121 [Google Scholar]
 Lissauer, J. J., Ragozzine, D., Fabrycky, D. C., et al. 2011, ApJS, 197, 8 [Google Scholar]
 Lopez, E. D., & Fortney, J. J. 2014, ApJ, 792, 1 [Google Scholar]
 Lovis, C., Dumusque, X., Santos, N. C., et al. 2011, ArXiv eprints [arXiv:1107.5325] [Google Scholar]
 Luger, R., Sestovic, M., Kruse, E., et al. 2017, Nat. Astron., 1, 0129 [Google Scholar]
 Luhman, K. L. 2018, AJ, 156, 271 [Google Scholar]
 Maldonado, J., Micela, G., Baratella, M., et al. 2020, A&A, 644, A68 [EDP Sciences] [Google Scholar]
 Mandt, K., Mousis, O., & Chassefière, E. 2015, Icarus, 254, 259 [Google Scholar]
 Mann, A. W., Feiden, G. A., Gaidos, E., Boyajian, T., & von Braun, K. 2015, ApJ, 804, 64 [Google Scholar]
 Mann, A. W., Dupuy, T., Kraus, A. L., et al. 2019, ApJ, 871, 63 [Google Scholar]
 Marboeuf, U., Thiabaud, A., Alibert, Y., Cabral, N., & Benz, W. 2014, A&A, 570, A36 [Google Scholar]
 Marfil, E., Tabernero, H. M., Montes, D., et al. 2020, in Contributions to the XIV.0 Scientific Meeting (virtual) of the Spanish Astronomical Society, 156 [Google Scholar]
 Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20 [NASA ADS] [Google Scholar]
 McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference (Stéfan van der Walt and Jarrod Millman), 51 [Google Scholar]
 Millholland, S., Wang, S., & Laughlin, G. 2017, ApJ, 849, L33 [NASA ADS] [CrossRef] [Google Scholar]
 MiretRoig, N., Galli, P. A. B., Brandner, W., et al. 2020, A&A, 642, A179 [Google Scholar]
 Morales, J. L., & Nocedal, J. 2011, ACM Trans. Math. Softw., 38, 7:1 [Google Scholar]
 Morley, C. V., Kreidberg, L., Rustamkulov, Z., Robinson, T., & Fortney, J. J. 2017, ApJ, 850, 121 [Google Scholar]
 Mortier, A., Zapatero Osorio, M. R., Malavolta, L., et al. 2020, MNRAS, 499, 5004 [Google Scholar]
 Mousis, O., Deleuil, M., Aguichine, A., et al. 2020, ApJ, 896, L22 [Google Scholar]
 Nelder, J. A., & Mead, R. 1965, Comput. J., 7, 308 [Google Scholar]
 Neves, V., Bonfils, X., Santos, N. C., et al. 2012, A&A, 538, A25 [Google Scholar]
 Noyes, R. W., Weiss, N. O., & Vaughan, A. H. 1984, ApJ, 287, 769 [Google Scholar]
 Parviainen, H., & Aigrain, S. 2015, MNRAS, 453, 3821 [Google Scholar]
 Pecaut, M. J.,& Mamajek, E. E. 2013, ApJS, 208, 9 [Google Scholar]
 Penoyre, Z., Belokurov, V., Wyn Evans, N., Everall, A., & Koposov, S. E. 2020, MNRAS, 495, 321 [Google Scholar]
 Pepe, F., Cristiani, S., Rebolo, R., et al. 2021, A&A, 645, A96 [CrossRef] [EDP Sciences] [Google Scholar]
 Pérez, F., & Granger, B. E. 2007, Comput. Sci. Eng., 9, 21 [Google Scholar]
 Perrakis, K., Ntzoufras, I., & Tsionas, E. G. 2014, Comput. Stat. Data Anal., 77, 54 [Google Scholar]
 Plez, B. 2012, Astrophysics Source Code Library [record ascl:1205.004] [Google Scholar]
 Queloz, D., Eggenberger, A., Mayor, M., et al. 2000, A&A, 359, L13 [NASA ADS] [Google Scholar]
 Queloz, D., Henry, G. W., Sivan, J. P., et al. 2001, A&A, 379, 279 [Google Scholar]
 Quillen, A. C., & French, R. S. 2014, MNRAS, 445, 3959 [Google Scholar]
 Quirrenbach, A., Amado, P. J., Caballero, J. A., et al. 2014, Proceedings of the SPIE: Groundbased and Airborne Instrumentation for Astronomy V, 9147, 12 [Google Scholar]
 Rainer, M., Borsa, F., & Affer, L. 2020, Exp. Astron., 49, 73 [Google Scholar]
 Rein, H. 2012, MNRAS, 427, L21 [NASA ADS] [Google Scholar]
 Rein, H., & Liu, S.F. 2012, A&A, 537, A128 [Google Scholar]
 Rein, H., & Tamayo, D. 2015, MNRAS, 452, 376 [Google Scholar]
 Rice, K., Malavolta, L., Mayo, A., et al. 2019, MNRAS, 484, 3731 [Google Scholar]
 Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, J. Astron. Telesc. Instrum. Syst., 1, 014003 [Google Scholar]
 Ryabchikova, T., Piskunov, N., Kurucz, R. L., et al. 2015, Phys. Scr, 90, 054005 [Google Scholar]
 Schneider, J., Dedieu, C., Le Sidaner, P., Savalle, R., & Zolotukhin, I. 2011, A&A, 532, A79 [Google Scholar]
 Schweitzer, A., Passegger, V. M., Cifuentes, C., et al. 2019, A&A, 625, A68 [Google Scholar]
 Sirianni, M., Jee, M. J., Benítez, N., et al. 2005, PASP, 117, 1049 [Google Scholar]
 Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163 [Google Scholar]
 Smith, J. C., Stumpe, M. C., Van Cleve, J. E., et al. 2012, PASP, 124, 1000 [Google Scholar]
 Sotin, C., Grasset, O., & Mocquet, A. 2007, Icarus, 191, 337 [Google Scholar]
 Stassun, K. G., & Torres, G. 2016, AJ, 152, 180 [Google Scholar]
 Stassun, K. G., Collins, K. A., & Gaudi, B. S. 2017, AJ, 153, 136 [Google Scholar]
 Stumpe, M. C., Smith, J. C., Catanzarite, J. H., et al. 2014, PASP, 126, 100 [Google Scholar]
 Suárez Mascareño, A., Faria, J. P., Figueira, P., et al. 2020, A&A, 639, A77 [Google Scholar]
 Tabernero, H. M., Dorda, R., Negueruela, I., & GonzálezFernández, C. 2018, MNRAS, 476, 3106 [Google Scholar]
 Tabernero, H. M., Dorda, R., Negueruela, I., & Marfil, E. 2021, A&A, 646, A98 [EDP Sciences] [Google Scholar]
 Tamayo, D., Silburt, A., Valencia, D., et al. 2016, ApJ, 832, L22 [Google Scholar]
 Tamayo, D., Cranmer, M., Hadden, S., et al. 2020, Proc. Nat. Acad. Sci., 117, 18194 [Google Scholar]
 Tamayo, D., Gilbertson, C., & ForemanMackey, D. 2021, MNRAS, 501, 4798 [Google Scholar]
 Teske, J., Xuesong Wang, S., Wolfgang, A., et al. 2020, ApJS, submitted [arXiv:2011.11560] [Google Scholar]
 Thiabaud, A., Marboeuf, U., Alibert, Y., et al. 2014, A&A, 562, A27 [Google Scholar]
 Tinetti, G., Drossart, P., Eccleston, P., et al. 2016, Proceedings of SPIE: Space Telescopes and Instrumentation 2016: Optical, Infrared, and Millimeter Wave, 9904, 99041X [Google Scholar]
 Turbet, M., Bolmont, E., Ehrenreich, D., et al. 2020, A&A, 638, A41 [Google Scholar]
 van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Comput. Sci. Eng., 13, 22 [Google Scholar]
 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
 Weiss, L. M., Marcy, G. W., Petigura, E. A., et al. 2018, AJ, 155, 48 [Google Scholar]
 Wolf, E. T. 2017, ApJ, 839, L1 [Google Scholar]
 Zacharias, N., Finch, C. T., Girard, T. M., et al. 2012, VizieR Online Data Catalog: I/322A [Google Scholar]
 Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577 [Google Scholar]
 Zechmeister, M., Dreizler, S., Ribas, I., et al. 2019, A&A, 627, A49 [NASA ADS] [EDP Sciences] [Google Scholar]
 Zeng, L., Sasselov, D. D., & Jacobsen, S. B. 2016, ApJ, 819, 127 [Google Scholar]
 Zhu, C., Byrd, R. H., Lu, P., & Nocedal, J. 1997, ACM Trans. Math. Softw., 23, 550 [Google Scholar]
HARPS stands for High Accuracy Radial velocity Planet Searcher and HARPSN for High Accuracy Radial velocity Planet Searcher North. CARMENES stands for Calar Alto highResolution search for M dwarfs with Exoearths with Nearinfrared and optical Échelle Spectrographs. ESPRESSO stands for Échelle SPectrograph for Rocky Exoplanets and Stable Spectroscopic Observations. TESS stands for Transiting Exoplanet Survey Satellite and TRAPPIST for TRAnsiting Planets and PlanetesImals Small Telescope.
A detailed description of the ESPRESSO DRS can be found in the ESPRESSO pipeline user manual available at espressopiperecipes
VOSA is publicly available online http://svo2.cab.intacsic.es/theory/vosa/
The modified version of batman is available at https://github.com/odemangeon/batman. It prevents the code to stay trapped in an infinite loop for highly eccentric orbits.
Several of the Python packages used for this work are publicly available on Github: radvel at https://github.com/CaliforniaPlanetSearch/radvel, george at https://github.com/dfm/george, batman at https://github.com/lkreidberg/batman, emcee at https://github.com/dfm/emcee, ldtk at https://github.com/hpparvi/ldtk.
bayev is available at https://github.com/exord/bayev
Confirmed planets with lower masses that can be found in exoplanet.eu and the NASA exoplanet archive were all measured through transittiming variations.
All Tables
ESPRESSO RV, FWHM, BIS, contrast, S_{index}, H_{α}, NaD, and BERV measurements for L 9859
All Figures
Fig. 1
Absolute magnitude (in the G Gaia bandpass) vs. color (magnitude difference between the Gaia bands G_{BP} and G_{RP}): L 9859 A (TOI175) is located in the Gaia colormagnitude diagram together with the mean sequences of young clusters and moving groups (Luhman 2018) and the main sequence of stars (Cifuentes et al. 2020). The error bars of L 9859 A are smaller than the symbol size. The gray area represents the 1 σ dispersion of field M dwarfs. 

In the text 
Fig. 2
GLSP of the RV and activity indicators from ESPRESSO (a) and HARPS (b) data. The last row for both instruments presents the window function. The vertical dotted lines indicate from right to left the orbital period of the planets b, c, d, and e, the planetary candidate 05, and half and the full stellar rotation period (assumed here to be 80 days). The horizontal lines indicate the amplitude levels corresponding to 10% (dashed line), 1% (dotdashed line) and 0.1% (dotted line) of the FAP. The amplitudes of the GLSP s are expressed using the ZK normalization described in Zechmeister & Kürster (2009, Eq. (5)). The FAP levels are computed using the analytical relation also described in Zechmeister & Kürster (2009) for this normalization. We display the GLSP of the BIS for completeness and comparison with C19, but we caution that the reliability of BIS measurements from CCF s for M dwarfs is uncertain (Rainer et al. 2020). 

In the text 
Fig. 3
GLSP of the TESS LC. The format of the this figure is identical to that in Fig. 2. In particular, the power of the GLSP is normalized using the ZK normalization. The highest peak in this periodogram is for a period of 93 days. 

In the text 
Fig. 4
Outcome of the fit of the four planets model: top left: RV time series along with the best model (solid green line), which includes the planetary signals and best prediction from the GP stellar activity model. The 1 σ uncertainties from the GP prediction are also displayed (shaded green area). For this plot, we subtracted the systemic velocity and the instruments offsets from the RV data (see values in Table H.1). Bottom left: Time series of the residuals of the best model. Right: Zoom on a small portion of the time series to better visualize the shorttimescale variations. 

In the text 
Fig. 5
Outcome of the fit of the fourplanet model: GLSP s of the RV time series (top) and of the planetary (second) and stellar activity (third) models sampled at the same time as the RV data, GLSP of the time series of the residuals (fourth) andthe window function (bottom). The vertical lines on the GLSP s correspond to the orbital periods of planets b, c, d, and e, half and the full rotation period (estimated at 80 days) from right to left. 

In the text 
Fig. 6
Phasefolded HARPS and ESPRESSO RV s, best model (top) and residuals (bottom) for the four planets. The HARPS data, presented in Sect. 2.1.1, are displayed with empty blue circles, and the ESPRESSO data, presented in Sect. 2.1.2, are displayed with orange circles. The filled orange circles are for the data taken before the fiber change of ESPRESSO. The empty orange circles are for data taken after the change. For clarity, the error bars of the HARPS and ESPRESSO data points are not displayed. For this plot, the stellar activity model has been subtracted from each data point. The points with error bars in red correspond to averages of the data within evenly spaced bins in orbital phase, whose size is 0.07 orbital period. The best model is shown with a green line. Before the subtraction of the stellar activity model the RMS of the RV data is 3.5, 3.4, and 3.2 m s^{−1} for HARPS, ESPRESSO_{pre}, and ESPRESSO_{post}, respectively. After the subtraction of the stellar activity model, it is 2.9, 2.5, and 2.3 m s^{−1} for HARPS, ESPRESSO_{pre}, and ESPRESSO_{post}, respectively. Finally, after subtraction of the planetary model, the RMS of the residuals is 1.8, 1.2, and 0.7 m s^{−1} for HARPS, ESPRESSO_{pre}, and ESPRESSO_{post}, respectively. 

In the text 
Fig. 7
Phasefolded TESS LC, best model (top) and residuals (bottom) for the three transiting planets. The data presented in Sect. 2.2 are displayed in black. For clarity, the error bars are not displayed. The points with error bars in red correspond to averages of the data within evenly spaced bins in orbital phase whose size corresponds to 5 min. The best model is shown with a black line. The standard deviation of the raw and binned residuals is indicated above each residual plot. 

In the text 
Fig. 8
Massradius diagram of the smallplanet population. Each point represents a confirmed exoplanet for which the mass and radius are measured with a relative precision better than 50%. These data have been extracted from exoplanet.eu (Schneider et al. 2011). The shape of the points indicates the technique used to measure the mass of the planet: circles for RV, and squares for transit timing variations. The color of the point reflects the intensity of the incident flux. The level of transparency of the error bars indicates the relative precision of the planetary bulk density. The better the precision, the more opaque the error bars. The three transiting planets in the L 9859 system are labeled and appear circled in black. The labeled blue stars indicate the Solar System planets. The colored dashed lines are the massradius models from Zeng et al. (2016). The gray region indicates the maximum collision stripping of the mantle. The shaded horizontal blue line represent the radius gap (Fulton et al. 2017). L 9859 b is in a sparsely populated region of the parameter space and currently the lowestmass planet with a mass measured through RV. Lower planetary masses have all been measured through transittiming variation, e.g., for Trappist1 h (Gillon et al. 2017) to the left of L 9859 b. This plot has been produced using the code available at https://github.com/odemangeon/massradius_diagram. 

In the text 
Fig. 9
Ternary diagrams showing the internal composition (mass fractions of the gas (H and He), the volatile (water) and the refractory elements) for the three transiting planets in system L 9859. 

In the text 
Fig. 10
Transmission spectrum metric vs. planetary radius diagram of the smallplanet population. Each point represents a confirmed exoplanets with mass and radius measured with a relative precision better than 50%. These data have been extracted from the exoplanet archive. The shape of the points indicates the technique used to measure the mass of the planet: circles for RV, and squares for transit timing variations. The color of the point reflects the equilibrium temperature of the planet. The level of transparency of the error bars indicates the relative precision of the planetary bulk density. The better the precision, the more opaque the error bars. The three transiting planets in system L 9859 are labeled and appear circled in black. We also display the names of the other planets with the highest transmission spectrum metrics. This plot has been produced using the code available at https://github.com/odemangeon/massradius_diagram. 

In the text 
Fig. C.1
posterior distributions of the main hyperparameters of the rotational kernel (Equation ((C.1))) 

In the text 
Fig. E.1
Phasefolded TESS LC assuming the best model ephemerides of planet e (a) and planetary candidate 5 (b). The black points are the TESS data points at the original cadence. The red line is the data binned in phase using bins of 15 min. The dashed pink and brown lines are the expected transit signal assuming that the planets have the same radius as planet d (see Table H.1). 

In the text 
Fig. E.2
Periodogram provided by the BLS search in the TESS data. The dashed vertical pink and brown links indicate the orbital period of planet e and planetary candidate, respectively. There is no significant power at these periods. 

In the text 
Fig. F.1
Outcome of the fit of the threeplanet model: The format of this figure is identical to that in Fig. 4, but isdescribed again here for convenience. (Top left) RV time series along with the best model (solid green line)that includes the planetary signals and best prediction from the GP stellar activity model. The 1 σ uncertainties from the GP prediction are also displayed (shaded green area). For this plot, we subtracted the systemic velocity and the instruments offsets from the RV data (see values in Table H.1). (Bottom left) Time series of the residuals of the best model. (Right) Zoom on a small portion of the time series to better visualize the shorttimescale variations. 

In the text 
Fig. F.2
Outcome of the fit of the threeplanet model: The format of this figure is identical to that in Fig. 5, but is described again here for convenience. GLSP s of the RV time series (top) and of the planetary (second) and stellar activity (third) models sampled at the same times as the RV data. GLSP of the time series of the residuals (fourth) and the window function (bottom). The vertical lines on the GLSP s correspond to the orbital periods of planets b, c, and d, half and thefull rotation period (estimated at 80 days) from right to left. 

In the text 
Fig. F.3
Outcome of the fit of the fourplanet model regarding the FWHM: The structure of this figure is similar to that of Fig. 4 or Fig. F.1, but the FWHM data and model are displayed instead of the RV data and model. 

In the text 
Fig. F.4
Outcome of the fit of the fourplanet model regarding the FWHM: The structure of this figure is similar to that of Fig. 5 or Fig. F.2, but the FWHM data and model are displayed instead of the RV data and model. 

In the text 
Fig. F.5
Iterative GLSP for the four planets model: GLSP of the RV data (top) and the window function (bottom). The GLSP s of the data shown in the previous row minus the model for planets c, e, and d, the stellar activity model, and the model for planet b are displayed in the second, third, fourth, fifth, and sixth row, respectively. 

In the text 
Fig. F.6
Outcome of the fit of the fiveplanet model: The format of this figure is identical to that in Figs. 4 and F.1. 

In the text 
Fig. F.7
Outcome of the fit of the fiveplanet model: The format of this figure is identical to that in Figs. 5 and F.2. 

In the text 
Fig. G.1
Corner plot showing the main internal structure parameters of L 9859 b. We show the mass fraction of the inner core, the mass fraction of water, the Si and Mg mole fraction in the mantle, the Fe mole fraction in the inner core, and the mass of gas (log scale). The values at the top of each column are the mean and 5% and 95% quantiles. 

In the text 
Fig. G.2
Same as Fig. G.1 for L 9859 c. 

In the text 
Fig. G.3
Same as Fig. G.1 for L 9859 d. 

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.