© ESO 2021

1 Introduction

In the pastyears, radial velocity (RV) instruments such as HARPS¹ (Mayor et al. 2003), HARPS-N¹ (Cosentino et al. 2012), and more recently, CARMENES¹ (Quirrenbach et al. 2014) and ESPRESSO¹ (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., Astudillo-Defru 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). Space-based transit surveys such as Kepler/K2 (Borucki et al. 2010; Howell et al. 2014) and TESS¹ (Ricker et al. 2015) and even ground-based surveys such as TRAPPIST¹ (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: Planet-planet 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 98-59 system, also known as the TESS Object of Interest 175 (TOI-175) 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; Bailer-Jones et al. 2018) M-dwarf 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⁻³ 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 follow-up RV campaign with ESPRESSO (Pepe et al. 2021) aimed at refining the mass of the planets in the L 98-59 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 98-59 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 High-resolution 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.C-0838, 1102.C-0339, and 0102.C-0525) between October 17, 2018 (barycentric Julian date, BJD = 2458408.5), and April 28, 2019 (BJD = 2458601.5). HARPS is a fiber-fed cross-dispersed echelle spectrograph operating in a temperature- and pressure-regulated 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 Lomb-Scargle 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 signal-to-noise 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 (Astudillo-Defru et al. 2017b). Their median precision (1 σ uncertainty) is 2.08 m s⁻¹.

