Issue 
A&A
Volume 649, May 2021



Article Number  A26  
Number of page(s)  29  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202039767  
Published online  06 May 2021 
Six transiting planets and a chain of Laplace resonances in TOI178
^{1}
Observatoire Astronomique de l’Université de Genève,
Chemin Pegasi 51,
Versoix, Switzerland
email: adrien.leleu@unige.ch
^{2}
Physikalisches Institut, University of Bern,
Gesellsschaftstrasse 6,
3012
Bern, Switzerland
^{3}
Centre for Exoplanet Science, SUPA School of Physics and Astronomy, University of St Andrews,
North Haugh,
St Andrews
KY16 9SS, UK
^{4}
IMCCE, UMR8028 CNRS, Observatoire de Paris, PSL Univ., Sorbonne Univ.,
77 av. DenfertRochereau,
75014
Paris, France
^{5}
Aix Marseille Univ, CNRS, CNES, LAM,
Marseille, France
^{6}
Department of Physics, University of Warwick,
Gibbet Hill Road,
Coventry
CV4 7AL, UK
^{7}
Centre for Exoplanets and Habitability, University of Warwick,
Gibbet Hill Road,
Coventry
CV4 7AL, UK
^{8}
Astrobiology Research Unit, Université de Liège,
Allée du 6 Août 19C,
4000
Liège, Belgium
^{9}
Institute of Planetary Research, German Aerospace Center (DLR),
Rutherfordstrasse 2,
12489
Berlin, Germany
^{10}
Space sciences, Technologies and Astrophysics Research (STAR) Institute, Université de Liège,
Allée du 6 Août 19C,
4000
Liège, Belgium
^{11}
School of Physics and Astronomy, University of Leicester,
Leicester
LE1 7RH, UK
^{12}
Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto,
CAUP, Rua das Estrelas,
4150762
Porto, Portugal
^{13}
Centro de Astrofísica da Universidade do Porto,
Rua das Estrelas,
4150762
Porto, Portugal
^{14}
Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre,
4169007
Porto, Portugal
^{15}
Instituto de Astrofísica de Canarias,
38200
La Laguna,
Tenerife, Spain
^{16}
Departamento de Astrofísica, Universidad de La Laguna,
38206
La Laguna,
Tenerife, Spain
^{17}
Camino El Observatorio 1515,
Las Condes,
Santiago, Chile
^{18}
ETH Zürich, Institute for Particle Physics and Astrophysics,
Zurich, Switzerland
^{19}
Institut de Ciències de l’Espai (ICE, CSIC),
Campus UAB, Can Magrans s/n,
08193
Bellaterra, Spain
^{20}
Institut d’Estudis Espacials de Catalunya (IEEC),
08034
Barcelona, Spain
^{21}
ESTEC, European Space Agency,
2201AZ,
Noordwijk, The Netherlands
^{22}
Departmento de Astrofísica, Centro de Astrobiologia (CSICINTA), ESAC campus,
28692
Villanueva de la Cãda (Madrid), Spain
^{23}
Space Research Institute, Austrian Academy of Sciences,
Schmiedlstrasse 6,
8042
Graz, Austria
^{24}
Center for Space and Habitability,
Gesellsschaftstrasse 6,
3012
Bern, Switzerland
^{25}
Université Grenoble Alpes, CNRS, IPAG,
38000
Grenoble, France
^{26}
Department of Astronomy, Stockholm University, AlbaNova University Center,
10691
Stockholm, Sweden
^{27}
Institute of Optical Sensor Systems, German Aerospace Center (DLR),
Rutherfordstrasse 2,
12489
Berlin, Germany
^{28}
Department of Earth, Atmospheric and Planetary Science, Massachusetts Institute of Technology,
77 Massachusetts Avenue,
Cambridge,
MA
02139, USA
^{29}
Admatis,
Miskok, Hungary
^{30}
Université de Paris, Institut de physique du globe de Paris, CNRS,
75005
Paris, France
^{31}
CFisUC, Department of Physics, University of Coimbra,
3004516
Coimbra, Portugal
^{32}
INAF – Osservatorio Astronomico di Trieste,
via G. B. Tiepolo 11,
34143
Trieste, Italy
^{33}
INAF – Osservatorio Astrofisico di Torino,
via Osservatorio 20,
10025
Pino Torinese, Italy
^{34}
Lund Observatory, Dept. of Astronomy and Theoretical Physics, Lund University,
Box 43,
22100
Lund, Sweden
^{35}
INAF, Istituto di Astrofisica e Planetologia Spaziali,
via del Fosso del Cavaliere 100,
00133
Roma, Italy
^{36}
European Southern Observatory,
Alonso de Coórdova 3107,
Vitacura,
Región Metropolitana, Chile
^{37}
Leiden Observatory, University of Leiden,
PO Box 9513,
2300 RA
Leiden, The Netherlands
^{38}
Department of Space, Earth and Environment, Chalmers University of Technology, Onsala Space Observatory,
43992
Onsala, Sweden
^{39}
Dipartimento di Fisica, Università degli Studi di Torino,
via Pietro Giuria 1,
10125
Torino, Italy
^{40}
Center for Astronomy and Astrophysics, Technical University Berlin,
Hardenberstrasse 36,
10623
Berlin, Germany
^{41}
Astronomy Unit, Queen Mary University of London,
Mile End Road,
London
E1 4NS, UK
^{42}
Cavendish Laboratory, JJ Thomson Avenue,
Cambridge
CB3 0HE, UK
^{43}
University of Vienna, Department of Astrophysics,
Türkenschanzstrasse 17,
1180
Vienna, Austria
^{44}
Department of Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology,
Cambridge,
MA
02139, USA
^{45}
Departamento de Astronomía, Universidad de Chile,
Camino El Observatorio 1515,
Las Condes,
Santiago, Chile
^{46}
Centro de Astrofísica y Tecnologías Afines (CATA),
Casilla 36D,
Santiago, Chile
^{47}
Facultad de Ingeniería y Ciencias, Universidad Adolfo Ibáñez,
Av. Diagonal las Torres 2640,
Peñalolén,
Santiago, Chile
^{48}
Millennium Institute for Astrophysics, Chile
^{49}
Konkoly Observatory, Research Centre for Astronomy and Earth Sciences,
1121
Budapest,
Konkoly Thege Miklós út 1517, Hungary
^{50}
Brorfelde Observatory,
Observator Gyldenkernes Vej 7,
4340
Tølløse, Denmark
^{51}
DTU Space, National Space Institute, Technical University of Denmark,
Elektrovej 327,
2800
Lyngby, Denmark
^{52}
Institut d’astrophysique de Paris, UMR7095 CNRS, Université Pierre & Marie Curie, 98bis blvd. Arago,
75014
Paris, France
^{53}
INAF, Osservatorio Astronomico di Padova,
Vicolo dell’Osservatorio 5,
35122
Padova, Italy
^{54}
Astrophysics Group, Keele University,
Staffordshire,
ST5 5BG, UK
^{55}
Department of Physics, University of Warwick,
xsCoventry, UK
^{56}
INAF – Osservatorio Astronomico di Palermo,
Piazza del Parlamento 1,
90134
Palermo, Italy
^{57}
IFPU,
Via Beirut 2,
34151
Grignano Trieste, Italy
^{58}
Instituto de Astronomía, Universidad Católica del Norte,
Angamos 0610,
Antofagasta, Chile
^{59}
Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências da Universidade de Lisboa,
Campo Grande,
1749016
Lisboa, Portugal
^{60}
Department of Astrophysics, University of Vienna,
Tuerkenschanzstrasse 17,
1180
Vienna, Austria
^{61}
INAF, Osservatorio Astrofisico di Catania,
Via S. Sofia 78,
95123
Catania, Italy
^{62}
Dipartimento di Fisica e Astronomia “Galileo Galilei”, Università degli Studi di Padova,
Vicolo dell’Osservatorio 3,
35122
Padova, Italy
^{63}
Space Science Data Center, ASI, via del Politecnico snc,
00133
Roma, Italy
^{64}
Departmentof Physics, University of Warwick,
Gibbet Hill Road,
Coventry
CV4 7AL, UK
^{65}
Fundación G. Galilei – INAF (Telescopio Nazionale Galileo),
Rambla J. A. Fernández Pérez 7,
38712
Breña Baja,
La Palma, Spain
^{66}
INAF – Osservatorio Astronomico di Brera,
Via E. Bianchi 46,
23807
Merate, Italy
^{67}
Institut für Geologische Wissenschaften, Freie Universität Berlin,
12249
Berlin, Germany
^{68}
Paris Observatory, LUTh UMR 8102,
92190
Meudon, France
^{69}
School of Physics & Astronomy, University of Birmingham,
Edgbaston,
Birmingham,
B15 2TT, UK
^{70}
ELTE Eötvös Loránd University, Gothard Astrophysical Observatory,
9700
Szombathely,
Szent Imre h. u. 112, Hungary
^{71}
MTAELTE Exoplanet Research Group,
9700
Szombathely,
Szent Imre h. u. 112, Hungary
^{72}
INAF, Osservatorio Astronomico di Roma,
via Frascati 33,
00078
Monte Porzio Catone,
Roma, Italy
^{73}
Institute of Astronomy, University of Cambridge,
Madingley Road,
Cambridge,
CB3 0HA, UK
^{74}
Centro de Astrobiología (CSICINTA),
Crta. Ajalvir km 4,
28850
Torrejón de Ardoz,
Madrid, Spain
Received:
27
October
2020
Accepted:
14
December
2020
Determining the architecture of multiplanetary 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, TOI178 has been the subject of particular attention since the first TESS observations hinted at the possible presence of a near 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 TOI178 harbours at least six planets in the superEarth to miniNeptune 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.23}^{+0.28} to 0.177_{−0.061}^{+0.055} 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 TOI178 (H = 8.76 mag, J = 9.37 mag, V = 11.95 mag) 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 superEarths and miniNeptunes.
Key words: techniques: photometric / techniques: spectroscopic / celestial mechanics / planets and satellites: detection / planets and satellites: dynamical evolution and stability
© A. Leleu et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Since the discovery of the first exoplanet orbiting a Sunlike 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 additional candidates to known systems and by better constraining the orbital architecture, radii, and masses of the known planets.
In particular, chains of planets in meanmotion 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 quasicircular orbits leads to a high probability of capture in firstorder 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 repeat itself, trapping the planets in a chain of MMRs and leading to very closely packed configurations (Lissauer et al. 2011). However, resonant configurations are not the most common orbital arrangements (Fabrycky et al. 2014). As the protoplanetary disc dissipates, the eccentricity damping lessens, which can lead to instabilities in packed systems (see, for example, Terquem & Papaloizou 2007; Pu & Wu 2015; Izidoro et al. 2017).
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 nearresonant systems, such as HD 158259, the tides seem to have pulled the configuration entirely out of resonance (Hara et al. 2020). However, through gentle tidal evolution it is possible to retain a resonant state even with null eccentricities through threebody resonances (Morbidelli 2002; Papaloizou 2015). Such systems are too finetuned 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 threebody 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), Kepler60 (Goździewski et al. 2016), Kepler80 (MacDonald et al. 2016), Kepler223 (Mills et al. 2016), Trappist1 (Gillon et al. 2017; Luger et al. 2017), and K2138 (Christiansen et al. 2018; Lopez et al. 2019). Out of these six systems, only K2138 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. 2021 for the case of Trappist1).
In this study, we present photometric and RV observations of TOI178, a V = 11.95 mag, Ktype 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 followup 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.
2 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 TOI178 with periods of 6.55 d, 10.35 d, and 9.96 d. Based on these data, TOI178 was identified as a potential coorbital system (Leleu et al. 2019) with two planets oscillating around the same period. This prompted ESPRESSO RV measurements and two sequences of simultaneous groundbased photometricobservations with NGTS and SPECULOOS. From the latter, no transit was observed for TOI178.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 TOI178.02 transit combined with the three transits observed by TESS at high signaltonoise ratio (S/N) (above 10) was interpreted as an additional sign of the strong TTVs expected in a coorbital 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 coorbital 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 coorbital 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 TOI178.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 groundbased 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 midtransit of the 20.7day 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 P = 15.2318 d (see Appendix C, Fig. C.1), 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 day.
As we will detail in the following sections, we can confirm the detection of a 6.56day and a 9.96day 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).
3 TOI178 stellar characterisation
Fortysix ESPRESSO observations (see Sect. 4.2) of TOI178, a V = 11.95 mag Kdwarf, 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 coadded 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 stellareffective 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 macroturbulent 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 SpecMatchEmp (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 TOI178 using the infrared flux method (IRFM; Blackwell & Shallis 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 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 TOI178 using stellar evolutionary models, using as inputs T_{eff}, 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 heliumtometal 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 TOI178 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 and t_{⋆} = 7.1 Gyr. All stellar parameters are shown in Table 1.
Stellar properties of TOI178, including the methods used to derive them.
Fig. 1 Light curves from TESS Sector 2 described in Sects. 2 and 4.1.1. Unbinned data are shown as grey points, and data in 30 min bins are shown as black circles. The best fitting transit model for the system is shown in black; the associated parameter values are shown in Tables 3 and 4. The positions of the transits are marked with lines coloured according to the legend. The photometry before and after the midsector gap are shown in the top and bottom panels, respectively. As the first transit of planet f (thick teal line) occurred precisely between the two transits of the similarly sized planet g (period of 20.71 d  red lines), the three transits were originally thought to have arisen due to a single planet, which was originally designated TOI178.02 with a period of 10.35 d (see Sect. 2). 
4 Data
4.1 Photometric data
In order to determine the orbital configuration of the TOI178 planetary system, we obtained photometric time series observations from multiple telescopes, as detailed below.
4.1.1 TESS
Listed as TIC 251 848 941 in the TESS Input Catalog (TIC; Stassun et al. 2018, 2019), TOI178 was observed by TESS in Sector 2, by camera 2, from August 22, 2018 to September 20, 2018. The individual frames were processed into 2minute 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 ‘NotaNumber’ flux or flux error values. After these quality cuts, the TESS light curve of TOI178 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.
4.1.2 CHEOPS
CHEOPS, the first ESA smallclass 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^{−1} for a V = 9 magnitude star. The CHEOPS instrument is composed of an f/8 RitcheyChretien onaxis telescope (~30 cm diameter) equipped with a single frametransfer backside illuminated chargecoupled device (CCD) detector. The satellite was successfully launched from Kourou (French Guiana) into a ~700 km altitude Sunsynchronous orbit on December 18, 2019. CHEOPS took its first image on February 7, 2020, and, after it passed the inorbit commissioning (IOC) phase, routine operations started on March 25, 2020. More details on the mission can be found in Benz et al. (2021), and the first results have recently been presented in Lendl et al. (2020).
The versatility of the CHEOPS mission allows for spacebased followup photometry of planetary systems identified by the TESS mission. This is particularly useful for completing the inventory of multiplanetary systems whose outer transiting planets have periods beyond ~10 days.
We obtained four observation runs (or visits) of TOI178 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 nearcontinuous 11day photometrictime 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 lowEarth orbit of CHEOPS, the spacecrafttarget 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 TOI178, 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 sizefixed 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 TOI178, was set at R = 28.5′′. The DRP estimates the level of contamination by simulating the fieldofview of TOI178 using the Gaia star catalogue (Gaia Collaboration 2018) to determinate the location and flux of the stars throughout the duration of the visit. In the case of TOI178, 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 detrending (see Sect. 5.1). For the second visit, careful removal of one ‘telegraphic’ pixel (a pixel with a nonstable 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 DEFAULT 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, nonastrophysical shortterm 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).
Log of CHEOPS observations of TOI178.
4.1.3 NGTS
The NGTS (Wheatley et al. 2018) facility consists of twelve 20cm diameter robotic telescopes and is situated at the ESO Paranal Observatory in Chile. The individual NGTS telescopes have a wide fieldofview of 8 squaredegrees and a plate scale of 5′′ pixel^{−1}. The DONUTS autoguiding algorithm (McCormac et al. 2013) affords the NGTS telescopes subpixel 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).
TOI178 was observed using NGTS multitelescope observations on two nights. On UT September 10, 2019, TOI178 was observed using six NGTS telescopes during the predicted transit event of TESS candidate TOI178.02. However, the NGTS data for this night rule out a transit occurring during the observations. A second predicted transit event of TOI178.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 TOI178 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 TOI178 were identified using the second Gaia data release (DR2; Gaia Collaboration 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.
Fig. 2 Light curves from simultaneous observations of TOI178 by NGTS (green stars) and SPECULOOSSouth (purple triangles), described in Sects. 4.1.3 and 4.1.4, respectively. Unbinned data are shown as grey points,and data in 15minute bins are shown as green stars (NGTS) and purple triangles (SPECULOOS). The observationsoccurred on September 11, 2019 (top panel) and October 12, 2019 (bottom panel). The transit model is shown in black. The position of the odd transit of candidate TOI178.02 is shown with a dashed red line, the transit of planet g (which corresponds to an even transit of TOI178.02) is shown with the solid red line, and the transits of planet b are shown with purple lines. 
4.1.4 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 1m 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 deepdepletion CCD camera whose 13.5 μm pixel size corresponds to 0.35″ on the sky (fieldofview = 12′× 12′). Observations were performed on the nights of October 10, 2019 (for ≃8 h) and November 11, 2019 (for ≃8 h) 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 defocus 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 timevarying telluric water vapour.
4.2 ESPRESSO data
The RV datawe 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 readout (HR 21). The source on fibre B is the FabryPerot 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^{−1}. We also included the time series of Hα measurements, the full width half maximum (FWHM), and the Sindex 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. 2021). The RV time series with nominal error bars is shown in Fig. 5.
Fig. 3 Extraction of 80 × 80 arcsec of the CHEOPS fieldofview for two different data frames at the beginning (top) and the end (bottom) of the second visit (see Table 2). The TOI178 PSF is shown at the centre, with the DEFAULT DRP photometricaperture 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. 
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.
Fig. 4 Similar to Fig. 1, but instead displaying data from the four CHEOPS visits described in Sect. 4.1.2. Top panel: 11 d observation. Transits of planets c and d at ~ 2459075 BJD occur too close to each other for their corresponding lines to be individually visible in the figure. Bottomleft panel: subsequent observation scheduled to confirm the presence of a planet with a 20.7 d period (planet g), which overlaps with a transit of planet e. Bottomright panel: final observation scheduled to confirm the presence of a planet fitting in the Laplace resonance (planet f, with a period of ~ 15.23 d), which overlaps with a transit of planet b. 
Fig. 5 ESPRESSO RV data of TOI178. 
5.1 Analysis of the photometry
5.1.1 Identification of the solution
The first release of candidates from the TESS alerts of Sector 2 included three planet candidates in TOI178 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 boxedleastsquare (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 misattribution 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 thecorresponding 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 predetrended 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), and the second iteration of the BLS (shown in the middle panel) yielded a1.91 d, 3.23 d, 6.55 d, 9.96 d, and 20.71 d solution, which explained all the significant dips observed in the NGTS/SPECULOOS data (Fig. 2) and the first visit of CHEOPS (Fig. 4 – top). This solution was later confirmed by the predicted double transit observed by the third visit of CHEOPS (Fig. 4 – bottom left).
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 resonant structure of the system (see Sect. 6). The 15.23 d candidate was confirmed by a fourth CHEOPS visit (Fig. 4 – bottom right). In the next section, we develop the characterisation of this sixplanet solution: 1.91 d, 3.23 d, 6.55 d, 9.96 d, 15.23 d, and 20.71 d.
Fig. 6 Successive use of the BLS algorithm to identify the new candidates. Top panel: green curve has its highest power at 9.96 d once the three highest peaks (i = 3) at 18.35 d (black), 10.37 d (blue), and 9.17 d (orange) are ignored. Middle panel: BLS after the removal of the 9.96 d signal in the light curve, with no peaks ignored (j = 0). Bottom panel: BLS after the removal of the 20.71 d signal as well. 
5.1.2 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 SPECULOOS 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 2004, 2006; Higson et al. 2019; Buchner 2016, 2017) for estimating Bayesian posteriors and evidence. In our analysis, the fitted parameters for each planet were: the planettostar radius ratio R_{p} ∕R_{⋆}, the impact parameter b, the orbital period P, and the midtransit 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 semimajor 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 limbdarkening coefficients, which were parametrised with the (q_{1}, q_{2}) sampling scheme of Kipping (2013). Normal priors were placed on these limbdarkening 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 offirst to fourthorder polynomials in the recorded external parameters (most importantly: time, background level, position of the pointspreadfunction (PSF) centroid, spacecraft roll angle, and contamination), as well as GPs (celerite ForemanMackey et al. 2017 and george Ambikasaran et al. 2015) 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 rollangledependent photometric variations (cf. Sect. 4.1.2). First to secondorder 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 phasefolded transits in the CHEOPS data presented in Fig. 7.
All planets appear to be in the superEarth to miniNeptune 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 Agol et al. (2021) in the case of TRAPPIST1, 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 TRAPPIST1 system (Agol et al. 2021).
5.2 Analysis of the radial velocities
5.2.1 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 asdefined in Hara et al. (2017). This method outputs a figure that has a similar aspect as a regular 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 (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 quasiperiodic 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 offset and smoothed ancillary indicators (Hα, FWHM, neither, or both), where the smoothing is done with a GP regression with a Gaussian 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 quasiperiodic 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 crossvalidation (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 16day 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.2day 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 oneday 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 TOI178 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 phasefolded 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 withthe most significant RV signals, are the clearest. As a final remark, the signals corresponding to the transiting planets havebeen fitted, and the strongest periodic signals occur at 38 and 16.3 days, which is compatible with the activity periods seenin the ancillary indicators.
Fig. 7 Detrended CHEOPS light curves phasefolded to the periods of each of the planets, with signals of the other planets removed. Unbinned data are shown with grey points, data in 15min bins are shown with coloured circles, and samples drawn from the posterior distribution of the global fit are shown with coloured lines. 
Fig. 8 Impact parameter (top) and inclination (bottom) of the six planets of the TOI178 system. In the top panel, the solid line shows the evolution of the impact parameter as a function of the period assuming that all planets are in the same plane (obtained by linear fit). The dashed line shows that an outer planet in the same plane will not transit if its period is above 26.4 d. 
5.2.2 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 lowfrequency 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 quasiperiodic 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, we give the 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 M_{⊕} and ρ_{⊕}.
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 modeldependent. A longer baseline would be suitable, so that the difference between 15.2 and 16 would be greater than the frequency resolution.
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. 
Fig. 10 Phasefolded RV. Error bars correspond to nominal errors. 
Instantaneous distances from resonances for the six planets of TOI178, expressed in terms of nearresonant angles ϕ.
6 Dynamics
6.1 A Laplace resonant chain
Meanmotion 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 TOI178, candidates b and c are close to a secondorder MMR (q = 2): P_{c} ∕P_{b} = 1.6915 ≈ 5∕3, while c, d, e, f, and g are close, pairwise, to firstorder 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 (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 ‘superperiod’: (1)
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 superperiod, 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 superperiods of all four 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: (2)
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 2 458 350.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. (3) 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 bcd 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 longterm stability of this system is discussed in Sect. 6.3.
Estimated values of the Laplace angles of planets c to g of TOI178 towards the beginning of TESS Sector 2 at 2 458 350.0 BJD (August 2018).
6.2 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 Kepler223 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 TOI178. The five external planets (c to g) orbiting TOI178 are involved in a 2:4:6:9:12 resonant chain. All consecutive pairs of planets in the chain are close to firstorder MMRs (1:2, 2:3, 2:3, 3:4). Moreover, as in the Kepler223 system, there are also strong interactions between nonconsecutive planets. Indeed, planet e and planet g (which are nonconsecutive) 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 Kepler223, 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 twoyear span. A more precise determination of the evolution of Laplace angles (withmore 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.
Fig. 11 Example of the evolution of the Laplace angles over 100 yr, starting from TESS observation of Sector 2, using the masses and orbital parameters from Tables 3 and 4. 
Fig. 12 Stability indicator (defined in Eq. (3)) for TOI178 as a function of the periods and eccentricity of planet f (top). The red areas show the initial conditions of unstable trajectories, while black shows stable (quasiperiodic) trajectories. Bottom panel: same stability criterion with respect to the initial periods of planets f and g. The dashed white lines show the observed periods reported in Table 4. 
6.3 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, with the exception of those of planet f, for which the initial 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, 1993): (3)
where and 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 yr (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 quasivertical 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), highorder twobody MMRs between planets f and g (horizontal strips), and highorder twobody 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 xaxis 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 TOI178. The bottom panel indicates whether the Laplace angles ψ_{1}, ψ_{2}, and ψ_{3} librate or not. More precisely, the colour code corresponds to the number of angles that librate during the first 200 yr 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 quasiperiodic, 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 (quasiperiodic) 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 highorder 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.
Fig. 13 Same stability indicator as in Fig. 12 (top) and the number of librating Laplace angles (bottom) for TOI178 as a function of the period and mass of planet f. The dashed white line shows the observed period of planet f reported in Table 4. 
6.4 Expected TTVs
Since planets c, d, e, f, and g are, pairwise, close to firstorder MMRs, we expect TTVs to occur over the superperiod (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 yr, 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 superperiods, the sixyear evolution of the TTVs shows hints of the longterm 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 260day 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.
7 Internal structure
7.1 Minimum mass of the protoplanetary disc
The total mass of the planets detected in TOI178 is approximately . Assuming a mass fraction of H and He of maximum ~20% (similar to those of Uranus and Neptune), the amount of heavy elements in the planets is at least 14.85 M_{⊕} (16th quantile). This number can be compared with the mass of heavy elements one would expect in a protoplanetary disc similar to the Minimum Mass Solar Nebula (MMSN) but around TOI178. Assuming a disc mass equal to 1% of the stellar mass, and using the metallicity of TOI178 ([Fe/H] = −0.23, corresponding to a metal content of Z ~ 0.0061), the mass of heavy elements in such a disc would equal ~13 M_{⊕}, 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 TOI178, 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 TOI178 system would point towards a formation channel similar to the one envisioned for the Trappist1 system (Schoonenberg et al. 2019).
7.2 Massradius relation
The six planets we detected in the TOI178 system are in the superEarth to miniNeptune range, with radii ranging from to 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 and , 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 .
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: Trappist1, K2138, Kepler60, Kepler80, and Kepler223. The diversity of planetary composition in TOI178 is clearly visible in the diagram: The two inner planets have a radius compatible with a gasfree structure, whereas the others contain water and/or gas. This is similar to the Kepler80 system, where the two innermost planets are compatible with a gasfree structure and the two outermost ones likely contain gas. Planets in Kepler223, Kepler60, and Trappist1 seem to have a more homogeneous structure: All planets in Kepler223 have a gas envelope, and the planets in Kepler60 and Trappist1 have small gas envelopes (below ~.01% of the total mass).
Considering the TOI178 system in greater detail, planets d, e, and g are located above the pure water line and definitely contain a nonzero gas mass fraction. Planet d and, depending on its mass, planet g, must contain a large gas fraction.
Figure 16 is a diagram showing the location of the TOI178 planets in stellar insolation versus the planetary radius. The colour code illustrates the density of exoplanets. This diagrams clearly shows the socalled 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.
Fig. 14 Example of TTVs of the six planets of TOI178, starting from the TESS observation of Sector 2 and spanning 6 yr, using the masses and orbital parameters from Tables 3 and 4. 
7.3 Comparison with other systems in Laplace resonance
The difference between the TOI178 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 Kepler60, Kepler80, and Kepler223, 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 K2138 are also 1σ consistent with such behaviour. In Trappist1, the density of planets is always higher than 4 g cm^{−3} and increases (with the exception of Trappist1 f) with decreasing equilibrium temperature. This is likely due to the presence of more ices in planets far from their stars (Agol et al. 2021). Contrary to the three Kepler systems, in the TOI178 system the density of the planets is not a growing function of the equilibrium temperature. Indeed, TOI178f has a density higher than that of planet e, and TOI178d 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 TOI178 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.
7.4 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, 2017) and has already been used in Mortier et al. (2020) and Delrez et al. (2021). 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 ironsulfur 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 150 000 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 semiamplitude 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 semiamplitude, which depend on the stellar radius and mass. In order to take this correlation into account, it is therefore important to fit the planetary system at once, 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 they add up to one) for the solid part. We also assumed: the mass fraction of water to be smaller than 50% (Thiabaud et al. 2014; 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 TOI178 with other systems with a Laplace resonance, we also performed a Bayesian analysis for the K2138, Kepler60, Kepler80, and Kepler223 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 Trappist1 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. 2021).
The posteriors distributions of the two most important parameters (mass fractions and composition of the mantle) of each planet in TOI178 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 Kepler223d, which is larger than that of planet e in the same system. Kepler223d 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 K2138, the gas fraction is compatible (within 1sigma error bars) with a monotonous function of the equilibrium temperature. In the case of TOI178, 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_{⊕} 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.
Fig. 15 TOI178 planets compared to other known transiting exoplanets with radius and mass uncertainties less than 40% (grey) and other systems known to harbour a Laplace resonance. Data on known exoplanets were taken from the NASA Exoplanet Archive on 18 September 2020. The dashed lines show theoretical massradius curves for some idealised compositions (Zeng et al. 2019). The six planets orbiting around TOI178 are indicated; the colour of the points and error bars give the equilibrium temperature. The seven planets orbiting Trappist1 are shown with diamonds, and the parameters are taken from Agol et al. (2021). The three planets orbiting Kepler60 are shown with ‘X’ marks, and the parameters are taken from JontofHutter et al. (2016). The six planets orbiting K2138 are shown with squares, and the parameters are taken from Lopez et al. (2019). The four planets orbiting Kepler80 are shown with inverse triangles, and the parameters are taken from MacDonald et al. (2016). The four planets orbiting Kepler223 are shown with regular triangles, and the parameters are taken from Mills et al. (2016). 
Fig. 16 Diagram of the position of the TOI178 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 highdensity and/or highly populated regions. 
Fig. 17 Densities of planets in the TOI178, Trappist1, K2138, Kepler60, Kepler80, and Kepler223 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. 
8 Discussion and conclusions
In this study, we presented new observations of TOI178 by CHEOPS, ESPRESSO, NGTS, and SPECULOOS. Thanks to this followup 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 reattribute 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 followup observations. In total, we therefore announce six planets in the superEarth to miniNeptune 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 quasicoplanar, with projected mutual inclinations between the outer planets estimated at about 0.1 deg, a feature that is also visible in other systems with threebody resonances, such as Trappist1 (Agol et al. 2021) and the Galilean moons. Current ephemeris and mass estimations indicate that the TOI178 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 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 firstorder (q = 1) and secondorder (q = 2) MMRs. These periods result from the equations detailed in Appendix C. In the TOI178 system, as well as in similar systems in Laplace resonance, planet pairs are virtually all nearly firstorder MMRs. The most likely of the periods shown in Table 7 are therefore the firstorder 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 TOI178, 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 TOI178 allows for further characterisation of the system by photometric measurements, radial velocities, and transit spectroscopy. These measurements will be essential 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, hydrogendominated atmospheres through escape (Kubyshkina et al. 2018, 2019), whereas all the other planets may have retained part of their primordial gas envelope. In this respect, the different planets of the TOI178 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 EELT 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 TOI178, 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 TOI178 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 thehigh 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.
Followup transit observations should also unveil TTVs for all but the innermost planet of the TOI178 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 twoplanet 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 TOI178 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.
Fig. 18 Gas mass in the planets of the TOI178, K2138, Kepler60, Kepler80, and Kepler223 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. 
Acknowledgements
The authors acknowledge support from the Swiss NCCR PlanetS and the Swiss National Science Foundation. Y.A. and M.J.H. acknowledge the support of the Swiss National Fund under grant 200020_172746. A.C.C. and T.W. acknowledge support from STFC consolidated grant number ST/M001296/1. This work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR10EQPX2901) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche. SH acknowledges CNES funding through the grant 837319. Based on data collected under the NGTS project at the ESO La Silla Paranal Observatory. The NGTS facility is operated by the consortium institutes with support from the UK Science and Technology Facilities Council (STFC) project ST/M001962/1. The Belgian participation to CHEOPS has been supported by the Belgian Federal Science Policy Office (BELSPO) in the framework of the PRODEX Program, and by the University of Liège through an ARC grant for Concerted Research Actions financed by the WalloniaBrussels Federation. V.A. acknowledges the support from FCT through Investigador FCT contract nr. IF/00650/2015/CP1273/CT0001. We acknowledge support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grants ESP201680435C21R, ESP201680435C22R, PGC2018098153BC33, PGC2018098153BC31, ESP201787676C51R, MDM20170737 Unidad de Excelencia “María de Maeztu” Centro de Astrobiología (INTACSIC), as well as the support of the Generalitat de Catalunya/CERCA programme. The MOC activities have been supported by the ESA contract No. 4000124370. S.C.C.B. acknowledges support from FCT through FCT contracts nr. IF/01312/2014/CP1215/CT0004. X.B., S.C., D.G., M.F. and J.L. acknowledge their role as ESAappointed CHEOPS science team members. ABr was supported by the SNSA. A.C. acknowledges support by CFisUC projects (UIDB/04564/2020 and UIDP/04564/2020), GRAVITY (PTDC/FISAST/7002/2020), ENGAGE SKA (POCI010145FEDER022217), and PHOBOS (POCI010145FEDER029932), funded by COMPETE 2020 and FCT, Portugal. This work was supported by FCT  Fundaçãopara a Ciência e a Tecnologia through national funds and by FEDER through COMPETE2020  Programa OperacionalCompetitividade e Internacionalização by these grants: UID/FIS/04434/2019; UIDB/04434/2020; UIDP/04434/2020; PTDC/FISAST/32113/2017 and POCI010145FEDER 032113; PTDC/FISAST/28953/2017 and POCI010145FEDER028953; PTDC/FISAST/28987/2017 and POCI010145FEDER028987. O.D.S.D. is supported in the form of work contract (DL 57/2016/CP1364/CT0004) funded by national funds through FCT. B.O.D. acknowledges support from the Swiss National Science Foundation (PP00P2190080). M.F. and C.M.P. gratefully acknowledgethe support of the Swedish National Space Agency (DNR 65/19, 174/18). D.G. gratefully acknowledges financial support from the CRT foundation under Grant No. 2018.2323 “Gaseousor rocky? Unveiling the nature of small worlds”. E.G. gratefully acknowledges support from the David and Claudia Harding Foundation in the form of a WintonExoplanet Fellowship. M.G. is an F.R.S.FNRS Senior Research Associate. J.I.G.H. acknowledges financial support from Spanish Ministry of Science and Innovation (MICINN) under the 2013 Ramón y Cajal programme RYC201314875. J.I.G.H., A.S.M., R.R., and C.A.P. acknowledge financial support from the Spanish MICINN AYA201786389P. A.S.M. acknowledges financial support from the Spanish Ministry of Science and Innovation (MICINN) under the 2019 Juan de la Cierva Programme. MNG ackowledges support from the MIT Kavli Institute as a Juan Carlos Torres Fellow. J.H. acknowledges the support of the Swiss National Fund under grant 200020_172746. KGI is the ESA CHEOPS Project Scientist and is responsible for the ESA CHEOPS Guest Observers Programme.She does not participate in, or contribute to, the definition of the Guaranteed Time Programme of the CHEOPS mission through which observations described in this paper have been taken, nor to any aspect of target selection forthe programme. J.S.J. acknowledges support by FONDECYT grant 1201371, and partial support from CONICYT project Basal AFB170002. A.J. acknowledges support from ANID  Millennium Science Initiative  ICN12_009 and from FONDECYT project 1171208. P.M. acknowledges support from STFC research grant number ST/M001040/1. N.J.N is supported by the contract and exploratory project IF/00852/2015, and projects UID/FIS/04434/2019, PTDC/FISOUT/29048/2017, COMPETE2020: POCI010145FEDER028987 & FCT: PTDC/FISAST/28987/2017. N.J.N is supported by the contract and exploratory project IF/00852/2015, and project PTDC/FISOUT/29048/2017. This work was also partially supported by a grant from the Simons Foundation (PI Queloz, grant number 327127). Acknowledges support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grant PGC2018098153B C33, as well as the support of the Generalitat de Catalunya/CERCA programme. S.G.S. acknowledges support from FCT through FCT contract nr. CEECIND/00826/2018 and POPH/FSE (EC). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This project has been supported by the Hungarian National Research, Development and Innovation Office (NKFIH) grants GINOP2.3.215201600003, K119517, K125015, and the City of Szombathely under Agreement No. 67.17721/2016. This research received funding from the MERAC foundation, from the European Research Council under the European Union’s Horizon 2020 research and innovation programme (grant agreement no 803193/ BEBOP, and from the Science and Technology Facilites Council (STFC, grant no ST/S00193X/1). V.V.G. is an F.R.SFNRS Research Associate. J.I.V. acknowledges support of CONICYTPFCHA/Doctorado Nacional21191829. M. R. Z. O. acknowledges financial support from projects AYA201679425C32P and PID2019109522GBC51 from the Spanish Ministry of Science, Innovation and Universities.
Appendix A Inspection of the CHEOPS data
The four CHEOPS visits were automatically processed through the DRP with individual frames undergoing various calibrations and corrections; aperture photometry was subsequently conducted for four aperture radii, as highlighted in Sect. 4.1.2 and detailed in Hoyer et al. (2020). The light curves produced for all runs in this study – which are often referred to as ‘raw’ in order to indicate no postprocessing detrending has taken place and which were obtained with the DEFAULT aperture – are shown in Fig. A.1. For the first, third, and fourth visits, standard data processing within the DRP was performed; however, for the second run, careful treatment of telegraphic pixels was needed.
Fig. A.1 DRPproduced light curves (DEFAULT aperture) of the four CHEOPS visits to TOI178 presented in this work. Threeσ outliers have been removed for better visualisation. 
In the CHEOPS CCD, there is a large number of hot pixels (see, for example, Fig. 3). Moreover, the behaviour ofsome normal pixels can change to an abnormal state within the duration of a visit. For example, a pixel can become ‘hot’ after an SAA crossing of the satellite. These pixels are called telegraphic due to their unstable response during the observations, and they can disturb the photometry if located inside the photometric aperture. To rule out the possibility that the detected transit events in the light curves correspond to the effect of telegraphic pixels, the data frames were carefully inspected and compared to the detection map of hot pixels delivered by the CHEOPS DRP (see details in Hoyer et al. 2020). By doing this, one telegraphic pixel was detected inside the DEFAULT aperture at the end of the second TOI178 visit. The exact CCD location of this abnormal pixel is shown in Fig. 3. The effect of this pixel in the photometry is shown in the form of a jump in flux in the light curve of the visit (top panel, Fig. A.2) at BJD ~2 459 076.5, which corresponds to the flux increase of the pixel (middle panel, Fig. A.2). After correcting the data by simply cancelling the flux of this pixel through the full observation and repeating the photometry extraction with the same aperture (R = 25′′), we removed the flux jump in the light curve (bottom panel, Fig. A.2). No telegraphic pixels were detected in the other TOI178 visits.
Fig. A.2 TOI178 normalised light curve of the second visit (grey symbols) shown in the top panel with its 10 min smoothed version overplotted (blue symbols). The flux jumps produced by the appearance of a new hot pixel are marked with the dashed vertical lines. The light curve of the detected telegraphic pixel is presented in the middle panel, showing its anomalous behaviour at the end of the visit. The light curve extracted from the corrected data is shown in the bottom panel. For better visualisation, the light curves are presented after a 3σ clipping andcorrected by a secondorder polynomial in time. 
Appendix B Analysis of the radial velocity data
B.1 Method
In this appendix, we describe the analysis of the RV data alone. To search for planet detections, we computed the ℓ_{1} periodogram of the RV, as defined in Hara et al. (2017). This tool is based on a sparse recovery technique called the basis pursuit algorithm (Chen et al. 1998). It aims to find a representation of the RV time series as a sum of a small number of sinusoids whose frequencies are in the input grid. The ℓ_{1} periodogram has the advantage, over a regular periodogram, of searching for several periodic components at the same time, therefore drastically reducing the impact of aliasing (Hara et al. 2017).
The ℓ_{1}periodogram has three inputs: a frequency grid, onto which the signal will be decomposed in the Fourier domain; a noise model in the form of a covariance matrix; and a base model. The base model represents offsets, trends, or activity models and can be understood as follows. The principle behind the ℓ_{1}periodogram is to consider the signal in the Fourier domain and to minimise the sum of the absolute value of the Fourier coefficients on a discretised frequency grid (their ℓ_{1} norm) whileensuring that the inverse Fourier transform is close enough to the model in a certain, precise sense. Due to the ℓ_{1} norm penalty, frequencies ‘compete’ against one another to have nonzero coefficients. However, one might assume that certain 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 ℓ_{1} periodogram 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.
Fig. B.1 Periodograms of the ancillary indicator time series (Hα, FWHM, bisector span, and shown in blue, orange, green, and red, respectively). Top: periodograms of the raw ancillary indicator time series. Bottom: periodograms of the same time series when the signal corresponding to the maximum peak is injected in the base model. 
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 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 36day and 16day 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 . 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 quasiperiodic term, as given by the formula (B.1)
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, Sindex, 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 Sindex 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^{−1}, τ = 0, 2, 4, 6 d, P_{act} = 36.5 d, σ_{QP} = 0,1,2,3,4 m s^{−1}, 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 FWHMtime 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.
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^{−1}, σ_{R} = 1.5 m s^{−1}, τ = 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 span, the frequency resolution is not good enough to distinguish between 36 and 45 days.
We find that signals at 36, 16, and 3.2 days are consistently included in the model. The 3.2day signal presents a median FAP of 0.002 and can therefore be confidently detected. When the base model consists of only an offset, 36 and 16 days are systematically significant with FAP < 5%, but their significance decreases when activity indicators are included in the model. Since these periodicities also appear in ancillary indicators, we conclude that they are due to activity. Signals at 6.5 and 9.9 days appear in ≈15% of the models, but with low significance. Although detections cannot be claimed from the RV data only, we note that there are signals at 2.08 d (alias of 1.91 d) and 9.9 d. We find that they are in phase with the photometric signals within 1.5σ, which further strengthens the detection of transiting planets at these periods. We finally note a small hint at 1.2 days (which can also appear at its alias, 5.6 days, in certain models).
The period of planet f is 15.2 days, and an activity signal in that region seems to be present (≈ 16 d, depending on ancillary indicators). The periods are such that 1∕(1∕15.23−1∕16) = 316 days, which is greater than the observation time span. To check whether our activity model allows us to detect planet f in the RV (and yield a meaningful mass estimate), we added to the base model all the known transiting planets, except the 15.2day one, as well as the smoothed Hα indicators and sine functions at 15.7 (closest to 15.2 days) and 36 days. From that we obtain Fig. B.3, where a signal at 15.2 days appears with an FAP of ≈20%.
The 20.7day transiting planet could correspond to the peak at 21.6 d in the ℓ_{1} periodogram (Fig. B.2, top). When restricting the frequency grid to 0 to 0.55 cycles per day, the best CV model yields Fig. 9, where a signal appears at 20.6 days. However, signals close to 20.7 days seem to disappear when the stellar activity model is changed. The planetary signature might be hidden in the RV due to stellar effects. Further observations would allow us to better disentangle stellar and planetary signals.
B.4 Detections: conclusion
In conclusion, we can claim an independent detection of the 3.2 d planet with RV. We find significant signatures at 36 and 15–16 days, which we attribute to activity. We find signals at 1.91 d (or its alias 2.08 d), 6.5 d, and 9.9 d in phase with the detected transits. Wesee a signal at 21 d for some models, which could correspond to the 20.7 d planet, but it seems that there are not enough points to disentangle a planetary signal from activity at this period. We note that modelling an activity signal at 15.7 days leaves a signal at 15.2 days, which likely stems from the presence of a planet at this period.
Finally, we checked the consistency in the phase of the signals fitted onto the photometric and RV data. We performed an MCMC computation of the orbital elements on the RVs with exactly the same priors as model 2, except that we set a flat prior on the phases. We find that the uncertainties on the phase corresponding to planets b, c, d, e, f, and g correspond respectively to 7, 2, 3.3, 4, and 100% of the period. The phases from transits are included, respectively, in the 1.5, 1, 2, 1, and 1 σ intervals derived from the RVs, such that we deem the phases derived from RVs consistent with the transits.
Fig. B.2 ℓ_{1}periodogram of the ESPRESSO RV corresponding to the best CV score with different linear base models for the data: Top: only one offset. Middle: offset, smoothed Hα. Bottom: smoothed Hα and smoothed FWHM. 
Fig. B.3 ℓ_{1}periodogram when all transiting planets but the 15.2d one are added to the base models, as well as a smoothed Hα indicator and sinusoids at 15.7 and 36 days. 
B.5 Mass and density estimates
The RV measurements allowed us to obtain mass estimates of the planets. These estimates can depend on the model of stellar 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 semiperiodic kernel. In both cases, we added the Hα time series smoothed with a GP, as described in Appendix 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 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 quasiperiodic kernel k of the form (B.2)
We imposed a Gaussian truncated prior on with σ = 100m^{2} s^{−2}. We imposed a flat prior on ν between 2π∕50 and 2π∕30 rad day^{−1} and a loguniform 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.
Inclusion in the top 20% best models of different periodicities.
Fig. B.4 Histograms of the weighted residuals of the maximum likelihood model (top) and their variogram (bottom) for activity modelled as two sinusoidal functions (model 1). 
We ran an adaptive MCMC algorithm as described in Delisle et al. (2018), implementing spleaf (Delisle et al. 2020)^{10}, 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 (Delisle et al. 2018). The posterior medians as well as the 1 σ credible intervals arereported in Table B.2 for both activity models. The kernel (B.2) is close to the stochastic harmonic oscillator (SHO), as defined in ForemanMackey 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.
Mass estimation with activity model as two sinusoids, Hα model, and correlated noise.
Modelling errors might leave a trace in the residuals. In Hara et al. (2019), it is shown that if the model used for theanalysis 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 ShapiroWilk normality test (Shapiro & Wilk 1965) and found a pvalue 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 nonnormality 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
TOI178 is in a configuration where successive pairs of planets are at the same distance to the exact neighbouring firstorder 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 resonance (k_{2} + q_{2}) : k_{2}, where k_{i} and q_{i} are integers. We hence write the nearresonant angles: (C.1)
where λ_{i} is the mean longitude of planet i. The associated distances to the resonances read: (C.2)
where n_{i} is the mean motion of planet i. A Laplace relation exists between these three planets if: (C.3)
where j_{1} and j_{2} are integers. In addition, the invariance by rotation of the Laplace angle (D’Alembert relation) gives: (C.4)
As a result, the Laplace relation requires: (C.5)
which translates, for the period of the third planet, as: (C.6)
where P_{1,2} is the superperiod associated with Δ_{1}. From this we can compute the periods of potential additional planets that would continue the Laplace resonant chain of TOI178. 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: (C.7)
where P_{f,g} is the superperiod between the near firstorder 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).
Fig. C.1 Superperiod of planet f with respect to e (in red) and g (in black) for a set of firstorder MMRs. The purple line indicates the value of the superperiod (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 nearintersection of the righthand side of a red curve, the lefthand side of a black curve, and the horizontal purple line. 
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 millionpoint) 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 TOI178 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 by the gas envelope onto the solid part. Given the mass of the gas envelope in all planets, this approximation is justified.
Fig. D.1 Histogram of the prediction error from our DNN on the test set. The yaxis is in log scale, and the xaxis covers an error from −1% to 1%. 
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.
Fig. D.2 Corner plot showing the main parameters of the internal structure of planet b. The parameters are: the core mass fraction; the mantle mass fraction; the water mass fraction (all relative to the solid planet); the molar fractions of Fe, Si, and Mg in the mantle; the molar fraction of Fe inthe core; and the mass of gas (log scale). The dashed lines give the positions of the 16 and 84% quantiles, and the number at the top of each column gives the median and the 16 and 84% quantiles. 
References
 Agol, E., Dorn, C., Grimm, S. L., et al. 2021, Planet. Sci. J., 2, 1 [Google Scholar]
 Alibert, Y., & Venturini, J. 2019, A&A, 626, A21 [CrossRef] [EDP Sciences] [Google Scholar]
 Alibert, Y., Carron, F., Fortier, A., et al. 2013, A&A, 558, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ambikasaran, S., ForemanMackey, D., Greengard, L., Hogg, D. W., & O’Neil, M. 2015, IEEE Transactions on Pattern Analysis and Machine Intelligence, 38, 2 [Google Scholar]
 Baluev, R. V. 2008, MNRAS, 385, 1279 [NASA ADS] [CrossRef] [Google Scholar]
 Barbary, K. 2016, J. Open Source Softw., 1, 58 [Google Scholar]
 Benz, W., Broeg, C., Fortier, A., et al. 2021, Exp. Astron, 51, 109 [Google Scholar]
 Bertin, E., & Arnouts, S. 1996, A&AS, 117, 393 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blackwell, D. E., & Shallis, M. J. 1977, MNRAS, 180, 177 [NASA ADS] [CrossRef] [Google Scholar]
 Bonfanti, A., Delrez, L., Hooton, M. J., et al. 2021, A&A, 646, A157 [EDP Sciences] [Google Scholar]
 Bruntt, H., Bedding, T. R., Quirion, P. O., et al. 2010, MNRAS, 405, 1907 [Google Scholar]
 Bryant, E. M., Bayliss, D., McCormac, J., et al. 2020, MNRAS, 494, 5872 [Google Scholar]
 Buchner, J. 2016, Stat. Comput., 26, 383 [CrossRef] [Google Scholar]
 Buchner, J. 2017, ArXiv eprints [arXiv:1707.04476] [Google Scholar]
 Burdanov, A., Delrez, L., Gillon, M., & Jehin, E. 2018, in Handbook of Exoplanets, (Berlin: Springer International Publishing AG), 130, 130 [Google Scholar]
 Cabrera, J., Csizmadia, S., Erikson, A., Rauer, H., & Kirste, S. 2012, A&A, 548, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Castelli, F., & Kurucz, R. L. 2003, IAU Symp., 210, A20 [Google Scholar]
 Chen, S. S., Donoho, D. L., & Saunders, M. A. 1998, SIAM J. Sci. Comput., 20, 33 [Google Scholar]
 Christiansen, J. L., Crossfield, I. J. M., Barentsen, G., et al. 2018, AJ, 155, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Claret, A., & Bloemen, S. 2011, A&A, 529, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Correia, A. C. M., Delisle, J.B., & Laskar, J. 2018, Planets in MeanMotion Resonances and the System Around HD45364, eds. H. J. Deeg, & J. A. Belmonte (USA: Nasa), 12 [Google Scholar]
 da CostaLuis, C., L., S., Mary, H., et al. 2018, https://doi.org/10.5281/zenodo.1468033 [Google Scholar]
 Dawson, R. I., & Fabrycky, D. C. 2010, ApJ, 722, 937 [NASA ADS] [CrossRef] [Google Scholar]
 Delisle, J. B. 2017, A&A, 605, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delisle, J.B., Laskar, J., Correia, A. C. M., & Boué, G. 2012, A&A, 546, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delisle, J. B., Ségransan, D., Dumusque, X., et al. 2018, A&A, 614, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delisle, J. B., Hara, N., & Ségransan, D. 2020, A&A, 638, A95 [EDP Sciences] [Google Scholar]
 Delrez, L., Gillon, M., Queloz, D., et al. 2018, Proc. SPIE 10700, 107001I [Google Scholar]
 Delrez, L., Ehrenreich, D., Alibert, Y., et al. 2021, Nat. Astron., submitted [Google Scholar]
 Dorn, C., Khan, A., Heng, K., et al. 2015, A&A, 577, A83 [EDP Sciences] [Google Scholar]
 Dorn, C., Venturini, J., Khan, A., et al. 2017, A&A, 597, A37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Emsenhuber, A., Mordasini, C., Burn, R., et al. 2020, A&A, submitted [arXiv:2007.05561] [Google Scholar]
 Espinoza, N., Kossakowski, D., & Brahm, R. 2019, MNRAS, 490, 2262 [NASA ADS] [CrossRef] [Google Scholar]
 Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ, 790, 146 [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
 ForemanMackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220 [Google Scholar]
 Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, AJ, 154, 109 [Google Scholar]
 Futyan, D., Fortier, A., Beck, M., et al. 2020, A&A, 635, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gillon, M. 2018, Nat. Astron., 2, 344 [NASA ADS] [CrossRef] [Google Scholar]
 Gillon, M., Triaud, A. H. M. J., Demory, B.O., et al. 2017, Nature, 542, 456 [Google Scholar]
 Goździewski,K., Migaszewski, C., Panichi, F., & Szuszkiewicz, E. 2016, MNRAS, 455, L104 [Google Scholar]
 Gupta, A., & Schlichting, H. E. 2019, MNRAS, 487, 24 [Google Scholar]
 Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hakim, K., Rivoldini, A., Van Hoolst, T., et al. 2018, Icarus, 313, 61 [NASA ADS] [CrossRef] [Google Scholar]
 Haldemann, J., Alibert, Y., Mordasini, C., & Benz, W. 2020, A&A, 643, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hara, N. C., Boué, G., Laskar, J., & Correia, A. C. M. 2017, MNRAS, 464, 1220 [NASA ADS] [CrossRef] [Google Scholar]
 Hara, N. C., Boué, G., Laskar, J., Delisle, J. B., & Unger, N. 2019, MNRAS, 489, 738 [NASA ADS] [CrossRef] [Google Scholar]
 Hara, N. C., Bouchy, F., Stalport, M., et al. 2020, A&A, 636, L6 [CrossRef] [EDP Sciences] [Google Scholar]
 Hayashi, C. 1981, Prog. Theore. Phys. Suppl., 70, 35 [Google Scholar]
 Haywood, R. D., Collier Cameron, A., Queloz, D., et al. 2014, MNRAS, 443, 2517 [NASA ADS] [CrossRef] [Google Scholar]
 Henrard, J., & Lemaitre, A. 1983, Celest. Mech., 30, 197 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Higson, E., Handley, W., Hobson, M., & Lasenby, A. 2019, Stat. Comput., 29, 891 [CrossRef] [Google Scholar]
 Høg, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27 [Google Scholar]
 Hoyer, S., Guterman, P., Demangeon, O., et al. 2020, A&A, 635, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Irwin, M. J., Lewis, J., Hodgkin, S., et al. 2004, Proc. SPIE, 5493, 411 [NASA ADS] [CrossRef] [Google Scholar]
 Izidoro, A., Ogihara, M., Raymond, S. N., et al. 2017, MNRAS, 470, 1750 [Google Scholar]
 Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, SPIE Conf. Ser., 9913, 99133E [Google Scholar]
 JontofHutter, D., Ford, E. B., Rowe, J. F., et al. 2016, ApJ, 820, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D. M. 2013, MNRAS, 435, 2152 [Google Scholar]
 Kovács, G., Zucker, S., & Mazeh, T. 2002, A&A, 391, 369 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kreidberg, L. 2015, PASP, 127, 1161 [NASA ADS] [CrossRef] [Google Scholar]
 Kubyshkina, D., Fossati, L., Erkaev, N. V., et al. 2018, A&A, 619, A151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kubyshkina, D., Cubillos, P. E., Fossati, L., et al. 2019, ApJ, 879, 26 [Google Scholar]
 Kurucz, R. L. 2013, Astrophysics Source Code Library [record ascl:1303.024] [Google Scholar]
 Laskar, J. 1990, Icarus, 88, 266 [Google Scholar]
 Laskar, J. 1993, Phys. D Nonlinear Phen., 67, 257 [Google Scholar]
 Lee, M. H., & Peale, S. J. 2002, ApJ, 567, 596 [Google Scholar]
 Leleu, A., Robutel, P., & Correia, A. C. M. 2015, A&A, 581, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leleu, A., LilloBox, J., Sestovic, M., et al. 2019, A&A, 624, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lendl, M., Csizmadia, S., Deline, A., et al. 2020, A&A, 643, A94 [EDP Sciences] [Google Scholar]
 Lissauer, J. J., Fabrycky, D. C., Ford, E. B., et al. 2011, Nature, 470, 53 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Lithwick, Y., Xie, J., & Wu, Y. 2012, ApJ, 761, 122 [Google Scholar]
 Lopez, E. D., & Fortney, J. J. 2014, ApJ, 792, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Lopez, T. A., Barros, S. C. C., Santerne, A., et al. 2019, A&A, 631, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Luger, R., Sestovic, M., Kruse, E., et al. 2017, Nat. Astron., 1, 0129 [Google Scholar]
 MacDonald, M. G., Ragozzine, D., Fabrycky, D. C., et al. 2016, AJ, 152, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Marboeuf, U., Thiabaud, A., Alibert, Y., Cabral, N., & Benz, W. 2014, A&A, 570, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marigo, P., Girardi, L., Bressan, A., et al. 2017, ApJ, 835, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Maxted, P. F. L. 2016, A&A, 591, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mayor, M., & Queloz, D. 1995, Nature, 378, 355 [Google Scholar]
 McCormac, J., Pollacco, D., Skillen, I., et al. 2013, PASP, 125, 548 [NASA ADS] [CrossRef] [Google Scholar]
 Mills, S. M., Fabrycky, D. C., Migaszewski, C., et al. 2016, Nature, 533, 509 [Google Scholar]
 Morbidelli, A. 2002, Modern Celestial Mechanics : Aspects of Solar System Dynamics (London: Taylor & Francis) [Google Scholar]
 Mordasini, C., Alibert, Y., Benz, W., Klahr, H., & Henning, T. 2012, A&A, 541, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mortier, A., Zapatero Osorio, M. R., Malavolta, L., et al. 2020, MNRAS, 499, 5004 [Google Scholar]
 Murray, C. A., Delrez, L., Pedersen, P. P., et al. 2020, MNRAS, 495, 2446 [CrossRef] [Google Scholar]
 Newville, M., Stensitzki, T., Allen, D. B., & Ingargiola, A. 2014, https://zenodo.org/record/11813#.YC5nHXnjKUk [Google Scholar]
 Noyes, R. W. 1984, in Space Research in Stellar Activity and Variability, eds. A. Mangeney, & F. Praderie (Paris: Meudon, Observatoire), 113 [Google Scholar]
 Papaloizou, J. C. B. 2015, Int. J. Astrobiol., 14, 291 [CrossRef] [Google Scholar]
 Papaloizou, J. C. B., & Terquem, C. 2010, MNRAS, 405, 573 [NASA ADS] [Google Scholar]
 Pepe, F., Cristiani, S., Rebolo, R., et al. 2021, A&A, 645, A96 [CrossRef] [EDP Sciences] [Google Scholar]
 Petit, A. C., Laskar, J., & Boué, G. 2018, A&A, 617, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Piskunov, N., & Valenti, J. A. 2017, A&A, 597, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pu, B., & Wu, Y. 2015, ApJ, 807, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Queloz, D., Henry, G. W., Sivan, J. P., et al. 2001, A&A, 379, 279 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rein, H., & Liu, S. F. 2012, A&A, 537, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rein, H., & Spiegel, D. S. 2015, MNRAS, 446, 1424 [NASA ADS] [CrossRef] [Google Scholar]
 Rivera, E. J., Laughlin, G., Butler, R. P., et al. 2010, ApJ, 719, 890 [Google Scholar]
 Samuel, B., Leconte, J., Rouan, D., et al. 2014, A&A, 563, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schoonenberg, D., Liu, B., Ormel, C. W., & Dorn, C. 2019, A&A, 627, A149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Scuflaire, R., Théado, S., Montalbán, J., et al. 2008, Ap&SS, 316, 83 [Google Scholar]
 Shapiro, S. S., & Wilk, M. B. 1965, Biometrika, 52, 591 [Google Scholar]
 Skilling, J. 2004, AIP Conf. Ser., 735, 395 [Google Scholar]
 Skilling, J. 2006, Bayesian Anal., 1, 833 [Google Scholar]
 Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163 [Google Scholar]
 Smith, J. C., Stumpe, M. C., Van Cleve, J. E., et al. 2012, PASP, 124, 1000 [Google Scholar]
 Smith, A. M. S., Eigmüller, P., Gurumoorthy, R., et al. 2020, Astron. Nachr., 341, 273 [Google Scholar]
 Sneden, C. 1973, PhD thesis, University of Texas, Austin, USA [Google Scholar]
 Sotin, C., Grasset, O., & Mocquet, A. 2007, Icarus, 191, 337 [Google Scholar]
 Sousa, S. G. 2014, ARES + MOOG: A Practical Overview of an Equivalent Width (EW) Method to Derive Stellar Parameters (USA: Nasa), 297 [Google Scholar]
 Sousa, S. G., Santos, N. C., Adibekyan, V., DelgadoMena, E., & Israelian, G. 2015, A&A, 577, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Speagle, J. S. 2020, MNRAS, 493, 3132 [Google Scholar]
 Stassun, K. G., Oelkers, R. J., Pepper, J., et al. 2018, AJ, 156, 102 [NASA ADS] [CrossRef] [Google Scholar]
 Stassun, K. G., Oelkers, R. J., Paegert, M., et al. 2019, AJ, 158, 138 [Google Scholar]
 Stumpe, M. C., Smith, J. C., Catanzarite, J. H., et al. 2014, PASP, 126, 100 [Google Scholar]
 Terquem, C., & Papaloizou, J. C. B. 2007, ApJ, 654, 1110 [NASA ADS] [CrossRef] [Google Scholar]
 Thiabaud, A., Marboeuf, U., Alibert, Y., et al. 2014, 562, A27 [Google Scholar]
 Valenti, J. A.,& Piskunov, N. 1996, A&AS, 118, 595 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Venturini, J., Guilera, O. M., Haldemann, J., Ronco, M. P., & Mordasini, C. 2020, A&A, 643, L1 [CrossRef] [EDP Sciences] [Google Scholar]
 Wheatley, P. J., West, R. G., Goad, M. R., et al. 2018, MNRAS, 475, 4476 [NASA ADS] [CrossRef] [Google Scholar]
 Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868 [Google Scholar]
 Yee, S. W., Petigura, E. A., & von Braun, K. 2017, ApJ, 836, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Zeng, L., Jacobsen, S. B., Sasselov, D. D., et al. 2019, Proc. Natl. Acad. Sci., 116, 9723 [Google Scholar]
Padova and Trieste Stellar Evolutionary Code http://stev.oapd.inaf.it/cgibin/cmd.
We note that if the presence of a valley seems robust, it could be due to effects that are not related to evaporation, e.g. core cooling (Gupta & Schlichting 2019), or from combined formation and evolution effects (Venturini et al. 2020).
All Tables
Instantaneous distances from resonances for the six planets of TOI178, expressed in terms of nearresonant angles ϕ.
Estimated values of the Laplace angles of planets c to g of TOI178 towards the beginning of TESS Sector 2 at 2 458 350.0 BJD (August 2018).
Mass estimation with activity model as two sinusoids, Hα model, and correlated noise.
All Figures
Fig. 1 Light curves from TESS Sector 2 described in Sects. 2 and 4.1.1. Unbinned data are shown as grey points, and data in 30 min bins are shown as black circles. The best fitting transit model for the system is shown in black; the associated parameter values are shown in Tables 3 and 4. The positions of the transits are marked with lines coloured according to the legend. The photometry before and after the midsector gap are shown in the top and bottom panels, respectively. As the first transit of planet f (thick teal line) occurred precisely between the two transits of the similarly sized planet g (period of 20.71 d  red lines), the three transits were originally thought to have arisen due to a single planet, which was originally designated TOI178.02 with a period of 10.35 d (see Sect. 2). 

In the text 
Fig. 2 Light curves from simultaneous observations of TOI178 by NGTS (green stars) and SPECULOOSSouth (purple triangles), described in Sects. 4.1.3 and 4.1.4, respectively. Unbinned data are shown as grey points,and data in 15minute bins are shown as green stars (NGTS) and purple triangles (SPECULOOS). The observationsoccurred on September 11, 2019 (top panel) and October 12, 2019 (bottom panel). The transit model is shown in black. The position of the odd transit of candidate TOI178.02 is shown with a dashed red line, the transit of planet g (which corresponds to an even transit of TOI178.02) is shown with the solid red line, and the transits of planet b are shown with purple lines. 

In the text 
Fig. 3 Extraction of 80 × 80 arcsec of the CHEOPS fieldofview for two different data frames at the beginning (top) and the end (bottom) of the second visit (see Table 2). The TOI178 PSF is shown at the centre, with the DEFAULT DRP photometricaperture 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. 

In the text 
Fig. 4 Similar to Fig. 1, but instead displaying data from the four CHEOPS visits described in Sect. 4.1.2. Top panel: 11 d observation. Transits of planets c and d at ~ 2459075 BJD occur too close to each other for their corresponding lines to be individually visible in the figure. Bottomleft panel: subsequent observation scheduled to confirm the presence of a planet with a 20.7 d period (planet g), which overlaps with a transit of planet e. Bottomright panel: final observation scheduled to confirm the presence of a planet fitting in the Laplace resonance (planet f, with a period of ~ 15.23 d), which overlaps with a transit of planet b. 

In the text 
Fig. 5 ESPRESSO RV data of TOI178. 

In the text 
Fig. 6 Successive use of the BLS algorithm to identify the new candidates. Top panel: green curve has its highest power at 9.96 d once the three highest peaks (i = 3) at 18.35 d (black), 10.37 d (blue), and 9.17 d (orange) are ignored. Middle panel: BLS after the removal of the 9.96 d signal in the light curve, with no peaks ignored (j = 0). Bottom panel: BLS after the removal of the 20.71 d signal as well. 

In the text 
Fig. 7 Detrended CHEOPS light curves phasefolded to the periods of each of the planets, with signals of the other planets removed. Unbinned data are shown with grey points, data in 15min bins are shown with coloured circles, and samples drawn from the posterior distribution of the global fit are shown with coloured lines. 

In the text 
Fig. 8 Impact parameter (top) and inclination (bottom) of the six planets of the TOI178 system. In the top panel, the solid line shows the evolution of the impact parameter as a function of the period assuming that all planets are in the same plane (obtained by linear fit). The dashed line shows that an outer planet in the same plane will not transit if its period is above 26.4 d. 

In the text 
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. 

In the text 
Fig. 10 Phasefolded RV. Error bars correspond to nominal errors. 

In the text 
Fig. 11 Example of the evolution of the Laplace angles over 100 yr, starting from TESS observation of Sector 2, using the masses and orbital parameters from Tables 3 and 4. 

In the text 
Fig. 12 Stability indicator (defined in Eq. (3)) for TOI178 as a function of the periods and eccentricity of planet f (top). The red areas show the initial conditions of unstable trajectories, while black shows stable (quasiperiodic) trajectories. Bottom panel: same stability criterion with respect to the initial periods of planets f and g. The dashed white lines show the observed periods reported in Table 4. 

In the text 
Fig. 13 Same stability indicator as in Fig. 12 (top) and the number of librating Laplace angles (bottom) for TOI178 as a function of the period and mass of planet f. The dashed white line shows the observed period of planet f reported in Table 4. 

In the text 
Fig. 14 Example of TTVs of the six planets of TOI178, starting from the TESS observation of Sector 2 and spanning 6 yr, using the masses and orbital parameters from Tables 3 and 4. 

In the text 
Fig. 15 TOI178 planets compared to other known transiting exoplanets with radius and mass uncertainties less than 40% (grey) and other systems known to harbour a Laplace resonance. Data on known exoplanets were taken from the NASA Exoplanet Archive on 18 September 2020. The dashed lines show theoretical massradius curves for some idealised compositions (Zeng et al. 2019). The six planets orbiting around TOI178 are indicated; the colour of the points and error bars give the equilibrium temperature. The seven planets orbiting Trappist1 are shown with diamonds, and the parameters are taken from Agol et al. (2021). The three planets orbiting Kepler60 are shown with ‘X’ marks, and the parameters are taken from JontofHutter et al. (2016). The six planets orbiting K2138 are shown with squares, and the parameters are taken from Lopez et al. (2019). The four planets orbiting Kepler80 are shown with inverse triangles, and the parameters are taken from MacDonald et al. (2016). The four planets orbiting Kepler223 are shown with regular triangles, and the parameters are taken from Mills et al. (2016). 

In the text 
Fig. 16 Diagram of the position of the TOI178 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 highdensity and/or highly populated regions. 

In the text 
Fig. 17 Densities of planets in the TOI178, Trappist1, K2138, Kepler60, Kepler80, and Kepler223 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. 

In the text 
Fig. 18 Gas mass in the planets of the TOI178, K2138, Kepler60, Kepler80, and Kepler223 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. 

In the text 
Fig. A.1 DRPproduced light curves (DEFAULT aperture) of the four CHEOPS visits to TOI178 presented in this work. Threeσ outliers have been removed for better visualisation. 

In the text 
Fig. A.2 TOI178 normalised light curve of the second visit (grey symbols) shown in the top panel with its 10 min smoothed version overplotted (blue symbols). The flux jumps produced by the appearance of a new hot pixel are marked with the dashed vertical lines. The light curve of the detected telegraphic pixel is presented in the middle panel, showing its anomalous behaviour at the end of the visit. The light curve extracted from the corrected data is shown in the bottom panel. For better visualisation, the light curves are presented after a 3σ clipping andcorrected by a secondorder polynomial in time. 

In the text 
Fig. B.1 Periodograms of the ancillary indicator time series (Hα, FWHM, bisector span, and shown in blue, orange, green, and red, respectively). Top: periodograms of the raw ancillary indicator time series. Bottom: periodograms of the same time series when the signal corresponding to the maximum peak is injected in the base model. 

In the text 
Fig. B.2 ℓ_{1}periodogram of the ESPRESSO RV corresponding to the best CV score with different linear base models for the data: Top: only one offset. Middle: offset, smoothed Hα. Bottom: smoothed Hα and smoothed FWHM. 

In the text 
Fig. B.3 ℓ_{1}periodogram when all transiting planets but the 15.2d one are added to the base models, as well as a smoothed Hα indicator and sinusoids at 15.7 and 36 days. 

In the text 
Fig. B.4 Histograms of the weighted residuals of the maximum likelihood model (top) and their variogram (bottom) for activity modelled as two sinusoidal functions (model 1). 

In the text 
Fig. C.1 Superperiod of planet f with respect to e (in red) and g (in black) for a set of firstorder MMRs. The purple line indicates the value of the superperiod (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 nearintersection of the righthand side of a red curve, the lefthand side of a black curve, and the horizontal purple line. 

In the text 
Fig. D.1 Histogram of the prediction error from our DNN on the test set. The yaxis is in log scale, and the xaxis covers an error from −1% to 1%. 

In the text 
Fig. D.2 Corner plot showing the main parameters of the internal structure of planet b. The parameters are: the core mass fraction; the mantle mass fraction; the water mass fraction (all relative to the solid planet); the molar fractions of Fe, Si, and Mg in the mantle; the molar fraction of Fe inthe core; and the mass of gas (log scale). The dashed lines give the positions of the 16 and 84% quantiles, and the number at the top of each column gives the median and the 16 and 84% quantiles. 

In the text 
Fig. D.3 Same as Fig. D.2 but for planet c. 

In the text 
Fig. D.4 Same as Fig. D.2 but for planet d. 

In the text 
Fig. D.5 Same as Fig. D.2 but for planet e. 

In the text 
Fig. D.6 Same as Fig. D.2 but for planet f. 

In the text 
Fig. D.7 Same as Fig. D.2 but for planet g. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.