© A. Leleu et al. 2021

1 Introduction

Since the discovery of the first exoplanet orbiting a Sun-like star by Mayor & Queloz (1995), the diversity of observed planetary systems has continued to challenge our understanding of their formation and evolution. As an ongoing effort to understand these physical processes, observational facilities strive to get a full picture of exoplanetary systems by looking for 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 mean-motion resonances (MMRs) are ‘Rosetta Stones’ of the formation and evolution of planetary systems. Indeed, our current understanding of planetary system formation theory implies that such configurations are a common outcome of protoplanetary discs: Slow convergent migration of a pair of planets in quasi-circular orbits leads to a high probability of capture in first-order MMRs – the period ratio of the two planets is equal to (k + 1)∕k, with k an integer (Lee & Peale 2002; Correia et al. 2018). As the disc strongly damps the eccentricities of the protoplanets, this mechanism can 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 near-resonant 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 three-body resonances (Morbidelli 2002; Papaloizou 2015). Such systems are too fine-tuned to result from scattering events and hence can be used to constrain the outcome of protoplanetary discs (Mills et al. 2016).

A Laplace resonance, in reference to the configuration of Io, Europa, and Ganymede, is a three-body resonance where each consecutive pair of bodies is in, or close to, two planet MMRs. To date, only a few systems have been observed in a chain of Laplace resonances: GJ 876 (Rivera et al. 2010), Kepler-60 (Goździewski et al. 2016), Kepler-80 (MacDonald et al. 2016), Kepler-223 (Mills et al. 2016), Trappist-1 (Gillon et al. 2017; Luger et al. 2017), and K2-138 (Christiansen et al. 2018; Lopez et al. 2019). Out of these six systems, only K2-138 has so far been observed by both transit and radial velocity (RV), mainly due to the relative faintness of the other host stars in the visible (V magnitudes greater than ~14). However, TTVs could also be used to estimate the mass of their planets (see, for example, Agol et al. 2021 for the case of Trappist-1).

In this study, we present photometric and RV observations of TOI-178, a V = 11.95 mag, K-type star that was first observed by TESS¹ in its Sector 2. We jointly analyse the photometric data of TESS, two nights of NGTS² and SPECULOOS³ data, and 285 hours of CHEOPS⁴ observations, along with 46 ESPRESSO⁵ RV points. This follow-up effort allows us to decipher the architecture of the system and demonstrate the presence of a chain of Laplace resonances between the five outer planets.

We begin in Sect. 2 by presenting the rationale that led to the CHEOPS observation sequence (two visits totalling 11 d followed by two shorter visits). In Sect. 3, we describe the parameters of the star. In Sect. 4, we present the photometric and RV data we use in the paper. In Sect. 5, we show how these data led us to the detection of six planets – the outer five of which are in a chain of Laplace resonances – and to constrain their parameters. In Sect. 6, we explain the resonant state of the system, discuss its stability, and describe the transit timing variations (TTVs) that this system could potentially exhibit in the coming years. Finally, we discuss the internal structure of the planets in Sect. 7, and conclusions are presented in Sect. 8.

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 TOI-178 with periods of 6.55 d, 10.35 d, and 9.96 d. Based on these data, TOI-178 was identified as a potential co-orbital system (Leleu et al. 2019) with two planets oscillating around the same period. This prompted ESPRESSO RV measurements and two sequences of simultaneous ground-based photometricobservations with NGTS and SPECULOOS. From the latter, no transit was observed for TOI-178.02 (P = 10.35 d) in September 2019; however, a transit ascribed to this candidate was detected one month later by NGTS and SPECULOOS. The abovementioned absence of a TOI-178.02 transit combined with the three transits observed by TESS at high signal-to-noise ratio (S/N) (above 10) was interpreted as an additional sign of the strong TTVs expected in a co-orbital configuration. This solution was supported by RV data that were consistent with the horseshoe orbits of objects with similar masses (Leleu et al. 2015).

A continuous 11 d CHEOPS observation (split into two visits for scheduling reasons) was therefore performed in August 2020 in order to confirm the orbital configuration of the system; as the instantaneous period of both members of the co-orbital pair will always be smaller than 11 d, at least one transit of both targets should thus be detectable. Analysis of this light curve led to the confirmation of the presence of two of the planets already discovered by TESS (in this study denoted as planets d and e, with periods of 6.55 d and 9.96 d, respectively) and the detection of two new inner transiting planets (denoted planets b and c, with periods of 1.9 d and 3.2 d, respectively). However, one of the planets belonging to the proposed co-orbital pair (with a period of 10.35 d) was not apparent in the light curve. A potential hypothesis was that the first and third transit of the TOI-178.02 candidate (P = 10.35 d) during TESS Sector 2 belonged to a planet of twice the period (20.7 d), while the second transit belonged to another planet of unknown period. This scenario was supported by the two ground-based observations mentioned above since a P = 20.7 d planet would not transit during the September 2019 observation window. Using the ephemerides from fitting the TESS, NGTS, and SPECULOOS data, we predicted the mid-transit of the 20.7-day candidate to occur between UTC 14:06 and 14:23 on September 7, 2020 with 98% certainty. A third visit of CHEOPS observed the system around this epoch for 13.36 h and confirmed the presence of this new planet (g) at the predicted time with consistent transit parameters. Further analysis of all available photometric data found two possibilities for the unknown period of the additional planet: ~12.9 d or ~ 15.24 d.

Careful analysis of the whole system additionally revealed that planets c, d, e, and g were in a Laplace resonance. In order to fit in the resonant chain, the unknown period of the additional planet could have only two values: P = 13.4527 d or 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 ${15.231915}_{-0.000095}^{+0.000115}$ day.

As we will detail in the following sections, we can confirm the detection of a 6.56-day and a 9.96-day period planet by TESS using the new observations presented in this paper. Furthermore, we can announce the detection of planets with periods of 1.91, 3.24, 15.23, and 20.71 d (see Tables 3 and 4 for the complete parameters of the system).

3 TOI-178 stellar characterisation

Forty-six ESPRESSO observations (see Sect. 4.2) of TOI-178, a V = 11.95 mag K-dwarf, have been used to determine the stellar spectral parameters. These observations were first shifted and stacked to produce a combined spectrum. We then used the publicly available spectral analysis package SME (Spectroscopy Made Easy; Valenti & Piskunov 1996; Piskunov & Valenti 2017) version 5.22 to model the co-added ESPRESSO spectrum. We selected the ATLAS12 model atmosphere grids (Kurucz 2013) and atomic line data from VALD to compute synthetic spectra, which were fitted to the observations using a χ² -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 macro-turbulent velocity was modelled and found to be 1.2 ± 0.9 km s⁻¹, and micro-turbulent velocity was fixed to 0.91 km s⁻¹ following the formulation in Bruntt et al. (2010). The surface gravity log g was constrained from the line wings of the Ca I triplet (6102, 6122, and 6162 Å) and the Ca I 6439 Å line with a fixed T_eff and Ca abundance.

We checked our model with the Na I doublet that is sensitive to both T_eff and log g, and, finally, we also tested the MARCS 2012 (Gustafsson et al. 2008) model atmosphere grids. The measured parameters are listed in Table 1. The SME results are in agreement with the empirical SpecMatch-Emp (Yee et al. 2017)code, which fits the stellar optical spectra to a spectral library of 404 M5 to F1 stars, resulting in T_eff = 4316 ± 70 K, log g = 4.45 ± 0.15, and [Fe∕H] = −0.29 ± 0.05 dex. We also used ARES+MOOG (Sousa 2014; Sousa et al. 2015; Sneden 1973) to conduct the spectroscopic analysis on the same combined ESPRESSO spectra, and, although we derived consistent parameters (T_eff = 4500 ± 230 K, log g = 4.38 ± 0.62, and [Fe∕H] = −0.34 ± 0.10), the large errors are indicative of the difficulties in using equivalent width methods with colder stars: Spectral lines are more crowded in the spectrum.

