Six transiting planets and a chain of Laplace resonances in TOI-178

Determining the architecture of multi-planetary systems is one of the cornerstones of understanding planet formation and evolution. Resonant systems are especially important as the fragility of their orbital configuration ensures that no significant scattering or collisional event has taken place since the earliest formation phase when the parent protoplanetary disc was still present. In this context, TOI-178 has been the subject of particular attention since the first TESS observations hinted at a 2:3:3 resonant chain. Here we report the results of observations from CHEOPS, ESPRESSO, NGTS, and SPECULOOS with the aim of deciphering the peculiar orbital architecture of the system. We show that TOI-178 harbours at least six planets in the super-Earth to mini-Neptune regimes, with radii ranging from 1.152(-0.070/+0.073) to 2.87(-0.13/+0.14) Earth radii and periods of 1.91, 3.24, 6.56, 9.96, 15.23, and 20.71 days. All planets but the innermost one form a 2:4:6:9:12 chain of Laplace resonances, and the planetary densities show important variations from planet to planet, jumping from 1.02(+0.28/-0.23) to 0.177(+0.055/-0.061) times the Earth's density between planets c and d. Using Bayesian interior structure retrieval models, we show that the amount of gas in the planets does not vary in a monotonous way, contrary to what one would expect from simple formation and evolution models and unlike other known systems in a chain of Laplace resonances. The brightness of TOI-178 allows for a precise characterisation of its orbital architecture as well as of the physical nature of the six presently known transiting planets it harbours. The peculiar orbital configuration and the diversity in average density among the planets in the system will enable the study of interior planetary structures and atmospheric evolution, providing important clues on the formation of super-Earths and mini-Neptunes.


Introduction
Since the discovery of the first exoplanet orbiting a Sun-like star by Mayor & Queloz (1995), the diversity of observed planetary systems has continued to challenge our understanding of their formation and evolution. As an ongoing effort to understand these physical processes, observational facilities strive to get a full picture of exoplanetary systems by looking for addi-CHEOPS Fellow Winton Fellow Juan Carlos Torres Fellow tional candidates to known systems and by better constraining the orbital architecture, radii, and masses of the known planets.
In particular, chains of planets in mean-motion resonances (MMRs) are 'Rosetta Stones' of the formation and evolution of planetary systems. Indeed, our current understanding of planetary system formation theory implies that such configurations are a common outcome of protoplanetary discs: Slow convergent migration of a pair of planets in quasi-circular orbits leads to a high probability of capture in first-order MMRs -the period ratio of the two planets is equal to (k + 1)/k, with k an integer (Lee & Peale 2002;Correia et al. 2018). As the disc strongly damps the eccentricities of the protoplanets, this mechanism can For planets that remained in resonance and are close enough to their host stars, tides become the dominant force that affects the architecture of the systems, which can then lead to a departure of the period ratio from the exact resonance (Henrard & Lemaitre 1983;Papaloizou & Terquem 2010;Delisle et al. 2012). In some near-resonant systems, such as HD 158259, the tides seem to have pulled the configuration entirely out of resonance . However, through gentle tidal evolution it is possible to retain a resonant state even with null eccentricities through three-body resonances (Morbidelli 2002;Papaloizou 2015). Such systems are too fine-tuned to result from scattering events and hence can be used to constrain the outcome of protoplanetary discs (Mills et al. 2016).
A Laplace resonance, in reference to the configuration of Io, Europa, and Ganymede, is a three-body resonance where each consecutive pair of bodies is in, or close to, two planet MMRs. To date, only a few systems have been observed in a chain of Laplace resonances: GJ 876 (Rivera et al. 2010), Kepler-60 (Goździewski et al. 2016), Kepler-80 (MacDonald et al. 2016), Kepler-223 (Mills et al. 2016), Trappist-1 (Gillon et al. 2017;Luger et al. 2017), and K2-138 (Christiansen et al. 2018;Lopez et al. 2019). Out of these six systems, only K2-138 has so far been observed by both transit and radial velocity (RV), mainly due to the relative faintness of the other host stars in the visible (V magnitudes greater than ∼14). However, TTVs could also be used to estimate the mass of their planets (see, for example, Agol et al. (2020) for the case of Trappist-1).
In this study, we present photometric and RV observations of TOI-178, a V = 11.95 mag, K-type star that was first observed by TESS 1 in its Sector 2. We jointly analyse the photometric data of TESS, two nights of NGTS 2 and SPECULOOS 3 data, and 285 hours of CHEOPS 4 observations, along with 46 ESPRESSO 5 RV points. This follow-up effort allows us to decipher the architecture of the system and demonstrate the presence of a chain of Laplace resonances between the five outer planets.
We begin in Sect. 2 by presenting the rationale that led to the CHEOPS observation sequence (two visits totalling 11 d followed by two shorter visits). In Sect. 3, we describe the parameters of the star. In Sect. 4, we present the photometric and RV data we use in the paper. In Sect. 5, we show how these data led us to the detection of six planets -the outer five of which are in a chain of Laplace resonances -and to constrain their parameters. In Sect. 6, we explain the resonant state of the system, discuss its stability, and describe the transit timing variations (TTVs) that this system could potentially exhibit in the coming years. Finally, we discuss the internal structure of the planets in Sect. 7, and conclusions are presented in Sect. 8.

CHEOPS observation strategy
The CHEOPS observation consisted of one long double visit (11 days) followed by two short visits (a few hours each) at precise dates. We explain in this section how we came up with this particular observation strategy. Details on all the data used and acquired, as well as their analysis, are presented in the following sections.
The first release of candidates from the TESS alerts of Sector 2 included three planet candidates in TOI-178 with periods of 6.55 d, 10.35 d, and 9.96 d. Based on these data, TOI-178 was identified as a potential co-orbital system (Leleu et al. 2019) with two planets oscillating around the same period. This prompted ESPRESSO RV measurements and two sequences of simultaneous ground-based photometric observations with NGTS and SPECULOOS. From the latter, no transit was observed for TOI-178.02 (P = 10.35 d) in September 2019; however, a transit ascribed to this candidate was detected one month later by NGTS and SPECULOOS. The abovementioned absence of a TOI-178.02 transit combined with the three transits observed by TESS at high S/N (above 10) was interpreted as an additional sign of the strong TTVs expected in a co-orbital configuration. This solution was supported by RV data that were consistent with the horseshoe orbits of objects with similar masses (Leleu et al. 2015).
A continuous 11 d CHEOPS observation (split into two visits for scheduling reasons) was therefore performed in August 2020 in order to confirm the orbital configuration of the system; as the instantaneous period of both members of the co-orbital pair will always be smaller than 11 d, at least one transit of both targets should thus be detectable. Analysis of this light curve led to the confirmation of the presence of two of the planets already discovered by TESS (in this study denoted as planets d and e, with periods of 6.55 d and 9.96 d, respectively) and the detection of two new inner transiting planets (denoted planets b and c, with periods of 1.9 d and 3.2 d, respectively). However, one of the planets belonging to the proposed co-orbital pair (with a period of 10.35 d) was not apparent in the light curve. A potential hypothesis was that the first and third transit of the TOI-178.02 candidate (P = 10.35 d) during TESS Sector 2 belonged to a planet of twice the period (20.7 d), while the second transit belonged to another planet of unknown period. This scenario was supported by the two ground-based observations mentioned above since a P = 20.7 d planet would not transit during the September 2019 observation window. Using the ephemerides from fitting the TESS, NGTS, and SPECULOOS data, we predicted the mid-transit of the 20.7-day candidate to occur between UTC 14:06 and 14:23 on September 7, 2020 with 98% certainty. A third visit of CHEOPS observed the system around this epoch for 13.36 h and confirmed the presence of this new planet (g) at the predicted time with consistent transit parameters. Further analysis of all available photometric data found two possibilities for the unknown period of the additional planet: ∼ 12.9 d or ∼ 15.24 d.
Careful analysis of the whole system additionally revealed that planets c, d, e, and g were in a Laplace resonance. In order to fit in the resonant chain, the unknown period of the additional planet could have only two values: P = 13.4527 d or , the latter value being more consistent with the RV data. A fourth CHEOPS visit was therefore scheduled for 3 October, and it detected the transit for the additional planet at a period of 15.231915 +0.000115 −0.000095 day. As we will detail in the following sections, we can confirm the detection of a 6.56-day and a 9.96-day period planet by TESS using the new observations presented in this paper. Furthermore, we can announce the detection of planets with periods of 1. 91, 3.24, 15.23, and 20.71 d (see Tables 3 and 4 for the complete parameters of the system).