In addition to the RV measurements, C19 provided the measurement of several stellar activity indicators: the full width at half maximum (FWHM) of the cross-correlation 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 S-index 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.C-0744, 1102.C-0958, and 1104.C-0350). ESPRESSO (Pepe et al. 2021) is also a fibre-fed high-resolution echelle spectrograph operating in a temperature- and pressure-regulated vacuum chamber. It covers wavelengths from 380 to 788 nm with an average spectral resolution of R = 140 000 in its single UT high-resolution mode (HR21, slow-readout 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 Data-Reduction Software (DRS)². It computes the CCF of the sky-subtracted 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 photo-noise-limited 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 S-index (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⁻¹ (a factor 2.6better than the HARPS RV s). At about the middle of our RV campaign, in June 2019, the fiber-link 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 High-precision photometry with TESS

L 98-59 (TIC 307 210 830, TOI-175) 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 pre-search data-conditioned 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 (Foreman-Mackey et al. 2017; Foreman-Mackey 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 $$S(\omega )=\sqrt{\frac{2}{\pi }}\frac{S_0{\omega }_0^4}{{({\omega }^2-{\omega }_0^2)}^2+{\omega }^2{\omega }_0^2/Q^2},$$(1)

where Q, the quality factor, was fixed to $\frac{1}{\sqrt{2}}$, S₀ is the amplitude, and ω₀ is the angular frequency corresponding to the break point in the power spectral density of this kernel.

The fit was performed using an affine-invariant ensemble sampler for MCMC (Goodman & Weare 2010) implemented in the Python package emcee (Foreman-Mackey 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₀, and we obtained a posterior providing an estimate of $S_0={82.38}_{-5.77}^{+6.59}$ ppm, using the median and the 68% confidence interval. For ω₀, we used a uniform prior between −20 and 15 on lnω₀, and we obtained an estimate of $\ln {\omega }_0={1.17}_{-0.09}^{+0.08}$ (in lnday⁻¹). 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 burn-in. The initial positions for this first run were drawn from the prior for S₀ and ω₀ 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 mono-modal, 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 mid-transit time. This reduced the number of data points and the computation time.

3 Characterization of the M dwarf L 98-59 A

According to K19, L 98-59 A is an M3V star. The derivation of accurate stellar properties through high-resolution spectroscopy for M stars is complicated because blended lines prevail. We thus used several approaches to characterize L 98-59 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 98-59 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 Gaia-DR2 parallax by Bailer-Jones et al. 2018) and the well-sampled photometric spectral energy distribution (SED, see Appendix A.1.3), we can derive a reliable estimate of the absolute bolometric luminosity of L 98-59: 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 98-59 A to be ${0.303}_{-0.023}^{+0.026}{R}_{\odot }$ using the Stefan-Boltzmann law. This agrees well (better than 1 σ) with the literature value derived by K19 from the mass-radius relations for M and K dwarfs of Boyajian et al. (2012). We derived the mass of L 98-59 A using the Virtual Observatory SED Analyzer online tools³ (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 BT-Settl evolutionary tracks (Allard et al. 2012). Finally, we determined the age of L 98-59 A using the photometry and distance provided by Gaia (see Appendix A.2 for more details). We compared the location of L 98-59 A in the color-magnitude diagram (see Fig. 1) to mean sequences of stellar members of the β Pictoris moving group (~20 Myr, Miret-Roig et al. 2020), the Tucana-Horologium 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 98-59 has an age consistent with that of the field. This age estimate is confirmed by our kinematics analysis, which indicates that L 98-59 A is a thin-disk 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 98-59 A are provided in Table H.1.

Fig. 1 Absolute magnitude (in the G Gaia bandpass) vs. color (magnitude difference between the Gaia bands GBP and GRP): L 98-59 A (TOI-175) is located in the Gaia color-magnitude 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 98-59 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 98-59 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 98-59 A is a member of the Galactic thin-disk population (see Table A.4), only stars belonging to the thin-disk population were selected. The selection of the thin-disk 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 thin-disk 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 (star-to-star 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 98-59 A is ${80.9}_{-5.3}^{+5.0}$ 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 S-index all show peaks at this period with a false-alarm 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 S-index, 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 79-day periodicity is a confirmation of the 80-day period identified in the GLSP s of the spectroscopic time series presented in Fig. 2. However, the 93- and 115-day 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 98-59. 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% (dot-dashed 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 light-curve modeling

4.1 Search for additional planets in the L 98-59 system

K19 and C19 confirmed the presence of three transiting planets in the L 98-59 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. 2-b) 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. 2-a) 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 high-resolution spectroscopy data. The main objective of this second step is to assess the presence of additional planets in system L 98-59 (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⁴ of the Python package batman⁵ (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 limb-darkening 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 affine-invariant ensemble sampler for MCMC implemented in the Python package emcee (Foreman-Mackey 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 pre-maximization 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 pre-minimization, 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⁻¹ for the stellar activity signal compared to ≲2 m s⁻¹ for the semi-amplitude 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 quasi-periodic kernel for the stellar activity. The mathematical expression of the kernel of this GP is $$K_{\text{RV}}(t_i,t_j)=A_{\text{RV}}{}^2\exp \left[-\frac{{(t_i-t_j)}^2}{2{\tau }_{\text{decay}}^2}-\frac{{\sin }^2\left(\frac{\pi }{P_{\text{rot}}}\text{|}{t}_i-t_j\text{|}\right)}{2{\gamma }^2}\right],$$(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⁵ (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 quasi-periodic 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 Lillo-Box 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 quasi-periodic 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₀) and one Keplerian function per planet in the system. The parameters of each Keplerian are the semi-amplitude (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⁵ (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 fiber-link change of ESPRESSO: HARPS, ESPRESSO before (pre), and ESPRESSO after the intervention (post). We used ESPRESSO_pre as RV reference, meaning that v₀ 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 pre-minimization 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 98-59 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 98-59 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 four-planet 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 semi-amplitude 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 ($\mathcal{Z}$) of our models using the Perrakis algorithm (Perrakis et al. 2014) using the Python implementation bayev⁶ (Díaz et al. 2014). We computed the logarithm of $\mathcal{Z}$ 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 four-planet 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 five-planet models. The BIC values of the five-planet model is significantly higher (ΔBIC = 3 for the RV+FWHM analysis), while the Bayesian evidence of these two models is very similar ($\Delta \ln \mathcal{Z}=0.4$).

We thus conclude that our additional ESPRESSO RV campaign allows us to identify one additional planet in the L 98-59 system: a fourth planet, hereafter planet L 98-59 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 98-59 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 long-term 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 short-timescale variations.

Fig. 5 Outcome of the fit of the four-planet 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.

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 RV-only 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 98-59 system

In compact multiplanetary systems such as L 98-59, the assumption of long-term stability of the system can bring strong constraints on the planetary masses and orbital properties. Both K19 and C19 performed N-body 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, N-body simulations (10⁴ 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 machine-learning algorithm that estimates the probability that the simulated system is stable on the long term (typically 10⁹ 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 N-body simulation (within 10⁴ 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⁵ versions of the L 98-59 systems described by the last 10⁵ 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 N-body 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 long-term stability estimated for the last 10⁵ 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 long-term stability, and they are reported in Table H.1. The phase-folded data (RV and photometry) and the best model are displayed in Figs. 6 and 7. The main impact of the long-term dynamical stability condition is on the eccentricity of planet c, which decreases from ${0.147}_{-0.048}^{+0.044}$ to ${0.103}_{-0.058}^{+0.045}$. 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 N-body simulations.

Finally, in order to assess whether planets c and d are in mean motion resonance, we performed an additional N-body simulation for each iteration of the system with a probability of long-term 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⁻⁴ year∕2π (which corresponds to ~ 10³ 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)) $${\theta }_i=2{\lambda }_d-{\lambda }_c-{\omega }_i,i\in [\text{c},\text{d}],$$

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 Phase-folded 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, ESPRESSOpre, and ESPRESSOpost, respectively. After the subtraction of the stellar activity model, it is 2.9, 2.5, and 2.3 m s−1 for HARPS, ESPRESSOpre, and ESPRESSOpost, 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, ESPRESSOpre, and ESPRESSOpost, 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 ${88.21}_{-0.27}^{+0.35}$ degrees and predict the impact parameter of planet e (${1.47}_{-0.30}^{+0.27}$) and candidate 5 (${2.27}_{-0.43}^{+0.46}$). 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 low-mass 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 quasi-periodic 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 98-59, 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 Phase-folded 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 Four-planet 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 semi-amplitude of ${0.46}_{-0.17}^{+0.20},\text{m}{\text{s}}^{-1}$ and a mass of ${0.40}_{-0.15}^{+0.16}$ M_⊕ (half the mass of Venus), L 98-59 b currently is the lowest-mass exoplanet measured through RV⁷. 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⁻¹ in multiplanetary systems even when there is stellar activity. The relative precision of the RV semi-amplitude 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 super-Earths around M dwarfs that can currently be achieved (Suárez Mascareño et al. 2020; Lillo-Box 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 mass-radius 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 ${3.06}_{-0.37}^{+0.33}$ and ${2.46}_{-0.82}^{+0.66}$ 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 ${285}_{-17}^{+18}$ K, the planetary candidate 5, if confirmed, would orbit in the habitable zone of its parent star.

Fig. 8 Mass-radius diagram of the small-planet 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 98-59 system are labeled and appear circled in black. The labeled blue stars indicate the Solar System planets. The colored dashed lines are the mass-radius 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 98-59 b is in a sparsely populated region of the parameter space and currently the lowest-mass planet with a mass measured through RV. Lower planetary masses have all been measured through transit-timing variation, e.g., for Trappist-1 h (Gillon et al. 2017) to the left of L 98-59 b. This plot has been produced using the code available at https://github.com/odemangeon/mass-radius_diagram.

5.3 Internal composition of the three transiting super-Earths

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 98-59. 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 98-59 d according to the Bayesian analysis are substantially different: the mode of the water-mass fraction distribution is at ~ 0.3, whereas the mode of the gas mass peaks at ~10⁻⁶ 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 98-59 d than in L 98-59 b and L 98-59 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 helium-free 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 water-dominated 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 water-dominated 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 TRAPPIST-1 system (Gillon et al. 2017; Luger et al. 2017) for comparison. Turbet et al. (2020) stressed the impact of irradiation on a water-dominated 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 TRAPPIST-1 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 TRAPPIST-1 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 TRAPPIST-1 b and c could still harbor significant amounts of water.

L 98-59 b is similar in mass and radius to TRAPPIST-1 d. However, it is significantly more irradiated (T_eq = 288 ± 5.6 K for TRAPPIST-1 d Gillon et al. 2017). L 98-59 b might thus undergo or have undergone efficient hydrodynamic escape. L 98-59 c and d are more massive than any of the TRAPPIST-1 planets, but they are also more irradiated (T_eq = 400.1 ± 7.7 K for TRAPPIST-1 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 98-59 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 98-59.

6 Conclusion: L 98-59, a benchmark system for super-Earth 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⁸ (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 so-called “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 98-59 system along with the distributions of these metrics derived by the authors from their sample. From Table 2, we conclude that system L 98-59 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 98-59 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 98-59 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⁻¹, the semi-amplitude of the expected Rossiter–McLaughling effect (e.g., Queloz et al. 2000) is 40 cm s⁻¹, 1 m s⁻¹, and 37 m s⁻¹ for planets b, c,and d, respectively. These amplitudes might at least for planets c and d be within the reach of high-resolution spectrographs such as ESPRESSO. This would give us access to the spin-orbital angle in this system, which would further constrain itsarchitecture and the possible mechanisms of its formation and migration.

The fact that L 98-59 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 98-59) around M dwarfs for which the planetary masses and radius of at least two planets have been measured. The other six are TRAPPIST-1 (Gillon et al. 2017), LTT-3780 (Cloutier et al. 2020), TOI-1266 (Demory et al. 2020), LHS-1140 (Lillo-Box et al. 2020), K2-146 (Hamann et al. 2019), and Kepler-138 (Jontof-Hutter et al. 2015). With a V magnitude of 11.7 and a distance of 10.6 pc, L 98-59 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 98-59 are comfortably above the thresholds proposed by Kempton et al. (2018) for a super-Earth 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 well-characterized small-planet population. L 98-59 b and c are the two planets with the highest TSM value below 1.5 R_⊕, and L 98-59 d has the second highest above this value. These three planets thus belong to the most favorable warm to temperate (T_eq < 650 K) super-Earths (R_p < 1.5 R_⊕) for an atmospheric characterization. Furthermore, L 98-59 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).

Fig. 10 Transmission spectrum metric vs. planetary radius diagram of the small-planet 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 98-59 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/mass-radius_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 AyA2017-84089. 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 (EXOPLANETS-A). 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/FIS-AST/32113/2017 and POCI-01-0145-FEDER-032113; PTDC/FIS-AST/28953/2017 and POCI-01-0145-FEDER-028953; PTDC/FIS-AST/28987/2017 and POCI-01-0145-FEDER-028987; PTDC/FIS-AST/30389/2017 and POCI-01-0145-FEDER-030389. 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 RYC-2013-14 875. J.I.G.H., A.S.M., R.R., and C.A.P. acknowledge financial support from the Spanish MICINN AYA2017-86 389-P. 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/FIS-PAR/0037/2019, PTDC/FIS-OUT/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 2018-24-HH.0. The financial contribution from the agreement ASI-INAF n.2018-16-HH.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 98-59 A

A.1 Atmospheric parameters of L 98-59 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 98-59 A.

A.1.1 Spectral synthesis with STEPARSYN

We employed the BT-Settl model grid (Allard et al. 2012), the radiative transfer code turbospectrum (Plez 2012), and a VALD3-based 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 98-59 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 (Foreman-Mackey et al. 2013), a Markov chain Monte Carlo (MCMC) method used to fully sample the underlying distribution of the stellar parameters of L 98-59. 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): $\sqrt{{\zeta }^2+{(v\sin i)}^2}=3.78\pm 0.44$ km s⁻¹.

A.1.2 Machine-learning regression with odusseas

The ODUSSEAS software (Antoniadis-Karnavas et al. 2020) receives a 1D spectrum and its resolution as input. The pseudo-equivalent widths are measured and used as input for a supervised machine-learning algorithm (ridge regression model) that is used to derive the spectroscopic parameters T_eff and [Fe∕H]. The implementation relies on the machine-learning 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 machine-learning models at this resolution, after taking the intrinsic uncertainties of the reference sample parameters during the machine-learning 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 BT-Settl model (Allard et al. 2012) for its treatment of dust and clouds, which is important for low-mass stars. Because of the small distance of 10.6194 pc (inferred from GAIA parallaxes, Bailer-Jones 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 alpha-element enhancement ([α/Fe]) of $-{0.03}_{-0.13}^{+0.16}$ dex and a luminosity of L = 0.01128 ± 0.00042 L_⊙.

A.1.4 K19 approach

K19 estimated T_eff and log g from two mostly independent derivations. T_eff was derived using the Stefan-Boltzman law. The required bolometric luminosity was estimated from V- and K-band 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 mass-luminosity relation of Benedict et al. (2016) and the mass-radius 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 98-59 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 machine-learning estimates, which rely on the high-resolution 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 machine-learning uncertainties appear to be underestimated. For [Fe∕H], the four estimates are compatible within 1.6 σ. As expected, the spectral synthesis and machine-learning 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.

In this paper, which focuses on the characterization of the planets in system L 98-59, we need to conclude with one final set of T_eff, [Fe∕H], and log g estimates. In order to keep a physically self-consistent 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 high-resolution 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 machine-learning 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 98-59 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 ${0.241}_{-0.069}^{+0.097}{\text{M}}_{\odot }$. 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 ${0.311}_{-0.081}^{+0.10}{\text{M}}_{\odot }$. Our third method relies on the mass-luminosity 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 (Bailer-Jones 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 low-mass 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 (BT-Settl 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 mass-luminosity relation for M dwarfs of Benedict et al. (2016). They derived the luminosity from K-band observations.

Table A.3 gathers all these estimates of the radius and mass of L 98-59 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 98-59 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.

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 Stefan-Boltzmann 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 98-59 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 98-59 A, we used the accurate photometry and distance provided by Gaia. We constructed the color-magnitude diagram shown in Fig. 1, where we also depict the well-known empirically determined mean sequences of stellar members of the β Pictoris moving group (~20 Myr, Miret-Roig et al. 2020), the Tucana-Horologium 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 98-59 A is feasible without any systematic effect. From its location in the Gaia color-magnitude diagram, we infer that L 98-59 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 98-59 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 98-59 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 98-59 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 98-59 A (see Table A.4). L 98-59 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.

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.

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 quasi-periodic signal. Its equation is $$k(\tau )=\frac{B}{2+C}{e}^{-\tau /L}\left[\cos \left(\frac{2\pi \tau }{P_{\text{rot}}}\right)+(1+C)\right],$$(C.1)

and it is taken from Foreman-Mackey 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 L-BFGS-B 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 (${190}_{-134}^{+189}$ days). It is also worth noting that the retrieved amplitude and timescales are low and short: ${0.11}_{-0.01}^{+0.05}$ ppm for the amplitude, and ${6.3}_{-1.0}^{+2.6}$ days for the timescales. In particular, $\frac{L}{P_{rot}}$ 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 model-independent approach to determining the rotational modulation in the TESS LC. TESS is designed for high-precision relative photometry (as opposed to high-precision 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.

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⁻³ 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 limb-darkening coefficients, we used Gaussian PDF s whose first two moments were defined using the Python package ldtk¹¹ (Parviainen & Aigrain 2015). Using a library of synthetic stellar spectra, it computes the limb-darkening 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 limb-darkening 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₀) 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 quasi-periodic kernels. The prior PDF s of the two amplitudes (A_RV, A_FWHM) are uniform between zero and the maximum of the peak-to-peak 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 98-59, 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 quasi-periodic 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 peak-to-peak 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 transit-timing variations. We then performed a box least-squares (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. Phase-folding 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 phase-folded 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 Phase-folded 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 98-59 system

As mentioned in Sect. 4.1.3, in order to assess the presence of additional planets in the L 98-59 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 semi-amplitude 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 23-day signal.

Consequently, we performed other analyses with five planets. The fits converge toward a significant detection of the semi-amplitude 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 four-planet analyses (Fig. 5).

Fig. F.1 Outcome of the fit of the three-planet 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 short-timescale variations.

Fig. F.2 Outcome of the fit of the three-planet 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 four-planet 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 four-planet 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 five-planet 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 five-planet model: The format of this figure is identical to that in Figs. 5 and F.2.

Appendix G Internal composition of three transiting super-earths

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 semi-amplitude 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 semi-amplitude, 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 98-59 are shown in Figs. G.1 to G.3.

Fig. G.1 Corner plot showing the main internal structure parameters of L 98-59 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.

Fig. G.2 Same as Fig. G.1 for L 98-59 c.

Fig. G.3 Same as Fig. G.1 for L 98-59 d.

Appendix H Additional table

References

All Tables

All Figures

Fig. 1 Absolute magnitude (in the G Gaia bandpass) vs. color (magnitude difference between the Gaia bands GBP and GRP): L 98-59 A (TOI-175) is located in the Gaia color-magnitude 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 98-59 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% (dot-dashed 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 short-timescale variations.

In the text

Fig. 5 Outcome of the fit of the four-planet 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 Phase-folded 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, ESPRESSOpre, and ESPRESSOpost, respectively. After the subtraction of the stellar activity model, it is 2.9, 2.5, and 2.3 m s−1 for HARPS, ESPRESSOpre, and ESPRESSOpost, 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, ESPRESSOpre, and ESPRESSOpost, respectively.

In the text

Fig. 7 Phase-folded 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 Mass-radius diagram of the small-planet 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 98-59 system are labeled and appear circled in black. The labeled blue stars indicate the Solar System planets. The colored dashed lines are the mass-radius 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 98-59 b is in a sparsely populated region of the parameter space and currently the lowest-mass planet with a mass measured through RV. Lower planetary masses have all been measured through transit-timing variation, e.g., for Trappist-1 h (Gillon et al. 2017) to the left of L 98-59 b. This plot has been produced using the code available at https://github.com/odemangeon/mass-radius_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 98-59.
In the text

Fig. 10 Transmission spectrum metric vs. planetary radius diagram of the small-planet 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 98-59 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/mass-radius_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 Phase-folded 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 three-planet 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 short-timescale variations.

In the text

Fig. F.2 Outcome of the fit of the three-planet 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 four-planet 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 four-planet 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 five-planet 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 five-planet 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 98-59 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 98-59 c.
In the text

	Fig. G.3 Same as Fig. G.1 for L 98-59 d.
In the text