Planets with radii between that of the Earth and Neptune hereafter referred to as subNeptunes are found in closein orbits around more than half of all Sunlike starsSUP12SUP. However their composition formation and evolution remain poorly understoodSUP3SUP. The study of multiplanetary systems offers an opportunity to investigate the outcomes of planet formation and evolution while controlling for initial conditions and environment. Those in resonance with their orbital periods related by a ratio of small integers are particularly valuable because they imply a system architecture practically unchanged since its birth. Here we present the observations of six transiting planets around the bright nearby star HD 110067. We find that the planets follow a chain of resonant orbits. A dynamical study of the innermost planet triplet allowed the prediction and later confirmation of the orbits of the rest of the planets in the system. The six planets are found to be subNeptunes with radii ranging from 1.94RSUBSUB to 2.85RSUBSUB. Three of the planets have measured masses yielding low bulk densities that suggest the presence of large hydrogendominated atmospheres.
['Astrophysics - Earth and Planetary Astrophysics']
A resonant sextuplet of subNeptunes transiting the bright star HD 110067
{'A resonant sextuplet of sub-Neptunes transiting the bright star HD 110067': "R. Luque 1 , ∗ , H. P. Osborn 2 , 3 , † , A. Leleu 4 , 2 , † , E. Pall'e 5 , 6 , † , A. Bonfanti 7 , O. Barrag'an 8 , T. G. Wilson 9 , 10 , 11 , C. Broeg 2 , 12 , A. Collier Cameron 9 , M. Lendl 4 , P. F. L. Maxted 13 , Y. Alibert 12 , 2 , D. Gandolfi 14 , J.-B. Delisle 4 , M. J. Hooton 15 , J. A. Egger 2 , G. Nowak 16 , 5 , 6 , M. Lafarga 10 , 11 , D. Rapetti 17 , 18 , J. D. Twicken 17 , 19 , J. C. Morales 20 , 21 , I. Carleo 5 , 22 , J. Orell-Miquel 5 , 6 , V. Adibekyan 23 , 24 , R. Alonso 5 , 6 , A. Alqasim 25 , P. J. Amado 26 , D. R. Anderson 10 , 11 , G. Anglada-Escud'e 20 , 21 , T. Bandy 27 , T. B'arczy 28 , D. Barrado Navascues 29 , S. C. C. Barros 30 , 31 , W. Baumjohann 7 , D. Bayliss 10 , J. L. Bean 1 , M. Beck 4 , T. Beck 2 , W. Benz 2 , 12 , N. Billot 4 , X. Bonfils 32 , L. Borsato 33 , A. W. Boyle 34 , A. Brandeker 35 , E. M. Bryant 25 , 10 , J. Cabrera 36 , S. Carrazco Gaxiola 37 , 38 , 39 , D. Charbonneau 40 , S. Charnoz 41 , D. R. Ciardi 34 , W. D. Cochran 42 , K. A. Collins 40 , I. J. M. Crossfield 43 , Sz. Csizmadia 36 , P. E. Cubillos 22 , 7 , F. Dai 44 , 34 , M. B. Davies 45 , H. J. Deeg 5 , 6 , M. Deleuil 46 , A. Deline 4 , L. Delrez 47 , 48 , O. D. S. Demangeon 30 , 31 , B.-O. Demory 12 , 2 , D. Ehrenreich 4 , 49 , A. Erikson 36 , E. Esparza-Borges 5 , 6 , B. Falk 50 , A. Fortier 2 , 12 , L. Fossati 7 , M. Fridlund 51 , 52 , A. Fukui 53 , 5 , J. Garcia-Mejia 40 , S. Gill 10 , M. Gillon 47 , E. Go ff o 14 , 54 , Y. Gomez Maqueo Chew 37 , M. Gudel 55 , E. W. Guenther 54 , M. N. Gunther 27 , A. P. Hatzes 54 , Ch. Helling 7 , K. M. Hesse 3 , S. B. Howell 17 , S. Hoyer 46 , K. Ikuta 56 , K. G. Isaak 27 , J. M. Jenkins 17 , T. Kagetani 56 , L. L. Kiss 57 , 58 , T. Kodama 53 , J. Korth 59 , K. W. F. Lam 36 , J. Laskar 60 , D. W. Latham 40 , A. Lecavelier des Etangs 61 , J. P. D. Leon 56 , J. H. Livingston 62 , 63 , 64 , D. Magrin 33 , R. A. Matson 65 , E. C. Matthews 66 , C. Mordasini 2 , 12 , M. Mori 56 , M. Moyano 67 , M. Munari 68 , F. Murgas 5 , 6 , N. Narita 53 , 62 , 5 , V. Nascimbeni 33 , G. Olofsson 35 , H. L. M. Osborne 25 , R. Ottensamer 55 , I. Pagano 68 , H.", 'The HD 110067 planetary system': "layer (H and He). The equation of state (EOS) of water is the one of ref.[144], the core EOS is the one of ref.[145], and we use for EOS of ref.[146] for the silicate mantle. The thickness of the gas envelope, which depends on the planetary age, mass, etc., is derived from ref.[147]. Note that the influence of the gas layer on the innermost planet (compression and thermal e ff ect) is not included in our model, as the mass of gas layer for the six planets is small (see below). The planetary Si / Mg / Fe molar ratio in all planets is assumed to be equal to the stellar one. The prior distribution of the mass fractions of the three innermost layers (core, mantle, and water layer) is assumed to be uniform on the simplex - the surface defined by the sum of the three mass fractions equal to one. In addition, the mass fraction of the water layer is assumed to be 50% at most [148, 149], and for the gas mass, we use a uniform log prior.  \nThe results from this analysis are shown in Extended Data Fig. 6. Our model shows that the gas mass content of all planets is of the order of 10 -3 M ⊕ to 10 -1 M ⊕ (median value, see Table S5), with the notable exception of HD 110067 e (median value of ∼ 10 -7 M ⊕ using the masses from M ethod I, ∼ 10 -3 M ⊕ using M ethod II). The apparent lack of an atmosphere of planet e (located just outside of planet d, which is the most gas-rich of the system, according to the internal structure models) is puzzling. If confirmed by future better determination of its density, the origin of the peculiar internal structure of planet e will have to be understood in the context of the very fragile architecture of the whole HD 110067 system. On the other hand, the water fraction for all planets is essentially unconstrained, due to the still large uncertainty in the planetary masses. However, according to simulations of combined planetary formation and evolution, independently of the accretion mechanism (planetesimal- or pebble-based) all the planets in the system have masses and radii consistent with a formation beyond the ice line [150-152]. Therefore, it is possible that even though the water content is unconstrained in our model, the planets' cores are rich in volatiles. JWST observations of some atmospheric trace gases (particularly ammonia, methane, and / or methanol) could be used as a proxy for the presence of a deep or shallow surface that could break the degeneracies from internal composition models using bulk density measurements alone [153, 154].", 'Bern, Switzerland': "- 3 Department of Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA\n- 4 Observatoire Astronomique de l'Universit'e de Gen'eve, Chemin Pegasi 51, 1290 Versoix, Switzerland\n- 5 Instituto de Astrofisica de Canarias, Via Lactea s / n, 38200 La Laguna, Tenerife, Spain\n- 6 Departamento de Astrofisica, Universidad de La Laguna, Astrof'ısico Francisco Sanchez s / n, 38206 La", 'Laguna, Tenerife, Spain': '- 9 Centre for Exoplanet Science, SUPA School of Physics and Astronomy, University of St Andrews, North', 'Haugh, St Andrews KY16 9SS, UK': '- 10 Department of Physics, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK\n- 11 Centre for Exoplanets and Habitability, University of Warwick, Coventry, CV4 7AL, UK\n- 12 Center for Space and Habitability, University of Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland\n- 13 Astrophysics Group, Lennard Jones Building, Keele University, Sta ff ordshire, ST5 5BG, United Kingdom  \n- 15 Cavendish Laboratory, JJ Thomson Avenue, Cambridge CB3 0HE, UK  \n16  \nInstitute of Astronomy, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University,', 'Grudzicadzka 5, 87-100 Torun, Poland': '- 17 NASA Ames Research Center, Mo ff ett Field, CA 94035, USA\n- 18 Research Institute for Advanced Computer Science, Universities Space Research Association, Washington,', 'DC 20024, USA': "19 SETI Institute, Mountain View, CA 94043, USA  \n- 20 Institut de Ciencies de l'Espai (ICE, CSIC), Campus UAB, Can Magrans s / n, 08193 Bellaterra, Spain\n- 21 Institut d'Estudis Espacials de Catalunya (IEEC), Gran Capit'a 2-4, 08034 Barcelona, Spain", '4 The HD 110067 planetary system': 'Extended Data Fig. 4 Results from the ground-based campaign to detect HD 110067 f . a, ∆ WAIC for each of the constrained period bins when compared to a transit-free model. b,c, Best-fit decorrelated photometry with ( b ) and without ( c ) a transit model. Each light curve from each telescope has been o ff set for clarity. Error bars represent 1 σ uncertainties.  \n<!-- image -->  \nExtended Data Fig. 5 Results from the two radial velocity analyses to measure the mass of each of the planets in the HD 110067 system . Each histogram represents the posterior density function (pdf) of the radial velocity semi-amplitudes as inferred from M ethod I (red) and M ethod II (blue). The area underneath each histogram is normalized to unity.  \n<!-- image -->  \nExtended Data Fig. 6 Gas mass fraction of the HD 110067 planets as a function of their equilibrium temperature . We infer two values per planet by assuming the di ff erent planetary masses from our M ethod I (red) and M ethod II (blue) radial velocity analyses. The boxes, orange lines, green triangles, and red stars represent respectively the 25% and 75% percentiles, medians, means, and modes of the posterior distributions. The opacity of the vertical lines is proportional to the posterior distribution.  \n<!-- image -->  \nExtended Data Table 1 CHEOPS observing log. Filler observations aim to catch transits serendipitously in between time-critical observations with higher priority. Boldface notes indicate that a transit event was detected in the data.', 'References': '- [27] Stassun, K. G. et al. The TESS Input Catalog and Candidate Target List. Astron. J. 156 (3), 102 (2018).\n- [28] Stumpe, M. C. et al. Kepler Presearch Data Conditioning I-Architecture and Algorithms for Error Correction in Kepler Light Curves. Proc. Acad. Sci. Pacific 124 (919), 985 (2012).\n- [29] Stumpe, M. C. et al. Multiscale Systematic Error Correction via WaveletBased Bandsplitting in Kepler Data. Proc. Acad. Sci. Pacific 126 , 100 (2014).\n- [30] Smith, J. C. et al. Kepler Presearch Data Conditioning II - A Bayesian Approach to Systematic Error Correction. Proc. Acad. Sci. Pacific 124 , 1000 (2012).\n- [31] Jenkins, J. M. The Impact of Solar-like Variability on the Detectability of Transiting Terrestrial Planets. Astrophys. J. 575 , 493-505 (2002).\n- [32] Jenkins, J. M. et al. Radziwill, N. M. & Bridger, A. (eds) Transiting planet search in the Kepler pipeline . (eds Radziwill, N. M. & Bridger, A.) Software and Cyberinfrastructure for Astronomy , Vol. 7740 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series , 77400D (2010).\n- [33] Jenkins, J. M. et al. Kepler Data Processing Handbook: Transiting Planet Search. Kepler Science Document KSCI-19081-003 (2020).\n- [34] Twicken, J. D. et al. Kepler Data Validation I - Architecture, Diagnostic Tests, and Data Products for Vetting Transiting Planet Candidates. Proc. Acad. Sci. Pacific 130 (6), 064502 (2018).\n- [35] Li, J. et al. Kepler Data Validation II-Transit Model Fitting and Multipleplanet Search. Proc. Acad. Sci. Pacific 131 (996), 024506 (2019).\n- [36] Guerrero, N. M. et al. The TESS Objects of Interest Catalog from the TESS Prime Mission. Astrophys. J. Suppl. 254 (2), 39 (2021).\n- [37] Fausnaugh, M. M., Burke, C. J., Ricker, G. R. & Vanderspek, R. Calibrated Full-frame Images for the TESS Quick Look Pipeline. Res. Not. Amer. Astron. Soc. 4 (12), 251 (2020).\n- [38] Hedges, C. et al. TOI-2076 and TOI-1807: Two Young, Comoving Planetary Systems within 50 pc Identified by TESS that are Ideal Candidates for Further Follow Up. Astron. J. 162 (2), 54 (2021).\n- [39] Osborn, H. et al. Two Warm Neptunes transiting HIP 9618 revealed by TESS &Cheops. Mon. Not. R. Astron. Soc. 523 (2), 3069-3089 (2023).', '14 The HD 110067 planetary system': "- [15] Quirrenbach, A. et al. Evans, C. J., Bryant, J. J. & Motohara, K. (eds) The CARMENES M-dwarf planet survey . (eds Evans, C. J., Bryant, J. J. & Motohara, K.) Ground-based and Airborne Instrumentation for Astronomy VIII , Vol. 11447, 114473C. International Society for Optics and Photonics (SPIE, 2020).\n- [16] Cosentino, R. et al. McLean, I. S., Ramsay, S. K. & Takami, H. (eds) Harps-N: the new planet hunter at TNG . (eds McLean, I. S., Ramsay, S. K. & Takami, H.) Ground-based and Airborne Instrumentation for Astronomy IV , Vol. 8446, 84461V. International Society for Optics and Photonics (SPIE, 2012).\n- [17] Holman, M. J. & Murray, N. W. The Use of Transit Timing to Detect Terrestrial-Mass Extrasolar Planets. Science 307 (5713), 1288-1291 (2005).\n- [18] Fulton, B. J. et al. The California-Kepler Survey. III. A Gap in the Radius Distribution of Small Planets. Astron. J. 154 , 109 (2017).\n- [19] Van Eylen, V. et al. An asteroseismic view of the radius valley: stripped cores, not born rocky. Mon. Not. R. Astron. Soc. 479 , 4786-4795 (2018).\n- [20] Kasting, J. F., Whitmire, D. P. & Reynolds, R. T. Habitable Zones around Main Sequence Stars. Icarus 101 (1), 108-128 (1993).\n- [21] Kopparapu, R. K. et al. Habitable Zones around Main-sequence Stars: Dependence on Planetary Mass. Astrophys. J. 787 , L29 (2014).\n- [22] Izidoro, A. et al. Formation of planetary systems by pebble accretion and migration. Hot super-Earth systems from breaking compact resonant chains. Astron. Astrophys. 650 , A152 (2021).\n- [23] Fabrycky, D. C. et al. Architecture of Kepler's Multi-transiting Systems. II. New Investigations with Twice as Many Candidates. Astrophys. J. 790 (2), 146 (2014).\n- [24] Zeng, L. et al. Growth model interpretation of planet size distribution. Proc. Nat. Acad. Sci. USA 116 (20), 9723-9728 (2019).\n- [25] Kempton, E. M. R. et al. A Framework for Prioritizing the TESS Planetary Candidates Most Amenable to Atmospheric Characterization. Proc. Acad. Sci. Pacific 130 (993), 114401 (2018).\n- [26] Otegi, J. F., Bouchy, F. & Helled, R. Revisited mass-radius relations for exoplanets below 120 M ⊕ . Astron. Astrophys. 634 , A43 (2020).", '0.1.1 TESS photometry': "HD 110067 was observed in sectors S23 and S49 (2020 March 18 to 2020 April 16 and 2022 February 26 to 2022 March 26, respectively) of the Transiting Exoplanet Survey Satellite (TESS) [4]. The star was included on the TESS Candidate Target List (CTL) [27] and therefore the target was observed at 120 s cadence in both sectors. The target images are processed by the TESS Science Processing Operations Center (SPOC) pipeline at NASA Ames [5], which calibrates the pixels, performs simple aperture photometry (SAP), flags poor-quality data, and removes systematic trends to create the so-called 'Presearch Data Conditioning' light curve (PDCSAP) [2830]. Finally, the SPOC pipeline runs a wavelet-based transiting planet search for periodic exoplanets [31-33] which in the case of the S23 data revealed two thresholdcrossing events (TCEs) - i.e. candidate planets - which passed Data Validation checks [34, 35]. Manual vetting of these two candidates resulted in the assignment of two TESS Object of Interest (TOI), namely TOI-1835.01 and TOI-1835.02 [36].  \nThe first TCE was alerted by SPOC as TOI-1835.01 with a period of 5.641 d. Although two of the three transits associated with this ephemeris appeared to be of similar depth and duration, the third did not appear to be associated with a clear transit. However, this could have been due to proximity to a systematic dip caused by a momentum dump. A second TCE, TOI-1835.02, was alerted as a single transit at 1948 . 98 TJD TESS Julian Date (TJD ≡ BJD -2457000), where BJD is the Barycentric Julian Date in units of days, as opposed to the period of 11.107 d proposed by SPOC due to an apparent discrepancy in the transit depth between the two transits purportedly linked by the transiting planet search. On the other hand, a peak in the background flux ruled out the planetary origin of an apparent transit feature at TJD = 1940.5. Both the depth and shape of this feature were strongly dependent on the detrending method used in the lightcurve, which is indicative of its spurious nature.  \nHD 110067 was later reobserved by TESS in S49. In order to confirm the TCEs and plan immediate follow-up, we downloaded the TESS Image CAlibrator Full Frame Images (TICA FFI) [37] only a week after downlink for each of the first and then second orbits. We computed a light curve from the 10-minute cadence TICA FFI cut-outs using SAP and clipped regions of high brightness due to the Earth and Moon, as well as parts of the light curve a ff ected by systematics such as momentum dumps. The first orbit alone revealed at least 5 new clear transit features. These all appeared to have varying depths and durations, and none were compatible with the 5.64 d period implied by TOI-1835.01. The second orbit showed three more clear transit events, making a total of eight in S49 and five in S23. A 20-second cadence target pixel file was later made available for HD 110067, which resulted in higher-precision photometry.", '0.1.2 CHEOPS photometry': "The CHaracterising ExOPlanets Satellite (CHEOPS) mission is a European Space Agency small-class mission dedicated to studying bright, nearby exoplanet host stars for the purpose of making high-precision photometric observations of transiting planets [6]. We collected 19 separate visits of HD 110067 with CHEOPS between 2022 April 11 and 2022 May 17 under Guaranteed Time Observing programs ID-048  \nand ID-031. The goal of these observations is (1) to confirm the true orbital period of single- and duo-transiting planet candidates and (2) to improve the planetary radius precision and ephemeris of confirmed planets. This has been done for large planets producing deep eclipses from the ground [47, 48], and for small planets from space [39, 49, 50]. An observing log summarizing the duration of each visit, its average observing e ffi ciency (considering the gaps produced by Earth occultations or passages over the South Atlantic Anomaly along the spacecraft's low-Earth orbit), and photometric precision are presented in Extended Data Table 1.  \nTo provide the highest quality photometric precision, we opted to perform custom photometric extraction of the CHEOPS imagettes using point-spread-function (PSF) photometry as implemented by the PIPE package [51, 52]. For bright targets such as HD 110067, light curves generated with PIPE exhibit lower median absolute di ff erences than those generated by the CHEOPS Data Reduction Pipeline [53]. The shorter cadence of the CHEOPS imagettes allows a higher cadence light curve, and PSF detrending is also better at removing trends due to systematic factors and background stars. As various PSF models have already been generated and vary as a function of stellar temperature, we opted to use a PSF generated using the star HD 189733, with a similar spectral type of HD 110067. In order to preserve intervisit flux di ff erences, we normalized the entire CHEOPS data together instead of individually. This revealed clear visit-to-visit flux di ff erences due to stellar rotation with an amplitude larger than that of TESS (as stellar activity is typically more pronounced at bluer bandpasses). The final light curves used in our subsequent analyses are shown in Fig. 1 and Fig. S1 in the Supplementary.", '0.1.3 Ground-based photometric campaign': 'We carried out a campaign on the night of May 23rd, 2022 to attempt to confirm the 41.05-day period orbit of HD 110067 f as predicted by our resonance chain analysis. Photometric observations were taken using 14 telescopes using seven di ff erent filters, which observed from various locations to continuously cover a temporal baseline of more than 11 hours (between 2022-05-23UT22:52:55 and 2022-05-24UT10:01:33). This window is long enough to catch the 5-hour transit expected from 02:52 to 07:12UT. However, no single location was able to cover both ingress and egress. A summary of the observations is shown in Extended Data Table 2. Details from each individual observation are shown below. Extended Data Fig. 4 shows the data and best-fit models as discussed in Sect. 0.3.4.', 'Teide Observatory': "We observed HD 110067 on May 23rd, 2022 using the MuSCAT2 instrument installed at the 1.5-m Telescopio Carlos S'anchez (TCS) located at the Teide Observatory, Spain [54]. The images were taken simultaneously in g , r , i , and zs filters with the telescope heavily defocused and with short exposure times of 3 to 5 sec, depending on the band, to avoid saturation. Relative light curves for each band and instrument of HD 110067 were extracted by aperture photometry using a custom pipeline [55] with optimal aperture radii of 8 . '' 1 to 11 . '' 3 depending on the band. Note  \nthat there was a technical problem on the dome of the TCS between BJD-2459723 = 0.488 and 0.526; we discarded the data taken during this period.  \nWe also observed HD 110067 on May 23rd, 2022 with one of the 1-m telescopes from Las Cumbres Observatory (LCO) global network located at the Teide Observatory, Spain [56]. The observations were obtained through Director's Discretionary Time program 2022A-005 (PI: Wilson). We collected 181 frames with an exposure time of 20 s, covering 2.5 h, using the 4096 × 4096 pix SINISTRO camera. The images were calibrated by the standard LCO BANZAI pipeline [57]. Di ff erential photometric data were extracted using AstroImageJ (AIJ) [58].", 'Paranal Observatory': "We observed HD 110067 on May 23rd, 2022 using the Next Generation Transit Survey (NGTS) facility located at ESO's Paranal Observatory in Chile [59]. NGTS consists of twelve 20-cm, f / 2.8 telescopes with Andor cameras and red-sensitive (600-900 nm) deep-depletion e2v CCDs. Nine NGTS telescopes observed from 23:14 to 04:35UT, covering a predicted transit ingress of HD 110067 f, and spanning an airmass range of 1.7-2.5. Two telescopes started observing two hours late due to a technical issue. All nine telescopes were defocused to avoid saturating the bright target star during the 10-second exposures. The NGTS camera shutters were not functional and so were kept open during the entire observing block. That caused the stars to streak during the 1.5-second readout sequences but without any apparent detrimental e ff ect on the photometry. Observing without using the shutters is now the standard operation mode of NGTS. We performed standard di ff erential aperture photometry, using large aperture radii of 6.5-8.0 pixels, and carefully selecting comparison stars to avoid those that exhibited variability. The light curve of each telescope was normalized individually and no detrending was performed.", 'F. L. Whipple Observatory': 'We observed HD 110067 on May 24th, 2022 using the Tierras instrument installed at the refurbished 1.3-m telescope located at the F. L. Whipple Observatory atop Mount Hopkins, Arizona, United States. The instrument is designed to regularly achieve a photometric precision of 250 ppm on a time scale of both 10 min and a complete observing season. The design choices that permit this precision include a four-lens focal reducer and field-flattener that increase the field-of-view of the telescope, a custom narrow bandpass filter centered around 863.5 nm to minimize precipitable water vapor errors, and a fully automated mode of operation [60]. A total of 1262 4-second exposures were gathered with Tierras for HD 110067. Astrometric calibrations were done in real-time during data gathering and were stored in WCS headers in the FITS files. The FITS files were then passed through the Tierras image reduction pipeline to perform bias corrections and image stitching (the CCD chip is read out through separate amplifiers). AIJ was used for photometric extraction. These observations were gathered shortly after Tierras started science operations, and the data were not flatfielded since knowledge of the flat-field was incomplete at the time. The RMS of the  \n15-minute binned data is 323 ppm. The photometric precision on this target is ultimately limited by scintillation, as the target was observed down to an airmass of 2.37. The observations were mildly a ff ected by cirrus.', "San Pedro M'artir Observatory": "Weobserved HD 110067 on May 24th, 2022 with the 1-m SAINT-EX telescope at the Observatorio Astron'omico Nacional de la Sierra de San Pedro M'artir in Baja California, Mexico [61]. SAINT-EX is equipped with a deep-depleted and back-illuminated Andor IKON CCD and a filter wheel. The observations were defocused and acquired in the 'zcut' filter, a custom filter optimized to reduce the systematic uncertainties in the light curves of red stars due to precipitable water vapor, with an exposure time of 10 s. The data were reduced with AIJ using the standard corrections for bias, flatfielding, and dark current. AIJ was also utilized to do the aperture photometry of the time series, producing the light curves and relevant meta-data. The observations were mildly a ff ected by high-altitude cirrus.", 'Haleakala Observatory': "We observed HD 110067 on May 24th, 2022 using the MuSCAT3 instrument mounted on the 2-m Faulkes Telescope North (FTN) at Haleakala Observatory on Maui, Hawaii, United States [62]. The images were taken simultaneously in g , r , i , and zs filters with the telescope heavily defocused and with short exposure times of 3 to 5 sec, depending on the band, to avoid saturation. Relative light curves for each band and instrument were extracted by aperture photometry using a custom pipeline [63] with optimal aperture radii of 8 . '' 1 to 11 . '' 3 depending on the band. There was a guiding issue on the FTN around BJD-2459723 = 0.795, which caused a large shift of the stellar positions on the detectors; we treated the MuSCAT3 data as two independent datasets separated by that time.", '0.1.4 High-resolution imaging': "As part of our standard process for validating transiting exoplanets, and to assess the possible contamination of bound or unbound companions on the derived planetary radii [64], we observed HD 110067 with near-infrared (NIR) adaptive optics (AO) imaging at Palomar Observatory and with optical speckle imaging at Gemini North. Gaia DR3 is also used to provide additional constraints on the presence of undetected stellar companions and wide companions. No close-in ( ≲ 1 '' ) stellar companions were detected by either the NIR adaptive optics or optical speckle imaging.", 'Palomar Observatory': "The Palomar Observatory observations of HD 110067 were made with the PHARO instrument [65] behind the natural guide star AO system P3K [66] on 2020 Jan 08 in a standard 5-point quincunx dither pattern with steps of 5 '' in the narrow-band Brγ filter. Each dither position was observed three times, o ff set in position from each other by 0.5 '' for a total of 15 frames; with an integration time of 1.4 seconds per frame, respectively for total on-source times of 21 seconds. PHARO has a pixel scale  \nof 0 . 025 '' per pixel for a total field of view of ∼ 25 '' . The sensitivities of the final combined AO image were determined by injecting simulated sources azimuthally around the primary target every 20 · at separations of integer multiples of the central source's FWHM [67]. The Palomar data have a sensitivity ∆ mag = 2 at 0.1 '' and ∆ mag = 9 at 1 '' ; the final sensitivity curve is shown in Fig. S2 of the Supplementary.", 'Gemini Observatory': "Weobserved HD 110067 with the 'Alopeke speckle imaging camera at Gemini North on 2020 June 10 [68]. We obtained five sets of 1000 frames, each frame having an integration time of 60 ms, obtaining images in each of the instrument's two bands (centered at 562 nm and 832 nm). The observations were reduced using our standard software pipeline [69] and reached a 5 σ sensitivity of ∆ mag = 7 (blue channel) and ∆ mag = 6 . 8 (red channel) at separations of 0.5 '' . The reconstructed speckle images show no evidence of additional nearby point sources. The final sensitivity curve is shown in Fig. S2 of the Supplementary.", 'Gaia Space Observatory': 'In addition to the high-resolution imaging, we have utilized Gaia to identify any wide stellar companions that may be bound members of the system [70, 71]. There are no additional widely separated companions identified by Gaia that have the same distance and proper motion as HD 110067. Additionally, the Gaia DR3 astrometry provides additional information on the possibility of inner companions that may have gone undetected by either Gaia or the high-resolution imaging data. The Gaia Renormalised Unit Weight Error (RUWE) is a metric, similar to a reduced chi-square, where values that are ≲ 1 . 4 indicate that the Gaia astrometric solution is consistent with the star being single whereas RUWE values ≳ 1 . 4 may indicate an astrometric excess noise, possibly caused the presence of an unseen companion [e.g., 72]. HD 110067 has a Gaia EDR3 RUWE value of 0.94 indicating that the astrometric fit is consistent with a single-star model.', 'Calar Alto Observatory': "We observed HD 110067 using the CARMENES instrument [15] installed at the 3.5m telescope of Calar Alto Observatory in Almer'ıa, Spain, between 3 July 2020 and 4 July 2021. We collected 39 high-resolution spectra under the observing programs F20-3.5-011 (PI: Nowak) and H20-3.5-013 (PI: Luque). Radial velocities and additional spectral indicators were derived using raccoon [73] and serval [74]. While the mean internal precision of the template matching serval RVs is 3 . 1 m s -1 , the precision of the cross-correlation method raccoon RVs is 2 . 9 m s -1 , so we used the latter in our analyses.", 'Roque de los Muchachos Observatory': "We observed HD 110067 with the HARPS-N spectrograph mounted at the 3.6 m Telescopio Nazionale Galileo [16] of Roque de los Muchachos observatory in  \nLa Palma, Spain, between 30 May 2020 and 4 May 2022. We collected 72 high-resolution spectra under the observing programs CAT19A 162 (PI: Nowak), CAT21A 119 (PI: Nowak) and ITP19 1 (PI: Pall'e) that were used to measure the photospheric properties of the star and precise radial velocities. Radial velocities and additional spectral indicators were derived using an online version of the DRS pipeline [75], the YABI tool, and serval [74]. Both the YABI- and serval -derived radial velocities have a median internal precision of 1 . 0 m s -1 , but we used the YABI ones (based on the cross-correlation method) in our final analyses for consistency with the CARMENES dataset.", '0.2.1 Photospheric parameters and abundances': 'To properly characterize the planetary system around HD 110067, we first conduct a series of analyses to determine the properties of the host star. We derive the stellar spectral parameters by applying the widely used ARES+MOOG tools to our co-added HARPS-N spectra [76-78]. ARES [79, 80] measures the equivalent widths of iron lines in the spectrum that are converted into stellar atmospheric parameters using the MOOG radiative transfer code [81] applied to Kurucz model atmospheres [82]. In Extended Data Table 3 we report the e ff ective temperature T e ff , surface gravity log g , and metallicity [Fe / H], obtained upon convergence of ionization and excitation equilibria within this method. Additionally, we measure the stellar v sin i from the HARPS-N spectra using ZASPE [83].  \nWe further study the photospheric parameters by conducting a classical curveof-growth analysis on our co-added HARPS-N spectrum using our aforementioned spectral parameters in order to obtain [Mg / H] and [Si / H] abundances for HD 110067. Utilizing the ARES+MOOG framework detailed above, we obtain the equivalent widths [80] for these elements, which are converted to abundances assuming local thermodynamic equilibrium [81, 82]. The specific details of this analysis are beyond the scope of this paper and can be found in refs.[84, 85]. We report the stellar abundances in Extended Data Table 3.', '0.2.2 Physical parameters': "Using our spectral parameters and the ATLAS [82, 86] and PHOENIX [87] catalogs, we build spectral energy distributions of HD 110067 that we compare to optical and infrared broadband photometry of the star (see Extended Data Table 3) to derive the stellar angular diameter and e ff ective temperature via the infrared flux method [88]. This is conducted in an MCMC approach [89, 90] within which we convert the angular diameter to the stellar radius using the Gaia EDR3 o ff set-corrected parallax [91] with model uncertainties accounted for using a Bayesian modeling averaging. We report the stellar radius R ⋆ in Extended Data Table 3.  \nLast, we complete our stellar characterization by determining the mass and age of HD 110067. We constrain two sets of stellar evolutionary models with help of our derived values for T e ff , log g , and R ⋆ [92]. On the one hand, we  \nuse an isochrone placement algorithm [93, 94] and interpolate over pre-computed grids of PARSECv1.2S [95] isochrones. On the other hand, we use the Code Li'egeois d' ' Evolution Stellaire [96] combined with a Levenberg-Marquadt minimization scheme [97] to optimize the best-fitting evolutionary track. The results from the two methods are combined to determine the mass and age of the star that is reported in Extended Data Table 3.  \nThe [Fe / H] and age of HD 110067 indicate that this star could belong either to the galactic thick disk stellar population or be an older member of the galactic thin disk. The values of [Mg / H] and [Si / H], being within 1 σ of [Fe / H], show that the star is not enhanced in α -capture elements and are indicative of a typical thin disk chemical composition. We determined the kinematic properties of HD 110067 by using the Gaia EDR3 astrometry to compute the Local Standard of Rest space velocities of this star following ref.[98]. From these velocities, we compute that the probability of kinematic membership in the galactic thin disk is 0 . 9911 ± 0 . 0029. Thus, we conclude that HD 110067 is on the older, more metal-poor end of the distribution of the galactic thin disk stellar population.", '0.3.1 Space-based photometry modeling': 'We performed simultaneous modeling of the space-based photometry. We used the quaternion-detrended TESS data combined with the PLD-detrended data for the missing S23 gaps, and the PIPE -detrended CHEOPS data for the three visits containing transits. We built transit models for the six planets with exoplanet [99]. Due to its nature as a rotating telescope on a near-Earth orbit, even PSF-detrended CHEOPS photometry can include systematic trends. However, these typically correlate with other measurements, for example, roll angle, background, or contaminant flux. In order to not bias the transit model and to better propagate uncertainties on the derived parameters, we performed CHEOPS decorrelation alongside our photometric transit modeling. We first fitted each CHEOPS transit individually alongside multiple possible decorrelation factors, allowing us to assess which decorrelation factors are most useful. This also enabled us to test whether such decorrelation is shared among all CHEOPS visits or individual to a single light curve. From this analysis, we included the following parameters in the linear correlation: position centroids, the second harmonic of the cosine of the roll-angle, cos 2 Φ , the change in telescope temperature, and quadratic trends with the x-y centroids. CHEOPS data have also been known to contain flux trends that vary stochastically as a function of roll angle over shorter frequencies [see e.g. 100]. These are not well removed using simple trigonometric functions, hence we also modeled a flexible spline shared between all visits to model shorter-timescale variation. To incorporate stellar variability, a floating mean and flux trend were also fitted to each CHEOPS visit, as well as an individual jitter term.  \nInformative priors were used on limb darkening parameters using the theoretical quadratic limb darkening parameters for TESS [101] and CHEOPS [102], with uncertainties inflated to 0.1 in all cases to guard against systematic o ff sets. The impact parameter and radius ratio are fitted from a broad uniform and log-normal prior,  \nrespectively, while the period and mid-transit epoch are fitted using broad normal priors from the transits identified and modeled above. Stellar parameters from Extended Data Table 3 were used as inputs to the model with Gaussian priors. Orbits were assumed circular in all cases which is a good approximation for planetary systems with multiple transiting planets [103-105]. The prior and posterior distributions of each parameter in the model are shown in Table S1 of the Supplementary.', '0.3.2 Properties of the unmatched transits': "Our first modeling of the TESS space-based photometry was able to account for a total of 5 transits of planet b (two in TESS S23 and three in TESS S49) and 4 transits of planet c (two each in TESS S23 and S49). However, this analysis left six 'unmatched' transits in the original TESS light curves. In order to pair the transits, we fitted each transit individually using a purely shape-based transit model agnostic to the orbital period using MonoTools [106]. From this analysis, we then compared each transit in duration-depth space, allowing us to clearly see that both transits from S23 shared unique regions of this parameter space with two more transits seen in S49 (duo-transits), while the two longest-duration transits seen only in S49 were solitary (mono-transits). Extended Data Fig. 1 shows this result.  \n.  \nWe then modeled both duo- and single-transits using MonoTools fitting. This allows long-period planets to be modeled in a way that the transit model is agnostic of the orbital period with the implied period distributions being manipulated using priors. This technique works for transits with single- or duo-transits. In the case of two transit events separated by a long gap, the planetary transit is fitted leaving the orbital period open, and the implied transit shape is used to calculate the probability for each of the possible period aliases. For single transits, potential orbital period windows are computed. In both cases, the period probability distribution comes from a combination of a simple period prior (longer period planets are geometrically disfavored) [107], an eccentricity prior (eccentric orbits are disfavored in multi-transiting systems) [108], and a stability prior using the orbits of other planets in the system (orbit-crossing is disallowed) [further details in 39, 109]. The resulting marginalized period predictions for planets HD 110067 d, e, f, and g are shown in Fig. S3, with posterior values of 21 . 6 + 2 . 9 -1 . 6 , 29 . 9 + 4 . 6 -3 . 3 , 40 . 1 + 7 . 1 -5 1 , and 47 . 0 ± 8 . 0 days, respectively.", '0.3.3 Continuing the resonant chain': 'In this section, we expand the analyses that led to the prediction of the orbits of planets HD 110067 e, f, and g based on the generalized Laplace resonant configuration of the three inner transiting planets in the system. We assume that all events mentioned in the previous section are transits that belong to planets that continue the resonant chain.  \nFor transiting systems, generalized three-body Laplace angles can be estimated in 0th order in eccentricity, defined as Ψ e = 0, from the times of mid-transit and the orbital period of the planets [see e.g., 110]. This estimation di ff ers from the actual generalized three-body Laplace angle proportionally to the eccentricities [eq. 15 of 111]. Interestingly, for known systems with a chain of three-body resonances, all  \nΨ e = 0 lie close to an equilibrium of the chain, as seen in Extended Data Fig. 2. The largest distance is ∼ 43 degrees for the inner triplet of K2-138 [112]. For HD 110067, the estimated angle Ψ e = 0 , bcd is also at about ∼ 44 degrees from its theorized 180degree equilibrium. Through the study of transit timing variations over several years, one can get constraints on the underlying generalized three-body Laplace angles. In known cases, one can see that these angles oscillate with amplitudes of a few tens of degrees at most around their equilibrium value, see Fig. 2 of [12] for Kepler-60, Fig. 25 of [13] for TRAPPIST-1.  \nAs shown above, the two events at 2646.088 TJD and 1937.851 TJD have fully consistent shapes. Among the probable periods computed with MonoTools (Fig. S3), Pe = 30 . 7931 days is the only one that continues the resonant chain, with Pe / Pd = 1 . 5007, landing inside the common 3:2 MMR (see Fig. S4 in the Supplementary). We compute the observed value of the associated generalized three-body Laplace angle Ψ e = 0 , cde = 169 . 995 deg, which is at only 10 deg from the expected 180-degree equilibrium. We hence predict a period of 30 . 7931 days for planet HD 110067 e if it is in the resonant chain.  \nFor the remaining two mono-transits, we try a set of first-order MMRs (2 / 1, 3 / 2, 4 / 3, 5 / 4, 6 / 5) between planet #4 and #5 and the same between planet #5 and #6 (hence 25 combinations). Each of these combinations has to be tested assuming that the transit at TJD = 2641.5778 belongs to the 5th planet, and TJD = 2656.0944 belongs to the 6th planet (case A), and vice-versa (case B). Fortunately, many of these 50 possibilities are excluded by existing data. We end up with 4 possibilities for case A and 9 for case B. As seen in Fig. 2, all known chains of Laplace resonances have either their estimated generalized three-body Laplace angle Ψ e = 0, or their actual generalized three-body Laplace angle Ψ close to an equilibrium of the chain. We will hence favor the configurations that are closest to an equilibrium of the chain. For each case, the distance of each estimated angle to its closest equilibrium ∆Ψ = | Ψ e = 0 -Ψ eq | is given in Extended Data Table 4. The case A2, with Pf / Pe = 4 / 3 and P g/ Pf = 4 / 3, comes out as a favorite, with the three outer generalized three-body Laplace angles at less than 20 deg from the closest equilibrium. In addition, one can note that 4 / 3 MMRs are relatively common in resonant chains (see Fig. S4 in the Supplementary).  \nFor completeness, we study the role that the eccentricity of the orbit plays in the prediction. To estimate the generalized three-body Laplace angle  \nΨ = l λ 1 -( l + m ) λ 2 + m λ 3 , (1)  \nat a given epoch, we estimate the value of the λ j as follow. Transits occur when the true longitude of the planet is equal to l 0 = -π/ 2. At first order in the eccentricity,  \nλ 0 = l 0 -2 e sin( l 0 -ϖ ) = -π 2 + 2 e cos( ϖ ) . (2)  \nWe then assume the planet to be in a circular, unperturbed orbit to compute the value of its mean longitude at the time of transit t 0, λ 0 = -π/ 2. We hence obtain  \nλ ( t ) = -π 2 + ( t -t 0) P 2 π. (3)  \nThe error on Ψ made by assuming zero eccentricity is hence, at first order [111]:  \n| Ψ -Ψ e = 0 | = | 2 le 1 cos ϖ 1 -2( l + m ) e 2 cos ϖ 2 + 2 me 3 cos ϖ 3 | . (4)  \nThis error can thus be substantial (several tens of degrees) if the eccentricities are of the order of several parts per hundred, as it is the case for Kepler-223 [110]. Therefore, we study if a given combination of the eccentricities and longitudes of periastron can make Ψ e = 0 closer to the equilibrium than Ψ actually is, or vice-versa. We check this for the cases presented in Extended Data Table 4.  \nEach case sets the orbital period of the planets and their mid-transit time. We estimate the planetary masses using the mass-radius relation from [26]. Then, varying the remaining parameters ki = ei cos ϖ i and hi = ei sin ϖ i , we minimize the cost function  \nC = A ( Ψ bcd ) + A ( Ψ cde ) + A ( Ψ e f g ) + A ( Ψ f g h ) (5)  \nover 200 years, where A ( Ψ ) = 2 π if Ψ circulates, and A ( Ψ ) is the peak-to-peak amplitude of libration of Ψ otherwise. For each case, 40 MCMC runs are conducted to minimize C , using REBOUND [113] for the N -body integration and samsam [114] for the MCMC. For each run, the ki and hi parameters are randomly initialized in the [ -0 . 05 , 0 . 05] range, which are also their boundaries during the MCMC runs. This allows eccentricities that are comparable to those of Kepler-60 [115] and Kepler-223 [110], which are other known chains for which the inner planets are far enough from the star to not have their eccentricities damped by tides.  \nThe best solution of each fit is shown in Extended Data Fig. 3. Case A2 is the only one for which the best solutions consistently have a peak-to-peak amplitude of the generalized three-body Laplace angles below 50 degrees on average across the four angles. In all other cases, we were not able to find values of the ki and kj parameters below 85 deg of amplitude on average, with the exception of case A0 for which an average of ≈ 66 deg was reached. The best solutions found across all MCMC runs for the A2 and A0 cases integrated for 1000 years of evolution are shown in Fig. S5 of the Supplementary. This analysis shows that the A2 case remains the one with the highest potential of being close to an equilibrium, while showing that all other cases cannot have an amplitude of libration smaller than 66 deg on average across their generalized three-body Laplace angles, regardless of the values of the ki and kj parameters. The case A2, with Pf / Pe = 4 / 3 and P g/ Pf = 4 / 3, is hence our prediction for the outer architecture of the HD 110067 system.', 'Recovering the missing cadences of TESS S23 observations': 'Based on the dynamical analysis presented above, the likely orbital periods associated with the two mono-transits observed in TESS S49 are approximately 41.05 and 54.74 days, respectively. According to this prediction, both planets transited their host star during TESS S23 observations, but at a time when the photometry was highly a ff ected by scattered light and sky background contamination. The Earth was a significant source of scattered light at the beginning of both S23 orbits (TJD = 1928.09 and TJD = 1941.83) and the Moon was a significant source of scattered light for a few days after the beginning of the second orbit (between 1942 and 1947 TJD). The cadences a ff ected were flagged by SPOC, thus not leaving enough valid data to derive cotrending basis vectors and missing in the PDCSAP light curve.  \nOur custom extraction using the PLD method from Sect. 0.1.1 was able to recover the missing data, showing two mono-transits at 1943.6 and 1944.1 TJD (Fig. 1 and also Fig. S6 in the Supplementary). Using MonoTools , we confirmed that the transits were consistent in duration-depth space with the two mono-transits from TESS S49 and separated by an integer number of orbits that matched the orbital periods predicted in our dynamical analysis for planets f and g. With this data reduction, all six planets in the HD 110067 system have been detected in transit at least twice, allowing a precise orbital period determination if we impose priors based on the hypothesis that all planets are trapped in a chain of first-order MMRs. Additionally, we recovered an additional transit of planet b at the beginning of TESS S23 observations.', 'Modeling of the ground-based photometric campaign': 'Targeted observations of HD 110067 were carried out on the night of May 23rd, 2022 to attempt to confirm the 41.05-day orbit of HD 110067 f as predicted via our resonance chain analysis. In order to reveal whether a transit was present in the combined dataset, we built a combined photometric model using all ground-based observations. In order to remove spurious systematic trends in a way that does not bias any transit fit, we opted to perform simultaneous linear decorrelation of each photometric data set using the various meta-data time series available. In all cases, for example, we included an airmass term in the decorrelation as well as a measure of the FWHM width. We also included two position-centroid terms (for MuSCAT-2 and -3), information on comparison star total counts and FWHM width (for LCO, Tierras, and SAINT-EX), and interpolated color time series derived from the relative shift in flux across bands in the MuSCAT-2 and -3 filters [as used in 49]. In all cases, the metadata were normalized to a time series with µ = 0, σ = 1 and modeled using a single scaling parameter with a normal prior of µ = 0, σ = 0 . 5. Quadratic limb darkening parameters were also constrained using normal priors dictated by theoretical limb darkening parameters as computed for each of the nine passbands using LDTK [116] and with inflated uncertainties following the methodology of the space-based photometric analysis. Each of the four time series observed by MuSCAT-3 was split into two around an observing gap that occurred due to the star passing close to the zenith at 2022-05-24UT06:57. Individual time series were used for each of the nine NGTS  \ntelescopes, which were decorrelated independently. The final result is 24 individual photometric time series. An o ff set was also applied to each light curve, as well as a single global slope parameter to include the possibility of stellar activity.  \nThe transit parameters were constrained based on those found in a fit of the TESS S49 mono-transit. The predicted period used was 4 / 3 × Pe = 41 . 051 ± 0 . 1 d, with the uncertainty implying a divergence from the perfect integer period of 2 . 4 × 10 -3 -larger than those values found for the inner three planets. We limited the period to 41 . 0 ± 0 . 2 days to ensure a transit fit that could be explored with the temporal baseline of the photometry. Due to the non-continuous nature of the photometry, the probability density function of the observed transit time is likely to be asymmetric and could potentially have multiple minima. Therefore, analyses using classical sampling techniques (Markov Chain Monte Carlo, Hamiltonian Monte Carlo, etc.) may not reveal the full picture. In order to initially test this, we kept all other parameters equal but split the range of periods covered by the time series into 36 bins across the anticipated period range (40 . 8 < P < 41 . 2) and fitted a constrained model for each. This would allow us to see the variation in the goodness-of-the-fit as a function of transit epoch. The transit model was built with exoplanet [99] and optimized using pymc3 [117], specifically with the pymc3-ext sampling which enabled correlated parameters for each time series to be grouped together, speeding up the computation. In order to assess whether or not a transit model was justified over a flat model, we used the Watanabe-Akaike information criterion (WAIC) [118, 119] as implemented in arviz [120].  \nOur results show a preference for a transit at the expected period P = 41 . 04 ± 0 . 01 days, with a ∆ WAIC of 9.5 over a transit-free model, as can be seen in Extended Data Fig. 4. The majority of instruments showed a weak preference for a peak at P ∼ 41 . 05d, with the exception of LCO (which observed no in-transit data) and SAINT-EX (which is the most a ff ected by cirrus). This is only equivalent to moderate evidence for a ∼ 41 . 051d period of HD110067 f. The lower two panels of Extended Data Fig. 4 show that both models (with or without transit) fit reasonably well. This is in part because systematic e ff ects dominate over astrophysical signals for transits with depths below 1000 ppm, especially when the target star is observed at a low airmass. A further peak in WAIC is seen at ∼ 40 . 9 d, but this hypothetical transit is covered only by the initial 1 h of MuSCAT-2 photometry and does not fit our predicted period, hence we consider it spurious. This campaign shows that such transits are at the very limit of what is possible with ground-based observations. However, observations during a more favorable observing season and without technical issues such as the meridian flip of MuSCAT-3 (which unluckily coincided with the expected egress) may have constrained better the presence of a transit.', 'Modeling of the radial velocity data': "We carried out an initial frequency-based exploration of the CARMENES and HARPS-N spectroscopic data sets to see which significant signals are present and those related to stellar activity using periodograms [121]. Figures S7 and S8 in the Supplementary show that the dominant signal in the generalized Lomb-Scargle periodograms of both the radial velocities and main activity indicators (CCF-FWHM,  \ndi ff erential line width, Mount Wilon's S-index, H α emission) is attributable to the rotational period of the star, measured photometrically to be approximately 20 days using TESS and CHEOPS data. As it is well known, stellar activity induces spurious radial velocity signals [e.g., 122-125], which should be properly removed to unveil induced Keplerian motions in the star. We followed two independent approaches to model the data and minimize the impact of stellar activity e ff ects on the detection and mass determination of the planets in the system.", 'M ethod I: SNfit and breakpoint algorithm': 'Spurious radial velocity signals induced by stellar activity come from the line shape variations in stellar spectra. Those can be quantified through the full-width-at-halfmaximum (FWHM) and the asymmetry of the cross-correlation function (CCF) computed from the spectra [e.g., 126-128]. Following [129], we first fit Skew Normal (SN) functions to the CCFs available from HARPS-N and CARMENES. An SN function is not only characterized by a location and a scale parameter (which are the counterparts of the mean and standard deviation of a Gaussian), but it has a further free parameter that expresses its skewness (hereafter denoted with γ ). For each observation, through the SN-fit we were able to retrieve the stellar radial velocity ( RV , quantified through the SN median), the FWHMSN, the contrast A , and the asymmetry γ . The errors σ RV of the RV measurements were inferred using a bootstrap approach. Denoting with fCCF the flux of a CCF data point, each point was perturbed by sampling values from a Normal distribution whose standard deviation is equal to p f CCF, since the errors a ff ecting the CCF data points are expected to be Poissonian.  \nAfter that, we applied the breakpoint ( bp ) method [130] to both the HARPS-N and CARMENES RV time series. The algorithm has been designed to detect those locations along the RV time series where the correlation changes against the vector [FWHMSN, A , γ ] are statistically significant. The goal is to then detrend the RV time series by applying a piece-wise interpolation to each segment found by the bp algorithm rather than performing an overall correction to the whole time series. In this way we are able to better correct for the contamination of stellar variability as shown by ref.[130] and ref.[131]. Finally, we jointly analyzed the RV time series using the MCMCI code [132], where we switched o ff the interaction with stellar evolutionary models to speed up the computations. We set up the detrending function on each piece-wise stationary segment found by the bp algorithm as a polynomial of the following form  \n<!-- image -->  \nwhere ( kt , kF , kA , k γ , kR ) is the vector of the polynomial orders whose optimal value has been established by launching several MCMC preliminary runs and selecting that  \ncombination which produces the minimum Bayesian Information Criterion (BIC) [133].  \nAfter performing a longer MCMCI run made of 4 independent runs (300 000 steps each), which successfully converged as checked through the Gelman-Rubin test [134], we retrieved the posterior distributions of the system parameters. Their median values along with their error bars at the 1 σ level are reported in Tables S2 and S3 of the Supplementary. The full RV time series and the phase-folded RVs of those planets whose detection is above the 3 σ level (planets d and f) are shown in Fig. S9 of the Supplementary.  \nM ethod II: multi -dimensional GP On the other hand, we also perform a multidimensional Gaussian Process (GP) approach to characterize the stellar and planetary signals in our RV time series as in [135, 136]. This approach has been proven useful to disentangle stellar and planetary signals in multi-planet systems [e.g., 137, 138]. We create N -dimensional GP models, including N time-series A i , as  \nA 1 = A 1 G ( t ) + B 1 ˙ G ( t ) . . . A N = ANG ( t ) + BN ˙ G ( t ) , (7)  \nwhere the variables A 1, B 1, · · · , AN , BN , are free parameters which relate the individual time-series to G ( t ) and ˙ G ( t ). In this approach, G ( t ) is assumed to be a latent (unobserved) variable that represents the projected area of the visible stellar disk that is covered by active regions as a function of time.  \nWe model the stellar signal using a GP whose covariance between two times ti and t j is given by  \nγ QP , i , j = exp       -sin 2 [ π ( ti -t j ) / P GP] 2 λ 2 P -( ti -t j ) 2 2 λ 2 e       , (8)  \nwhere γ QP , i , j is the Quasi-Periodic (QP) kernel, whose hyperparameters are, P GP, the GP characteristic period, λ p, the inverse of the harmonic complexity, and λ e, the long term evolution timescale.  \nWe perform a two-dimensional GP model between the RVs and FWHM. We note that these quantities are equivalent in the HARPS-N and CARMENES data. The multidimensional covariance matrix was created using the kernel given in Eq. (8) and its derivatives [135, 136]. We assume that RVs can be described as A i = AiG ( t ) + Bi ˙ G ( t ), while the FWHM time series is described as A i = AiG ( t ). The planetary signals were included in the model as the mean function of the RV time series. We use N Keplerian signals (where N is the number of planetary signals); each one of them depends on the time of minimum conjunction t 0, orbital period P , and Doppler semi-amplitude, K . All orbits are fixed to be circular, so the eccentricity and angle of periastron are fixed. For the FWHM the mean function was treated as an o ff set,  \nnoting that we include a di ff erent o ff set per instrument. We also include a jitter term per time series and per instrument to account for unaccounted systematic errors.  \nWe perform MCMC samplings of the parameter space using the code pyaneti [136, 139]. We sample the parameter space with 250 walkers and create the posterior distributions with the last 5000 iterations of converged chains with a thin factor of 10. This leads to posterior distributions of 125 000 points for each sampled parameter. Figure 2 shows the spectroscopic time series resulting from this joint analysis. Median values along with their 1 σ uncertainties are reported in Table S4 of the Supplementary.  \nModeling techniques employing GP are particularly subject to overfitting giving their flexibility to reproduce the data [see e.g., 140, 141]. To test the robustness of our GP model, we carried out a cross-validation analysis. We repeat the twodimensional GP model described in this section but applied only to the HARPS-N data. Then, we create a predictive model with the inferred parameters and overlay the CARMENESdata with median o ff sets subtracted (similar to a training / evaluation set for machine learning algorithms). Figure S10 shows this analysis, zoomed in to the 2021 observing campaign. The plot shows that the RV data is in agreement with the predictive model, suggesting that our assumption that the stellar signal imprinted in the RVs can be described with a two-dimensional GP is valid for the time span of our observations. For the CCF FWHM CARMENES data, the correlation with the RV measurements is not as strong as for the HARPS-N one, thus the prediction is less accurate in this case.  \nBoth the two methods clearly detect planets HD 110067 d and f, M ethod II also detects planet b, however, the detection levels slightly di ff er. Recalling that the two di ff erent techniques are based on di ff erent RV extraction methodologies and on a different treatment of stellar activity, on the one hand, the slight output tension suggests that the RV data alone do not strongly constrain all the six Keplerian signals. We tested additional stellar mitigation approaches, such as sinusoid-fitting at the stellar rotation period and its harmonics or GP decorrelation as a function of time only, but they all turned unsuccessful at constraining the masses of any of the planets (only the time-dependent GP model could recover the signal of planet f, but with a much larger uncertainty, K f = 2 . 0 ± 1 . 0 m s -1 ). On the other hand, the RV semi-amplitudes inferred from the two methods shown in the manuscript are compatible within ∼ 1.5 σ , with the statistical tension ∆ I -II below ∼ 1 σ for planets b, c, d, and g. Extended Data Fig. 5 displays the pairs of posterior density functions for the RV semi-amplitudes of each planet for comparison.  \nFor HD 110067 f, after imposing a Gaussian prior centered around 41.05 days in both M ethod I and M ethod II RV models, we indeed recover a significant RV signal with a detection level of ∼ 3 σ . The planet transits only twice in the TESS data and the ground-based photometric campaign hints at moderate evidence for a planetary transit compatible with this value. Therefore, to secure an independent detection of planet f from spectroscopy, we performed an RV-only analysis imposing a uniform unbounded prior (between 30 to 54 days) to the orbital period. The MCMC converged and detected a clear Keplerian signal with a period Pf , uni = 40 . 2 ± 0 . 2 d. Thus, the RV  \ndata independently suggest the presence of a putative planet having a period close to 41.05 d. The tension at the ∼ 4 σ level with the predictions from the resonant chain model and the transit observations prove that the current RV data set cannot fully constrain the entire architecture of the planetary system.  \nWefinally checked whether there are further Keplerian signals within the RV time series. In particular, given that both planets e and g were not detected via our previous RV analyses where model-dependent values of the orbital periods were imposed as priors, we investigated the presence of potential planetary signals at di ff erent periods that could be attributed to planet e or g. To this end, we performed an MCMC run, where we modeled planets b and c (which are clearly confirmed by the transit events) along with planets d and f (which are clearly detected also in the RV time series). As a result, we produced the Generalized Lomb Scargle periodograms [121] of the residuals, obtained after subtracting the Keplerian signals of all four planets from the activity-cleaned time series (Fig. S11). The high false alarm probability level of the highest peak in both the HARPS-N (18%) and CARMENES (10%) residuals suggests that there are no signals left in the RV data that could be associated with other planets or a misidentification of the orbital periods of planets e and g.', '0.3.5 Final model': 'We computed a final model of the photometric and spectroscopic data sets of the HD 110067 system. Based on the analyses above, neither the light curves nor the radial velocities are precise enough to constrain the eccentricity of the planets. Assuming circular orbits, the photometry and radial velocity thus only constrain jointly the period and phase of a given planet in the system. However, the transit data dominate the precision of these two quantities (by several orders of magnitude). Therefore, for our final model, we opt to perform an independent analysis of the photometry and radial velocity datasets, where priors inform the planet periods and phases in the radial velocity model based on the posterior distributions of the photometry-only fit. Besides, the large number of free parameters in each of the models makes it computationally expensive to run a joint fit, not to mention the complications for numerical samplers to explore the vast multi-dimensional parameter space. Table 1 shows the most relevant planetary parameters of the system based on the photometric fit from Table S1, the radial velocity fit using M ethod II from Table S4, and the stellar parameters from Extended Data Table 3. A corner plot with the posterior distribution of the fitted transit parameters is shown in Fig. S12. The resulting best-fit models and corresponding credibility bands are presented in Fig. 1 for the TESS and CHEOPS photometry and Fig. 2 for the radial velocities.', '0.3.6 Planetary internal structures': 'Using a Bayesian analysis [142, 143], we computed the possible internal structures of the six planets of the system, using the results provided in Tables 1, Extended Data Table 3 and the planetary masses from M ethod I and M ethod II. The forward model used to compute the likelihood is based on a four-layer structure: a central core (iron and sulfur), a silicate mantle (containing Si, Mg, and Fe), a water layer, and a gas', '34 The HD 110067 planetary system': "- [40] Vanderburg, A. et al. TESS Spots a Compact System of Super-Earths around the Naked-eye Star HR 858. Astrophys. J. 881 (1), L19 (2019).\n- [41] Deming, D. et al. Spitzer Secondary Eclipses of the Dense, Modestlyirradiated, Giant Exoplanet HAT-P-20b Using Pixel-level Decorrelation. Astrophys. J. 805 (2), 132 (2015).\n- [42] Luger, R. et al. EVEREST: Pixel Level Decorrelation of K2 Light Curves. Astron. J. 152 (4), 100 (2016).\n- [43] Luger, R. et al. starry: Analytic Occultation Light Curves. Astron. J. 157 , 64 (2019).\n- [44] Lightkurve Collaboration et al. Lightkurve: Kepler and TESS time series analysis in Python. Astrophysics Source Code Library (2018).\n- [45] Gilliland, R. L. et al. Kepler Mission Stellar and Instrument Noise Properties. Astrophys. J. Suppl. 197 (1), 6 (2011).\n- [46] Van Cleve, J. E. et al. That's How We Roll: The NASA K2 Mission Science Products and Their Performance Metrics. Proc. Acad. Sci. Pacific 128 (965), 075002 (2016).\n- [47] Schanche, N. et al. TOI-2257 b: A highly eccentric long-period sub-Neptune transiting a nearby M dwarf. Astron. Astrophys. 657 , A45 (2022).\n- [48] Ulmer-Moll, S. et al. Two long-period transiting exoplanets on eccentric orbits: NGTS-20 b (TOI-5152 b) and TOI-5153 b. Astron. Astrophys. 666 , A46 (2022).\n- [49] Osborn, A. et al. TOI-431 / HIP 26013: a super-Earth and a sub-Neptune transiting a bright, early K dwarf, with a third RV planet. Mon. Not. R. Astron. Soc. 507 (2), 2782-2803 (2021).\n- [50] Tuson, A. et al. TESS and CHEOPS discover two warm sub-Neptunes transiting the bright K-dwarf HD 15906. Mon. Not. R. Astron. Soc. 523 (2), 3090-3118 (2023).\n- [51] Szab'o, G. M. et al. The changing face of AU Mic b: stellar spots, spinorbit commensurability, and transit timing variations as seen by CHEOPS and TESS. Astron. Astrophys. 654 , A159 (2021).\n- [52] Morris, B. M. et al. CHEOPS precision phase curve of the Super-Earth 55 Cancri e. Astron. Astrophys. 653 , A173 (2021).\n- [53] Hoyer, S. et al. Expected performances of the Characterising Exoplanet Satellite (CHEOPS). III. Data reduction pipeline: architecture and simulated performances. Astron. Astrophys. 635 , A24 (2020).  \n- [54] Narita, N. et al. MuSCAT2: four-color simultaneous camera for the 1.52-m Telescopio Carlos S'anchez. J. Astron. Telesc. Instrum. Syst. 5 (1), 015001 (2019).\n- [55] Parviainen, H. et al. MuSCAT2 multicolour validation of TESS candidates: an ultra-short-period substellar object around an M dwarf. Astron. Astrophys. 633 , A28 (2020).\n- [56] Brown, T. M. et al. Las Cumbres Observatory Global Telescope Network. Proc. Acad. Sci. Pacific 125 (931), 1031 (2013).\n- [57] McCully, C. et al. Real-time processing of the imaging data from the network of Las Cumbres Observatory Telescopes using BANZAI. Proc. SPIE 10707 , 107070K (2018).\n- [58] Collins, K. A., Kielkopf, J. F., Stassun, K. G. & Hessman, F. V. AstroImageJ: Image Processing and Photometric Extraction for Ultra-precise Astronomical Light Curves. Astron. J. 153 , 77 (2017).\n- [59] Wheatley, P. J. et al. The Next Generation Transit Survey (NGTS). Mon. Not. R. Astron. Soc. 475 (4), 4476-4493 (2018).\n- [60] Garcia-Mejia, J. et al. Marshall, H. K., Spyromilio, J. & Usuda, T. (eds) The Tierras Observatory: An ultra-precise photometer to characterize nearby terrestrial exoplanets . (eds Marshall, H. K., Spyromilio, J. & Usuda, T.) Groundbased and Airborne Telescopes VIII , Vol. 11445 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series , 114457R (2020).\n- [61] Demory, B. O. et al. A super-Earth and a sub-Neptune orbiting the bright, quiet M3 dwarf TOI-1266. Astron. Astrophys. 642 , A49 (2020).\n- [62] Narita, N. et al. Evans, C. J., Bryant, J. J. & Motohara, K. (eds) MuSCAT3: a 4-color simultaneous camera for the 2m Faulkes Telescope North . (eds Evans, C. J., Bryant, J. J. & Motohara, K.) Ground-based and Airborne Instrumentation for Astronomy VIII , Vol. 11447, 114475K. International Society for Optics and Photonics (SPIE, 2020).\n- [63] Fukui, A. et al. Measurements of Transit Timing Variations for WASP-5b. Pub. Astron. Soc. Japan 63 , 287-300 (2011).\n- [64] Ciardi, D. R., Beichman, C. A., Horch, E. P. & Howell, S. B. Understanding the E ff ects of Stellar Multiplicity on the Derived Planet Radii from Transit Surveys: Implications for Kepler, K2, and TESS. Astrophys. J. 805 (1), 16 (2015).\n- [65] Hayward, T. L. et al. PHARO: A Near-Infrared Camera for the Palomar Adaptive Optics System. Proc. Acad. Sci. Pacific 113 (779), 105-118 (2001).", '36 The HD 110067 planetary system': "- [66] Dekany, R. et al. PALM-3000: Exoplanet Adaptive Optics for the 5 m Hale Telescope. Astrophys. J. 776 (2), 130 (2013).\n- [67] Furlan, E. et al. The Kepler Follow-up Observation Program. I. A Catalog of Companions to Kepler Stars from High-Resolution Imaging. Astron. J. 153 (2), 71 (2017).\n- [68] Scott, N. J. et al. Twin High-resolution, High-speed Imagers for the Gemini Telescopes: Instrument description and science verification results. Frontiers in Astronomy and Space Sciences 8 , 138 (2021).\n- [69] Howell, S. B., Everett, M. E., Sherry, W., Horch, E. & Ciardi, D. R. Speckle Camera Observations for the NASA Kepler Mission Follow-up Program. Astron. J. 142 (1), 19 (2011).\n- [70] Mugrauer, M. & Michel, K.-U. Gaia search for stellar companions of TESS Objects of Interest. Astronomische Nachrichten 341 (10), 996-1030 (2020).\n- [71] Mugrauer, M. & Michel, K.-U. Gaia search for stellar companions of TESS Objects of Interest II. Astronomische Nachrichten 342 (6), 840-864 (2021).\n- [72] Ziegler, C. et al. SOARTESS Survey. I. Sculpting of TESS Planetary Systems by Stellar Companions. Astron. J. 159 (1), 19 (2020).\n- [73] Lafarga, M. et al. The CARMENES search for exoplanets around M dwarfs. Radial velocities and activity indicators from cross-correlation functions with weighted binary masks. Astron. Astrophys. 636 , A36 (2020).\n- [74] Zechmeister, M. et al. Spectrum radial velocity analyser (SERVAL). Highprecision radial velocities and two alternative spectral indicators. Astron. Astrophys. 609 , A12 (2018).\n- [75] Cosentino, R. et al. Ramsay, S. K., McLean, I. S. & Takami, H. (eds) HARPSN @ TNG, two year harvesting data: performances and results . (eds Ramsay, S. K., McLean, I. S. & Takami, H.) Ground-based and Airborne Instrumentation for Astronomy V , Vol. 9147, 91478C. International Society for Optics and Photonics (SPIE, 2014).\n- [76] Santos, N. C. et al. SWEET-Cat: A catalogue of parameters for Stars With ExoplanETs. I. New atmospheric parameters and masses for 48 stars with planets. Astron. Astrophys. 556 , A150 (2013).\n- [77] Sousa, S. G. ARES + MOOG: A Practical Overview of an Equivalent Width (EW) Method to Derive Stellar Parameters , 297-310 (Springer, 2014).\n- [78] Sousa, S. G. et al. SWEET-Cat 2.0: The Cat just got SWEETer. Higher quality spectra and precise parallaxes from Gaia eDR3. Astron. Astrophys. 656 , A53  \n(2021).  \n- [79] Sousa, S. G., Santos, N. C., Israelian, G., Mayor, M. & Monteiro, M. J. P. F. G. Anewcode for automatic determination of equivalent widths: Automatic Routine for line Equivalent widths in stellar Spectra (ARES). Astron. Astrophys. 469 (2), 783-791 (2007).\n- [80] Sousa, S. G., Santos, N. C., Adibekyan, V., Delgado-Mena, E. & Israelian, G. ARES v2: new features and improved performance. Astron. Astrophys. 577 , A67 (2015).\n- [81] Sneden, C. A. Carbon and Nitrogen Abundances in Metal-Poor Stars. Ph.D. thesis, THE UNIVERSITY OF TEXAS AT AUSTIN. (1973).\n- [82] Kurucz, R. L. SYNTHE spectrum synthesis programs and line data (Astrophysics Source Code Library, 1993).\n- [83] Brahm, R., Jord'an, A., Hartman, J. & Bakos, G. ZASPE: A Code to Measure Stellar Atmospheric Parameters and their Covariance from Spectra. Mon. Not. R. Astron. Soc. 467 (1), 971-984 (2017).\n- [84] Adibekyan, V. Z. et al. Chemical abundances of 1111 FGK stars from the HARPS GTO planet search program. Galactic stellar populations and planets. Astron. Astrophys. 545 , A32 (2012).\n- [85] Adibekyan, V. et al. Identifying the best iron-peak and α -capture elements for chemical tagging: The impact of the number of lines on measured scatter. Astron. Astrophys. 583 , A94 (2015).\n- [86] Castelli, F. & Kurucz, R. L. Piskunov, N., Weiss, W. W. & Gray, D. F. (eds) New Grids of ATLAS9 Model Atmospheres . (eds Piskunov, N., Weiss, W. W. &Gray, D. F.) Modelling of Stellar Atmospheres , Vol. 210 of IAU Symposium , A20 (2003).\n- [87] Allard, F. Booth, M., Matthews, B. C. & Graham, J. R. (eds) The BT-Settl Model Atmospheres for Stars, Brown Dwarfs and Planets . (eds Booth, M., Matthews, B. C. & Graham, J. R.) Exploring the Formation and Evolution of Planetary Systems , Vol. 299, 271-272 (2014).\n- [88] Blackwell, D. E. & Shallis, M. J. Stellar angular diameters from infrared photometry. Application to Arcturus and other stars; with e ff ective temperatures. Mon. Not. R. Astron. Soc. 180 , 177-191 (1977).\n- [89] Schanche, N. et al. WASP-186 and WASP-187: two hot Jupiters discovered by SuperWASP and SOPHIE with additional observations by TESS. Mon. Not. R. Astron. Soc. 499 (1), 428-440 (2020).  \n- [133] Schwarz, G. Estimating the Dimension of a Model. Ann. Stat. 6 (2), 461-464 (1978) .\n- [134] Gelman, A. & Rubin, D. B. Inference from Iterative Simulation Using Multiple Sequences. Stat. Sci. 7 , 457-472 (1992).\n- [135] Rajpaul, V., Aigrain, S., Osborne, M. A., Reece, S. & Roberts, S. A Gaussian process framework for modelling stellar activity signals in radial velocity data. Mon. Not. R. Astron. Soc. 452 (3), 2269-2291 (2015).\n- [136] Barrag'an, O., Aigrain, S., Rajpaul, V. M. & Zicher, N. PYANETI - II. A multidimensional Gaussian process approach to analysing spectroscopic timeseries. Mon. Not. R. Astron. Soc. 509 (1), 866-883 (2022). .\n- [137] Barrag'an, O. et al. The young HD 73583 (TOI-560) planetary system: two 10M ⊕ mini-Neptunes transiting a 500-Myr-old, bright, and active K dwarf. Mon. Not. R. Astron. Soc. 514 (2), 1606-1627 (2022).\n- [138] Zicher, N. et al. One year of AU Mic with HARPS - I. Measuring the masses of the two transiting planets. Mon. Not. R. Astron. Soc. 512 (2), 3060-3078 (2022).\n- [139] Barrag'an, O., Gandolfi, D. & Antoniciello, G. PYANETI: a fast and powerful software suite for multiplanet radial velocity and transit fitting. Mon. Not. R. Astron. Soc. 482 , 1017-1030 (2019).\n- [140] Cale, B. L. et al. Diving Beneath the Sea of Stellar Activity: Chromatic Radial Velocities of the Young AU Mic Planetary System. Astron. J. 162 (6), 295 (2021).\n- [141] Blunt, S. et al. Overfitting A ff ects the Reliability of Radial Velocity Mass Estimates of the V1298 Tau Planets. arXiv e-prints arXiv:2306.08145 (2023).\n- [142] Dorn, C. et al. Can we constrain the interior structure of rocky exoplanets from mass and radius measurements? Astron. Astrophys. 577 , A83 (2015).\n- [143] Dorn, C. et al. A generalized Bayesian inference method for constraining the interiors of super Earths and sub-Neptunes. Astron. Astrophys. 597 , A37 (2017).\n- [144] Haldemann, J., Alibert, Y., Mordasini, C. & Benz, W. AQUA: a collection of H2O equations of state for planetary models. Astron. Astrophys. 643 , A105 (2020).\n- [145] Hakim, K. et al. A new ab initio equation of state of hcp-Fe and its implication on the interior structure and mass-radius relations of rocky super-Earths. Icarus 313 , 61-78 (2018).  \n(2021).  \n- [160] Cannon, A. J. & Pickering, E. C. The Henry Draper catalogue 0h, 1h, 2h, and 3h. Annals of Harvard College Observatory 91 , 1-290 (1918).\n- [161] Gaia Collaboration et al. Gaia Early Data Release 3. Summary of the contents and survey properties. Astron. Astrophys. 649 , A1 (2021).\n- [162] Yoss, K. M. & Gri ffi n, R. F. Radial Velocities and DDO, BV Photometry of Henry Draper G5-M Stars Near the North Galactic Pole. Journal of Astrophysics and Astronomy 18 , 161 (1997).\n- [163] Skrutskie, M. F. et al. The Two Micron All Sky Survey (2MASS). Astron. J. 131 , 1163-1183 (2006).\n- [164] Delisle, J. B. Analytical model of multi-planetary resonant chains and constraints on migration scenarios. Astron. Astrophys. 605 , A96 (2017).", 'Data Availability': 'The TESS observations used in this study are publicly available at the Mikulski Archive for Space Telescopes (https: // archive.stsci.edu / missions-and-data / tess). The CHEOPS observations used in this study are available at the CHEOPS mission archive (https: // cheops-archive.astro.unige.ch / archive browser / ). The groundbased photometry and high-resolution imaging observations are uploaded to ExoFOP (https: // exofop.ipac.caltech.edu / tess / target.php?id = 347332255) and are publicly available. CARMENES and HARPS-N reduced spectra, together with the derived CCF-based radial velocities and spectral indicators are available in Zenodo (https: // doi.org / 10.5281 / zenodo.8211589). All reduced transit photometry and radial velocity measurements used in this work are also provided in Zenodo (https: // doi.org / 10.5281 / zenodo.8211589).', 'Code Availability': 'We used the following publicly available codes, resources and Python packages to reduce, analyze and interpret our observations of HD 110067: numpy [155], matplotlib [156], astropy [157], lightkurve [44], PIPE [51, 52], AstroImageJ [58], raccoon [73], serval [74], ARES [79, 80], MOOG [81], ZASPE [83], emcee [158], CLES [96], exoplanet [99], MonoTools [106], pymc3 [117], ArviZ [120], GLS [121], MCMCI [132], and pyaneti [136, 139]. We can share the code used in the data reduction or data analysis on request.', 'Acknowledgements': "We acknowledge the use of public TESS data from pipelines at the TESS Science O ffi ce and at the TESS Science Processing Operations Center. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center for the production of the SPOC data products. CHEOPS is an ESA mission in partnership with Switzerland with important contributions to the payload and the ground segment from Austria, Belgium, France, Germany, Hungary, Italy, Portugal, Spain, Sweden, and the United Kingdom. The CHEOPS Consortium would like to gratefully acknowledge the support received by all the agencies, o ffi ces, universities, and industries involved. Their flexibility and willingness to explore new approaches were essential to the success of this mission. CARMENES acknowledges financial support from the Agencia Estatal de Investigaci'on of the Ministerio de Ciencia e Innovaci'on MCIN / AEI / 10.13039 / 501100011033 and the ERDF 'A way of making Europe' through projects PID2019-107061GB-C61, PID2019107061GB-C66, PID2021-125627OB-C31, and PID2021-125627OB-C32, from the Centre of Excellence 'Severo Ochoa' award to the Instituto de Astrof'ısica de Canarias (CEX2019-000920-S), from the Centre of Excellence 'Mar'ıa de Maeztu' award to the Institut de Ci'encies de l'Espai (CEX2020-001058-M), and from the Generalitat de Catalunya / CERCA programme. Based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated on the island of La Palma by the Fundaci'on Galileo Galilei of the INAF (Istituto Nazionale di Astrofisica) at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias. This article is based on observations made with the MuSCAT2 instrument, developed by ABC, at Telescopio Carlos S'anchez operated on the island of Tenerife by the IAC in the Spanish Observatorio del Teide. This paper is based on observations made with the MuSCAT3 instrument, developed by the Astrobiology Center and under financial supports by JSPS KAKENHI (JP18H05439) and JST PRESTO (JPMJPR1775), at Faulkes Telescope North on Maui, HI, operated by the Las Cumbres Observatory. Tierras is supported by grants from the John Templeton Foundation and the Harvard Origins of Life Initiative. The opinions expressed in this publication are those of the authors and do not necessarily reflect the views of the John Templeton Foundation. The NGTS facility is operated by the consortium institutes with support from the UK Science and Technology Facilities Council (STFC) under projects ST / M001962 / 1 and ST / S002642 / 1. Some of the observations presented in this paper were carried out at the Observatorio Astron'omico Nacional on the Sierra de San Pedro M'artir (OAN-SPM), Baja California, M'exico. This work makes use of observations from the Las Cumbres Observatory global telescope network. Some of the observations in this paper made use of the High-Resolution Imaging instrument Alopeke and were obtained under Gemini LLP Proposal Number GN-S-2021A-LP105. Alopeke was funded by the NASA Exoplanet Exploration Program and built at the NASA Ames Research Center by Steve B. Howell, Nic Scott, Elliott P. Horch, and Emmett Quigley. Alopeke was mounted on the Gemini North telescope of the international Gemini Observatory, a program of NSF OIR Lab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative  \nagreement with the National Science Foundation. on behalf of the Gemini partnership: the National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigaci'on y Desarrollo (Chile), Ministerio de Ciencia, Tecnolog'ıa e Innovaci'on (Argentina), Minist'erio da Ciˆencia, Tecnologia, Inovac¸ ˜oes e Comunicac¸ ˜oes (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea). This work was supported by the KESPRINT collaboration, an international consortium devoted to the characterization and research of exoplanets discovered with space-based missions. R.Lu. thanks Prof. Daniel Fabrycky for helpful discussions regarding the orbital dynamics of the HD 110067 system. R.Lu. acknowledges funding from University of La Laguna through the Margarita Salas Fellowship from the Spanish Ministry of Universities ref. UNI / 551 / 2021-May 26, and under the EU Next Generation funds. This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40 182901 and 51NF40 205606. A.C.Ca. and T.G.Wi. acknowledge support from STFC consolidated grant numbers ST / R000824 / 1 and ST / V000861 / 1, and UKSA grant number ST / R003203 / 1. O.Ba. acknowledges that has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 865624) M.Le. acknowledges support of the Swiss National Science Foundation under grant number PCEFP2 194576. P.F.L.Ma. acknowledges support from STFC research grant number ST / M001040 / 1. Y.Al. acknowledges support from the Swiss National Science Foundation (SNSF) under grant 200020 192038. D.Ga. gratefully acknowledges financial support from the CRT foundation under Grant No. 2018.2323 'Gaseous or rocky? Unveiling the nature of small worlds'. J.A.Eg. acknowledges support from the Swiss National Science Foundation (SNSF) under grant 200020 192038. G.No. thanks for the research funding from the Ministry of Education and Science programme the Excellence Initiative - Research University conducted at the Centre of Excellence in Astrophysics and Astrochemistry of the Nicolaus Copernicus University in Torun, Poland. D.Ra. was supported by NASA under award number NNA16BD14C for NASA Academic Mission Services. M.La. acknowledges funding from a UKRI Future Leader Fellowship, grant number MR / S035214 / 1. V.Ad. is supported by FCT through national funds by the following grants UIDB / 04434 / 2020, UIDP / 04434 / 2020, and 2022.06962.PTDC. P.J.Am. acknowledges financial support from the grants CEX2021-001131-S and PID2019-109522GB-C52, both funded by MCIN / AEI / 10.13039 / 501100011033 and by 'ERDF: A way of making Europe'. S.C.C.Ba. acknowledges support from FCT through FCT contracts nr. IF / 01312 / 2014 / CP1215 / CT0004. X.Bo., S.Ch., D.Ga., M.Fr. and J.La. acknowledge their role as ESA-appointed CHEOPS science team members. L.Bo., G.Br., V.Na., I.Pa., G.Pi., R.Ra., G.Sc., V.Si., and T.Zi. acknowledge support from CHEOPS ASI-INAF agreement n. 2019-29-HH.0. A.Br. was supported by the SNSA. Contributions at the Mullard Space Science Laboratory by E.M.Br. were supported by STFC through the consolidated grant ST / W001136 / 1. S.C.Ga. acknowledges support from UNAM PAPIIT-IG101321. D.Ch. and J.G-M. thank the sta ff at the F. L. Whipple Observatory for their assistance  \nin the refurbishment and maintenance of the 1.3-m telescope. W.D.Co. acknowledges support from NASA grant 80NSSC23K0429. This is University of Texas Center for Planetary Systems Habitability Contribution 0063. K.A.Co. acknowledges support from the TESS mission via subaward s3449 from MIT. H.J.De. acknowledges support from the Spanish Research Agency of the Ministry of Science and Innovation (AEI-MICINN) under grant PID2019-107061GB-C66, DOI: 10.13039 / 501100011033. This project was supported by the CNES. The Belgian participation to CHEOPS has been supported by the Belgian Federal Science Policy O ffi ce (BELSPO) in the framework of the PRODEX Program, and by the University of Li'ege through an ARC grant for Concerted Research Actions financed by the Wallonia-Brussels Federation. L.De. is an F.R.S.-FNRS Postdoctoral Researcher. This work was supported by FCT - Fundac¸˜ao para a Ciˆencia e a Tecnologia through national funds and by FEDER through COMPETE2020 - Programa Operacional Competitividade e Internacionalizac˜ao by these grants: UID / FIS / 04434 / 2019, UIDB / 04434 / 2020, UIDP / 04434 / 2020, PTDC / FIS-AST / 32113 / 2017 & POCI-010145-FEDER- 032113, PTDC / FIS-AST / 28953 / 2017 & POCI-01-0145-FEDER028953, PTDC / FIS-AST / 28987 / 2017 & POCI-01-0145-FEDER-028987, O.D.S.De. is supported in the form of work contract (DL 57 / 2016 / CP1364 / CT0004) funded by national funds through FCT. B.-O.De. acknowledges support from the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract number MB22.00046. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (project F our A ces grant agreement No 724427). It has also been carried out in the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). D.Eh. acknowledges financial support from the Swiss National Science Foundation for project 200021 200726. E.EBo. acknowledges financial support from the European Union and the State Agency of Investigation of the Spanish Ministry of Science and Innovation (MICINN) under the grant PRE2020-093107 of the Pre-Doc Program for the Training of Doctors (FPISO) through FSE funds. M.Fr. gratefully acknowledges the support of the Swedish National Space Agency (DNR 65 / 19, 174 / 18). J.G-M. acknowledges support by the National Science Foundation through a Graduate Research Fellowship under grant No. DGE1745303 and by the Ford Foundation through a Ford Foundation Predoctoral Fellowship, administered by the National Academies of Sciences, Engineering, and Medicine. The contributions at the University of Warwick by S.Gi. have been supported by STFC through consolidated grants ST / L000733 / 1 and ST / P000495 / 1. M.Gi. is F.R.S.-FNRS Research Director. Y.G.M.Ch. acknowledges support from UNAMPAPIIT-IG101321. E.Go. acknowledge the support by the Thueringer Ministerium fuer Wirtschaft, Wissenschaft und Digitale Gesellschaft. M.N.Gu. is the ESA CHEOPS Project Scientist and Mission Representative, and as such also responsible for the Guest Observers (GO) Programme. M.N.Gu. does not relay proprietary information between the GO and Guaranteed Time Observation (GTO) Programmes, and does not decide on the definition and target selection of the GTO Programme. A.P.Ha. acknowledges support by DFG grant HA 3279 / 12-1 within the DFG Schwerpunkt  \nSPP 1992. Ch.He. acknowledges support from the European Union H2020-MSCAITN-2019 under Grant Agreement no. 860470 (CHAMELEON). S.Ho. gratefully acknowledges CNES funding through the grant 837319. This work is partly supported by JST CREST Grant Number JPMJCR1761. K.G.Is. 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 for the programme. J.Ko. gratefully acknowledges the support of the Swedish National Space Agency (SNSA, DNR 2020-00104) and of the Swedish Research Council (VR: Etableringsbidrag 2017-04945). K.W.F.La. was supported by Deutsche Forschungsgemeinschaft grants RA714 / 14-1 within the DFG Schwerpunkt SPP 1992, Exploring the Diversity of Extrasolar Planets. 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. A.L-E. acknowledges support from the CNES (Centre national d''etudes spatiales, France). This work is partly supported by Astrobiology Center SATELLITE Research project AB022006. This work is partly supported by JSPS KAKENHI Grant Number JP18H05439 and JST CREST Grant Number JPMJCR1761. H.L.M.Os. acknowledges funding support by STFC through a PhD studentship. H.Pa. acknowledges the support by the Spanish Ministry of Science and Innovation with the Ramon y Cajal fellowship number RYC2021-031798-I. This work was also partially supported by a grant from the Simons Foundation (PI Queloz, grant number 327127). S.N.Qu. acknowledges support from the TESS mission via subaward s3449 from MIT. S.N.Qu. acknowledges support from the TESS GI Program under award 80NSSC21K1056 (G03268). L.Sa. acknowledges support from UNAM PAPIIT project IN110122. N.C.Sa. acknowledges funding by the European Union (ERC, FIERCE, 101052347). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. N.Sc. acknowledges support from the Swiss National Science Foundation (PP00P2-163967 and PP00P2-190080) and NASA under award number 80GSFC21M0002. S.G.So. acknowledge support from FCT through FCT contract nr. CEECIND / 00826 / 2018 and POPH / FSE (EC). Gy.M.Sz. acknowledges the support of the Hungarian National Research, Development and Innovation O ffi ce (NKFIH) grant K-125015, a PRODEX Experiment Agreement No. 4000137122, the Lendulet LP2018-7 / 2021 grant of the Hungarian Academy of Science and the support of the city of Szombathely. A.Tu. acknowledges funding support from the STFC via a PhD studentship. V.V.Ey. acknowledges support by STFC through the consolidated grant ST / W001136 / 1. V.V.Gr. is an F.R.SFNRS Research Associate. J.Ve. acknowledges support from the Swiss National Science Foundation (SNSF) under grant PZ00P2 208945. N.A.Wa. acknowledges UKSA grant ST / R004838 / 1. N.Wa. is partly supported by JSPS KAKENHI Grant Number JP21K20376.  \nThe HD 110067 planetary system", 'Authors and A ffi liations': 'Author Contributions R.Lu, H.P.Os, A.Le., E.Pa., A.Bo., O.Ba., and T.G.Wi. conceived the project and contributed notably to the writing of this manuscript. R.Lu. and H.P.Os. led the analysis of the photometric data. A.Le. led the dynamical analysis of the system and developed the method with J.-B.De. to predict the orbits of the planets based on their resonant state within the chain. R.Lu., A.Bo. and O.Ba. led the analysis of the radial velocity data and the stellar activity mitigation. T.G.Wi. led the stellar characterization with the help of V.Ad., S.G.So., A.Bo., V.V.Gr., S.Sa., and W.D.Co. Y.Al. and J.A.Eg. led the analysis of the internal structures, while L.Fo. and A.Bo. performed the atmospheric evolution simulations. D.Ra., J.D.Tw., and J.M.Je. improved the TESS data reduction to recover the missing cadences a ff ected by reflected light and high background. R.Lu., E.Pa., and G.No. planned and obtained the time for the observations with CARMENES and HARPS-N. CARMENES observations were made possible by M.La., J.C.Mo., P.J.Am., A.Qu., and I.Ri. HARPS-N observations were made possible by I.Ca., J.O-M., F.Mu., H.J.De., J.Ko., D.Ga., J.H.Li., W.D.Co., E.W.Gu., V.V.Ey., H.L.M.Os., S.Re., E.Go., F.Da., and K.W.F.La. High-resolution imaging observations from Palomar and Gemini North were made possible by A.W.Bo., D.R.Ci, I.J.M.Cr., S.B.Ho., E.Ma., and J.E.Sc. Ground-based photometric observations to catch the transit of planet f were made possible by the MuSCAT2 (R.Lu., E.Pa., N.Na., J.H.Li., K.Ik., E.E-B., J.O-M., N.Wa., F.Mu., G.No., A.Fu., H.Pa., M.Mo., T.Ka., J.P.Le., and T.Ko.), LCO (T.G.Wi., R.Lu., H.P.Os., E.Pa., A.Le., A.Tu., M.J.Ho., Y.Al., and D.Ga.), NGTS (H.P.Os., S.Gi., D.Ba., D.R.An., M.Mo., A.M.S.Sm., E.M.Br., and S.Ud.), Tierras (J.G-M. and D.Ch.), SAINT-EX (N.Sc., Y.G-M-C., L.Sa., S.C-G., and B.-O.De.), and MuSCAT3 (N.Na., J.H.Li., K.Ik., N.Wa., A.Fu., M.Mo., T.Ka., J.P.Le., and T.Ko.) instruments. The remaining authors provided key contributions to the development of the TESS and CHEOPS mission. All authors read and commented on the manuscript, and helped with its revision.  \nCorresponding author Correspondence to Rafael Luque (rluque@uchicago.edu).  \nCompeting interests The authors declare no competing interests.  \nThe HD 110067 planetary system  \nExtended Data Figures and Tables  \nExtended Data Fig. 1 Transit duration versus transit depth for all unassigned transits in TESS data . TESS Sector 23 and Sector 49 are shown as di ff erent colors. The numbers above each transit denote the mid-transit time in TJD. Contours represent percentile levels, the innermost one corresponding to the 50th percentile and the outermost to the 99th percentile by increments of 10%. The transit of planet f in PLD photometry is marked with ∗ to indicate that its properties are heavily a ff ected by pre-transit systematic noise.  \n<!-- image -->  \nExtended Data Fig. 2 Generalized three-body Laplace angles for known systems in resonant chains. Included are the Galilean satellites, Kepler-60 [12, 115], Kepler-80 [159], K2-138 [112], Kepler-223 [110], TRAPPIST-1 [13], and TOI-178 [10]. Measurements belonging to the same system are marked with the same color. The line marks the observed distance to the theorized equilibrium (marked with a circle). The distances are estimated at 0th order in eccentricity [110, 111]. For most systems, a single estimation of the generalized Laplace angle is made, while [110] made an estimation for each Kepler quarter.  \n<!-- image -->  \nExtended Data Fig. 3 Observed distance from the equilibrium for all the simulated scenarios in which planets f and g continue the resonant chain . The y-axis is converted to the mean peak-to-peak amplitude from the generalized three-body Laplace angle using the following expression: mean( A ( Ψ i )) = C / 4. Case A2 remains the one that has the potential to be the closest to an equilibrium.  \n<!-- image -->', 'Extended Data Table 2 Ground-based photometric campaign observing log.': 'Extended Data Table 3 Stellar parameters of HD 110067. Error bars represent 1 σ uncertainties.  \nExtended Data Table 4 Distance of the estimated generalized three-body Laplace angle Ψ e = 0 to the closest equilibrium for all period ratios that are not excluded by available observations . Case A assumes that the mono-transit at 2641.5778 TJD belongs to the 5th planet, and 2656.0944 TJD belongs to the 6th planet. Case B assumes the opposite. The flag column ∗ = 1 indicates that the position of the equilibria varies with the masses of the planets, in which case the equilibrium is recomputed using [164], with masses computed using a mass-radius relation for sub-Neptunes [26]. For ∗ = 0, the equilibrium is 180 degrees.', '10 The HD 110067 planetary system': 'Supplementary Information Supplementary Information is available for this paper.  \nRights and permissions Reprints and permissions information is available at www.nature.com / reprints.'}