TOI-178 stellar characterisation
Forty-six ESPRESSO observations (see Sect. 4.2) of TOI-178, a V = 11.95 mag K-dwarf, have been used to determine the stellar spectral parameters. These observations were first shifted and stacked to produce a combined spectrum. We then used the publicly available spectral analysis package SME (Spectroscopy Made Easy; Valenti & Piskunov 1996;Piskunov & Valenti 2017) version 5.22 to model the co-added ESPRESSO spectrum. We selected the ATLAS12 model atmosphere grids (Kurucz 2013) and atomic line data from VALD to compute synthetic spectra, which were fitted to the observations using a χ 2 -minimising procedure. We modelled different spectral lines to obtain different photospheric parameters starting with the line wings of Hα, which are particularly sensitive to the stellar effective temperature T eff . We then proceeded with the metal abundances and the projected rotational velocity v sin i , which were modelled with narrow lines between 5900 and 6500 Å. We found similar values for [Fe/H], [Ca/H], and [Na/H]. The macro-turbulent velocity was modelled and found to be 1.2±0.9 km s −1 , and microturbulent velocity was fixed to 0.91 km s −1 following the formulation in Bruntt et al. 2010. The surface gravity log g was constrained from the line wings of the Ca i triplet (6102, 6122, and 6162 Å) and the Ca i 6439 Å line with a fixed T eff and Ca abundance.
We checked our model with the Na i doublet that is sensitive to both T eff and log g, and, finally, we also tested the MARCS 2012 (Gustafsson et al. 2008) model atmosphere grids. The measured parameters are listed in Table 1. The SME results are in agreement with the empirical SpecMatch-Emp (Yee et al. 2017) code, which fits the stellar optical spectra to a spectral library of 404 M5 to F1 stars, resulting in T eff = 4316 ± 70 K, log g = 4.45 ± 0.15, and [Fe/H] = −0.29 ± 0.05 dex. We also used ARES+MOOG (Sousa 2014;Sousa et al. 2015;Sneden 1973) to conduct the spectroscopic analysis on the same combined ESPRESSO spectra, and, although we derived consistent parameters (T eff = 4500 ± 230 K, log g = 4.38 ± 0.62, and [Fe/H] = −0.34 ± 0.10), the large errors are indicative of the difficulties in using equivalent width methods with colder stars: Spectral lines are more crowded in the spectrum.
Using these precise spectral parameters as priors on stellar atmospheric model selection, we determined the radius of TOI-178 using the infrared flux method (IRFM; Blackwell & Shal- Wright et al. (2010) lis 1977) in a Markov chain Monte Carlo (MCMC) approach. The IRFM computes the stellar angular diameter and effective temperature by comparing observed broadband optical and infrared fluxes and synthetic photometry obtained from convolution of the considered filter throughputs, using the known zeropoint magnitudes, with the stellar atmospheric model, with the stellar radius then calculated using the parallax of the star. For this study, we retrieved the Gaia G, G BP , and G RP , 2MASS J, H, and K, and WISE W1 and W2 fluxes and relative uncertainties from the most recent data releases (Skrutskie et al. 2006;Wright et al. 2010;Gaia Collaboration et al. 2018, respectively), and utilised the stellar atmospheric models from the atlas Catalogues (Castelli & Kurucz 2003), to obtain R = 0.651 ± 0.011 R , and T eff = 4352 ± 52 K, in agreement with the spectroscopic T eff used as a prior. We inferred the mass and age of TOI-178 using stellar evolutionary models, using as inputs T e f f , R and Fe/H, with evolutionary tracks and isochrones generated by two grids of models separately, the PARSEC 6 v1.2S code (Marigo et al. 2017) and the CLES code (Code Liégeois d'Évolution Stellaire; Scuflaire et al. 2008), with the reported values representing a careful combination of results from both sets of models. This was done because the sets of models differ slightly in their approaches (reaction rates, opacity and overshoot treatment, and helium-to-metal enrichment ratio), and thus by comparing masses and ages derived from both grids it is possible to include systematic uncertainties within the modelling of the position of TOI-178 on evolutionary tracks and isochrones. A detailed discussion of combining the PARSEC and CLES models to determine masses and ages can be found in Bonfanti et al. (2021). For this study, we derived M = 0.647 +0.035 −0.032 M and t = 7.1 +6.2 −5.4 Gyr. All stellar parameters are shown in Table 1.

Photometric data
In order to determine the orbital configuration of the TOI-178 planetary system, we obtained photometric time series observations from multiple telescopes, as detailed below.

TESS
Listed as TIC 251848941 in the TESS Input Catalog (TIC; Stassun et al. 2018Stassun et al. , 2019, TOI-178 was observed by TESS in Sector 2, by camera 2, from August 22, 2018 to September 20, 2018. The individual frames were processed into 2-minute cadence observations and reduced by the Science Processing Operations Center (SPOC; Jenkins et al. 2016) into light curves made publicly available at the Mikulski Archive for Space Telescopes (MAST). For our analysis, we retrieved the Presearch Data Conditioning Single Aperture Photometry (PDCSAP) light curve data, using the default quality bitmask, which have undergone known systematic correction (Smith et al. 2012;Stumpe et al. 2014). Lastly, we rejected data points flagged as of bad quality by SPOC (QUALITY > 0) and those with 'Not-a-Number' flux or flux error values. After these quality cuts, the TESS light curve of TOI-178 contained 18,316 data points spanning 25.95 d. The full dataset with the transits of the six identified planets is shown in Fig. 1.

CHEOPS
CHEOPS, the first ESA small-class mission, is dedicated to the observation of bright stars (V 12 mag) that are known to host planets and performs ultra high precision photometry, with the precision limited by stellar photon noise of 150 ppm/min for a V = 9 magnitude star. The CHEOPS instrument is composed of an f/8 Ritchey-Chretien on-axis telescope (∼30 cm diameter) equipped with a single frame-transfer back-side illuminated charge-coupled device (CCD) detector. The satellite was successfully launched from Kourou (French Guiana) into a ∼700 km altitude Sun-synchronous orbit on December 18, 2019. CHEOPS took its first image on February 7, 2020, and, after it passed the in-orbit commissioning (IOC) phase, routine operations started on March 25, 2020. More details on the mission can be found in Benz et al. (2020), and the first results have recently been presented in Lendl et al. (2020).
The versatility of the CHEOPS mission allows for spacebased follow-up photometry of planetary systems identified by the TESS mission. This is particularly useful for completing the inventory of multi-planetary systems whose outer transiting planets have periods beyond ∼10 days.
We obtained four observation runs (or visits) of TOI-178 with CHEOPS between 4 August 2020 and 3 October 2020 as part of the Guaranteed Time, totalling 11.88 days on target with the observing log shown in Table 2. The majority of this time was spent during the first two visits (with lengths of 99.78 h and 164.06 h), which were conducted sequentially beginning on August 4, 2020 and ending on August 15, 2020; as such, we achieved a near-continuous 11-day photometric time series. The runs were split due to scheduling constraints, with a 0.84 h gap between visits. The third and fourth visits were conducted to confirm the additional planets predicted in the scenario presented in Sect. 2. They took place on September 7, 2020 and October 3, 2020 and lasted for 13.36 h and 8.00 h, respectively.
Due to the low-Earth orbit of CHEOPS, the spacecraft-target line of sight was interrupted by Earth occultations and passages through the South Atlantic Anomaly (SAA), where no data were downlinked. This resulted in gaps in the photometry on CHEOPS orbit timescales (around 100 min). For our observations of TOI-178, this resulted in light curve efficiencies of 51%, 54%, 65%, and 86%. For all four visits, we used an exposure time of 60 s.
These data were automatically processed with the latest version of the CHEOPS data reduction pipeline (DRP v12;Hoyer et al. 2020). This includes image calibration (bias, gain, nonlinearity, dark current, and flat field) and instrumental and environmental corrections (cosmic rays and smearing trails of field stars and background). The DRP also performs aperture photometry for a set of three size-fixed apertures -R=22.5 (RINF), 25 (DEFAULT), and 30 (RSUP) -plus one extra aperture that minimises the contamination from nearby field stars (ROPT), which, in the case of TOI-178, was set at R=28.5 . The DRP estimates the level of contamination by simulating the field-ofview of TOI-178 using the GAIA star catalogue (Gaia Collaboration et al. 2018) to determinate the location and flux of the stars throughout the duration of the visit. In the case of TOI-178, the mean contamination level was below 0.1% and was mostly modulated by the rotation around the target of a nearby star of Gaia G=13.3 mag at a projected sky distance of 60.8 from the target. The contamination as a function of time (or, equivalently, as a function of the roll angle of the satellite) is provided as a product of the DRP for further de-trending (see Sect.5.1). For the second visit, careful removal of one 'telegraphic' pixel (a pixel with a  Table 2). The TOI-178 PSF is shown at the centre, with the DEFAULT DRP photometric aperture represented by the dashed black circles. The telegraphic pixel location that appeared close to the end of the observation is marked by the red circle. non-stable abnormal behaviour during the visit) within the photometric aperture was needed. The location of this telegraphic pixel is shown in red in Fig.3, and details regarding the detection and correction are described in Appendix A. Following the reductions, we found that the light curves obtained using the DE-FAULT aperture (R=25 ) yielded the lowest RMS for all visits and, as such, are used for this study (see Appendix A, Fig. A.1).
Lastly, it became apparent that -due to the nature of the CHEOPS orbit and the rotation of the CHEOPS field around the target -photometric, non-astrophysical short-term trends -either from a varying background, a nearby contaminating source, or other sources -can be found in the data on orbital timescales. These effects can be efficiently modelled by using a Gaussian process (GP) with roll angle as input (e.g. Lendl et al. 2020;Bonfanti et al. 2021), as we will discuss in our data analysis (Sect. 5.1.2). Following this, the average noise over a 3 h sliding window for the four visits of the G = 11.15 mag target was found to be 63.9, 64.2, 66.3, and 75.8 ppm. In all cases, this marginally improved upon the precision of the light curves that we had previously simulated for these observation windows using the CHEOPSim tool (Futyan et al. 2020).

NGTS
The NGTS (Wheatley et al. 2018) facility consists of twelve 20-cm diameter robotic telescopes and is situated at the ESO Paranal Observatory in Chile. The individual NGTS telescopes have a wide field-of-view of 8 square-degrees and a plate scale of 5 pixel −1 . The DONUTS auto-guiding algorithm (McCormac et al. 2013) affords the NGTS telescopes sub-pixel level guiding. Simultaneous observations using multiple NGTS telescopes have been shown to yield ultra high precision light curves of exoplanet transits (Bryant et al. 2020;Smith et al. 2020).
TOI-178 was observed using NGTS multi-telescope observations on two nights. On UT September 10, 2019, TOI-178 was observed using six NGTS telescopes during the predicted transit event of TESS candidate TOI-178.02. However, the NGTS data for this night rule out a transit occurring during the observations. A second predicted transit event of TOI-178.02 was observed on the night of UT October 11, 2019 using seven NGTS telescopes, and on this night the transit event was robustly detected by NGTS. A total of 13,991 images were obtained on the first night, and 12,854 were obtained on the second. For both nights, the images were taken using the custom NGTS filter (520 -890 nm) with an exposure time of 10 s. All NGTS observations of TOI-178 were performed with airmass < 2 and under photometric sky conditions.
The NGTS images were reduced using a custom pipeline that uses the SEP library to perform source extraction and aperture photometry (Bertin & Arnouts 1996;Barbary 2016). A selection of comparison stars with brightnesses, colours, and CCD positions similar to those of TOI-178 were identified using the second GAIA data release (DR2) (Gaia Collaboration et al. 2018). More details on the photometry pipeline are provided in Bryant et al. (2020).
The NGTS light curves are presented in Fig. 2. They show transit events for planet b and planet g on the nights of September 11, 2019 and October 1, 2019, respectively.

SPECULOOS
The SPECULOOS Southern Observatory (SSO; Gillon 2018;Burdanov et al. 2018;Delrez et al. 2018) is located at ESO's Paranal Observatory in Chile and is part of the SPECULOOS project. The facility consists of a network of four robotic 1-m telescopes (Callisto, Europa, Ganymede, and Io). Each robotic SSO telescope has a primary aperture of 1 m and a focal length of 8 m, and is equipped with a 2k×2k deep-depletion CCD camera whose 13.5 µm pixel size corresponds to 0.35" on the sky (fieldof-view = 12 x12 ). Observations were performed on the nights of October 10, 2019 (for 8 hours) and November 11, 2019 (for 8 hours) with three SPECULOOS telescopes on sky simultaneously (SSO/Io, SSO/Europa, and SSO/Ganymede). These observations were carried out in a Sloan i' filter with exposure times of 10 s. A small de-focus was applied to avoid saturation as the target was too bright for SSO. Light curves were extracted using the SSO pipeline (Murray et al. 2020) and are shown in purple in Fig. 2. For each observing night, the SSO pipeline used the casutools software (Irwin et al. 2004) to perform automated differential photometry and to correct for systematics caused by time-varying telluric water vapour.

ESPRESSO data
The RV data we analyse consist of 46 ESPRESSO points 7 . Each measurement was taken in high resolution (HR) mode with an integration time of 20 min using a single telescope (UT) and slow read-out (HR 21). The source on fibre B is the Fabry-Perot interferometer. Observations were made with a maximum airmass of 1.8 and a minimum 30 • separation from the Moon.
The measurements span from September 29, 2019 to January 20, 2020 and have an average nominal error bar of 93 cm/s. We also included the time series of Hα measurements, the full width half maximum (FWHM), and the S -index in our analysis. The velocity and ancillary indicators are extracted from the spectra with the standard ESPRESSO pipeline v 2.0.0 (Pepe et al. 2020). The RV time series with nominal error bars is shown in Fig. 5.

Detections and parameter estimations
In this section, we analyse the photometric and spectral data in order to derive the orbital and planetary parameters of the six planets in the system.

Identification of the solution
The first release of candidates from the TESS alerts of Sector 2 included three planet candidates in TOI-178 with periods of 6.55 d, 10.35 d, and 9.96 d, which transited four, three, and two times, respectively. In addition, our analysis of this dataset with the DST (Détection Spécialisée de Transits Cabrera et al. 2012) yielded two additional candidates: a clear 3.23 d signal and a fainter 1.91 d signal. Upon receiving visits 1 and 2 from CHEOPS (Table 2), a study of the CHEOPS data alone with successive applications of the boxed-least-square (BLS) algorithm (Kovács et al. 2002) retrieved the 6.55 d, 3.23 d, and 1.91 d signals, in phase agreement with the TESS data. An additional dip, consistent in epoch with a transit of the 9.9 d candidate, was also identified; however, it could also marginally correspond to a transit of the 10.35 d candidate. The transit of one of these candidates was hence missing, pointing to a mis-attribution of transits in the TESS Sector 2 data. In order to identify new possible solutions for the available data (TESS Sector 2, NGTS visits in September and October 2019, and CHEOPS visits 1 and 2), we individually predetrended each light curve and subtracted the signal of the 6.55 d, 3.23 d, and 1.91 d candidates. We then applied the BLS algorithm to this dataset for the first time. The resulting periodogram (Fig. 6, top panel) had: several peaks of similar power due to the existence of multiple transits; dips of similar depths and durations; and, overall, a low number of transits per candidate spread along a long baseline. We hence explored different solutions by successively ignoring some of the highest peaks of the BLS periodogram, proceeding as follows: On the first periodogram, we saved the most likely candidate of period P 1 and removed the corresponding signal in the light curve; we then applied the BLS a second time to obtain a second candidate of period P 2 . That created a first pair of candidates, c i, j = c 0,0 , where i and j are the number of peaks that have been ignored in the first or second iteration of the BLS, respectively (in this case, no peaks have been ignored). We then repeated this process, but ignoring the result of the BLS for periods less than 0.2 d away from the highest peak   of the periodogram (in black in the top panel of Fig. 6). Since we ignored the black peak, the second highest is the blue one, which we assigned to period P 1 , and we removed the associated signal in the light curve. For this candidate, i = 1 since we ignored one peak in the first BLS. For the second iteration of the BLS, we did not ignore any peaks and hence take the largest one ( j = 0), leading to the pair of candidates c 1,0 . We repeated this process 25 times, yielding 25 potential pairs of candidates c 0≤i≤4,0≤ j≤4 . For all of these potential solutions, we modelled the transits of the five candidates -1.91 d, 3.23 d, 6.55 d, and c i, j -using the batman package (Kreidberg 2015) and ran an MCMC on the pre-detrended light curve to estimate the relative likelihood of the different c i, j . The c 3,0 pair was favoured (its first BLS is shown with the green curve in the top panel of Fig. 6 Applying the BLS algorithm to the residuals of the available photometric data (bottom panel of Fig. 6), two mutually exclusive peaks appeared, 12.9 d and 15.24 d, which shared the odd transit of the previous TOI178.02 candidate. The 12.9 d signal was slightly favoured by the BLS analysis; however, the global fit of the light curve favoured the ∼ 15.24 d signal. In addition, this solution was very close to the period that would fit the res-onant structure of the system (see Sect. 6). The 15.23 d candidate was confirmed by a fourth CHEOPS visit (Fig. 4bottom right). In the next section, we develop the characterisation of this six-planet solution: 1.91 d, 3.23 d, 6.55 d, 9.96 d, 15.23 d, and 20.71 d.

Determination of radii and orbital parameters
To characterise the system, we performed a joint fit of the TESS, NGTS, and CHEOPS photometry. As the NGTS and SPECU-LOOS data presented in Fig. 2 cover the same epoch of observations, we only included the NGTS data in our fit as they had a smaller RMS photometric scatter.
The fit was performed with the juliet package (Espinoza et al. 2019), which uses batman (Kreidberg 2015) for the modelling of transits and the nested sampling dynesty algorithm (Speagle 2020;Skilling 2004Skilling , 2006aHigson et al. 2019;Buchner 2014Buchner , 2017Skilling 2006b) for estimating Bayesian posteriors and evidence. In our analysis, the fitted parameters for each planet were: the planet-to-star radius ratio R p /R , the impact parameter b, the orbital period P, and the mid-transit time T 0 . We also fitted for the stellar density ρ , which, together with the orbital period P of each planet, defines through Kepler's third law a value for the scaled semi-major axis a/R of each planet. We assumed a normal prior for the stellar density based on the value and uncertainty reported in Table 1 (Sect. 3) and wide uniform priors for the other transit parameters. The orbits were assumed to be circular, as justified in Sect. 6.3. For each bandpass (TESS, CHEOPS, and NGTS), we fitted two quadratic limb-darkening coefficients, which were parametrised with the (q 1 , q 2 ) sampling scheme of Kipping (2013). Normal priors were placed on these limb-darkening coefficients based on Claret & Bloemen (2011).
We modelled the correlated noise present in the light curves simultaneously with the planetary signals to ensure a full propagation of the uncertainties. We first performed individual analyses of each light curve in order to select for each of them the best correlated noise model based on Bayesian evidence. We explored a large range of models for the CHEOPS light curves, consisting of first-to fourth-order polynomials in the recorded external parameters (most importantly: time, background level, position of the point-spread-function (PSF) centroid, spacecraft roll angle, and contamination), as well as GPs (celerite Foreman-Mackey et al. 2017 and george Ambikasaran et al. 2014) against time, roll angle, or a combination of the two. We found that a Matérn 3/2 GP against roll angle was strongly favoured for all visits to account for the roll-angle-dependent photometric variations (cf. Sect. 4.1.2). First-to second-order polynomials in time and x − y centroid position were also needed for the first three visits. For the TESS light curve, we used a Matérn 3/2 GP against time, and we used a linear function of airmass for the NGTS light curves. For each light curve, we also fitted a jitter term, which was added quadratically to the error bars of the data points, to account for any underestimation of the uncertainties or any excess noise not captured by our modelling.
The results of our joint fit are displayed in Tables 3 (planets b, c, and d) and 4 (planets e, f , and g). The transits of all detected planets are shown in Figs. 1 (TESS), 2 (NGTS), and 4 (CHEOPS), with the phase-folded transits in the CHEOPS data presented in Fig. 7.
All planets appear to be in the super-Earth to mini-Neptune range, with the inner two planets falling to either side of the radius valley (Fulton et al. 2017). The inclination of the planets are also worth mentioning: The projected inclination of the four outer planets differ only by about 0.1 deg. As discussed in Table 3. Fitted and derived results for planets b, c, and d associated with the fits to the photometry and spectroscopy described in Sects. 5.1 and 5.2, respectively. a T eq = T eff √ R /a ( f (1 − A B )) 1/4 , assuming an efficient heat redistribution ( f = 1/4) and a null Bond albedo (A B = 0).  Table 4. Fitted and derived results for planets e, f , and g associated with the fits to the photometry and spectroscopy described in Sects. 5.1 and  Agol et al. (2020) in the case of TRAPPIST-1, it is unlikely that the underlying inclinations and ascending nodes are scattered given the clustering of the projected inclination; hence, the outer part of the system is probably quite flat. In addition, the uncertainty on the inclination of the inner planets allows for the nearcoplanarity of the entire system, a feature that was also found in the TRAPPIST-1 system (Agol et al. 2020).

Detections
In this section, we consider the 46 ESPRESSO data points only and look for potential planet detections. Our analysis follows the same steps as in Hara et al. (2020) and is described in detail in Appendix B. To search for potential periodicities, we computed the 1 -periodogram 8 of the RV as defined in Hara et al. (2017). This method outputs a figure that has a similar aspect as a regu- lar periodogram, but with fewer peaks due to aliasing. The peaks can be assigned a false alarm probability (FAP), whose interpretation is close to the FAP of a regular periodogram peak. A preliminary analysis of ancillary indicators Hα, FWHM, bisector span (Queloz et al. 2001), and log R HK (Noyes 1984) revealed that they exhibit statistically significant periodicities at ≈ 36 days and ≈ 16 days, such that stellar activity effects are to be expected in the RVs, especially at these periods. In our analysis, stellar activity has been taken into account both with a linear model constructed with activity indicators smoothed with a GP regression similarly to Haywood et al. (2014), which we call the base model, and with a Gaussian noise model with a white, correlated (Gaussian kernel), and quasi-periodic component.
Radial velocity signals found to be statistically significant might vary from one activity model to another. To confirm the robustness of our detections, we tested whether signal detections can be claimed for a variety of noise models, following Hara et al. (2020). This approach consists of three steps. First, we computed the 1 -periodogram of the data on a grid of models. The linear base models considered include an an offset and smoothed ancillary indicators (Hα, FWHM, neither, or both), where the smoothing is done with a GP regression with a Gaus- sian kernel. The noise models we considered are Gaussian with three components: white, red with a Gaussian kernel, and quasiperiodic. According to the analysis of the ancillary indicators, the quasi-periodic term of the noise is fixed at 36.5 days. We considered a grid of values for each noise component (amplitude and decay timescale) and computed the 1 -periodograms. Secondly, we ranked the noise models with a procedure based on cross-validation (CV). Finally, we examined the distribution of FAPs of each signal in the 20% highest ranked models.
Here we will succinctly describe the conclusions of our analysis and refer the reader to Appendix B for a more detailed presentation. The 1 -periodogram corresponding to the highest ranked model is represented in Fig. 9. The periods at which the peaks occur are shown in red and account for most of the signals that might appear with varying assumptions on the noise (or activity) model. More precisely, we find the following results. We note that these results stem from an analysis of over 1300 noise models and that they are not all portrayed in Fig. 9, which only shows the results from the highest ranked noise model.
First, our analysis yields a consistent, significant detection of signals close to 3.2, ≈ 36, and ≈ 16 days. Radial velocity then allows an independent detection of planet c. Signals at ≈ 36 and ≈ 16 d appear in ancillary indicators, and as such we attributed these apparent periodicities in RV to stellar activity. The 16-day signal is, however, very likely partly due to planet f . Indeed, when the activity signals close to 40 and 16 days are modelled and the signal of transiting planets is removed, one finds a residual signal at 15.1 or 15.2 days, even though the 16 and 15.2-day signals are very close.
Second, we find signals, though not statistically significant ones, at 6.5 and 9.8 days (consistent with planets d and e) and at 2.08 days, which is the one-day alias (see Dawson & Fabrycky 2010) of 1.91 days (planet b). Third, we do not consistently find a candidate near 20.7 days. However, a 20.6 d signal appears in the highest ranked model, and a stellar activity signal might hide the signal corresponding to TOI-178 g.
Fourth, the signal at 43.3 days appearing in Fig. 9 might be a residual effect of an imperfect correction of the activity. However, we have not strictly ruled out the possibility that it stems from a planetary companion at 45 days. As will be discussed in the conclusion, this period corresponds to one of the possible ways to continue the Laplace resonance beyond planet g. Finally, depending on the assumptions, hints at 1.2 or 5.6 days (aliases of each other) can appear.
For comparison, we performed the RV analysis with an iterative periodogram approach. This analysis is able to show signals corresponding to 15.2, 3.2, and 6.5 days; however, it is unable to clearly establish the significance of the 3.2 d signal and fails to unveil candidates at 1.91 and 9.9 days.
The photometric data allow us to independently detect six planets at 1. 91, 3.24, 6.56, 9.96, 15.23, and 20.71 days. As detailed in Appendix B, we find that the phases of RV and photometric signals are consistent within 2σ. We phase-folded the RV data at the periods given in Tables 3 and 4, which are shown in Fig. 10 with periods increasing from top to bottom. The variations at 3.24 and 15.2 days, corresponding to the planets with the most significant RV signals, are the clearest. As a final remark, the signals corresponding to the transiting planets have been fitted, and the strongest periodic signals occur at 38 and 16.3 days, which is compatible with the activity periods seen in the ancillary indicators.

Mass and density estimations
To estimate the planetary masses, we fitted circular orbits to the radial velocities. As shown in Sect. 6, for the system to be stable, eccentricities cannot be greater than a few percent. We set the posterior distributions obtained from the fits to photometric data from Tables 3 and 4 as priors on T c and period. This approach, as opposed to a joint fit, is justified by the fact that, here, RVs brings very little information to the parameters constrained by photometry, and vice versa. Activity signals are clearly present in the RV data, and, depending on the activity model used, mass estimates may vary. To assess the model dependency of mass estimates, we used two different activity models. Both include as a linear predictor the smoothed Hα time series, as described above. In the first model, we represented activity as a sum of two sine functions at 36 and 16 days and a correlated noise with an exponential kernel. This was motivated by the fact that both periodicities appear in the ancillary indicators, but with different phases. The correlated noise models low-frequency variations, which can also be anticipated from the analysis of the indicators. The second model of activity consists of a correlated Gaussian noise with a quasi-periodic kernel. We computed the posterior distributions of the orbital parameters with an adaptive MCMC, as in Delisle et al. (2018). This analysis is presented in detail in Appendix B.5.
The mass estimates for each model are given in Tables B.2, 3, and 4, respectively. The 1σ intervals obtained with the two methods all have a large overlap. In Tables 3 and 4 Fig. 9. 1 -periodogram corresponding to the best noise models in terms of the CV score, computed on a grid of frequencies from 0 to 0.55 cycles per day. mass and density intervals in a conservative manner: The lower and upper bounds are respectively taken as the minimum lower bound and maximum upper bound we obtained with the two estimation methods. The mass estimates are given as the mean of the estimates obtained with the two methods, which are posterior medians. We took this approach and did not select the error bars of one model or another since model comparisons depend heavily on the prior chosen, which, in our case, would be rather arbitrary. We find planet masses and densities in the ranges of 1.50 +0.39 −0.44 7.72 +1.67 −1.52 M ⊕ and 0.177 +0.055 −0.061 1.02 +0.28 −0.23 ρ ⊕ . It appears that the mass of planet f (at 15.24 days) is the highest. This result might seem untrustworthy since there are activity signals close to 16 days in the ancillary indicators and since 1/(1/15.23 − 1/16) = 316 days, which is greater than the observation time span. However, the two activity models considered here -one of which includes a signal at 16 days -yield a similar mass estimates, and the convergence of MCMC was ensured by computing the number of effective samples. As a consequence, we deemed the mass estimate appropriate. Nonetheless, more stellar activity models can be considered, and the mass intervals might still evolve as more data come along and the results become less model-dependent. A longer baseline would be suitable, so that the difference between 15.2 and 16 would be greater than the frequency resolution.

A Laplace resonant chain
Mean-motion resonances are orbital configurations where the period ratio of a pair of planets is equal to, or oscillates around, a rational number of the form (k + q)/k, where k and q are integers. In TOI-178, candidates b and c are close to a second-order MMR (q = 2): P c /P b = 1.6915 ≈ 5/3, while c, d, e, f , and g are close, pairwise, to first-order MMRs (q = 1): P d /P c = 2.0249 ≈ 2/1, P e /P d = 1.5191 ≈ 3/2, P f /P e = 1.5290 ≈ 3/2, and P g /P f = 1.3595 ≈ 4/3. Pairs of planets lying just outside MMRs are common occurrences in systems observed by transit Table 6. Estimated values of the Laplace angles of planets c to g of TOI-178 towards the beginning of TESS Sector 2 at 2458350.0 BJD (August 2018). The derivative of the angles are averaged values from between August 2018 and August 2020, based on the solution presented in Tables 3  and 4 (Fabrycky et al. 2014). To study such pairs of planets, a relevant quantity is the distance to the exact resonance. Taking the pair of c and d as an example, the distance to the 2/1 MMR in the frequency space is given by 1n c − 2n d , where n c = 2π/P c is the mean motion of planet c. Following Lithwick et al. (2012), we name the associated timescale the 'super-period': The values of these quantities are given in Table 5 for planets b to g, along with the expression of the associated angles φ i . Transit timing variations are expected over the super-period, with amplitudes depending on the distance to the resonance, the mass of the perturbing planet, and the eccentricities of the pair (Lithwick et al. 2012). The fact that the super-periods of all three of these pairs are close to the same value from planet c outwards has additional implications : The difference between the angles φ i is evolving very slowly. In other words, there is a Laplace relation between consecutive triplets: implying that the system is in a 2:4:6:9:12 Laplace resonant chain. The values of the Laplace angles ψ j and derivatives are given in Table 6. The values of the ψ j are instantaneous and computed at the date 2458350.0 BJD, which is towards the beginning of the observation of TESS Sector 2. As no significant TTVs were determined over the last two years, the derivatives of the ψ j are average values over that period. The Laplace relations described in Eq.
(2) do not extend towards the innermost triplet of the system: According to Eqs. (C.3) and (C.4), P b,c should be equal to half of P c,d for the b-c-d triplet to form a Laplace relation, which is not the case (see Table 5). To continue the chain, planet b would have needed a period of ∼ 1.95 d. Its current period of 1.91 d could indicate that it was previously in the chain but was pulled away, possibly by tidal forces. Figure 11 shows the evolution of the Laplace angles when integrating the nominal solution given in Tables 3 and 4, starting at the beginning of the observation of TESS Sector 2. The three angles librate over the integrated time for the selected initial conditions, with combinations of periods ranging from a few years to several decades. The exact periods and amplitudes of these variations depend on the masses and eccentricities of the involved planets. The theoretical equilibria of the resonant angles are discussed in Sect. 6.2, while the long-term stability of this system is discussed in Sect. 6.3.

Equilibria of the resonant chain
For a given resonant chain, there might exist several equilibrium values around which the Laplace angles could librate (Delisle 2017). For instance, the four planets known to orbit Kepler-223 are observed to librate around one of the six possible equilibria predicted by theory (Mills et al. 2016;Delisle 2017).
We used the method described in Delisle (2017) to determine the position of the possible equilibria for the Laplace angles of TOI-178. The five external planets (c to g) orbiting TOI-178 are involved in a 2:4:6:9:12 resonant chain. All consecutive pairs of planets in the chain are close to first-order MMRs (1:2, 2:3, 2:3, 3:4). Moreover, as in the Kepler-223 system, there are also strong interactions between non-consecutive planets. Indeed, planet e and planet g (which are non-consecutive) are also close to a 1:2 MMR. As explained in Delisle (2017), this breaks the symmetry of the equilibria (i.e. the equilibrium is not necessarily at 180 deg), and the position of the equilibria for the Laplace angles depends on the planets' masses.
We solved for the position of these equilibria using the masses given in Tables 3 and 4. We also propagated the errors to estimate the uncertainty on the Laplace angle equilibria. We find two possible equilibria for the system, which are symmetric with respect to 0 deg. We provide the values of the Laplace angles corresponding to the first equilibrium in Table 6 (the second is simply obtained by taking ψ j → −ψ j for each angle).
It should be noted that these values correspond to the position of the fixed point around which the system is expected to librate. Depending on the libration amplitude, the instantaneous values of the Laplace angles can significantly differ from the equilibrium. For instance, in the case of Kepler-223, the amplitude of libration could be determined and is about 15 deg for all Laplace angles (Mills et al. 2016). In Table 6, we observe that all instantaneous values of Laplace angles (as of August 2018) are also found within 15 deg of the expected equilibrium. We note that ψ 1 was moving away from the equilibrium throughout the two years of observations (see the averaged derivatives in Table 6). This would imply a final (as of September 2020) distance from equilibrium of about 23 deg. On the other hand, ψ 2 was moving closer to the equilibrium and ψ 3 was only slowly evolving during this two-year span. A more precise determination of the evolution of Laplace angles (with more data and the detection of TTVs) would be needed to estimate the amplitude of libration of each angle. However, it is unlikely to find values so close to the equilibrium for each of the three angles by chance alone. Therefore, these results provide strong evidence that the system is indeed librating around the Laplace equilibrium.

Stability
Planets c to g are embedded in a resonant chain, which seems to greatly stabilise this planetary system. The distance between the planets being quite small, their eccentricities may become a major source of instability. This point is illustrated by Fig. 12, which shows a section of the system's phase space. The initial conditions and masses of the planets are displayed in Tables 3  and 4 tial period and eccentricity vary; all other planets start on circular orbits. The colour code corresponds to a stability index based on the diffusion of the main frequencies of the system defined as (Laskar 1990(Laskar , 1993:  Fig. 11. Example of the evolution of the Laplace angles over 100 years, starting from TESS observation of Sector 2, using the masses and orbital parameters from Tables 3 and 4. where n (1) f and n (2) f are the proper mean motion of planet f in degree per year, computed over the first half and the second half of the integration, respectively, implying that the stability increases from red to black (for more details, see Sect. 4.1 of Petit et al. 2018). This map shows that the eccentricity of planet f must not exceed a few hundredths if the system is to remain stable. The same stability maps (not reproduced here) made for each of the planets of the system lead to the same result: The planetary eccentricities have to be small in order to guarantee the system's stability. This constrain is verified if the system starts with sufficiently small eccentricities. In particular, we verified that starting a set of numerical integrations on circular orbits with masses and initial conditions close to those of the nominal solution and integrating the system over 100 000 years (i.e. about 19 million orbits of planet b or 1.3 million orbits of planet g) does not excite the eccentricities above 1/100 in most cases.
All the P i versus e i stability maps, such as the one in the top panel of Fig. 12, show long quasi-vertical structures that contain very stable regions in their central parts and less stable regions at their edges. These are mainly MMRs between two or three planets. In particular, the black area crossed by the dashed white line in the top panel of Fig. 12 corresponds to the resonant chain where the nominal system is located.
The bottom panel of Fig. 12 shows a different section of the phase space where P f and P g vary and all eccentricities are initialised at 0. This reveals part of the resonant structure of the system. The black regions are stable, while the green to red areas mark the instability caused by the resonance web. This figure still highlights the stability island in which the planetary system is located. This narrow region is surrounded by resonances: the n e − 3n f + 2n g = 0 Laplace resonance (central diagonal strip), high-order two-body MMRs between planets f and g (horizontal strips), and high-order two-body MMRs between planets e and f (vertical lines).
The extent of this resonant chain versus the initial period and mass of planet f is shown in Fig. 13. On both panels, the x-axis corresponds to the orbital period of planet f , while the yaxis corresponds to the planetary mass. The figure presents two different, but complementary, aspects of the dynamics of TOI-178. The bottom panel indicates whether the Laplace angles ψ 1 , ψ 2 , and ψ 3 librate or not. More precisely, the colour code corre-sponds to the number of angles that librate during the first 200 years of integration. The top panel present the same stability indicator as Fig. 12. Three regions stand out clearly from this figure, each with a different dynamical regime.
The central region (yellow in the bottom panel), where the three Laplace angles librate simultaneously (see also Fig. 11), shows the heart of the 2:4:6:9:12 resonant chain, where the nominal system is located. On the stability map (top), its dark black colour reveals a very low diffusion rate and therefore a very longterm stability. This central region seems not to depend strongly on the mass of planet f .
On the other hand, the modification of the orbital period within a fraction of a day changes the dynamics of the system. Indeed, outside of the central region where the three Laplace angles librate, which is about 0.015 days wide, chaotic layers are present (red in the top panel panel, dark blue in the bottom) where none of the three Laplace angles librate. Here, the red colour of the stability index corresponds to a significant, but moderate, diffusion rate. Thus, although in this region the trajectories are not quasi-periodic, the chaos has limited consequences (it is bounded) and probably does not lead to the destruction of the system.
Outside of these layers lies a very stable (quasi-periodic) region. The blue on the Laplace angle map shows that only one Laplace angle librates, while the others circulate. Although the 2:4:6:9:12 chain is broken, planets c, d, e, and g remain inside the 2:4:6:12 resonant chain, independently of the orbital period of planet f . The stability map indicates a very strong regularity of the whole region. Nevertheless, one can notice the presence of some narrow zones where the diffusion is more important; they are induced by high-order orbital resonances but do not have any significant consequence on the stability of the system.
The study described in this section gives only a very partial picture of the structure of the phase space (parameter space) of the problem. Nevertheless, it can be seen that as long as the system is in the complete resonant chain, which is the case for the nominal derived parameters and for the bulk of the posterior given in Tables 3 and 4, it remains stable.

Expected TTVs
Since planets c, d, e, f , and g are, pairwise, close to first-order MMRs, we expect TTVs to occur over the super-period (Eq. 1, see also Lithwick et al. 2012), which is roughly 260 days for all these pairs (Table 5). The amplitude of these TTVs depends on the masses and eccentricities of each planet of a pair. Since the stability analysis concludes that the system is more stable with eccentricities close to 0, we integrated the initial conditions described in Tables 3 and 4 over 6 years, starting during the observation of TESS Sector 2 using the rebound package (Rein & Liu 2012;Rein & Spiegel 2015). The resulting TTVs are shown in Fig. 14. In addition to the terms coming from the super-periods, the six-year evolution of the TTVs shows hints of the long-term evolution of the Laplace angles.
As the exact shape of the TTVs depends mainly on the masses of the involved planets, the sampling of the 260-day signal over a few years, in addition to a longer monitoring of the evolution of the Laplace angles, can provide precise constrains on the orbital configuration and masses of this system, as well as the presence of other planets in the chain. The detailed study of the TTVs of this system will be the subject of a forthcoming paper.   Table 4. heavy elements one would expect in a protoplanetary disc similar to the Minimum Mass Solar Nebula (MMSN) but around TOI-178. Assuming a disc mass equal to 1% of the stellar mass, and using the metallicity of TOI-178 ([Fe/H] = −0.23, corresponding to a metal content of Z ∼ 0.0061), the mass of heavy elements in such a disc would equal ∼ 13M ⊕ , which is remarkably similar to the minimum mass of heavy elements in planets, as mentioned above. The concept of the MMSN, scaled appropriately to reflect the reduced mass and metallicity of TOI-178, seems therefore to be as applicable for this system as for the Solar System. We note, however, that in (Hayashi 1981) the MMSN model is assumed to be based on the in situ formation of planets, whereas here we compare the mass of planets with the whole mass of solids in the inner parts. Another implication of this comparison is that a problem in term of available mass would appear should the planetary masses be revised to higher values or should another massive planet be detected in the same system. In such a case, the TOI-178 system would point towards a for-mation channel similar to the one envisioned for the Trappist-1 system (Schoonenberg et al. 2019).

Mass-radius relation
The six planets we detected in the TOI-178 system are in the super-Earth to mini-Neptune range, with radii ranging from 1.152 +0.073 −0.070 to 2.87 +0.14 −0.13 R ⊕ . Although the mass determination is limited by the extent of the available spectroscopic dataset, planets b and c appear to have roughly terrestrial densities of 0.98 +0.35 −0.31 ρ ⊕ and 1.02 +0.28 −0.23 , respectively, where ρ ⊕ is the density of the Earth. The outer planets seem to have significantly lower densities; in particular, we estimate the density of planet d to be 0.177 +0.055 −0.061 ρ ⊕ . Figure 15 shows the position of the six planets in a massradius diagram, in comparison with planets with mass and radius uncertainties less than 40% (light grey). Planets belonging to five systems in Laplace resonance are indicated in the same diagram: Trappist-1, K2-138, Kepler-60, Kepler-80, and Kepler 223. The diversity of planetary composition in TOI-178 is clearly visible in the diagram: The two inner planets have a radius compatible with a gas-free structure, whereas the others contain water and/or gas. This is similar to the Kepler-80 system, where the two innermost planets are compatible with a gas-free structure and the two outermost ones likely contain gas. Planets in Kepler-223, Kepler-60, and Trappist-1 seem to have a more homogeneous structure: All planets in Kepler-223 have a gas envelope, and the planets in Kepler-60 and Trappist-1 have small gas envelopes (below ∼ .01% of the total mass).
Considering the TOI-178 system in greater detail, planets d, e, and g are located above the pure water line and definitely contain a non-zero gas mass fraction. Planet d and, depending on its mass, planet g are located in a part of the diagram where no other planets exist (at least no planets with mass uncertainties smaller than 40%) and must contain a large gas fraction. Figure 16 is a diagram showing the location of the TOI-178 planets in stellar insolation versus the planetary radius. The colour code illustrates the density of exoplanets. This diagrams clearly shows the so-called evaporation valley 9 . Planets b and c are located below the valley, and their high densities could result from the evaporation of a primordial envelope. The outer planets are located above the valley and have probably preserved (part of) their primordial gas envelope.

Comparison with other systems in Laplace resonance
The difference between the TOI-178 system and other systems in Laplace resonance can clearly be seen in Fig. 17, where we show the density of planets as a function of their equilibrium temperature for the same systems as in Fig. 15. In Kepler-60, Kepler-80, and Kepler-223, the density of the planets decreases when the equilibrium temperature decreases. This can be understood as the effect of evaporation that removes part of the primordial gaseous envelope, this effect being stronger for planets closer to their stars. The densities of K2-138 are also 1-σ consistent with such behaviour. In Trappist-1, the density of planets is always higher than 4 g/cm 3 and increases (with the exception of Trappist-1 f) with decreasing equilibrium temperature. This is likely due to the presence of more ices in planets far from their A&A proofs: manuscript no. TOI178 stars (Agol et al. 2020). Contrary to the three Kepler systems, in the TOI-178 system the density of the planets is not a growing function of the equilibrium temperature. Indeed, TOI-178f has a density higher than that of planet e, and TOI-178d has a density smaller than that of planet e.
Planet f is substantially more massive than all the other planets in the system. From a formation perspective, one would expect this planet to have a density smaller than the other planets, in particular planet e, at the end of the formation phase, similar to what we observe, for example, for Jupiter, Saturn, Uranus, and Neptune. Since this planet is farther away from the star compared to planet e, evaporation should have been less effective for planet f compared to planet e. The combined effect of formation and evolution should therefore lead to a smaller density for planet f compared to planet e. Similarly, planet d is smaller than planet e and is located closer to the star. Using the same arguments, the combined effect of formation and evolution should have led to planet d having a density larger than that of planet e. The TOI-178 system therefore seems at odds with the general understanding of planetary formation and evaporation, where one would expect the density to decrease when the distance to the star increases, or when the mass of the planet increases.

Internal structure modelling
We used a Bayesian analysis to compute the posterior distribution of the internal planetary structure parameters. The method we used closely follows the one in Dorn et al. (2015) and Dorn et al. (2017) and has already been used in Mortier et al. (2020) and Delrez (2020submitted). Here we review the main physical assumptions of the model.
The model is split into two parts. The first is the forward model, which provides the planetary radius as a function of the internal structure parameters, and the second is the Bayesian analysis, which provides the posterior distribution of the internal structure parameters, given the observed radii, masses, and stellar parameters (in particular their composition). For the forward model, we assume that each planet is composed of four layers: an iron-sulfur inner core, a mantle, a water layer, and a gas layer. We used the equation of state (EOS) for Hakim et al. (2018) for the core, the EOS from Sotin et al. (2007) for the silicate mantle, and the EOS from Haldemann et al. (2020) for the water. These three layers constitute the 'solid' part of the planets. The thickness of the gas layer (assumed to be made of pure H/He) is computed as a function of the stellar age, mass, and radius of the solid part as well as irradiation from the star, using the formulas in Lopez & Fortney (2014). The internal structure parameters of each planet are therefore the iron molar fraction in the core, the Si and Mg molar fraction in the mantle, the mass fraction of all layers (inner core, mantle, and water), the age of the planet (equal to the age of the star), and the irradiation from the star. More technical details regarding the calculation of the forward model are given in Appendix D.
In the Bayesian analysis part of the model, we proceeded in two steps. We first generated 150000 synthetic stars, taking their masses, radii, effective temperatures, and ages at random following the stellar parameters computed in Sect. 3. The Fe/Si/Mg bulk molar ratios in the star are assumed to be solar, with an uncertainty of 0.05 (uncertainty on [Fe/H], see Sect. 3). For each of these stars, we generated 1000 planetary systems, varying the internal structure parameters of all planets and assuming that the bulk Fe/Si/Mg molar ratios are equal to the stellar ones. We then computed the transit depth and RV semi-amplitude for each of the planets and retained the models that fitted the observed data within the error bars. By doing so, we included the fact that all synthetic planets orbit a star with exactly the same parameters. Indeed, 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, rather than fitting each planet independently.
For the Bayesian analysis, we assumed the following two priors: first, that the mass fraction of the gas envelope is uniform in log; and second, that the mass fraction of the inner core, mantle, and water layer are uniform on the simplex (the surface on which   16. Diagram of the position of the TOI-178 planets in a stellar light intensity (relative to Earth) versus planetary radius. The marker size is inversely proportional to the density, and the colour code gives the density of exoplanets, from yellow for empty regions of the diagram to violet for high-density and/or highly populated regions.
they add up to one) for the solid part. We also assumed: the mass fraction of water to be smaller than 50% Marboeuf et al. 2014); the molar fraction of iron in the inner core to be uniform between 0.5 and 1; and the molar fraction of Si, Mg, and Fe in the mantle to be uniform on the simplex (they also add up to one). In order to compare TOI-178 with other systems with a Laplace resonance, we also performed a Bayesian analysis for the  systems to compute the probability distribution of the gas mass in the different planets. The parameters for all systems are taken from the references mentioned above, and for all systems we have assumed that [Si/H]=[Mg/H]=[Fe/H]. We did not consider the Trappist-1 system in this comparison as it is likely that the variations in the densities of these planets result from variations in their ice contents (Agol et al. 2020).
The posteriors distributions of the two most important parameters (mass fractions and composition of the mantle) of each planet in TOI-178 are shown in Appendix D, Figs. D.2 to D.7. We focus here on the mass of gas in each planet, and we plot the mass of the gaseous envelope for each planet as a function of their equilibrium temperature in Fig. 18.
In all three Kepler systems, the gas mass in planets generally decreases when the equilibrium temperature decreases, as can be seen in Fig. 18. One exception to this tendency is the mass of the envelope of Kepler-223d, which is larger than that of planet e in the same system. Kepler-223d is, however, more massive than planets c and e in the same system, and one would expect from formation models that the mass of the primordial gas envelope is a growing function of the total planetary mass. In K2-138, the gas fraction is compatible (within 1-sigma error bars) with a monotonous function of the equilibrium temperature. In the case of TOI-178, the mass of gas also globally increases when the equilibrium temperature decreases, with the notable exception of planet d. Indeed, a linear interpolation would provide a gas mass for planet d of the order of 10 −6 − 10 −5 M ⊕ , whereas the interior structure modelling gives values of 9.66 × 10 −3 M ⊕ Fig. 17. Densities of planets in the TOI-178, Trappist-1, K2-138, Kepler-60, Kepler-80, and Kepler-223 systems as a function of their equilibrium temperatures. The error bars give the 16% and 84% quantiles, and the marker is located at the median of the computed distribution. The colour code gives the planetary masses in Earth masses. The parameters of the planets are taken from the references mentioned in Fig. 15. and 2.56 × 10 −2 M ⊕ for the 16% and 84% quantiles. Even more intriguing, our results show that the amount of gas in planet d is larger than in planet e, the latter being both more massive and at a larger distance from the star. Indeed, from the joint probability distribution of all planetary parameters as provided by the Bayesian model, the probability that planet d has more gas than planet e is 92%.
From a formation point of view, one would generally expect that the mass of gas is a growing function of the core mass. From an evolution perspective, one would also expect that evaporation is more effective for planets closer to the star. Both would point towards a gas mass in planet d that should be smaller than that of planet e. The large amount of gas in planet d and, more generally, the apparent irregularity in the planetary envelope masses are surprising in view of the apparent regularity of the orbital configuration of the system, which was presented in Sect. 6.

Discussion and conclusions
In this study, we presented new observations of TOI-178 by CHEOPS, ESPRESSO, NGTS, and SPECULOOS. Thanks to this follow-up effort, we were able to determine the architecture of the system: Out of the three previously announced candidates at 6.56 d, 9.96 d, and 10.3 d, we confirm the first two (6.56 d and 9.96 d) and re-attribute the transits of the third to 15.23 d and 20.7 d planets, in addition to the detection of two new inner planets at 1.91 d and 3.24 d; all of these planets were confirmed by follow-up observations. In total, we therefore announce six planets in the super-Earth to mini-Neptune range, with orbital periods from 1.9 d to 20.7 d, all of which (with the exception of the innermost planet) are in a 2:4:6:9:12 Laplace resonant chain. All the orbits appear to be quasi-coplanar, with projected mutual inclinations between the outer planets estimated at about 0.1 deg, a feature that is also visible in other systems with three-body resonances, such as Trappist-1 (Agol et al. 2020) and the Galilean moons. Current ephemeris and mass estimations indicate that the TOI-178 system is very stable, with Laplace angles librating over decades. In the TESS Sector 29 data, made available during the referee process of this paper, we recovered the transits of planets b, c, d, e, and f at the predicted dates. The transit Table 7. Potential periods that would continue the resonant chain in a near (k + q) : k MMR with planet g, taking the super-period at 260 days (see equations in Appendix C). Changing the super-period by a few days typically changes the resulting period by less than 0.1 day. of planet g was not observed as it transited during the gap and was observed during the third CHEOPS visit (see Fig. 4).
As there is no theoretical reason for the resonant chain to stop at 20.7 days, and the current limit probably comes from the duration of the available photometric and RV datasets, we give in Table 7 the periods that would continue the resonant chain for first-order (q = 1) and second-order (q = 2) MMRs. These periods result from the equations detailed in Appendix C. In the TOI-178 system, as well as in similar systems in Laplace resonance, planet pairs are virtually all nearly first-order MMRs. The most likely of the periods shown in Table 7 are therefore the first-order solutions with low k, hence 45.00, 32.35, or 28.36 days; however, additional planets with such periods would not be transiting if they were in the same orbital plane as the others (see Fig. 8). We note that, for a star such as TOI-178, the inner boundary of the habitable zone lies around 0.2 AU, or at a period of the order of 40 days. Additional planets in the Laplace resonance could therefore orbit inside, or very close to, the habitable zone.
The brightness of TOI-178 allows for further characterisation of the system by photometric measurements, radial velocities, and transit spectroscopy. These measurements will be es- sential for further constraining the system, not only on its orbital architecture but also for the physical characterisation of the different planets.
As discussed above, the current mass and radius determinations show significant differences between the different planets. It appears that the two innermost planets are likely to be rocky, which may be due to the fact that they have lost their primary, hydrogen-dominated atmospheres through escape (Kubyshkina et al. 2018(Kubyshkina et al. , 2019, whereas all the other planets may have retained part of their primordial gas envelope. In this respect, the different planets of the TOI-178 system lie on both sides of the radius valley (Fulton et al. 2017). Therefore, reconstructing the past orbital and atmospheric history of this planet may provide clues regarding the origin of the valley. The system is located at a declination in the sky that makes it observable by most groundbased observatories around the globe. Furthermore, the host star is bright enough and the radii of the outer planets large enough to make them potentially amenable to optical and infrared transmission spectroscopy observations from both ground and space, particularly by employing the upcoming E-ELT and James Webb Space Telescope (JWST) facilities. Indeed, it has been shown that the JWST has the capacity to perform transmission spectroscopy of planets with radii down to 1.5 R ⊕ (Samuel et al. 2014).
Planets d and f are particularly interesting as their densities are very different from those of their neighbours, and they depart from the general tendency of planetary density (i.e. decreasing for decreasing equilibrium temperatures). The densities of planets d and f , in the context of the general trend seen in TOI-178, are difficult to understand in terms of formation and evaporation process, and they could be difficult to reproduce with planetary system formation models (Mordasini et al. 2012;Alibert et al. 2013;Emsenhuber et al. 2020). We stress that even though two different analyses yielded similar estimates of masses and densities, they were performed with only 46 data points. As such, these estimates need to be confirmed by further RV measurements, which would provide, in particular, a better frequency resolution and confirm the mass estimate of planet f .
The orbital configuration of TOI-178 is too fragile to survive giant impacts, or even significant close encounters: Fig. 12 shows that a sudden change in period of one of the planets of less than a few ∼ .01 d can render the system chaotic, while Fig. 13 shows that modifying a single period axis can break the resonant structure of the entire chain. Understanding, in a single framework, the apparent disorder in terms of planetary density on one side and the high level of order seen in the orbital architecture on the other side will be a challenge for planetary system formation models. Additional observations with CHEOPS and RV facilities will allow further constraining of the internal structures of all planets in the system, in particular the (lack of) similarity between the water fraction and gas mass fraction between planets.
Follow-up transit observations should also unveil TTVs for all but the innermost planet of the TOI-178 system (see Fig. 14), with two timescales: a ∼ 260 d signal with an amplitude of several minutes to a few tens of minutes and a larger signal over several years or decades. Future observations of this system can hence be used to measure the planetary masses and estimate eccentricities directly from TTVs as well as compare the results with masses derived from RV measurements.
Finally, the innermost planet, b, lies just outside the 3:5 MMR with planet c; however, it is too far to be part of the Laplace chain, which would require a period of ∼ 1.95d. Since the formation of the Laplace resonant chain probably results from a slow drift from a chain of two-planet resonances due to tidal effects (Papaloizou & Terquem 2010;Delisle et al. 2012;Papaloizou 2015;MacDonald et al. 2016), the current state of the system might constrain the dissipative processes that tore apart the innermost link of the chain while the rest of the configuration survived.
The TOI-178 system, as revealed by the recent observations described in this paper, contains a number of very important features: Laplace resonances, variation in densities from planet to planet, and a stellar brightness that allows a number of followup observations (photometric, atmospheric, and spectroscopic). It is therefore likely to become one of the Rosetta Stones for understanding planet formation and evolution, even more so if additional planets continuing the chain of Laplace resonances is discovered orbiting inside the habitable zone.   frequencies, and more largely certain signals, are in the data by default and should not be penalised by the 1 norm. We can define a linear model whose column vectors are not penalised, which we call the base model. The signals found to be statistically significant might vary depending on the frequency grid, base, and noise models. To explore this aspect, as in Hara et al. (2020), we computed the 1periodogram of the data with different assumptions regarding the noise covariance. We then ranked the covariance models via CV. That is, we fixed a frequency grid. Secondly, for every choice of base and noise models, we recorded which detections are announced. We assessed the score of the detections+noise models via CV. The data were separated randomly into a training set and a test set that contain, respectively, 70 and 30% of the data. The model was fitted onto the training set, and we computed the likelihood of the data on the test set. This operation was repeated 250 times, and we attributed the median of the 250 scores to the triplet base model, covariance model, and signal detected. To determine if a signal at a given period was significant, we studied the distribution of its FAP among the highest ranked models.

Appendix B.2: Definition of the alternative models
To define the alternative RV models we, as a preliminary step, analysed the Hα, FWHM, bisector span, and log R HK time series as provided by the ESPRESSO pipeline. We computed the residual periodograms as described in (Baluev 2008). These periodograms allowed us to take into account general linear base models that are fitted along candidate frequencies. We computed the periodograms and iteratively added a sinusoidal function whose frequency corresponds to the maximum peak of the periodogram. Iterations 1 and 2 are shown in the top and bottom panels of Fig. B.1, respectively. The 36-day and 16-day positions (top and bottom panels, respectively) are marked with dotted black lines. Pursuing the iterations, we find signal detections with FAP <10 −3 of periods at 36, 115, and 15.9 days for Hα; 35.5, 20.8, and 145 days for the FWHM; 36 and 16 days for the bisector span; and 36.7 and 16.5 days for the log R HK . The 36day periodicities are always detected with FAP < 10 −6 . These results indicate that activity effects in the RV at ≈ 36 and ≈ 16 days are to be expected, along with low frequency effects. These signals likely stem from the rotation period of the star, which creates signals at the fundamental frequency and the first harmonic.
We now turn to the RV and define the alternative noise models we explored. These are Gaussian, with a white component, an exponential decay, and a quasi-periodic term, as given by the formula where V kl is the element of the covariance matrix at row k and column l; δ k,l is the Kronecker symbol; σ k is the nominal uncertainty on the measurement k; and σ W , σ R , τ R , σ QP , and P act are the parameters of our noise model. A preliminary analysis on the ancillary indicators (FWHM, S-index, and Hα) showed that they all exhibit statistically significant variations at ≈ 40 (36 days), as well as a periodogram peak at 16 d. The 36 and 16 d signals exhibit phases compatible with each other at 1 sigma, except for the 36 d signal in the FWHM which is 3σ away from the phase fitted on the S-index and Hα. This points to a stellar rotation period of 36 d. We considered all the possible combinations of values for σ R and σ W in 0.0,0.5 1.0,1.25, 1.5,1.75,2 m/s, τ = 0, 2, 4, 6 d, P act = 36.5 d, σ QP = 0,1,2,3,4 m/s, and τ QP = 18, 36, or 72 d.
The computation of the 1 -periodogram was made assuming a certain base model. By this, we mean a linear model that is assumed to be in the data by default and will thus automatically be fitted. This base model could represent, for instance, offsets, trends, or certain periodic signals. We tried the following base models: one offset, one offset and smoothed Hα, one offset and smoothed FWHM, and one offset and the smoothed Hα and FWHM time series. The smoothing of a given indicator is done via a GP. This process has a Gaussian kernel, whose parameters (timescale and amplitude) have been optimised to maximise the likelihood of the data.

Appendix B.3: Results
In Fig. B.2, we represent the 1 -periodogram obtained with the noise model with the maximum CV score with three different base models (from top to bottom: without any indicator, with smoothed Hα, and with both smoothed Hα and smoothed FWHM). In all cases, we find signals at 35-40 d, 15 d, and 3.24 d, as well as a peak at 6.55 d. We can also find peaks at 21 d, 9.84 d, 2.08 d, 1.92 1.44, or 1.21 days, depending on the models used. We note that 2.08 and 1.92 d are aliases of each other, so that the presence of one or the other might originate from the same signal. The model with the overall highest CV score includes only the smoothed Hα time series in the base model and with the notations of Eq. (B.1), σ W = 1.75 m/s, σ R = 1.5 m/s, τ = 2 days, and σ QP = 0, and it corresponds to the middle panel in Fig. B.2.
We computed the 1 -periodogram on a grid from 0 to 0.95 cycles per day, hence ignoring periods below one day, with all combinations of the base and noise models. As in Hara et al. (2020), the models are ranked by CV. We then considered the 20% highest ranked models (all noise and base models considered), which we denote with CV 20 , and computed the number of times a signal is included in the model. A signal is considered to be included if a peak has a frequency within 1/T obs of a reference frequency 1/P 0 , where T obs is the observation time span and has an FAP below 0.5. We report these values in Table B.1 for the reference periods P 0 that correspond to signals appearing at least once in Fig. B.2. We note that, due to the short time  activity used. To account for this, we estimated masses with two different stellar activity models: (1) Activity is modelled as two sinusoidal signals plus a correlated noise, and (2) activity is modelled as a correlated Gaussian noise with a semi-periodic kernel.
In both cases, we added the Hα time series smoothed with a GP, as described in Sect. B.3. In model 1, the two sinusoids have periods of ≈ 40 days and ≈ 16 days. This is motivated by the fact that these periodicities appeared systematically in ancillary indicators, though with different phases. We therefore allowed 2.93 · 10 −1 0.384% -3.24 1.33 · 10 −2 100.0% 2.20 · 10 −3 6.5 1.31 · 10 −1 15.57% -9.9 3.41 · 10 −1 16.63% -15.2 8.38 · 10 −2 89.32% 6.83 · 10 −2 20.7 1.00 0.0% -36 (45) 5.24 · 10 −1 95.48% 2.57 · 10 −1 the phase to vary freely in the RVs. The priors on the stellar activity periods are taken as Gaussians with mean 16 days, σ = 1 day and mean 36.8 days, σ = 8 days, according to the variability of the position of the peaks appearing in the spectroscopic ancillary indicators. We further added free noise components: a white component and a correlated component with an exponential kernel. The prior on the variance is a truncated Gaussian with σ = 100 m 2 /s 2 , and the prior on the noise timescale is loguniform between 1h and 30 days. The second model is a GP (or here, correlated Gaussian noise) with a quasi-periodic kernel k of the form We imposed a Gaussian truncated prior on σ 2 W , σ 2 R with σ = 100 m 2 /s 2 . We imposed a flat prior on ν between 2π/50 and 2π/30 rad/day and a log-uniform prior on τ on 1h to 1000 days. For planets b, c, d, g, e, and f, the priors on the periods and the times of conjunction are set as Gaussian with means and variances set according to the constraints obtained from the joint fit of the TESS and CHEOPS data. The densities are computed by combining the posterior samples of mass and the posterior samples of radii, assuming the mass and radius estimates are independent.
We ran an adaptive MCMC algorithm as described in Delisle et al. (2018), implementing spleaf (Delisle et al. 2020) 11 , which offers optimised routines to compute the inverse of covariance matrices. We checked the convergence of the algorithm by ascertaining that 600 effective samples had been obtained for each variable ). The posterior medians as well as the 1 σ credible intervals are reported in Table B.2 for both activity models. The kernel (B.2) is close to the stochastic harmonic oscillator (SHO), as defined in Foreman-Mackey et al. (2017). We also tried the SHO kernel and simply imposed that the quality factor Q be greater than one half, the other parameters having flat priors. We found very similar results, and, as such, they are not reported in Table B.2.
Modelling errors might leave a trace in the residuals. In Hara et al. (2019), it is shown that if the model used for the analysis is correct, the residuals of the maximum likelihood model, appropriately weighted, should follow a normal distribution and not exhibit correlations. In Fig. B.4, we show the histogram of the residuals (in blue) and the probability distribution function of a normal variable. The two appear to be in agreement. We further computed the Shapiro-Wilk normality test (Shapiro & Wilk 11 https://gitlab.unige.ch/Jean-Baptiste.Delisle/spleaf 1965) and found a p-value of 0.78, which is compatible with normality. The variogram does not exhibit signs of correlations. The same analysis was performed with model 2, for which the residuals also exhibit neither non-normality nor correlation (pvalue of 0.99 on the residuals). We conclude that both models are compatible with the data.

Appendix C: Continuation of a Laplace resonant chain
TOI-178 is in a configuration where successive pairs of planets are at the same distance to the exact neighbouring first-order MMR. Generalising the configuration a bit, we define fictive planets 1, 2, and 3 such that planets 1 and 2 are close to the resonance (k 1 + q 1 ) : k 1 and planets 2 and 3 are close to the reso- nance (k 2 + q 2 ) : k 2 , where k i and q i are integers. We hence write the near-resonant angles: ϕ 1 = k 1 λ 1 − (k 1 + q 1 )λ 2 , ϕ 2 = k 2 λ 2 − (k 2 + q 2 )λ 3 , (C.1) where λ i is the mean longitude of planet i. The associated distances to the resonances read: ∆ 1 = k 1 n 1 − (k 1 + q 1 )n 2 , ∆ 2 = k 2 n 2 − (k 2 + q 2 )n 3 , (C.2) where n i is the mean motion of planet i. A Laplace relation exists between these three planets if: where j 1 and j 2 are integers. In addition, the invariance by rotation of the Laplace angle (D'Alembert relation) gives: As a result, the Laplace relation requires: n 3 ≈ k 2 n 2 − q 2 q 1 ∆ 1 k 2 + q 2 , (C.5) which translates, for the period of the third planet, as: where P 1,2 is the super-period associated with ∆ 1 . From this we can compute the periods of potential additional planets that would continue the Laplace resonant chain of TOI-178. Taking planets f and g as planets 1 and 2, the formula to compute the possible period of a planet x that could continue the resonant chain becomes: 2/3 with e 3/4 with e 4/5 with e 3/2 with g 4/3 with g 5/4 with g Fig. C.1. Super-period of planet f with respect to e (in red) and g (in black) for a set of first-order MMRs. The purple line indicates the value of the super-period (260 day) between the other pairs of the chain. Only two possible periods allowed planet f to be part of the resonant chain, which correspond to the near-intersection of the right-hand side of a red curve, the left-hand side of a black curve, and the horizontal purple line.
where P f,g is the super-period between the near first-order resonances of the known chain defined by Eq. (1), here ∼ 260 days, and k and q are integers such that planet x and f are near a (k + q)/k MMR. Some of the relevant periods are displayed in Table 7. A similar computation allowed us to determine the possible period of planet f prior to its confirmation by CHEOPS (see Fig. C.1).

Appendix D: Internal structure
We provide here the posterior distributions of the interior planetary models, as well as some more details on their calculations. As mentioned in the main text, deriving the posterior distribution of the internal structure parameters requires computing the radius of planets millions of times for different sets of parameters. In order to speed up this calculation, we first computed a large (5 million-point) database of internal structure models, varying the different parameters. This database was split randomly into three sets: one training set (80% of the models), one validation set, and one test set (the last two each containing 10% of the whole database). We then, in a second time, fitted a deep neural network (DNN) in order to be able to compute the radius of a planet very rapidly with a given set of internal structure parameters. The architecture of the DNN we used is made of six layers of 2048 nodes each, and we used the classical rectified linear unit (ReLU) as an activation function (Alibert & Venturini 2019). The DNN is trained for a few hundred epochs, using a learning rate that is progressively reduced from 1.e − 2 down to 1.e − 4. Our DNN allowed us to compute these radii with an average error below 0.25% and with an increased rapidity of many orders of magnitude (a few thousand models computed per second). Figure D.1 shows the prediction error we reach on the test set (which was not used for training). The error on the predicted radius is lower than 0.4% -an error much lower than the uncertainty on the radii we obtained for the TOI-178 planets -in 99.9% of the cases. It is finally important to note that this model does not include the compression effect that would be generated Article number, page 27 of 29 A&A proofs: manuscript no. TOI178 Fig. D.1. Histogram of the prediction error from our DNN on the test set. The y-axis is in log scale, and the x-axis covers an error from -1% to 1%. by the gas envelope onto the solid part. Given the mass of the gas envelope in all planets, this approximation is justified.
The following plots show the posterior distribution of the interior structure parameters. The parameters are: the inner core, mantle, and water mass fraction relative to the mass of the solid planet; the Fe, Si, and Mg molar fraction in the mantle; the Fe molar fraction in the core; and the mass of gas (log scale). It is important to remember that since the core, mantle, and water mass fractions add up to one, they are not independent. This is also the case for the Si, Mg, and Fe molar fraction in the mantle.