Using these precise spectral parameters as priors on stellar atmospheric model selection, we determined the radius of TOI-178 using the infrared flux method (IRFM; Blackwell & 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 zero-point 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 TOI-178 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⁶ v1.2S code (Marigo et al. 2017) and the CLES code (Code Liégeois d’Évolution Stellaire; Scuflaire et al. 2008), with the reported values representing a careful combination of results from both sets of models. This was done because the sets of models differ slightly in their approaches (reaction rates, opacity and overshoot treatment, and helium-to-metal enrichment ratio), and thus by comparing masses and ages derived from both grids it is possible to include systematic uncertainties within the modelling of the position of TOI-178 on evolutionary tracks and isochrones. A detailed discussion of combining the PARSEC and CLES models to determine masses and ages can be found in Bonfanti et al. (2021). For this study, we derived M_⋆ = 0.647${}_{-0.032}^{+0.035}{M}_{\odot }$ and t_⋆ = 7.1${}_{-5.4}^{+6.2}$ Gyr. All stellar parameters are shown in Table 1.

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

4.1.2 CHEOPS

CHEOPS, the first ESA small-class mission, is dedicated to the observation of bright stars (V ≲ 12 mag) that are known to host planets and performs ultra high precision photometry, with the precision limited by stellar photon noise of 150 ppm min⁻¹ for a V = 9 magnitude star. The CHEOPS instrument is composed of an f/8 Ritchey-Chretien on-axis telescope (~30 cm diameter) equipped with a single frame-transfer back-side illuminated charge-coupled device (CCD) detector. The satellite was successfully launched from Kourou (French Guiana) into a ~700 km altitude Sun-synchronous orbit on December 18, 2019. CHEOPS took its first image on February 7, 2020, and, after it passed the in-orbit commissioning (IOC) phase, routine operations started on March 25, 2020. More details on the mission can be found in Benz et al. (2021), and the first results have recently been presented in Lendl et al. (2020).

The versatility of the CHEOPS mission allows for space-based follow-up photometry of planetary systems identified by the TESS mission. This is particularly useful for completing the inventory of multi-planetary systems whose outer transiting planets have periods beyond ~10 days.

We obtained four observation runs (or visits) of TOI-178 with CHEOPS between 4 August 2020 and 3 October 2020 as part of the Guaranteed Time, totalling 11.88 days on target with the observing log shown in Table 2. The majority of this time was spent during the first two visits (with lengths of 99.78 h and 164.06 h), which were conducted sequentially beginning on August 4, 2020 and ending on August 15, 2020; as such, we achieved a near-continuous 11-day 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 low-Earth orbit of CHEOPS, the spacecraft-target line of sight was interrupted by Earth occultations and passages through the South Atlantic Anomaly (SAA), where no data were downlinked. This resulted in gaps in the photometry on CHEOPS orbit timescales (around 100 min). For our observations of TOI-178, this resulted in light curve efficiencies of 51, 54, 65, and 86%. For all four visits, we used an exposure time of 60 s.

These data were automatically processed with the latest version of the CHEOPS data reduction pipeline (DRP v12; Hoyer et al. 2020). This includes image calibration (bias, gain, non-linearity, dark current, and flat field) and instrumental and environmental corrections (cosmic rays and smearing trails of field stars and background). The DRP also performs aperture photometry for a set of three size-fixed apertures – R = 22.5′′ (RINF), 25′′ (DEFAULT), and 30′′ (RSUP) – plus one extra aperture that minimises the contamination from nearby field stars (ROPT), which, in the case of TOI-178, was set at R = 28.5′′. The DRP estimates the level of contamination by simulating the field-of-view of TOI-178 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 TOI-178, the mean contamination level was below 0.1% and was mostly modulated by the rotation around the target of a nearby star of Gaia G = 13.3 mag at a projected sky distance of 60.8′′ from the target. The contamination as a function of time (or, equivalently, as a function of the roll angle of the satellite) is provided as a product of the DRP for further de-trending (see Sect. 5.1). For the second visit, careful removal of one ‘telegraphic’ pixel (a pixel with a non-stable abnormal behaviour during the visit) within the photometric aperture was needed. The location of this telegraphic pixel is shown in red in Fig. 3, and details regarding the detection and correction are described in Appendix A. Following the reductions, we found that the light curves obtained using the 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, non-astrophysical short-term trends – either from a varying background, a nearby contaminating source, or other sources – can be found in the data on orbital timescales. These effects can be efficiently modelled by using a Gaussian process (GP) with roll angle as input (e.g. Lendl et al. 2020; Bonfanti et al. 2021), as we will discuss in our data analysis (Sect. 5.1.2). Following this, the average noise over a 3 h sliding window for the four visits of the G = 11.15 mag target was found to be 63.9, 64.2, 66.3, and 75.8 ppm. In all cases, this marginally improved upon the precision of the light curves that we had previously simulated for these observation windows using the CHEOPSim tool (Futyan et al. 2020).

4.1.3 NGTS

The NGTS (Wheatley et al. 2018) facility consists of twelve 20-cm diameter robotic telescopes and is situated at the ESO Paranal Observatory in Chile. The individual NGTS telescopes have a wide field-of-view of 8 square-degrees and a plate scale of 5′′ pixel⁻¹. The DONUTS auto-guiding algorithm (McCormac et al. 2013) affords the NGTS telescopes sub-pixel level guiding. Simultaneous observations using multiple NGTS telescopes have been shown to yield ultra high precision light curves of exoplanet transits (Bryant et al. 2020; Smith et al. 2020).

TOI-178 was observed using NGTS multi-telescope observations on two nights. On UT September 10, 2019, TOI-178 was observed using six NGTS telescopes during the predicted transit event of TESS candidate TOI-178.02. However, the NGTS data for this night rule out a transit occurring during the observations. A second predicted transit event of TOI-178.02 was observed on the night of UT October 11, 2019 using seven NGTS telescopes, and on this night the transit event was robustly detected by NGTS. A total of 13,991 images were obtained on the first night, and 12,854 were obtained on the second. For both nights, the images were taken using the custom NGTS filter (520–890 nm) with an exposure time of 10 s. All NGTS observations of TOI-178 were performed with airmass < 2 and under photometric sky conditions.

The NGTS images were reduced using a custom pipeline that uses the SEP library to perform source extraction and aperture photometry (Bertin & Arnouts 1996; Barbary 2016). A selection of comparison stars with brightnesses, colours, and CCD positions similar to those of TOI-178 were identified using the second Gaia data release (DR2; Gaia Collaboration 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 TOI-178 by NGTS (green stars) and SPECULOOS-South (purple triangles), described in Sects. 4.1.3 and 4.1.4, respectively. Unbinned data are shown as grey points,and data in 15-minute 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 TOI-178.02 is shown with a dashed red line, the transit of planet g (which corresponds to an even transit of TOI-178.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 1-m telescopes (Callisto, Europa, Ganymede, and Io). Each robotic SSO telescope has a primary aperture of 1 m and a focal length of 8 m, and is equipped with a 2k × 2k deep-depletion CCD camera whose 13.5 μm pixel size corresponds to 0.35″ on the sky (field-of-view = 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 de-focus was applied to avoid saturation as the target was too bright for SSO. Light curves were extracted using the SSO pipeline (Murray et al. 2020) and are shown in purple in Fig. 2. For each observing night, the SSO pipeline used the casutools software (Irwin et al. 2004) to perform automated differential photometry and to correct for systematics caused by time-varying telluric water vapour.

4.2 ESPRESSO data

The RV datawe analyse consist of 46 ESPRESSO points⁷. Each measurement was taken in high resolution (HR) mode with an integration time of 20 min using a single telescope (UT) and slow read-out (HR 21). The source on fibre B is the Fabry-Perot interferometer. Observations were made with a maximum airmass of 1.8 and a minimum 30° separation from the Moon.

The measurements span from September 29, 2019 to January 20, 2020 and have an average nominal error bar of 93 cm s⁻¹. We also included the time series of Hα measurements, the full width half maximum (FWHM), and the S-index in our analysis. The velocity and ancillary indicators are extracted from the spectra with the standard ESPRESSO pipeline v 2.0.0 (Pepe et al. 2021). The RV time series with nominal error bars is shown in Fig. 5.

Fig. 3 Extraction of 80 × 80 arcsec of the CHEOPS field-of-view for two different data frames at the beginning (top) and the end (bottom) of the second visit (see Table 2). The TOI-178 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. Bottom-left 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. Bottom-right 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 TOI-178.

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 TOI-178 with periods of 6.55 d, 10.35 d, and 9.96 d, which transited four, three, and two times, respectively. In addition, our analysis of this dataset with the DST (Détection Spécialisée de Transits Cabrera et al. 2012) yielded two additional candidates: a clear 3.23 d signal and a fainter 1.91 d signal. Upon receiving visits 1 and 2 from CHEOPS (Table 2), a study of the CHEOPS data alone with successive applications of the boxed-least-square (BLS) algorithm (Kovács et al. 2002) retrieved the 6.55 d, 3.23 d, and 1.91 d signals, in phase agreement with the TESS data. An additional dip, consistent in epoch with a transit of the 9.9 d candidate, was also identified; however, it could also marginally correspond to a transit of the 10.35 d candidate. The transit of one of these candidates was hence missing, pointing to a mis-attribution of transits in the TESS Sector 2 data.

In order to identify new possible solutions for the available data (TESS Sector 2, NGTS visits in September and October 2019, and CHEOPS visits 1 and 2), we individually pre-detrended 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₁ and removed thecorresponding signal in the light curve; we then applied the BLS a second time to obtain a second candidate of period P₂. 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₁, and we removed the associated signal in the light curve. For this candidate, i = 1 since we ignored one peak in the first BLS. For the second iteration of the BLS, we did not ignore any peaks and hence take the largest one (j = 0), leading to the pair of candidates c_1,0. We repeated this process 25 times, yielding 25 potential pairs of candidates c_{0≤i≤4,0≤j≤4}. For all of these potential solutions, we modelled the transits of the five candidates – 1.91 d, 3.23 d, 6.55 d, and c_i,j – using the batman package (Kreidberg 2015) and ran an MCMC on the pre-detrended light curve to estimate the relative likelihood of the different c_i,j.

The c_3,0 pair was favoured (its first BLS is shown with the green curve in the top panel of Fig. 6), 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 six-planet 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 planet-to-star radius ratio R_p ∕R_⋆, the impact parameter b, the orbital period P, and the mid-transit time T₀. We also fitted for the stellar density ρ_⋆, which, together with the orbital period P of each planet, defines through Kepler’s third law a value for the scaled semi-major axis a∕R_⋆ of each planet. We assumed a normal prior for the stellar density based on the value and uncertainty reported in Table 1 (Sect. 3) and wide uniform priors for the other transit parameters. The orbits were assumed to be circular, as justified in Sect. 6.3. For each bandpass (TESS, CHEOPS, and NGTS), we fitted two quadratic limb-darkening coefficients, which were parametrised with the (q₁, q₂) sampling scheme of Kipping (2013). Normal priors were placed on these limb-darkening coefficients based on Claret & Bloemen (2011).

We modelled the correlated noise present in the light curves simultaneously with the planetary signals to ensure a full propagation of the uncertainties. We first performed individual analyses of each light curve in order to select for each of them the best correlated noise model based on Bayesian evidence. We explored a large range of models for the CHEOPS light curves, consisting offirst- to fourth-order polynomials in the recorded external parameters (most importantly: time, background level, position of the point-spread-function (PSF) centroid, spacecraft roll angle, and contamination), as well as GPs (celerite Foreman-Mackey et al. 2017 and george Ambikasaran et al. 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 roll-angle-dependent photometric variations (cf. Sect. 4.1.2). First- to second-order polynomials in time and x − y centroid position were also needed for the first three visits. For the TESS light curve, we used a Matérn 3/2 GP against time, and we used a linear function of airmass for the NGTS light curves. For each light curve, we also fitted a jitter term, which was added quadratically to the error bars of the data points, to account for any underestimation of the uncertainties or any excess noise not captured by our modelling.

The results of our joint fit are displayed in Tables 3 (planets b, c, and d) and 4 (planets e, f, and g). The transits of all detected planets are shown in Figs. 1 (TESS), 2 (NGTS), and 4 (CHEOPS), with the phase-folded transits in the CHEOPS data presented in Fig. 7.

All planets appear to be in the super-Earth to mini-Neptune range, with the inner two planets falling to either side of the radius valley (Fulton et al. 2017). The inclination of the planets are also worth mentioning: The projected inclination of the four outer planets differ only by about 0.1 deg. As discussed in Agol et al. (2021) in the case of TRAPPIST-1, it is unlikely that the underlying inclinations and ascending nodes are scattered given the clustering of the projected inclination; hence, the outer part of the system is probably quite flat. In addition, the uncertainty on the inclination of the inner planets allows for the near-coplanarity of the entire system, a feature that was also found in the TRAPPIST-1 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 ℓ₁ -periodogram⁸ 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 $\log {R^{\prime }}_{\text{HK}}$ (Noyes 1984) revealed that they exhibit statistically significant periodicities at ≈ 36 days and ≈16 days, such that stellar activity effects are to be expected in the RVs, especially at these periods. In our analysis, stellar activity has been taken into account both with a linear model constructed with activity indicators smoothed with a GP regression similarly to Haywood et al. (2014), which we call the base model, and with a Gaussian noise model with a white, correlated (Gaussian kernel), and quasi-periodic component.

Radial velocity signals found to be statistically significant might vary from one activity model to another. To confirm the robustness of our detections, we tested whether signal detections can be claimed for a variety of noise models, following Hara et al. (2020). This approach consists of three steps. First, we computed the ℓ₁ -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 quasi-periodic. According to the analysis of the ancillary indicators, the quasi-periodic term of the noise is fixed at 36.5 days. We considered a grid of values for each noise component (amplitude and decay timescale) and computed the ℓ₁ -periodograms. Secondly, we ranked the noise models with a procedure based on cross-validation (CV). Finally, we examined the distribution of FAPs of each signal in the 20% highest ranked models. Here we will succinctly describe the conclusions of our analysis and refer the reader to Appendix B for a more detailed presentation. The ℓ₁ -periodogram corresponding to the highest ranked model is represented in Fig. 9. The periods at which the peaks occur are shown in red and account for most of the signals that might appear with varying assumptions on the noise (or activity) model. More precisely, we find the following results. We note that these results stem from an analysis of over 1300 noise models and that they are not all portrayed in Fig. 9, which only shows the results from the highest ranked noise model.

First, our analysis yields a consistent, significant detection of signals close to 3.2, ≈36, and ≈16 days. Radial velocity then allows an independent detection of planet c. Signals at ≈36 and ≈16 d appear in ancillary indicators, and as such we attributed these apparent periodicities in RV to stellar activity. The 16-day signal is, however, very likely partly due to planet f. Indeed, when the activity signals close to 40 and 16 days are modelled and the signal of transiting planets is removed, one finds a residual signal at 15.1 or 15.2 days, even though the 16 and 15.2-day signals are very close.

Second, we find signals, though not statistically significant ones, at 6.5 and 9.8 days (consistent with planets d and e) and at 2.08 days, which is the one-day alias (see Dawson & Fabrycky 2010) of 1.91 days (planet b). Third, we do not consistently find a candidate near 20.7 days. However, a 20.6 d signal appears in the highest ranked model, and a stellar activity signal might hide the signal corresponding to TOI-178 g.

Fourth, the signal at 43.3 days appearing in Fig. 9 might be a residual effect of an imperfect correction of the activity. However, we have not strictly ruled out the possibility that it stems from a planetary companion at 45 days. As will be discussed in the conclusion, this period corresponds to one of the possible ways to continue the Laplace resonance beyond planet g. Finally, depending on the assumptions, hints at 1.2 or 5.6 days (aliases of each other) can appear.

For comparison, we performed the RV analysis with an iterative periodogram approach. This analysis is able to show signals corresponding to 15.2, 3.2, and 6.5 days; however, it is unable to clearly establish the significance of the 3.2 d signal and fails to unveil candidates at 1.91 and 9.9 days.

The photometric data allow us to independently detect six planets at 1.91, 3.24, 6.56, 9.96, 15.23, and 20.71 days. As detailed in Appendix B, we find that the phases of RV and photometric signals are consistent within 2σ. We phase-folded the RV data at the periods given in Tables 3 and 4, which are shown in Fig. 10 with periods increasing from top to bottom. The variations at 3.24 and 15.2 days, corresponding to the planets 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 phase-folded to the periods of each of the planets, with signals of the other planets removed. Unbinned data are shown with grey points, data in 15-min 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 TOI-178 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 low-frequency variations, which can also be anticipated from the analysis of the indicators. The second model of activity consists of a correlated Gaussian noise with a quasi-periodic kernel. We computed the posterior distributions of the orbital parameters with an adaptive MCMC, as in Delisle et al. (2018). This analysis is presented in detail in Appendix B.5.

The mass estimates for each model are given in Tables B.2, 3, and 4, respectively. The 1σ intervals obtained with the two methods all have a large overlap. In Tables 3 and 4, 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 ${1.50}_{-0.44}^{+0.39}$ ${7.72}_{-1.52}^{+1.67}$ M_⊕ and ${0.177}_{-0.061}^{+0.055}$ ${1.02}_{-0.23}^{+0.28}$ ρ_⊕.

It appears that the mass of planet f (at 15.24 days) is the highest. This result might seem untrustworthy since there are activity signals close to 16 days in the ancillary indicators and since 1∕(1∕15.23 − 1∕16) = 316 days, which is greater than the observation time span. However, the two activity models considered here – one of which includes a signal at 16 days –yield a similar mass estimates, and the convergence of MCMC was ensured by computing the number of effective samples. As a consequence, we deemed the mass estimate appropriate. Nonetheless, more stellar activity models can be considered, and the mass intervals might still evolve as more data come along and the results become less model-dependent. A longer baseline would be suitable, so that the difference between 15.2 and 16 would be greater than the frequency resolution.

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 Phase-folded RV. Error bars correspond to nominal errors.

6 Dynamics

6.1 A Laplace resonant chain

Mean-motion resonances are orbital configurations where the period ratio of a pair of planets is equal to, or oscillates around, a rational number of the form (k + q)∕k, where k and q are integers. In TOI-178, candidates b and c are close to a second-order MMR (q = 2): P_c ∕P_b = 1.6915 ≈ 5∕3, while c, d, e, f, and g are close, pairwise, to first-order MMRs (q = 1): P_d ∕P_c = 2.0249 ≈ 2∕1, P_e ∕P_d = 1.5191 ≈ 3∕2, P_f ∕P_e = 1.5290 ≈ 3∕2, and P_g ∕P_f = 1.3595 ≈ 4∕3. Pairs of planets lying just outside MMRs are common occurrences in systems observed by transit (Fabrycky et al. 2014). To study such pairs of planets, a relevant quantity is the distance to the exact resonance. Taking the pair of c and d as an example, the distance to the 2∕1 MMR in the frequency space is given by 1n_c − 2n_d, where n_c = 2π∕P_c is the mean motion of planet c. Following Lithwick et al. (2012), we name the associated timescale the ‘super-period’: $$P_{c,d}\equiv \frac{1}{|(k+q)/P_d-k/P_c|}.$$(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 super-period, with amplitudes depending on the distance to the resonance, the mass of the perturbing planet, and the eccentricities of the pair (Lithwick et al. 2012). The fact that the super-periods of all 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: $$\begin{array}{lll}\text{d}{\psi }_1/\text{d}t\hfill & =\hfill & \text{d}({\varphi }_1-{\varphi }_2)/\text{d}t=1n_c-4n_d+3n_e\approx 0,\hfill \\ \text{d}{\psi }_2/\text{d}t\hfill & =\hfill & \text{d}({\varphi }_2-{\varphi }_3)/\text{d}t=2n_d-5n_e+3n_f\approx 0,\hfill \\ \text{d}{\psi }_3/\text{d}t\hfill & =\hfill & \frac{1}{2}\text{d}({\varphi }_3-{\varphi }_4)/\text{d}t=1n_e-3n_f+2n_g\approx 0,\hfill \end{array}$$(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 b-c-d triplet to form a Laplace relation, which is not the case (see Table 5). To continue the chain, planet b would have needed a period of ~1.95 d. Its current period of 1.91 d could indicate that it was previously in the chain but was pulled away, possibly by tidal forces.

Figure 11 shows the evolution of the Laplace angles when integrating the nominal solution given in Tables 3 and 4, starting at the beginning of the observation of TESS Sector 2. The three angles librate over the integrated time for the selected initial conditions, with combinations of periods ranging from a few years to several decades. The exact periods and amplitudes of these variations depend on the masses and eccentricities of the involved planets. The theoretical equilibria of the resonant angles are discussed in Sect. 6.2, while the long-term stability of this system is discussed in Sect. 6.3.

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 Kepler-223 are observed to librate around one of the six possible equilibria predicted by theory (Mills et al. 2016; Delisle 2017).

We used the method described in Delisle (2017) to determine the position of the possible equilibria for the Laplace angles of TOI-178. The five external planets (c to g) orbiting TOI-178 are involved in a 2:4:6:9:12 resonant chain. All consecutive pairs of planets in the chain are close to first-order MMRs (1:2, 2:3, 2:3, 3:4). Moreover, as in the Kepler-223 system, there are also strong interactions between non-consecutive planets. Indeed, planet e and planet g (which are non-consecutive) are also close to a 1:2 MMR. As explained in Delisle (2017), this breaks the symmetry of the equilibria (i.e. the equilibrium is not necessarily at 180 deg), and the position of the equilibria for the Laplace angles depends on the planets’ masses.

We solved for the position of these equilibria using the masses given in Tables 3 and 4. We also propagated the errors to estimate the uncertainty on the Laplace angle equilibria. We find two possible equilibria for the system, which are symmetric with respect to 0 deg. We provide the values of the Laplace angles corresponding to the first equilibrium in Table 6 (the second is simply obtained by taking ψ_j →−ψ_j for each angle).

It should be noted that these values correspond to the position of the fixed point around which the system is expected to librate. Depending on the libration amplitude, the instantaneous values of the Laplace angles can significantly differ from the equilibrium.For instance, in the case of Kepler-223, the amplitude of libration could be determined and is about 15 deg for all Laplace angles (Mills et al. 2016). In Table 6, we observe that all instantaneous values of Laplace angles (as of August 2018) are also found within 15 deg of the expected equilibrium. We note that ψ₁ 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, ψ₂ was moving closer to the equilibrium and ψ₃ was only slowly evolving during this two-year 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 TOI-178 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 (quasi-periodic) 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): $${\log }_{10}\left|n_f^{(1)}-n_f^{(2)}\right|,$$(3)

where $n_f^{(1)}$ and $n_f^{(2)}$ 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 quasi-vertical structures that contain very stable regions in their central parts and less stable regions at their edges. These are mainly MMRs between two or three planets. In particular, the black area crossed by the dashed white line in the top panel of Fig. 12 corresponds to the resonant chain where the nominal system is located.

The bottom panel of Fig. 12 shows a different section of the phase space where P_f and P_g vary and all eccentricities are initialised at 0. This reveals part of the resonant structure of the system. The black regions are stable, while the green to red areas mark the instability caused by the resonance web. This figure still highlights the stability island in which the planetary system is located. This narrow region is surrounded by resonances: the n_e − 3n_f + 2n_g = 0 Laplace resonance (central diagonal strip), high-order two-body MMRs between planets f and g (horizontal strips), and high-order two-body MMRs between planets e and f (vertical lines).

The extent of this resonant chain versus the initial period and mass of planet f is shown in Fig. 13. On both panels, the x-axis corresponds to the orbital period of planet f, while the y-axis corresponds to the planetary mass. The figure presents two different, but complementary, aspects of the dynamics of TOI-178. The bottom panel indicates whether the Laplace angles ψ₁, ψ₂, and ψ₃ 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 long-term stability. This central region seems not to depend strongly on the mass of planet f.

On the other hand, the modification of the orbital period within a fraction of a day changes the dynamics of the system. Indeed, outside of the central region where the three Laplace angles librate, which is about 0.015 days wide, chaotic layers are present (red in the top panel panel, dark blue in the bottom) where none of the three Laplace angles librate. Here, the red colour of the stability index corresponds to a significant, but moderate, diffusion rate. Thus, although in this region the trajectories are not quasi-periodic, the chaos has limited consequences (it is bounded) and probably does not lead to the destruction of the system.

Outside of these layers lies a very stable (quasi-periodic) region. The blue on the Laplace angle map shows that only one Laplace angle librates, while the others circulate. Although the 2:4:6:9:12 chain is broken, planets c, d, e, and g remain inside the 2:4:6:12 resonant chain, independently of the orbital period of planet f. The stability map indicates a very strong regularity of the whole region. Nevertheless, one can notice the presence of some narrow zones where the diffusion is more important; they are induced by high-order orbital resonances but do not have any significant consequence on the stability of the system.

The study described in this section gives only a very partial picture of the structure of the phase space (parameter space) of the problem. Nevertheless, it can be seen that as long as the system is in the complete resonant chain, which is the case for the nominal derived parameters and for the bulk of the posterior given in Tables 3 and 4, it remains stable.

Fig. 13 Same stability indicator as in Fig. 12 (top) and the number of librating Laplace angles (bottom) for TOI-178 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 first-order MMRs, we expect TTVs to occur over the super-period (Eq. (1), see also Lithwick et al. 2012), which is roughly 260 days for all these pairs (Table 5). The amplitude of these TTVs depends on the masses and eccentricities of each planet of a pair. Since the stability analysis concludes that the system is more stable with eccentricities close to 0, we integrated the initial conditions described in Tables 3 and 4 over 6 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 super-periods, the six-year evolution of the TTVs shows hints of the long-term evolution of the Laplace angles.

As the exact shape of the TTVs depends mainly on the masses of the involved planets, the sampling of the 260-day signal over a few years, in addition to a longer monitoring of the evolution of the Laplace angles, can provide precise constrains on the orbital configuration and masses of this system, as well as the presence of other planets in the chain. The detailed study of the TTVs of this system will be the subject of a forthcoming paper.

7 Internal structure

7.1 Minimum mass of the protoplanetary disc

The total mass of the planets detected in TOI-178 is approximately ${24.8}_{-6.23}^{+5.96}{M}_{\oplus }$. 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 TOI-178. Assuming a disc mass equal to 1% of the stellar mass, and using the metallicity of TOI-178 ([Fe/H] = −0.23, corresponding to a metal content of Z ~ 0.0061), the mass of heavy elements in such a disc would equal ~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 TOI-178, seems therefore to be as applicable for this system as for the Solar System. We note, however, that in Hayashi (1981) the MMSN model is assumed to be based on the in situ formation of planets, whereas here we compare the mass of planets with the whole mass of solids in the inner parts. Another implication of this comparison is that a problem in term of available mass would appear should the planetary masses be revised to higher values or should another massive planet be detected in the same system. In such a case, the TOI-178 system would point towards a formation channel similar to the one envisioned for the Trappist-1 system (Schoonenberg et al. 2019).

7.2 Mass-radius relation

The six planets we detected in the TOI-178 system are in the super-Earth to mini-Neptune range, with radii ranging from ${1.152}_{-0.070}^{+0.073}$ to ${2.87}_{-0.13}^{+0.14}$ R_⊕. Although the mass determination is limited by the extent of the available spectroscopic dataset, planets b and c appear to have roughly terrestrial densities of ${0.98}_{-0.31}^{+0.35}{\rho }_{\oplus }$ and ${1.02}_{-0.23}^{+0.28}$, respectively, where ρ_⊕ is the density of the Earth. The outer planets seem to have significantly lower densities; in particular, we estimate the density of planet d to be ${0.177}_{-0.061}^{+0.055}{\rho }_{\oplus }$.

Figure 15 shows the position of the six planets in a mass-radius diagram, in comparison with planets with mass and radius uncertainties less than 40% (light grey). Planets belonging to five systems in Laplace resonance are indicated in the same diagram: Trappist-1, K2-138, Kepler-60, Kepler-80, and Kepler-223. The diversity of planetary composition in TOI-178 is clearly visible in the diagram: The two inner planets have a radius compatible with a gas-free structure, whereas the others contain water and/or gas. This is similar to the Kepler-80 system, where the two innermost planets are compatible with a gas-free structure and the two outermost ones likely contain gas. Planets in Kepler-223, Kepler-60, and Trappist-1 seem to have a more homogeneous structure: All planets in Kepler-223 have a gas envelope, and the planets in Kepler-60 and Trappist-1 have small gas envelopes (below ~.01% of the total mass).

Considering the TOI-178 system in greater detail, planets d, e, and g are located above the pure water line and definitely contain a non-zero gas mass fraction. Planet d and, depending on its mass, planet g, must contain a large gas fraction.

Figure 16 is a diagram showing the location of the TOI-178 planets in stellar insolation versus the planetary radius. The colour code illustrates the density of exoplanets. This diagrams clearly shows the so-called evaporation valley⁹. 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 TOI-178, 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 TOI-178 system and other systems in Laplace resonance can clearly be seen in Fig. 17, where we show the density of planets as a function of their equilibrium temperature for the same systems as in Fig. 15. In Kepler-60, Kepler-80, and Kepler-223, the density of the planets decreases when the equilibrium temperature decreases. This can be understood as the effect of evaporation that removes part of the primordial gaseous envelope, this effect being stronger for planets closer to their stars. The densities of K2-138 are also 1σ consistent with such behaviour. In Trappist-1, the density of planets is always higher than 4 g cm⁻³ and increases (with the exception of Trappist-1 f) with decreasing equilibrium temperature. This is likely due to the presence of more ices in planets far from their stars (Agol et al. 2021). Contrary to the three Kepler systems, in the TOI-178 system the density of the planets is not a growing function of the equilibrium temperature. Indeed, TOI-178f has a density higher than that of planet e, and TOI-178d has a density smaller than that of planet e.

Planet f is substantially more massive than all the other planets in the system. From a formation perspective, one would expect this planet to have a density smaller than the other planets, in particular planet e, at the end of the formation phase, similar to what we observe, for example, for Jupiter, Saturn, Uranus, and Neptune. Since this planet is farther away from the star compared to planet e, evaporation should have been less effective for planet f compared to planet e. The combined effect of formation and evolution should therefore lead to a smaller density for planet f compared to planet e. Similarly, planet d is smaller than planet e and is located closer to the star. Using the same arguments, the combined effect of formation and evolution should have led to planet d having a density larger than that of planet e. The TOI-178 system therefore seems at odds with the general understanding of planetary formation and evaporation, where one would expect the density to decrease when the distance to the star increases, or when the mass of the planet increases.

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 iron-sulfur inner core, a mantle, a water layer, and a gas layer. We used the equation of state (EOS) for Hakim et al. (2018) for the core, the EOS from Sotin et al. (2007) for the silicate mantle, and the EOS from Haldemann et al. (2020) for the water. These three layers constitute the ‘solid’ part of the planets. The thickness of the gas layer (assumed to be made of pure H/He) is computed as a function of the stellar age, mass, and radius of the solid part as well as irradiation from the star, using the formulas in Lopez & Fortney (2014). The internal structure parameters of each planet are therefore the iron molar fraction in the core, the Si and Mg molar fraction in the mantle, the mass fraction of all layers (inner core, mantle, and water), the age of the planet (equal to the age of the star), and the irradiation from the star. More technical details regarding the calculation of the forward model are given in Appendix D.

In the Bayesian analysis part of the model, we proceeded in two steps. We first generated 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 semi-amplitude for each of the planets and retained the models that fitted the observed data within the error bars. By doing so, we included the fact that all synthetic planets orbit a star with exactly the same parameters. Indeed, planetary masses and radii are correlated by the fact that the fitted quantities are the transit depth and RV semi-amplitude, which depend on the stellar radius and mass. In order to take this correlation into account, it is therefore important to fit the planetary system at once, rather than fitting each planet independently.

For the Bayesian analysis, we assumed the following two priors: first, that the mass fraction of the gas envelope is uniform in log; and second, that the mass fraction of the inner core, mantle, and water layer are uniform on the simplex (the surface on which 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 TOI-178 with other systems with a Laplace resonance, we also performed a Bayesian analysis for the K2-138, Kepler-60, Kepler-80, and Kepler-223 systems to compute the probability distribution of the gas mass in the different planets. The parameters for all systems are taken from the references mentioned above, and for all systems we have assumed that [Si/H] = [Mg/H] = [Fe/H]. We did not consider the Trappist-1 system in this comparison as it is likely that the variations in the densities of these planets result from variations in their ice contents (Agol et al. 2021).

The posteriors distributions of the two most important parameters (mass fractions and composition of the mantle) of each planet in TOI-178 are shown in Appendix D, Figs. D.2 to D.7. We focus here on the mass of gas in each planet, and we plot the mass of the gaseous envelope for each planet as a function of their equilibrium temperature in Fig. 18.

In all three Kepler systems, the gas mass in planets generally decreases when the equilibrium temperature decreases, as can be seen in Fig. 18. One exception to this tendency is the mass of the envelope of Kepler-223d, which is larger than that of planet e in the same system. Kepler-223d is, however, more massive than planets c and e in the same system, and one would expect from formation models that the mass of the primordial gas envelope is a growing function of the total planetary mass. In K2-138, the gas fraction is compatible (within 1-sigma error bars) with a monotonous function of the equilibrium temperature. In the case of TOI-178, the mass of gas also globally increases when the equilibrium temperature decreases, with the notable exception of planet d. Indeed, a linear interpolation would provide a gas mass for planet d of the order of 10⁻⁶−10⁻⁵M_⊕, whereas the interior structure modelling gives values of 9.66 × 10⁻³M_⊕ and 2.56 × 10⁻²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 TOI-178 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 mass-radius curves for some idealised compositions (Zeng et al. 2019). The six planets orbiting around TOI-178 are indicated; the colour of the points and error bars give the equilibrium temperature. The seven planets orbiting Trappist-1 are shown with diamonds, and the parameters are taken from Agol et al. (2021). The three planets orbiting Kepler-60 are shown with ‘X’ marks, and the parameters are taken from Jontof-Hutter et al. (2016). The six planets orbiting K2-138 are shown with squares, and the parameters are taken from Lopez et al. (2019). The four planets orbiting Kepler-80 are shown with inverse triangles, and the parameters are taken from MacDonald et al. (2016). The four planets orbiting Kepler-223 are shown with regular triangles, and the parameters are taken from Mills et al. (2016).

Fig. 16 Diagram of the position of the TOI-178 planets in a stellar light intensity (relative to Earth) versus planetary radius. The marker size is inversely proportional to the density, and the colour code gives the density of exoplanets, from yellow for empty regions of the diagram to violet for high-density and/or highly populated regions.

Fig. 17 Densities of planets in the TOI-178, Trappist-1, K2-138, Kepler-60, Kepler-80, and Kepler-223 systems as a function of their equilibrium temperatures. The error bars give the 16 and 84% quantiles, and the marker is located at the median of the computed distribution. The colour code gives the planetary masses in Earth masses. The parameters of the planets are taken from the references mentioned in Fig. 15.

8 Discussion and conclusions

In this study, we presented new observations of TOI-178 by CHEOPS, ESPRESSO, NGTS, and SPECULOOS. Thanks to this follow-up effort, we were able to determine the architecture of the system: Out of the three previously announced candidates at 6.56 d, 9.96 d, and 10.3 d, we confirm the first two (6.56 d and 9.96 d) and re-attribute the transits of the third to 15.23 d and 20.7 d planets, in addition to the detection of two new inner planets at 1.91 d and 3.24 d; all of these planets were confirmed by follow-up observations. In total, we therefore announce six planets in the super-Earth to mini-Neptune range, with orbital periods from 1.9 d to 20.7 d, all of which (with the exception of the innermost planet) are in a 2:4:6:9:12 Laplace resonant chain. All the orbits appear to be quasi-coplanar, with projected mutual inclinations between the outer planets estimated at about 0.1 deg, a feature that is also visible in other systems with three-body resonances, such as Trappist-1 (Agol et al. 2021) and the Galilean moons. Current ephemeris and mass estimations indicate that the TOI-178 system is very stable, with Laplace angles librating over decades. In the TESS Sector 29 data, made available during the referee process of this paper, we recovered the transits of planets b, c, d, e, and f at the predicted dates. The transit of planet g was not observed as it transited during the gap and was observed during the third CHEOPS visit (see Fig. 4).

As there is no theoretical reason for the resonant chain to stop at 20.7 days, and the current limit probably comes from the duration of the available photometric and RV datasets, we give in Table 7 the periods that would continue the resonant chain for first-order (q = 1) and second-order (q = 2) MMRs. These periods result from the equations detailed in Appendix C. In the TOI-178 system, as well as in similar systems in Laplace resonance, planet pairs are virtually all nearly first-order MMRs. The most likely of the periods shown in Table 7 are therefore the first-order solutions with low k, hence 45.00, 32.35, or 28.36 days; however, additional planets with such periods would not be transiting if they were in the same orbital plane as the others (see Fig. 8). We note that, for a star such as TOI-178, the inner boundary of the habitable zone lies around 0.2 AU, or at a period of the order of  40 days. Additional planets in the Laplace resonance could therefore orbit inside, or very close to, the habitable zone.

The brightness of TOI-178 allows for further characterisation of the system by photometric measurements, radial velocities, and transit spectroscopy. These measurements will be 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, hydrogen-dominated 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 TOI-178 system lie on both sides of the radius valley (Fulton et al. 2017). Therefore, reconstructing the past orbital and atmospheric history of this planet may provide clues regarding the origin of the valley. The system is located at a declination in the sky that makes it observable by most ground-based observatories around the globe. Furthermore, the host star is bright enough and the radii of the outer planets large enough to make them potentially amenable to optical and infrared transmission spectroscopy observations from both ground and space, particularly by employing the upcoming E-ELT and James Webb Space Telescope (JWST) facilities. Indeed, it has been shown that the JWST has the capacity to perform transmission spectroscopy of planets with radii down to 1.5 R_⊕ (Samuel et al. 2014).

Planets d and f are particularly interesting as their densities are very different from those of their neighbours, and they depart from the general tendency of planetary density (i.e. decreasing for decreasing equilibrium temperatures). The densities of planets d and f, in the context of the general trend seen in TOI-178, are difficult to understand in terms of formation and evaporation process, and they could be difficult to reproduce with planetary system formation models (Mordasini et al. 2012; Alibert et al. 2013; Emsenhuber et al. 2020). We stress that even though two different analyses yielded similar estimates of masses and densities, they were performed with only 46 data points. As such, these estimates need to be confirmed by further RV measurements, which would provide, in particular, a better frequency resolution and confirm the mass estimate of planet f.

The orbital configuration of TOI-178 is too fragile to survive giant impacts, or even significant close encounters: Fig. 12 shows that a sudden change in period of one of the planets of less than a few ~.01 d can render the system chaotic, while Fig. 13 shows that modifying a single period axis can break the resonant structure of the entire chain. Understanding, in a single framework, the apparent disorder in terms of planetary density on one side and 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.

Follow-up transit observations should also unveil TTVs for all but the innermost planet of the TOI-178 system (see Fig. 14), with two timescales: a ~ 260 d signal with an amplitude of several minutes to a few tens of minutes and a larger signal over several years or decades. Future observations of this system can hence be used to measure the planetary masses and estimate eccentricities directly from TTVs as well as compare the results with masses derived from RV measurements.

Finally, the innermost planet, b, lies just outside the 3:5 MMR with planet c; however, it is too far to be part of the Laplace chain, which would require a period of ~ 1.95d. Since the formation of the Laplace resonant chain probably results from a slow drift from a chain of two-planet resonances due to tidal effects (Papaloizou & Terquem 2010; Delisle et al. 2012; Papaloizou 2015; MacDonald et al. 2016), the current state of the system might constrain the dissipative processes that tore apart the innermost link of the chain while the rest of the configuration survived.

The TOI-178 system, as revealed by the recent observations described in this paper, contains a number of very important features: Laplace resonances, variation in densities from planet to planet, and a stellar brightness that allows a number of follow-up 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 TOI-178, K2-138, Kepler-60, Kepler-80, and Kepler-223 systems as a function of their equilibrium temperatures. The error bars give the 16 and 84% quantiles and the marker is located at the median of the computed distribution. The colour code gives the planetary masses in Earth masses. The parameters of the planets are taken from the references mentioned in Fig. 15.

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 ANR-10-EQPX-29-01) 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 Wallonia-Brussels 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 ESP2016-80435-C2-1-R, ESP2016-80435-C2-2-R, PGC2018-098153-B-C33, PGC2018-098153-B-C31, ESP2017-87676-C5-1-R, MDM-2017-0737 Unidad de Excelencia “María de Maeztu”- Centro de Astrobiología (INTA-CSIC), 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 ESA-appointed 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/FIS-AST/7002/2020), ENGAGE SKA (POCI-01-0145-FEDER-022217), and PHOBOS (POCI-01-0145-FEDER-029932), 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/FIS-AST/32113/2017 and POCI-01-0145-FEDER- 032113; PTDC/FIS-AST/28953/2017 and POCI-01-0145-FEDER-028953; PTDC/FIS-AST/28987/2017 and POCI-01-0145-FEDER-028987. 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 (PP00P2-190080). 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 RYC-2013-14875. J.I.G.H., A.S.M., R.R., and C.A.P. acknowledge financial support from the Spanish MICINN AYA2017-86389-P. 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 AFB-170002. 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/FIS-OUT/29048/2017, COMPETE2020: POCI-01-0145-FEDER-028987 & FCT: PTDC/FIS-AST/28987/2017. N.J.N is supported by the contract and exploratory project IF/00852/2015, and project PTDC/FIS-OUT/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 PGC2018-098153-B- 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 GINOP-2.3.2-15-2016-00003, K-119517, K-125015, and the City of Szombathely under Agreement No. 67.177-21/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.S-FNRS Research Associate. J.I.V. acknowledges support of CONICYT-PFCHA/Doctorado Nacional-21191829. M. R. Z. O. acknowledges financial support from projects AYA2016-79425-C3-2-P and PID2019-109522GB-C51 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 post-processing 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 DRP-produced light curves (DEFAULT aperture) of the four CHEOPS visits to TOI-178 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 TOI-178 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 TOI-178 visits.

Fig. A.2 TOI-178 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 second-order 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 ℓ₁ -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 ℓ₁ -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 ℓ₁-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 ℓ₁-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 ℓ₁ norm) whileensuring that the inverse Fourier transform is close enough to the model in a certain, precise sense. Due to the ℓ₁ norm penalty, frequencies ‘compete’ against one another to have non-zero 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 ℓ₁ 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 ℓ₁ -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 $\log {R^{\prime }}_{\text{HK}}$ 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 $\log {R^{\prime }}_{\text{HK}}$ time series as provided by the ESPRESSO pipeline. We computed the residual periodograms as described in Baluev (2008). These periodograms allowed us to take into account general linear base models that are fitted along candidate frequencies. We computed the periodograms and iteratively added a sinusoidal function whose frequency corresponds to the maximum peak of the periodogram. Iterations 1 and 2 are shown in the top and bottom panels of Fig. B.1, respectively. The 36-day and 16-day positions (top and bottom panels, respectively) are marked with dotted black lines. Pursuing the iterations, we find signal detections with FAP < 10⁻³ of periods at 36, 115, and 15.9 days for Hα; 35.5, 20.8, and 145 days for the FWHM; 36 and 16 days for the bisector span; and 36.7 and 16.5 days for the $\log {R^{\prime }}_{\text{HK}}$. The 36-day periodicities are always detected with FAP < 10⁻⁶. These results indicate that activity effects in the RV at ≈36 and ≈16 days are to be expected, along with low frequency effects. These signals likely stem from the rotation period of the star, which creates signals at the fundamental frequency and the first harmonic.

We now turn to the RV and define the alternative noise models we explored. These are Gaussian, with a white component, an exponential decay, and a quasi-periodic term, as given by the formula $$\begin{array}{ll}{V}_{kl}\hfill & ={\delta }_{k,l}({\sigma }_k^2+{\sigma }_W^2)+{\sigma }_R^2{\text{e}}^{-\frac{{(t_k-t_l)}^2}{2{\tau }_R^2}}+{\sigma }_{QP}^2{\text{e}}^{-\frac{{(t_k-t_l)}^2}{{\tau }_{QP}^2}{\sin }^2\left(\frac{t_k-t_l}{P_{\text{act}}}\right)}\hfill \end{array},$$(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, S-index, and Hα) showed that they all exhibit statistically significant variations at ≈40 (36 days), as well as a periodogram peak at 16 d. The 36 and 16 d signals exhibit phases compatible with each other at 1 sigma, except for the 36 d signal in the FWHM which is 3σ away from the phase fitted on the S-index and Hα. This points to a stellar rotation period of 36 d. We considered all the possible combinations of values for σ_R and σ_W in 0.0,0.5 1.0,1.25, 1.5,1.75,2 m s⁻¹, τ = 0, 2, 4, 6 d, P_act = 36.5 d, σ_QP = 0,1,2,3,4 m s⁻¹, and τ_QP = 18, 36, or 72 d.

The computation of the ℓ₁-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 ℓ₁ -periodogram obtained with the noise model with the maximum CV score with three different base models (from top to bottom: without any indicator, with smoothed Hα, and with both smoothed Hα and smoothed FWHM). In all cases, we find signals at 35–40 d, 15 d, and 3.24 d, as well as a peak at 6.55 d. We can also find peaks at 21 d, 9.84 d, 2.08 d, 1.92 1.44, or 1.21 days, depending on the models used. We note that 2.08 and 1.92 d are aliases of each other, so that the presence of one or the other might originate from the same signal. The model with the overall highest CV score includes only the smoothed Hα time series in the base model and with the notations of Eq. (B.1), σ_W = 1.75 m s⁻¹, σ_R = 1.5 m s⁻¹, τ = 2 days, and σ_QP = 0, and it corresponds to the middle panel in Fig. B.2.

We computed the ℓ₁-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₂₀, 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₀, 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₀ 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.2-day 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.2-day 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.7-day transiting planet could correspond to the peak at 21.6 d in the ℓ₁ -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 semi-periodic 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² s⁻², and the prior on the noise timescale is log-uniform between 1h and 30 days. The second model is a GP (or here, correlated Gaussian noise) with a quasi-periodic kernel k of the form $$k(t;{\sigma }_W,{\sigma }_R,\lambda ,\nu )={\sigma }_W^2+{\sigma }_R^2{\text{e}}^{-\frac{t}{\tau }}\cos (\nu t).$$(B.2)

We imposed a Gaussian truncated prior on ${\sigma }_W^2,{\sigma }_R^2$ with σ = 100m² s⁻². We imposed a flat prior on ν between 2π∕50 and 2π∕30 rad day⁻¹ and a log-uniform prior on τ on 1h to 1000 days. For planets b, c, d, g, e, and f, the priors on the periods and the times of conjunction are set as Gaussian with means and variances set according to the constraints obtained from the joint fit of the TESS and CHEOPS data. The densities are computed by combining the posterior samples of mass and the posterior samples of radii, assuming the mass and radius estimates are independent.

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)¹⁰, 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 Foreman-Mackey et al. (2017). We also tried the SHO kernel and simply imposed that the quality factor Q be greater than one half, the other parameters having flat priors. We found very similar results, and, as such, they are not reported in Table B.2.

Modelling errors might leave a trace in the residuals. In Hara et al. (2019), it is shown that if the model used for 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 Shapiro-Wilk normality test (Shapiro & Wilk 1965) and found a p-value of 0.78, which is compatible with normality. The variogram does not exhibit signs of correlations. The same analysis was performed with model 2, for which the residuals also exhibit neither non-normality nor correlation (p-value of 0.99 on the residuals). We conclude that both models are compatible with the data.

Appendix C Continuation of a Laplace resonant chain

TOI-178 is in a configuration where successive pairs of planets are at the same distance to the exact neighbouring first-order MMR. Generalising the configuration a bit, we define fictive planets 1, 2, and 3 such that planets 1 and 2 are close to the resonance (k₁ + q₁) : k₁ and planets 2 and 3 are close to the resonance (k₂ + q₂) : k₂, where k_i and q_i are integers. We hence write the near-resonant angles: $$\begin{array}{ccc}{\phi }_1& =& k_1{\lambda }_1-(k_1+q_1){\lambda }_2,\\ {\phi }_2& =& k_2{\lambda }_2-(k_2+q_2){\lambda }_3,\end{array}$$(C.1)

where λ_i is the mean longitude of planet i. The associated distances to the resonances read: $$\begin{array}{ccc}{\Delta }_1& =& k_1{n}_1-(k_1+q_1)n_2,\\ {\Delta }_2& =& k_2{n}_2-(k_2+q_2)n_3,\end{array}$$(C.2)

where n_i is the mean motion of planet i. A Laplace relation exists between these three planets if: $$j_1{\Delta }_1-j_2{\Delta }_2\approx 0,$$(C.3)

where j₁ and j₂ are integers. In addition, the invariance by rotation of the Laplace angle (D’Alembert relation) gives: $$j_1{q}_1=j_2{q}_2.$$(C.4)

As a result, the Laplace relation requires: $$n_3\approx \frac{k_2{n}_2-\frac{q_2}{q_1}{\Delta }_1}{k_2+q_2},$$(C.5)

which translates, for the period of the third planet, as: $$P_3\approx \frac{k_2+q_2}{\frac{k_2}{P_2}-\frac{q_2}{q_1{P}_{1,2}}},$$(C.6)

where P_1,2 is the super-period associated with Δ₁. From this we can compute the periods of potential additional planets that would continue the Laplace resonant chain of TOI-178. Taking planets f and g as planets 1 and 2, the formula to compute the possible period of a planet x that could continue the resonant chain becomes: $$P_x=\frac{k+q}{\frac{k}{P_f}-\frac{q}{P_{f,g}}},$$(C.7)

where P_f,g is the super-period between the near first-order resonances of the known chain defined by Eq. (1), here ~ 260 days, and k and q are integers such that planet x and f are near a (k + q)∕k MMR. Some of the relevant periods are displayed in Table 7. A similar computation allowed us to determine the possible period of planet f prior to its confirmation by CHEOPS (see Fig. C.1).

Fig. C.1 Super-period of planet f with respect to e (in red) and g (in black) for a set of first-order MMRs. The purple line indicates the value of the super-period (260 day) between the other pairs of the chain. Only two possible periods allowed planet f to be part of the resonant chain, which correspond to the near-intersection of the right-hand side of a red curve, the left-hand side of a black curve, and the horizontal purple line.

Appendix D Internal structure

We provide here the posterior distributions of the interior planetary models, as well as some more details on their calculations. As mentioned in the main text, deriving the posterior distribution of the internal structure parameters requires computing the radius of planets millions of times for different sets of parameters. In order to speed up this calculation, we first computed a large (5 million-point) database of internal structure models, varying the different parameters. This database was split randomly into three sets: one training set (80% of the models), one validation set, and one test set (the last two each containing 10% of the whole database). We then, in a second time, fitted a deep neural network (DNN) in order to be able to compute the radius of a planet very rapidly with a given set of internal structure parameters. The architecture of the DNN we used is made of six layers of 2048 nodes each, and we used the classical rectified linear unit (ReLU) as an activation function (Alibert & Venturini 2019). The DNN is trained for a few hundred epochs, using a learning rate that is progressively reduced from 1.e − 2 down to 1.e − 4. Our DNN allowed us to compute these radii with an average error below 0.25% and with an increased rapidity of many orders of magnitude (a few thousand models computed per second). Figure D.1 shows the prediction error we reach on the test set (which was not used for training). The error on the predicted radius is lower than 0.4% – an error much lower than the uncertainty on the radii we obtained for the TOI-178 planets – in 99.9% of the cases.It is finally important to note that this model does not include the compression effect that would be generated 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 y-axis is in log scale, and the x-axis 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.

Fig. D.3 Same as Fig. D.2 but for planet c.

Fig. D.4 Same as Fig. D.2 but for planet d.

Fig. D.5 Same as Fig. D.2 but for planet e.

Fig. D.6 Same as Fig. D.2 but for planet f.

Fig. D.7 Same as Fig. D.2 but for planet g.

References

All Tables

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 mid-sector 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 TOI-178.02 with a period of 10.35 d (see Sect. 2).

In the text

Fig. 2 Light curves from simultaneous observations of TOI-178 by NGTS (green stars) and SPECULOOS-South (purple triangles), described in Sects. 4.1.3 and 4.1.4, respectively. Unbinned data are shown as grey points,and data in 15-minute 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 TOI-178.02 is shown with a dashed red line, the transit of planet g (which corresponds to an even transit of TOI-178.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 field-of-view for two different data frames at the beginning (top) and the end (bottom) of the second visit (see Table 2). The TOI-178 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. Bottom-left 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. Bottom-right 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 TOI-178.
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 phase-folded to the periods of each of the planets, with signals of the other planets removed. Unbinned data are shown with grey points, data in 15-min 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 TOI-178 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 Phase-folded 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 TOI-178 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 (quasi-periodic) 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 TOI-178 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 TOI-178, 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 TOI-178 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 mass-radius curves for some idealised compositions (Zeng et al. 2019). The six planets orbiting around TOI-178 are indicated; the colour of the points and error bars give the equilibrium temperature. The seven planets orbiting Trappist-1 are shown with diamonds, and the parameters are taken from Agol et al. (2021). The three planets orbiting Kepler-60 are shown with ‘X’ marks, and the parameters are taken from Jontof-Hutter et al. (2016). The six planets orbiting K2-138 are shown with squares, and the parameters are taken from Lopez et al. (2019). The four planets orbiting Kepler-80 are shown with inverse triangles, and the parameters are taken from MacDonald et al. (2016). The four planets orbiting Kepler-223 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 TOI-178 planets in a stellar light intensity (relative to Earth) versus planetary radius. The marker size is inversely proportional to the density, and the colour code gives the density of exoplanets, from yellow for empty regions of the diagram to violet for high-density and/or highly populated regions.
In the text

Fig. 17 Densities of planets in the TOI-178, Trappist-1, K2-138, Kepler-60, Kepler-80, and Kepler-223 systems as a function of their equilibrium temperatures. The error bars give the 16 and 84% quantiles, and the marker is located at the median of the computed distribution. The colour code gives the planetary masses in Earth masses. The parameters of the planets are taken from the references mentioned in Fig. 15.

In the text

Fig. 18 Gas mass in the planets of the TOI-178, K2-138, Kepler-60, Kepler-80, and Kepler-223 systems as a function of their equilibrium temperatures. The error bars give the 16 and 84% quantiles and the marker is located at the median of the computed distribution. The colour code gives the planetary masses in Earth masses. The parameters of the planets are taken from the references mentioned in Fig. 15.

In the text

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

Fig. A.2 TOI-178 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 second-order polynomial in time.

In the text

	Fig. B.1 Periodograms of the ancillary indicator time series (Hα, FWHM, bisector span, and $\log {R^{\prime }}_{\text{HK}}$ 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 Super-period of planet f with respect to e (in red) and g (in black) for a set of first-order MMRs. The purple line indicates the value of the super-period (260 day) between the other pairs of the chain. Only two possible periods allowed planet f to be part of the resonant chain, which correspond to the near-intersection of the right-hand side of a red curve, the left-hand side of a black curve, and the horizontal purple line.

In the text

	Fig. D.1 Histogram of the prediction error from our DNN on the test set. The y-axis is in log scale, and the x-axis covers an error from −1% to 1%.
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