The HD 110067 planetary system]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é^5,6,†,
A. Bonfanti^7,
O. Barragán^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é^20,21,
T. Bandy^27,
T. Bárczy^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. Goffo^14,54,
Y. Gomez Maqueo Chew^37,
M. Güdel^55,
E. W. Guenther^54,
M. N. Günther^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. Parviainen^5,6,
G. Peter^69,
G. Piotto^33,70,
D. Pollacco^10,
D. Queloz^71,15,
S. N. Quinn^40,
A. Quirrenbach^72,
R. Ragazzoni^33,70,
N. Rando^27,
F. Ratti^27,
H. Rauer^36,73,74,
S. Redfield^75,
I. Ribas^20,21,
G. R. Ricker^3,
A. Rudat^3,
L. Sabin^76,
S. Salmon^4,
N. C. Santos^30,31,
G. Scandariato^68,
N. Schanche^12,77,
J. E. Schlieder^78,
S. Seager^3,79,80,
D. Ségransan^4,
A. Shporer^3,
A. E. Simon^2,
A. M. S. Smith^36,
S. G. Sousa^30,
M. Stalport^48,
Gy. M. Szabó^81,82,
N. Thomas^2,
A. Tuson^15,
S. Udry^4,
A. M. Vanderburg^3,
V. Van Eylen^25,
V. Van Grootel^48,
J. Venturini^4,
I. Walter^83,
N. A. Walton^84,
N. Watanabe^56,
J. N. Winn^85,
T. Zingales^70

^1 Department of Astronomy & Astrophysics, University of Chicago, Chicago, IL 60637, USA

^2 Space Research and Planetary Sciences, Physics Institute, University of Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland

^3 Department of Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

^4 Observatoire Astronomique de l'Université de Genève, Chemin Pegasi 51, 1290 Versoix, Switzerland

^5 Instituto de Astrofisica de Canarias, Via Lactea s/n, 38200 La Laguna, Tenerife, Spain

^6 Departamento de Astrofisica, Universidad de La Laguna, Astrofísico Francisco Sanchez s/n, 38206 La Laguna, Tenerife, Spain

^7 Space Research Institute, Austrian Academy of Sciences, Schmiedlstrasse 6, A-8042 Graz, Austria

^8 Sub-department of Astrophysics, Department of Physics, University
of Oxford, Oxford, OX1 3RH, UK

^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

^11 Centre for Exoplanets and Habitability, University of Warwick, Coventry, CV4 7AL, UK

^12 Center for Space and Habitability, University of Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland

^13 Astrophysics Group, Lennard Jones Building, Keele University, Staffordshire, ST5 5BG, United Kingdom

^14 Dipartimento di Fisica, Universita degli Studi di Torino, via Pietro Giuria 1, I-10125, Torino, Italy

^15 Cavendish Laboratory, JJ Thomson Avenue, Cambridge CB3 0HE, UK

^16 Institute of Astronomy, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University, Grudzicadzka 5, 87-100 Torun, Poland

^17 NASA Ames Research Center, Moffett Field, CA 94035, USA

^18 Research Institute for Advanced Computer Science, Universities Space Research Association, Washington, DC 20024, USA

^19 SETI Institute, Mountain View, CA 94043, USA

^20 Institut de Ciencies de l'Espai (ICE, CSIC), Campus UAB, Can Magrans s/n, 08193 Bellaterra, Spain

^21 Institut d’Estudis Espacials de Catalunya (IEEC), Gran Capità 2-4, 08034 Barcelona, Spain

^22 INAF, Osservatorio Astrofisico di Torino, Via Osservatorio, 20, I-10025 Pino Torinese To, Italy

^23 Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150-762 Porto, Portugal

^24 Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, 4169-007 Porto, Portugal

^25 Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey, RH5 6NT, UK

^26 Instituto de Astrofísica de Andalucía (IAA-CSIC), Glorieta de la Astronomía s/n, 18008 Granada, Spain

^27 European Space Agency (ESA), European Space Research and Technology Centre (ESTEC), Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands

^28 Admatis, 5. Kandó Kálmán Street, 3534 Miskolc, Hungary

^29 Depto. de Astrofisica, Centro de Astrobiologia (CSIC-INTA), ESAC campus, 28692 Villanueva de la Cañada (Madrid), Spain

^30 Instituto de Astrofisica e Ciencias do Espaco, Universidade do Porto, CAUP, Rua das Estrelas, 4150-762 Porto, Portugal

^31 Departamento de Fisica e Astronomia, Faculdade de Ciencias, Universidade do Porto, Rua do Campo Alegre, 4169-007 Porto, Portugal

^32 Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

^33 INAF, Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, 35122 Padova, Italy

^34 Department of Astronomy, California Institute of Technology, 1200 E. California Blvd, Pasadena, CA 91106

^35 Department of Astronomy, Stockholm University, AlbaNova University Center, 10691 Stockholm, Sweden

^36 Institute of Planetary Research, German Aerospace Center (DLR), Rutherfordstrasse 2, 12489 Berlin, Germany

^37 Universidad Nacional Autónoma de México, Instituto de Astronomía, AP 70-264, Ciudad de México,  04510, México

^38 Department of Physics and Astronomy, Georgia State University, Atlanta, GA 30302-4106, USA

^39 RECONS Institute, Chambersburg, PA 17201, USA

^40 Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

^41 Université de Paris Cité, Institut de physique du globe de Paris, CNRS, 1 Rue Jussieu, F-75005 Paris, France

^42 McDonald Observatory and Center for Planetary Systems Habitability, The University of Texas, Austin, Texas, USA

^43 Department of Physics and Astronomy, University of  Kansas, Lawrence, KS, USA

^44 Division of Geological and Planetary Sciences, 1200 E California Boulevard, Pasadena, CA 91125, USA

^45 Centre for Mathematical Sciences, Lund University, Box 118, 221 00 Lund, Sweden

^46 Aix Marseille Univ, CNRS, CNES, LAM, 38 rue Frédéric Joliot-Curie, 13388 Marseille, France

^47 Astrobiology Research Unit, Université de Liège, Allée du 6 Août 19C, B-4000 Liège, Belgium

^48 Space sciences, Technologies and Astrophysics Research (STAR) Institute, Université de Liège, Allée du 6 Août 19C, 4000 Liège, Belgium

^49 Centre Vie dans l’Univers, Faculté des sciences, Université de Genève, Quai Ernest-Ansermet 30, 1211 Genève 4, Switzerland

^50 Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD, 21218, USA

^51 Leiden Observatory, University of Leiden, PO Box 9513, 2300 RA Leiden, The Netherlands

^52 Department of Space, Earth and Environment, Chalmers University of Technology, Onsala Space Observatory, 439 92 Onsala, Sweden

^53 Komaba Institute for Science, The University of Tokyo, 3-8-1 Komaba, Meguro, Tokyo 153-8902, Japan

^54 Th uringer Landessternwarte Tautenburg, Sternwarte 5, D-07778 Tautenburg, Germany

^55 Department of Astrophysics, University of Vienna, Türkenschanzstrasse 17, 1180 Vienna, Austria

^56 Department of Multi-Disciplinary Sciences, Graduate School of Arts and Sciences, The University of Tokyo, 3-8-1 Komaba, Meguro, Tokyo 153-8902, Japan

^57 Konkoly Observatory, HUN-REN Research Centre for Astronomy and Earth Sciences, 1121 Budapest, Konkoly Thege Miklós út 15-17, Hungary

^58 ELTE Eötvös Loránd University, Institute of Physics, Pázmány Péter sétány 1/A, 1117 Budapest, Hungary

^59 Lund Observatory, Division of Astrophysics, Department of Physics, Lund University, Box 43, SE-221 00 Lund, Sweden

^60 IMCCE, UMR8028 CNRS, Observatoire de Paris, PSL Univ., Sorbonne Univ., 77 av. Denfert-Rochereau, 75014 Paris, France

^61 Institut d'astrophysique de Paris, UMR7095 CNRS, Université Pierre & Marie Curie, 98bis blvd. Arago, 75014 Paris, France

^62 Astrobiology Center, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

^63 National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan

^64 Department of Astronomical Science, The Graduated University for Advanced Studies, SOKENDAI, 2-21-1, Osawa, Mitaka, Tokyo, 181-8588, Japan

^65 U.S. Naval Observatory, Washington, D.C. 20392, USA

^66 Max Planck Institute for Astronomy, Heidelberg, Germany

^67 Instituto de Astronomía, Universidad Católica del Norte, Angamos 0610, 1270709, Antofagasta, Chile

^68 INAF, Osservatorio Astrofisico di Catania, Via S. Sofia 78, 95123 Catania, Italy

^69 Institute of Optical Sensor Systems, German Aerospace Center (DLR), Rutherfordstrasse 2, 12489 Berlin, Germany

^70 Dipartimento di Fisica e Astronomia "Galileo Galilei", Universita degli Studi di Padova, Vicolo dell'Osservatorio 3, 35122 Padova, Italy

^71 ETH Zurich, Department of Physics, Wolfgang-Pauli-Strasse 2, CH-8093 Zurich, Switzerland

^72 Landessternwarte, Zentrum für Astronomie der Universität Heidelberg, D-69117 Heidelberg, Germany

^73 Zentrum für Astronomie und Astrophysik, Technische Universität Berlin, Hardenbergstr. 36, D-10623 Berlin, Germany

^74 Institut fuer Geologische Wissenschaften, Freie Universitaet Berlin, Maltheserstrasse 74-100,12249 Berlin, Germany

^75 Astronomy Department and Van Vleck Observatory, Wesleyan University, Middletown, CT 06459, USA

^76 Instituto de Astronomía, Universidad Nacional Autónoma de México, Apdo. Postal 877, 22860 Ensenada, B.C., Mexico

^77 Department of Astronomy, University of Maryland, College Park, MD 20742, USA

^78 NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

^79 Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

^80 Department of Aeronautics and Astronautics, MIT, 77 Massachusetts Avenue, Cambridge, MA 02139, USA

^81 ELTE Eötvös Loránd University, Gothard Astrophysical Observatory, 9700 Szombathely, Szent Imre h. u. 112, Hungary

^82 HUN-REN-ELTE Exoplanet Research Group, 9700 Szombathely, Szent Imre h. u. 112, Hungary

^83 German Aerospace Center (DLR), Institute of Optical Sensor Systems, Rutherfordstraße 2, 12489 Berlin

^84 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK

^85 Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

^∗Corresponding author. E-mail<EMAIL_ADDRESS>
^†These authors contributed equally to this work.

Planets with radii between that of the Earth and Neptune (hereafter referred to as “sub-Neptunes”) are found in close-in orbits around more than half of all Sun-like stars [1, 2]. Yet, their composition, formation, and evolution remain poorly understood [3].
The study of multi-planetary 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 sub-Neptunes with radii ranging from 1.94 to 2.85 $R_{\oplus}$. Three of the planets have measured masses, yielding low bulk densities that suggest the presence of large hydrogen-dominated atmospheres.

HD 110067 (TIC 347332255) is a bright K0-type star in the constellation of Coma Berenices with mass and radius of approximately 80 percent of the Sun's. The Transiting Exoplanet Survey Satellite (TESS) monitored HD 110067 as part of its observations of Sector 23 [4]. The data, processed by the TESS Science Processing Operations Center (SPOC) [5], exhibited several dips that could be associated with transiting planets. SPOC reported two candidates: one with an orbital period of 5.642 days based on three dips with apparently similar depth and duration, and a second one with an unconstrained orbital period based on a single event. TESS re-observed the star two years later in Sector 49, revealing nine further transits incompatible with the previously announced candidates (Fig. <ref>).

Combining TESS Sectors 23 and 49 we could associate a fraction of the transits with two new planet candidates: HD 110067 b, a planet with an orbital period of 9.114 days, and HD 110067 c, a planet with an orbital period of 13.673 days. For the remaining unidentified transit events (two in Sector 23 and four in Sector 49), we attempted to match them by modeling them individually using a purely shape-based transit model agnostic to the orbital period (Methods and Extended Data Fig.<ref>). Then, we compared them in duration-depth space, allowing us to identify two “duo-transits” (planets seen to transit only once in each of the two widely-separated sectors) and two single-transits (solitary transit events seen only in Sector 49). Duo-transits have orbital periods limited to a finite number of harmonics or aliases which are constrained at the short-period end by the extent of continuous photometry observed before/after transit and at the long end by the distance between transits. Targeted observations with the CHaracterising ExOPlanets Satellite (CHEOPS) [6] allowed us to rule out many of these aliases and confirm a third planet in the system, HD 110067 d, with an orbital period of 20.519 days (Fig. <ref>).

The orbital periods of planets HD 110067 b, c, and d (9.114, 13.673, and 20.519 days, respectively) have ratios very close to 3/2 ($P_c/P_b=1.5003 $ and $P_d/P_c=1.5007$). Mean-motion resonances (MMRs) are orbital configurations where the period ratio of a pair of planets is oscillating near a rational number of the form $(k+q)/k$, where $k$ and $q$ are integers. First-order MMRs (where $q=1$) are the most common among planetary systems, as well as resonant chains. For the two innermost pairs of planets (bc and cd), we define two resonant angles $\phi_1= 2 \lambda_b -3\lambda_c + \varpi_c$ and  $\phi_2= 2 \lambda_c -3\lambda_d + \varpi_c$, where $\lambda$ is the mean longitude of the planets and $\varpi$ its longitude of periastron. Since $d\lambda_x/dt=2\pi/P_x$, and given the aforementioned period ratios between $b$, $c$, and $d$, a generalized Laplace relation links the triplet of planets via $2/P_b-5/P_c+3/P_d \approx 0$ [7, 8, 9].  This relation ensures that the associated Laplace angle $ \Psi_{bcd} = \phi_1-\phi_2 = 2\lambda_b-5\lambda_c+3\lambda_d $, evolves slowly ($d\Psi_{bcd}/dt \approx 0$), indicating that the three planets might indeed be trapped in a Laplace chain of 3/2 resonances.

We explored the possibility that the remaining “unmatched” dips correspond to planets that continue the first-order generalized Laplace resonant chain of HD 110067 b, c, and d. Using a Laplace relation similar to the one above, we can predict the possible orbital period of the planets, in a similar fashion as the discoveries of TOI-178 f and TRAPPIST-1 h [10, 11]. We assume that the next planet in the system must continue the generalized Laplace chain of first-order resonances. Considering the most common first-order MMRs (2/1, 3/2, 4/3, 5/4, 6/5), the only continuation of the chain that matches the spare duo-transit supports the presence of a fourth planet in the system, HD 110067 e, with an orbital period of 30.7931 days (in a 3/2 MMR with planet d).

We are then left with two mono-transits in TESS Sector 49, which we assume that they correspond to two separate planets f and g that continue the chain. With single transits, we cannot use the same method as above and must introduce a novel argument. In nature, all known three-body resonant chains are close to an equilibrium, i.e. an equilibrium value of their angles $\Psi$ [12, 10, 13, 14] (Extended Data Fig. <ref>). Assuming low eccentricities, only the period and a single transit of each planet are required to estimate the value of $\Psi$. We can therefore try out different periods for planets f and g, then use their single transits to estimate $\Psi_{def}$ and $\Psi_{efg}$, and see if it lands close to an equilibrium of the chain. Here, we try the same set of first-order MMRs (2/1, 3/2, 4/3, 5/4, 6/5) between e and f and between f and g (50 combinations in total, see Methods). Out of those, the only combination not excluded by existing data that can be close to an equilibrium of the chain is the one where $P_f/P_e=4/3$ and $P_g/P_f=4/3$, yielding the three outer generalized Laplace angles at less than 20 degrees from the closest equilibrium (Methods, Extended Data Table 4, and Extended Data Fig. <ref>).

According to this prediction, planets HD 110067 f and g would have orbital periods of 41.0575 and 54.7433 days, respectively. If so, both planets would have transited during TESS Sector 23 observations, but during a time when the effect of scattered light and the relative contributions of the Earth and Moon to the image backgrounds were highly significant. These frames are typically discarded, but a detailed reprocessing of the TESS Sector 23 observations triggered by this prediction showed two additional transit events at 1943.6 TJD and 1944.1 TJD matching exactly our model based on the system's resonant dynamics (Fig. <ref>). Additionally, a ground-based multi-instrument photometric campaign to catch a predicted transit of HD 110067 f (Methods and Extended Data Fig. <ref>) recovered a statistically significant detection ($\Delta$WAIC = 9.5) with a consistent depth and duration at the expected time, confirming its orbital period.

Also, we collected high-precision radial velocities of the star with the CARMENES [15] and HARPS-N [16] instruments. The dominant signals in the data correspond to the star's magnetic activity rather than the planetary companions but, by applying state-of-the-art analyses to model stellar activity (Methods), we could independently confirm the detection of HD 110067 f, measuring an orbital period and phase (using agnostic priors on both quantities) that matched the transits. We used the radial velocities of the system to measure precise masses for three of the planets (HD 110067 b, d, and f) and place upper limits on the remaining ones (Fig. <ref>). Although transit timing variations caused by the planets' mutual gravitational interactions are expected in resonant systems [17], our photometric analysis did not measure any significant deviation from the linear ephemeris (below 5 min for the inner triplet), likely due to the low number of individual transits for each planet. Further monitoring of the system will enable an independent measurement of the planetary masses using this technique.

The HD 110067 planetary system is thus comprised of at least six transiting planets orbiting in a chain of first-order MMRs (3/2—3/2—3/2—4/3—4/3). The planets have radii ranging between 1.94 and 2.85 times the radius of the Earth, orbital periods between 9 and 55 days, and equilibrium temperatures between 440 and 800 K (Fig. <ref>). Thanks to the constraints on the planetary masses and their location above the radius valley [18, 19], our internal composition modeling concludes that all the planets in the system (with the exception perhaps of planet e, which remains undetected in the radial velocity data) must possess large hydrogen-dominated atmospheres to explain their relatively low bulk densities (see Supplementary). A summary of the most relevant properties of the system is presented in Table 1. We note that for both planets e and g their orbital period measurement relies on a prediction based on the dynamical properties of the system, but an independent third transit observation confirming each orbit has not been obtained yet. Given the low mutual inclination of the system (below 1 deg), additional transiting planets may yet be found at periods longer than 70 days, which would correspond to orbits within or beyond the habitable zone of the star [20, 21].

From an observational point of view, HD 110067 is the brightest star found to host more than four transiting exoplanets. The current delicate configuration of the planetary orbits in HD 110067 rules out any violent event over the billion-year history of the system [22], making it a rare “fossil” [23] to study migration mechanisms and the properties of its protoplanetary disk in a pristine environment. The combination of host star brightness and the inferred presence of extended atmospheres in the majority of its planets makes HD 110067 the most favorable multi-planetary sub-Neptune system to be observed in transmission spectroscopy with JWST (Fig. <ref>). HD 110067 offers a chance to gain insight into the nature of sub-Neptune planets and where, how, and under what conditions resonant chains form and survive.

Summary of the planetary parameters of the HD 110067 system.

Parameter$^{(a)}$     HD 110067 b     HD 110067 c      HD 110067 d

4cFit parameters

$P$ (d)                              $9.113678 \pm 1\times10^{-5}$           $13.673694\pm2.4\times10^{-5}$       $20.519617\pm4\times10^{-5}$

$t_0$$^{(b)}$ (d)                          $2640.15797\pm0.00036$            $2657.45704\pm0.0007$        $2708.20282\pm0.0008$

$R_{\rm p}/R_\star$                     $0.0256\pm0.0002$               $0.0278\pm0.0003$           $0.0332\pm0.0003$

$b = (a/R_\star)\cos i_{\rm p}$         $0.355\pm0.033$                 $0.155\pm0.078$             $0.488\pm0.023$

$a/R_\star$                             $21.66\pm0.24$                   $28.37\pm0.30$                $37.21\pm0.40$

$K$ (m s$^{-1}$)                       $2.03^{+0.62}_{-0.65}$            $< 1.55$                      $2.32^{+0.89}_{-0.88}$

4cDerived physical parameters

$M_{\rm p}$ ($M_\oplus$)                $5.69_{-1.82}^{+1.78}$        $< 6.3$                      $8.52_{-3.25}^{+3.31}$

$R_{\rm p}$ ($R_\oplus$)                $2.200 \pm 0.030$             $2.388 \pm 0.036$             $2.852 \pm 0.039$

$a_{\rm p}$ (au)                        $0.0793 \pm 0.00096$          $0.1039 \pm 0.0013$          $0.1362 \pm 0.0017$

$T_\textnormal{eq}$ (K)$^{(c)}$         $800 \pm 10$                  $699 \pm 9$                  $602 \pm 8$

$t_T$ (d)$^{(d)}$                    $0.12895\pm0.00078$               $0.15435\pm0.00121$           $0.15907\pm0.00158$

$i$ (deg)      $89.061\pm0.099$       $89.687\pm0.163$      $89.248\pm0.046$

$\dagger$HD 110067 e     HD 110067 f     $\dagger$HD 110067 g

4cDerived fit parameters

$P$ (d)                              $30.793091\pm1.2\times10^{-5}$           $41.05854\pm1\times10^{-4}$            $54.76992\pm2\times10^{-4}$

$t_0$ (d)                                   $2646.0919\pm0.0011$          $2641.5763\pm0.001$          $2656.0921\pm0.002$

$R_{\rm p}/R_\star$                     $0.0226\pm0.0004$               $0.03026\pm0.00039$               $0.0303\pm0.00051$

$b = (a/R_\star)\cos i_{\rm p}$         $0.113\pm0.074$                  $0.337\pm0.043$                 $0.338\pm0.087$

$a/R_\star$                             $48.77\pm0.52$                    $59.08\pm0.52$                    $71.59\pm0.78$

$K$ (m s$^{-1}$)                       $< 0.80$                          $1.09^{+0.40}_{-0.42}$            $< 1.30$

4cDerived physical parameters

$M_{\rm p}$ ($M_\oplus$)                $< 3.9$                      $5.04_{-1.94}^{+1.89}$        $<8.4$

$R_{\rm p}$ ($R_\oplus$)                $1.940 \pm 0.040$             $2.601 \pm 0.042$             $2.607 \pm 0.052$

$a_{\rm p}$ (au)                        $0.1785 \pm 0.0022$           $0.2163 \pm 0.0026$          $0.2621 \pm 0.0032$

$T_\textnormal{eq}$ (K)                 $533 \pm 7$                   $489 \pm 6$                   $440 \pm 6$

$t_T$ (d)                            $0.20265\pm0.00244$               $0.2152\pm0.0022$               $0.2361\pm0.0052$

$i$ (deg)      $89.867\pm0.089$       $89.673\pm0.046$      $89.729\pm0.073$

(a) Error bars denote the $68\%$ posterior credibility intervals. Upper limits correspond to $99.7\%$ posterior credibility intervals.
(b) Mid-transit epoch is measured using the full baseline of transits observed with TESS, CHEOPS, and ground-based facilities. Units are expressed in terms of the TESS Julian Date, TJD = BJD$-$2457000, where BJD is the Barycentric Julian Date in units of days. Reference time system is TDB (Barycentric Dynamical Time).
(c) Equilibrium temperatures were calculated assuming zero Bond albedo and perfect energy redistribution.
(d) Transit duration $t_T$ is measured from the first to the last contact.

$\dagger$ The orbital period for planets e and g have been measured from their assigned duo-transits placing a prior based on our prediction from the resonant chain analysis.

Space-based photometry from TESS and CHEOPS of HD 110067. Transits identified for each planet after our analyses are each associated with a different color. a-f, Photometric time series of TESS Sector 23 (a,b, 2-minute cadence), TESS Sector 49 (c,d, 20-second cadence), and CHEOPS (e,f). Red points in a show the reprocessed data affected by scattered light and high levels of sky background. Time units are in TESS Julian Date $\mathrm{(TJD \equiv BJD-2457000)}$, where BJD is the Barycentric Julian Date in units of days.

Detrended light curves and the final model are shown in b,d,f. g-l, TESS phase-folded transits of each planet in the system. The transit model and binned photometry are color-coded following the same convention. For HD 110067 b (g) and d (i), CHEOPS photometry is shown atop the TESS photometry with an arbitrary offset for clarity. m, top-down view of the planetary system. Axes show the distance from the central star as a function of the semi-major axis and equilibrium temperature of the planet.

Radial velocity data from CARMENES and HARPS-N of HD 110067. Time units are the same as Fig. <ref>. a-e, radial velocity (RV) and full-width at half-maximum (FWHM) time-series after being corrected by inferred offsets. Each panel shows: RV data together with full, stellar, and planetary inferred models (a); RV data with the stellar model subtracted (b); RV residuals (c); FWHM data together with the inferred stellar model (d); and FWHM residuals (e). HARPS-N (blue) and CARMENES (orange) measurements are shown with solid circles with 1$\sigma$ error bars with a semi-transparent error bar extension accounting for the inferred jitter. The solid lines show the inferred full model coming from our multi-dimensional Gaussian process model (Methods), lightly shaded areas showing the 1 and 2$\sigma$ credibility intervals of the corresponding model. For the RV time series (a) we also show the inferred stellar (red) and planetary (green) recovered signals with an offset for better clarity. f-k, phase-folded RV signals for all the planets following the subtraction of the systemic velocities, stellar signal, and other planets. Nominal RV observations are shown as light gray points. Solid points show data binned to a tenth of the orbital phase. The inferred model is shown with a solid line following the color-coding of Fig. <ref>. The planets with mass measurement uncertainties smaller than 3$\sigma$ (f,h,j) are marked with thicker lines.

Properties of the HD 110067 system compared to the known sub-Neptune-sized planet population.  a, period-radius diagram. b, mass-radius diagram. Mass-radius relationships for different planet compositions are taken from ref. [24]. The 3$\sigma$ mass upper limits are shown for planets HD 110067 c, e, and g; 1$\sigma$ error bars for the rest. c, host star brightness in $J$-band magnitude as a function of planet equilibrium temperature showing systems with five or more transiting planets. Lower magnitudes indicate brighter stars. In all plots, point size is proportional to the planet radius while the color represents a proxy of the expected atmospheric scale height in transmission spectroscopy using JWST following the metric of ref. [25]. Where mass is not known, we use the empirical mass-radius relation of ref. [26] to compute this metric. Planet population properties were retrieved from the NASA Exoplanet Archive on May 2023.

[1]
authorHoward, A. W. et al.
titlePlanet Occurrence within 0.25 AU of Solar-type Stars
from Kepler.
journal volume201 (2),
pages15 (year2012).

[2]
authorFressin, F. et al.
titleThe False Positive Rate of Kepler and the Occurrence
of Planets.
journal volume766 (2),
pages81 (year2013).

[3]
authorBean, J. L., authorRaymond, S. N. &
authorOwen, J. E.
titleThe Nature and Origins of Sub-Neptune Size
journalJ. Geophys. Res. Planets
volume126 (1), pagese06639

[4]
authorRicker, G. R. et al.
titleTransiting Exoplanet Survey Satellite (TESS).
journalJournal of Astronomical Telescopes,
Instruments, and Systems volume1 (1),
pages014003 (year2015).

[5]
authorJenkins, J. M. et al.
editorChiozzi, G. & editorGuzman, J. C.
(eds) titleThe TESS science processing operations
(eds editorChiozzi, G. & editorGuzman,
J. C.) booktitleSoftware and Cyberinfrastructure for
Astronomy IV, Vol. volume9913, pages99133E.
organizationInternational Society for Optics and Photonics
(publisherSPIE, year2016).

[6]
authorBenz, W. et al.
titleThe CHEOPS mission.
journalExper. Astron.
volume51 (1), pages109–151

[7]
authorSinclair, A. T.
titleThe orbital resonance amongst the Galilean
satellites of Jupiter.
journal volume171,
pages59–72 (year1975).

[8]
authorMorbidelli, A.
titleModern celestial mechanics : aspects of solar
system dynamics  (publisherAdvances in Astronomy &
Astrophysics, year2002).

[9]
authorPapaloizou, J. C. B.
titleThree body resonances in close orbiting planetary
systems: tidal dissipation and orbital evolution.
journalInt. J. Astrobiology
volume14 (2), pages291–304

[10]
authorLeleu, A. et al.
titleSix transiting planets and a chain of Laplace
resonances in TOI-178.
journal volume649,
pagesA26 (year2021).

[11]
authorLuger, R. et al.
titleA seven-planet resonant chain in TRAPPIST-1.
journalNat. Astron.
volume1, pages0129 (year2017).

[12]
authorGoździewski, K., authorMigaszewski, C.,
authorPanichi, F. & authorSzuszkiewicz, E.
titleThe Laplace resonance in the Kepler-60 planetary
journal volume455 (1),
pagesL104–L108 (year2016).

[13]
authorAgol, E. et al.
titleRefining the Transit-timing and Photometric Analysis
of TRAPPIST-1: Masses, Radii, Densities, Dynamics, and Ephemerides.
journal volume2 (1),
pages1 (year2021).

[14]
authorDai, F. et al.
titleTOI-1136 is a Young, Coplanar, Aligned Planetary
System in a Pristine Resonant Chain.
journal volume165 (2),
pages33 (year2023).

[15]
authorQuirrenbach, A. et al.
editorEvans, C. J., editorBryant, J. J. &
editorMotohara, K. (eds) titleThe CARMENES
M-dwarf planet survey.
(eds editorEvans, C. J., editorBryant, J. J.
& editorMotohara, K.) booktitleGround-based and
Airborne Instrumentation for Astronomy VIII, Vol. volume11447,
pages114473C. organizationInternational Society for
Optics and Photonics (publisherSPIE, year2020).

[16]
authorCosentino, R. et al.
editorMcLean, I. S., editorRamsay, S. K. &
editorTakami, H. (eds) titleHarps-N: the new
planet hunter at TNG.
(eds editorMcLean, I. S., editorRamsay, S. K.
& editorTakami, H.) booktitleGround-based and
Airborne Instrumentation for Astronomy IV, Vol. volume8446,
pages84461V. organizationInternational Society for
Optics and Photonics (publisherSPIE, year2012).

[17]
authorHolman, M. J. & authorMurray, N. W.
titleThe Use of Transit Timing to Detect Terrestrial-Mass
Extrasolar Planets.
volume307 (5713), pages1288–1291

[18]
authorFulton, B. J. et al.
titleThe California-Kepler Survey. III. A Gap in the
Radius Distribution of Small Planets.
journal volume154,
pages109 (year2017).

[19]
authorVan Eylen, V. et al.
titleAn asteroseismic view of the radius valley: stripped
cores, not born rocky.
journal volume479,
pages4786–4795 (year2018).

[20]
authorKasting, J. F., authorWhitmire, D. P. &
authorReynolds, R. T.
titleHabitable Zones around Main Sequence Stars.
journal volume101 (1),
pages108–128 (year1993).

[21]
authorKopparapu, R. K. et al.
titleHabitable Zones around Main-sequence Stars:
Dependence on Planetary Mass.
journal volume787,
pagesL29 (year2014).

[22]
authorIzidoro, A. et al.
titleFormation of planetary systems by pebble accretion
and migration. Hot super-Earth systems from breaking compact resonant
journal volume650,
pagesA152 (year2021).

[23]
authorFabrycky, D. C. et al.
titleArchitecture of Kepler's Multi-transiting Systems.
II. New Investigations with Twice as Many Candidates.
journal volume790 (2),
pages146 (year2014).

[24]
authorZeng, L. et al.
titleGrowth model interpretation of planet size
journal volume116 (20),
pages9723–9728 (year2019).

[25]
authorKempton, E. M. R. et al.
titleA Framework for Prioritizing the TESS Planetary
Candidates Most Amenable to Atmospheric Characterization.
journal volume130 (993),
pages114401 (year2018).

[26]
authorOtegi, J. F., authorBouchy, F. &
authorHelled, R.
titleRevisited mass-radius relations for exoplanets below
120 M$_{{\ensuremath{\oplus}}}$.
journal volume634,
pagesA43 (year2020).

§ METHODS

§.§ Data

§.§.§ 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) [28, 29, 30]. Finally, the SPOC pipeline runs a wavelet-based transiting planet search for periodic exoplanets [31, 32, 33] which in the case of the S23 data revealed two threshold-crossing 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].

The 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 $\mathrm{(TJD \equiv 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.

HD 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 affected 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.

Manual vetting of both S23 and S49 PDCSAP light curves revealed unexpectedly large systematic uncertainties (Gini's mean difference, i.e., the average point-to-point absolute difference, of 215 ppm in 3-hour bins) which were absent from the uncorrected SAP flux. This has been seen in many bright stars which possess stellar variability <cit.>. In order to better correct these systematics, we performed a custom extraction of the TESS light curves for both sectors using the quaternion detrending technique against spacecraft motion developed by [40]. This involved fitting a model consisting of a linear combination of a basis spline (with breakpoints spaced every 1.5 d to model long-timescale stellar or instrumental variability) and decorrelation with parameters linked to systematic flux changes, namely the means and standard deviations of the spacecraft quaternion time series (and the squared time series) within each exposure. Using a linear least squares technique (matrix inversion), we solved for the best-fit coefficients of our free parameters while iteratively excluding $3\sigma$ outliers from the fit until convergence was reached. After calculating the best-fit model for each aperture's systematics, we then subtracted it from the uncorrected light curve and identified the aperture that produced the light curve with the lowest photometric scatter. The final light curves used in our subsequent analyses (Gini's mean difference of 130 ppm in 3-hour bins) are shown in Fig. <ref>.

Finally, as discussed above, a large portion of the data was missing from the PDCSAP light curves due to high levels of scattered light and sky background from the Earth and Moon.
The dates affected coincide with the potential transit events of planets f and g based on our dynamical model prediction. Therefore, in order to recover data affected by scattered light, we performed a custom extraction of the TESS light curves for both sectors using a pixel level decorrelation (PLD) method [41, 42, 43] implemented in the  class of the community Python package  [44]. This method employs (i) a spline polynomial fit to describe stellar variability, (ii) Principal Component Analysis (PCA) eigenmodes to model the background light, and (iii) the PLD technique to account for pointing and mechanical effects. Before applying the , we add the background flux and errors estimated by the TESS SPOC pipeline back onto the SAP light curve. Flux level, fraction, and crowding adjustments are then applied to the corrected light curve. To automatically optimize the selection of parameter values for the PLD corrector, we evaluate the resulting light curve using the Savitsky-Golay Combined Differential Photometric Precision (sgCDPP) proxy algorithm [45, 46] implemented in . For a grid of PLD corrector parameter values, we calculate the harmonic mean of these quantities and select the corrected light curve that minimizes it. We use this data for the cadences missing in the quaternion-detrended light curve in our final analyses (marked with a different color in Fig. <ref>).

§.§.§ 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 and 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 [49, 39, 50]. An observing log summarizing the duration of each visit, its average observing efficiency (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.

To 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  package [51, 52]. For bright targets such as HD 110067, light curves generated with  exhibit lower median absolute differences 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 inter-visit flux differences, we normalized the entire CHEOPS data together instead of individually. This revealed clear visit-to-visit flux differences 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. <ref> and Fig. S1 in the Supplementary.

§.§.§ 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 different 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. <ref> shows the data and best-fit models as discussed in Sect. <ref>.

Teide Observatory
We observed HD 110067 on May 23rd, 2022 using the MuSCAT2 instrument installed at the 1.5-m Telescopio Carlos Sánchez (TCS) located at the Teide Observatory, Spain [54]. The images were taken simultaneously in $g$, $r$, $i$, and $z_s$ 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.\!\!^{\prime\prime}1$ to $11.\!\!^{\prime\prime}3$ depending on the band. Note that 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.

We 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\times4096$ pix SINISTRO camera. The images were calibrated by the standard LCO BANZAI pipeline [57]. Differential 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 effect on the photometry. Observing without using the shutters is now the standard operation mode of NGTS. We performed standard differential 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 flat-fielded since knowledge of the flat-field was incomplete at the time. The RMS of the 15-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 affected by cirrus.

San Pedro Mártir Observatory
We observed HD 110067 on May 24th, 2022 with the 1-m SAINT-EX telescope at the Observatorio Astronómico Nacional de la Sierra de San Pedro Mártir 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, flat-fielding, 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 affected 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 $z_s$ 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.\!\!^{\prime\prime}1$ to $11.\!\!^{\prime\prime}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.

§.§.§ 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 ($\lesssim 1\arcsec$) 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-$\gamma$ filter. Each dither position was observed three times, offset 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 of $0.025\arcsec$ per pixel for a total field of view of $\sim25\arcsec$. The sensitivities of the final combined AO image were determined by injecting simulated sources azimuthally around the primary target every $20^\circ $ at separations of integer multiples of the central source's FWHM [67]. The Palomar data have a sensitivity $\Delta {\rm mag} = 2$ at 0.1 and $\Delta {\rm mag} = 9$ at 1; the final sensitivity curve is shown in Fig. S2 of the Supplementary.

Gemini Observatory
We observed 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$\sigma$ sensitivity of $\Delta {\rm mag} = 7$ (blue channel) and $\Delta {\rm 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 $\lesssim 1.4$ indicate that the Gaia astrometric solution is consistent with the star being single whereas RUWE values $\gtrsim 1.4$ may indicate an astrometric excess noise, possibly caused the presence of an unseen companion <cit.>.  HD 110067 has a Gaia EDR3 RUWE value of 0.94 indicating that the astrometric fit is consistent with a single-star model.

§.§.§ Radial velocity monitoring

Calar Alto Observatory
We observed HD 110067 using the CARMENES instrument [15] installed at the 3.5-m 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  [73] and  [74]. While the mean internal precision of the template matching  RVs is $3.1\,\mathrm{m\,s^{-1}}$, the precision of the cross-correlation method  RVs is $2.9\,\mathrm{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 La 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é) 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  [74]. Both the YABI- and -derived radial velocities have a median internal precision of $1.0\,\mathrm{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.

§.§ Stellar parameters

§.§.§ 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  tools to our co-added HARPS-N spectra [76, 77, 78].  [79, 80] measures the equivalent widths of iron lines in the spectrum that are converted into stellar atmospheric parameters using the  radiative transfer code [81] applied to Kurucz model atmospheres [82]. In Extended Data Table 3 we report the effective temperature $T_{\mathrm{eff}}$, 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  [83].

We further study the photospheric parameters by conducting a classical curve-of-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  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.

§.§.§ Physical parameters

Using our spectral parameters and the  [82, 86] and  [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 effective 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 offset-corrected parallax [91] with model uncertainties accounted for using a Bayesian modeling averaging. We report the stellar radius $R_\star$ in Extended Data Table 3.

Last, 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_{\mathrm{eff}}$, $\log g$, and $R_\star$ [92]. On the one hand, we use an isochrone placement algorithm [93, 94] and interpolate over pre-computed grids of  [95] isochrones. On the other hand, we use the Code Liègeois d’Évolution 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.

The [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\sigma$ of [Fe/H], show that the star is not enhanced in $\alpha$-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 \pm 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.

§.§ Analysis

§.§.§ 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 -detrended CHEOPS data for the three visits containing transits. We built transit models for the six planets with  [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\Phi}$, 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 <cit.>. 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.

Informative 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 offsets. The impact parameter and radius ratio are fitted from a broad uniform and log-normal prior, respectively, 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, 104, 105]. The prior and posterior distributions of each parameter in the model are shown in Table S1 of the Supplementary.

§.§.§ 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  [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. <ref> shows this result.

We then modeled both duo- and single-transits using  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) <cit.>. 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 \pm 8.0$ days, respectively.

§.§.§ 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.

For transiting systems, generalized three-body Laplace angles can be estimated in 0th order in eccentricity, defined as $\Psi_{e=0}$, from the times of mid-transit and the orbital period of the planets <cit.>. This estimation differs from the actual generalized three-body Laplace angle proportionally to the eccentricities <cit.>. Interestingly, for known systems with a chain of three-body resonances, all $\Psi_{e=0}$ lie close to an equilibrium of the chain, as seen in Extended Data Fig. <ref>. The largest distance is $\sim 43$ degrees for the inner triplet of K2-138 [112]. For HD 110067, the estimated angle $\Psi_{e=0,bcd}$ is also at about $\sim 44$ degrees from its theorized 180-degree 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.

As shown above, the two events at 2646.088 TJD and 1937.851 TJD have fully consistent shapes. Among the probable periods computed with  (Fig. S3),
$P_e=30.7931$ days is the only one that continues the resonant chain, with $P_e/P_d=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 $\Psi_{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.

For 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. <ref>, all known chains of Laplace resonances have either their estimated generalized three-body Laplace angle $\Psi_{e=0}$, or their actual generalized three-body Laplace angle $\Psi$ 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 $\Delta\Psi=\lvert \Psi_{e=0}-\Psi_{eq} \rvert$ is given in Extended Data Table 4. The case A2, with $P_f/P_e=4/3$ and $P_g/P_f=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).

For completeness, we study the role that the eccentricity of the orbit plays in the prediction. To estimate the generalized three-body Laplace angle

\begin{equation}
\begin{aligned}
\Psi = l \lambda_1 - (l+m) \lambda_2 + m \lambda_3 , \\
\end{aligned}
\label{eq:LaplaceAngle}
\end{equation}

at a given epoch, we estimate the value of the $\lambda_j$ as follow.
Transits occur when the true longitude of the planet is equal to $l_0=-\pi/2$. At first order in the eccentricity,

\begin{equation}
\begin{aligned}
\lambda_0=l_0-2e\sin (l_0-\varpi) = -\frac{\pi}{2} + 2 e\cos(\varpi) \, .
\end{aligned}
\label{eq:lambdafl}
\end{equation}

We 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$, $\lambda_0=-\pi/2$. We hence obtain

\begin{equation}
\begin{aligned}
\lambda (t)= -  \frac{\pi}{2} + \frac{(t-t_0)}{P} 2\pi\, .
\end{aligned}
\label{eq:lambdaft}
\end{equation}

The error on $\Psi$ made by assuming zero eccentricity is hence, at first order [111]:

\begin{equation}
\begin{aligned}
\lvert \Psi - \Psi_{e=0} \rvert = \lvert 2 l e_1\cos \varpi_1  - 2 (l+m)  e_2\cos \varpi_2 + 2 m  e_3\cos \varpi_3 \rvert.
\end{aligned}
\label{eq:LaplaceDist}
\end{equation}

This 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 $\Psi_{e=0}$ closer to the equilibrium than $\Psi$ actually is, or vice-versa. We check this for the cases presented in Extended Data Table 4.

Each 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 $k_i=e_i \cos \varpi_i$ and $h_i=e_i \sin \varpi_i$, we minimize the cost function

\begin{equation}
\begin{aligned}
C=\mathcal{{A}} (\Psi_{bcd}) + \mathcal{{A}} (\Psi_{cde}) + \mathcal{{A}} (\Psi_{efg}) + \mathcal{{A}} (\Psi_{fgh})
\end{aligned}
\label{eq:cost}
\end{equation}

over 200 years, where $\mathcal{{A}} (\Psi)=2\pi$ if $\Psi$ circulates, and $\mathcal{{A}} (\Psi)$ is the peak-to-peak amplitude of libration of $\Psi$ otherwise. For each case, 40 MCMC runs are conducted to minimize $C$, using  [113] for the $N$-body integration and  [114] for the MCMC. For each run, the $k_i$ and $h_i$ 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.

The best solution of each fit is shown in Extended Data Fig. <ref>. 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 $k_i$ and $k_j$ parameters below $85$ deg of amplitude on average, with the exception of case A0 for which an average of $\approx 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 $k_i$ and $k_j$ parameters. The case A2, with $P_f/P_e=4/3$ and $P_g/P_f=4/3$, is hence our prediction for the outer architecture of the HD 110067 system.

§.§.§ Confirming the predictions

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 affected 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 affected were flagged by SPOC, thus not leaving enough valid data to derive cotrending basis vectors and missing in the PDCSAP light curve.

Our custom extraction using the PLD method from Sect. <ref> was able to recover the missing data, showing two mono-transits at 1943.6 and 1944.1 TJD (Fig. <ref> and also Fig. S6 in the Supplementary).
Using , 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 <cit.>. In all cases, the metadata were normalized to a time series with $\mu=0$, $\sigma=1$ and modeled using a single scaling parameter with a normal prior of $\mu=0$, $\sigma=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  [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 telescopes, which were decorrelated independently. The final result is 24 individual photometric time series. An offset was also applied to each light curve, as well as a single global slope parameter to include the possibility of stellar activity.

The transit parameters were constrained based on those found in a fit of the TESS S49 mono-transit. The predicted period used was $4/3 \times P_e = 41.051\pm0.1$ d, with the uncertainty implying a divergence from the perfect integer period of $2.4\times10^{-3}$ — larger than those values found for the inner three planets. We limited the period to $41.0\pm0.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  [99] and optimized using  [117], specifically with the  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  [120].

Our results show a preference for a transit at the expected period $P=41.04\pm0.01$ days, with a $\Delta$WAIC of 9.5 over a transit-free model, as can be seen in Extended Data Fig. <ref>. The majority of instruments showed a weak preference for a peak at $P\sim41.05$d, with the exception of LCO (which observed no in-transit data) and SAINT-EX (which is the most affected by cirrus). This is only equivalent to moderate evidence for a $\sim41.051$d period of HD110067 f. The lower two panels of Extended Data Fig. <ref> show that both models (with or without transit) fit reasonably well. This is in part because systematic effects 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 $\sim40.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, differential line width, Mount Wilon's S-index, H$\alpha$ 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 <cit.>, 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 effects on the detection and mass determination of the planets in the system.

Method I: SN-fit 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-half-maximum (FWHM) and the asymmetry of the cross-correlation function (CCF) computed from the spectra <cit.>. 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 $\gamma$). For each observation, through the SN-fit we were able to retrieve the stellar radial velocity ($\overline{RV}$, quantified through the SN median), the $\mathrm{FWHM_{SN}}$, the contrast $A$, and the asymmetry $\gamma$. The errors $\sigma_{\mathrm{RV}}$ of the $\overline{RV}$ measurements were inferred using a bootstrap approach. Denoting with $f_{CCF}$ the flux of a CCF data point, each point was perturbed by sampling values from a Normal distribution whose standard deviation is equal to $\sqrt{f_{\mathrm{CCF}}}$, since the errors affecting the CCF data points are expected to be Poissonian.

After 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 [$\mathrm{FWHM_{SN}}$, $A$, $\gamma$] 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  code [132], where we switched off 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

\begin{multline}
RV_{\star} = \beta_0 + \displaystyle\sum_{k=1}^{k_t}{\beta_{k,t} t^k} + \displaystyle\sum_{k=1}^{k_F}{\beta_{k,F} \mathrm{FWHM_{SN}}^k} + \displaystyle\sum_{k=1}^{k_A}{\beta_{k,A} A^k} + \\
+ \displaystyle\sum_{k=1}^{k_{\gamma}}{\beta_{k,\gamma} \gamma^k} + \displaystyle\sum_{k=1}^{k_R}{\beta_{k,R} \log{R'_{\mathrm{HK}}}^k},
\label{eq:RVactivity}
\end{multline}

where ($k_t$, $k_F$, $k_A$, $k_{\gamma}$, $k_R$) is the vector of the polynomial orders whose optimal value has been established by launching several MCMC preliminary runs and selecting that combination which produces the minimum Bayesian Information Criterion (BIC) [133].

After performing a longer  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$\sigma$ 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$\sigma$ level (planets d and f) are shown in Fig. S9 of the Supplementary.

Method 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 <cit.>.
We create $N-$dimensional GP models, including $N$ time-series $\mathcal{A}_i$,  as

\begin{equation}
\begin{matrix}
\mathcal{A}_1 =  A_{1} G(t) + B_{1} \dot{G}(t) \\
\vdots \\
\mathcal{A}_N =  A_{N} G(t) + B_{N} \dot{G}(t), \\
\end{matrix}
\label{eq:gps}
\end{equation}

where the variables $A_{1}$, $B_{1}$, $\cdots$, $A_{N}$, $B_{N}$, are free parameters which relate the individual time-series to $G(t)$ and $\dot{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.

We model the stellar signal using a GP whose covariance between two times $t_i$ and $t_j$ is given by

\begin{equation}
\gamma_{{\rm QP},i,j} = \exp
\left[
- \frac{\sin^2[\pi(t_i - t_j)/P_{\rm GP}]}{2 \lambda_{\rm P}^2}
- \frac{(t_i - t_j)^2}{2\lambda_{\rm e}^2}
\right],
\label{eq:gamma}
\end{equation}

where $\gamma_{{\rm QP},i,j}$ is the Quasi-Periodic (QP) kernel, whose hyperparameters are, , the GP characteristic period, , the inverse of the harmonic complexity, and , the long term evolution timescale.

We 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. (<ref>) and its derivatives [135, 136].  We assume that RVs can be described as $\mathcal{A}_i =  A_{i} G(t) + B_{i} \dot{G}(t)$, while the FWHM time series is described as $\mathcal{A}_i =  A_{i} G(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 offset, noting that we include a different offset per instrument. We also include a jitter term per time series and per instrument to account for unaccounted systematic errors.

We perform MCMC samplings of the parameter space using the code  [139, 136]. 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 <ref>
shows the spectroscopic time series resulting from this joint analysis. Median values along with their 1$\sigma$ uncertainties are reported in Table S4 of the Supplementary.

Modeling techniques employing GP are particularly subject to overfitting giving their flexibility to reproduce the data <cit.>. To test the robustness of our GP model, we carried out a cross-validation analysis. We repeat the two-dimensional 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 CARMENES data with median offsets 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.

Both the two methods clearly detect planets HD 110067 d and f, Method II also detects planet b, however, the detection levels slightly differ. Recalling that the two different techniques are based on different 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_{\rm f} = 2.0\pm1.0\,\mathrm{m\,s^{-1}}$). On the other hand, the RV semi-amplitudes inferred from the two methods shown in the manuscript are compatible within $\sim$ 1.5 $\sigma$, with the statistical tension $\Delta_{\mathrm{I-II}}$ below $\sim$ 1$\sigma$ for planets b, c, d, and g. Extended Data Fig. <ref> displays the pairs of posterior density functions for the RV semi-amplitudes of each planet for comparison.

For HD 110067 f, after imposing a Gaussian prior centered around 41.05 days in both Method I and Method II RV models, we indeed recover a significant RV signal with a detection level of $\sim$ 3$\sigma$. 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 $P_{f,\mathrm{uni}}=40.2\pm0.2$ d. Thus, the RV data independently suggest the presence of a putative planet having a period close to 41.05 d. The tension at the $\sim$ 4$\sigma$ 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.

We finally 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 different 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.

§.§.§ 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 Method 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. <ref> for the TESS and CHEOPS photometry and Fig. <ref> for the radial velocities.

§.§.§ 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 Method I and Method 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 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 effect) 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.

The results from this analysis are shown in Extended Data Fig. <ref>. Our model shows that the gas mass content of all planets is of the order of $10^{-3} M_\oplus$ to $10^{-1} M_\oplus$ (median value, see Table S5),
with the notable exception of HD 110067 e (median value of $\sim 10^{-7} M_\oplus$ using the masses from Method I, $\sim 10^{-3} M_\oplus$ using Method 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, 151, 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].

[27]
authorStassun, K. G. et al.
titleThe TESS Input Catalog and Candidate Target List.
journal volume156 (3),
pages102 (year2018).

[28]
authorStumpe, M. C. et al.
titleKepler Presearch Data Conditioning
I—Architecture and Algorithms for Error Correction in Kepler
Light Curves.
journal volume124 (919),
pages985 (year2012).

[29]
authorStumpe, M. C. et al.
titleMultiscale Systematic Error Correction via
Wavelet-Based Bandsplitting in Kepler Data.
journal volume126,
pages100 (year2014).

[30]
authorSmith, J. C. et al.
titleKepler Presearch Data Conditioning II - A Bayesian
Approach to Systematic Error Correction.
journal volume124,
pages1000 (year2012).

[31]
authorJenkins, J. M.
titleThe Impact of Solar-like Variability on the
Detectability of Transiting Terrestrial Planets.
journal volume575,
pages493–505 (year2002).

[32]
authorJenkins, J. M. et al.
editorRadziwill, N. M. & editorBridger,
A. (eds) titleTransiting planet search in the Kepler
(eds editorRadziwill, N. M. &
editorBridger, A.) booktitleSoftware and
Cyberinfrastructure for Astronomy, Vol. volume7740 of
seriesSociety of Photo-Optical Instrumentation Engineers
(SPIE) Conference Series, pages77400D (year2010).

[33]
authorJenkins, J. M. et al.
titleKepler Data Processing Handbook: Transiting Planet
howpublishedKepler Science Document KSCI-19081-003

[34]
authorTwicken, J. D. et al.
titleKepler Data Validation I - Architecture, Diagnostic
Tests, and Data Products for Vetting Transiting Planet Candidates.
journal volume130 (6),
pages064502 (year2018).

[35]
authorLi, J. et al.
titleKepler Data Validation II-Transit Model Fitting and
Multiple-planet Search.
journal volume131 (996),
pages024506 (year2019).

[36]
authorGuerrero, N. M. et al.
titleThe TESS Objects of Interest Catalog from the TESS
Prime Mission.
journal volume254 (2),
pages39 (year2021).

[37]
authorFausnaugh, M. M., authorBurke, C. J.,
authorRicker, G. R. & authorVanderspek, R.
titleCalibrated Full-frame Images for the TESS Quick Look
journalRes. Not. Amer. Astron. Soc. volume4 (12), pages251

[38]
authorHedges, C. et al.
titleTOI-2076 and TOI-1807: Two Young, Comoving Planetary
Systems within 50 pc Identified by TESS that are Ideal Candidates for Further
Follow Up.
journal volume162 (2),
pages54 (year2021).

[39]
authorOsborn, H. et al.
titleTwo Warm Neptunes transiting HIP 9618 revealed by
TESS & Cheops.
journal volume523 (2),
pages3069–3089 (year2023).

[40]
authorVanderburg, A. et al.
titleTESS Spots a Compact System of Super-Earths around
the Naked-eye Star HR 858.
journal volume881 (1),
pagesL19 (year2019).

[41]
authorDeming, D. et al.
titleSpitzer Secondary Eclipses of the Dense,
Modestly-irradiated, Giant Exoplanet HAT-P-20b Using Pixel-level
journal volume805 (2),
pages132 (year2015).

[42]
authorLuger, R. et al.
titleEVEREST: Pixel Level Decorrelation of K2 Light
journal volume152 (4),
pages100 (year2016).

[43]
authorLuger, R. et al.
titlestarry: Analytic Occultation Light Curves.
journal volume157,
pages64 (year2019).

[44]
authorLightkurve Collaboration et al.
titleLightkurve: Kepler and TESS time series analysis in
howpublishedAstrophysics Source Code Library

[45]
authorGilliland, R. L. et al.
titleKepler Mission Stellar and Instrument Noise
journal volume197 (1),
pages6 (year2011).

[46]
authorVan Cleve, J. E. et al.
titleThat's How We Roll: The NASA K2 Mission Science
Products and Their Performance Metrics.
journal volume128 (965),
pages075002 (year2016).

[47]
authorSchanche, N. et al.
titleTOI-2257 b: A highly eccentric long-period
sub-Neptune transiting a nearby M dwarf.
journal volume657,
pagesA45 (year2022).

[48]
authorUlmer-Moll, S. et al.
titleTwo long-period transiting exoplanets on eccentric
orbits: NGTS-20 b (TOI-5152 b) and TOI-5153 b.
journal volume666,
pagesA46 (year2022).

[49]
authorOsborn, A. et al.
titleTOI-431/HIP 26013: a super-Earth and a sub-Neptune
transiting a bright, early K dwarf, with a third RV planet.
journal volume507 (2),
pages2782–2803 (year2021).

[50]
authorTuson, A. et al.
titleTESS and CHEOPS discover two warm sub-Neptunes transiting the bright K-dwarf HD 15906.
journal volume523 (2),
pages3090–3118 (year2023).

[51]
authorSzabó, G. M. et al.
titleThe changing face of AU Mic b: stellar spots,
spin-orbit commensurability, and transit timing variations as seen by CHEOPS
and TESS.
journal volume654,
pagesA159 (year2021).

[52]
authorMorris, B. M. et al.
titleCHEOPS precision phase curve of the Super-Earth 55
Cancri e.
journal volume653,
pagesA173 (year2021).

[53]
authorHoyer, S. et al.
titleExpected performances of the Characterising
Exoplanet Satellite (CHEOPS). III. Data reduction pipeline: architecture and
simulated performances.
journal volume635,
pagesA24 (year2020).

[54]
authorNarita, N. et al.
titleMuSCAT2: four-color simultaneous camera for the
1.52-m Telescopio Carlos Sánchez.
journalJ. Astron. Telesc.
Instrum. Syst. volume5 (1),
pages015001 (year2019).

[55]
authorParviainen, H. et al.
titleMuSCAT2 multicolour validation of TESS candidates:
an ultra-short-period substellar object around an M dwarf.
journal volume633,
pagesA28 (year2020).

[56]
authorBrown, T. M. et al.
titleLas Cumbres Observatory Global Telescope Network.
journal volume125 (931),
pages1031 (year2013).

[57]
authorMcCully, C. et al.
titleReal-time processing of the imaging data from the
network of Las Cumbres Observatory Telescopes using BANZAI.
journal volume10707,
pages107070K (year2018).

[58]
authorCollins, K. A., authorKielkopf, J. F.,
authorStassun, K. G. & authorHessman, F. V.
titleAstroImageJ: Image Processing and Photometric
Extraction for Ultra-precise Astronomical Light Curves.
journal volume153,
pages77 (year2017).

[59]
authorWheatley, P. J. et al.
titleThe Next Generation Transit Survey (NGTS).
journal volume475 (4),
pages4476–4493 (year2018).

[60]
authorGarcia-Mejia, J. et al.
editorMarshall, H. K., editorSpyromilio,
J. & editorUsuda, T. (eds) titleThe Tierras
Observatory: An ultra-precise photometer to characterize nearby terrestrial
(eds editorMarshall, H. K.,
editorSpyromilio, J. & editorUsuda, T.)
booktitleGround-based and Airborne Telescopes VIII, Vol.
volume11445 of seriesSociety of Photo-Optical
Instrumentation Engineers (SPIE) Conference Series,
pages114457R (year2020).

[61]
authorDemory, B. O. et al.
titleA super-Earth and a sub-Neptune orbiting the bright,
quiet M3 dwarf TOI-1266.
journal volume642,
pagesA49 (year2020).

[62]
authorNarita, N. et al.
editorEvans, C. J., editorBryant, J. J. &
editorMotohara, K. (eds) titleMuSCAT3: a
4-color simultaneous camera for the 2m Faulkes Telescope North.
(eds editorEvans, C. J., editorBryant, J. J.
& editorMotohara, K.) booktitleGround-based and
Airborne Instrumentation for Astronomy VIII, Vol. volume11447,
pages114475K. organizationInternational Society for
Optics and Photonics (publisherSPIE, year2020).

[63]
authorFukui, A. et al.
titleMeasurements of Transit Timing Variations for
journal volume63,
pages287–300 (year2011).

[64]
authorCiardi, D. R., authorBeichman, C. A.,
authorHorch, E. P. & authorHowell, S. B.
titleUnderstanding the Effects of Stellar Multiplicity on
the Derived Planet Radii from Transit Surveys: Implications for Kepler, K2,
and TESS.
journal volume805 (1),
pages16 (year2015).

[65]
authorHayward, T. L. et al.
titlePHARO: A Near-Infrared Camera for the Palomar
Adaptive Optics System.
journal volume113 (779),
pages105–118 (year2001).

[66]
authorDekany, R. et al.
titlePALM-3000: Exoplanet Adaptive Optics for the 5 m
Hale Telescope.
journal volume776 (2),
pages130 (year2013).

[67]
authorFurlan, E. et al.
titleThe Kepler Follow-up Observation Program. I. A
Catalog of Companions to Kepler Stars from High-Resolution Imaging.
journal volume153 (2),
pages71 (year2017).

[68]
authorScott, N. J. et al.
titleTwin High-resolution, High-speed Imagers for the
Gemini Telescopes: Instrument description and science verification results.
journalFrontiers in Astronomy and Space Sciences
volume8, pages138 (year2021).

[69]
authorHowell, S. B., authorEverett, M. E.,
authorSherry, W., authorHorch, E. &
authorCiardi, D. R.
titleSpeckle Camera Observations for the NASA Kepler
Mission Follow-up Program.
journal volume142 (1),
pages19 (year2011).

[70]
authorMugrauer, M. & authorMichel, K.-U.
titleGaia search for stellar companions of TESS Objects
of Interest.
journalAstronomische Nachrichten
volume341 (10), pages996–1030

[71]
authorMugrauer, M. & authorMichel, K.-U.
titleGaia search for stellar companions of TESS Objects
of Interest II.
journalAstronomische Nachrichten
volume342 (6), pages840–864

[72]
authorZiegler, C. et al.
titleSOAR TESS Survey. I. Sculpting of TESS Planetary
Systems by Stellar Companions.
journal volume159 (1),
pages19 (year2020).

[73]
authorLafarga, M. et al.
titleThe CARMENES search for exoplanets around M dwarfs.
Radial velocities and activity indicators from cross-correlation functions
with weighted binary masks.
journal volume636,
pagesA36 (year2020).

[74]
authorZechmeister, M. et al.
titleSpectrum radial velocity analyser (SERVAL).
High-precision radial velocities and two alternative spectral indicators.
journal volume609,
pagesA12 (year2018).

[75]
authorCosentino, R. et al.
editorRamsay, S. K., editorMcLean, I. S. &
editorTakami, H. (eds) titleHARPS-N @ TNG, two
year harvesting data: performances and results.
(eds editorRamsay, S. K., editorMcLean, I. S.
& editorTakami, H.) booktitleGround-based and
Airborne Instrumentation for Astronomy V, Vol. volume9147,
pages91478C. organizationInternational Society for
Optics and Photonics (publisherSPIE, year2014).

[76]
authorSantos, N. C. et al.
titleSWEET-Cat: A catalogue of parameters for Stars With
ExoplanETs. I. New atmospheric parameters and masses for 48 stars with
journal volume556,
pagesA150 (year2013).

[77]
authorSousa, S. G.
titleARES + MOOG: A Practical Overview of an
Equivalent Width (EW) Method to Derive Stellar Parameters,
pages297–310 (publisherSpringer,

[78]
authorSousa, S. G. et al.
titleSWEET-Cat 2.0: The Cat just got SWEETer. Higher
quality spectra and precise parallaxes from Gaia eDR3.
journal volume656,
pagesA53 (year2021).

[79]
authorSousa, S. G., authorSantos, N. C.,
authorIsraelian, G., authorMayor, M. &
authorMonteiro, M. J. P. F. G.
titleA new code for automatic determination of equivalent
widths: Automatic Routine for line Equivalent widths in stellar Spectra
journal volume469 (2),
pages783–791 (year2007).

[80]
authorSousa, S. G., authorSantos, N. C.,
authorAdibekyan, V., authorDelgado-Mena, E. &
authorIsraelian, G.
titleARES v2: new features and improved performance.
journal volume577,
pagesA67 (year2015).

[81]
authorSneden, C. A.
titleCarbon and Nitrogen Abundances in Metal-Poor
Ph.D. thesis, schoolTHE UNIVERSITY OF TEXAS AT AUSTIN.

[82]
authorKurucz, R. L.
titleSYNTHE spectrum synthesis programs and line
data  (publisherAstrophysics Source Code Library,

[83]
authorBrahm, R., authorJordán, A.,
authorHartman, J. & authorBakos, G.
titleZASPE: A Code to Measure Stellar Atmospheric
Parameters and their Covariance from Spectra.
journal volume467 (1),
pages971–984 (year2017).

[84]
authorAdibekyan, V. Z. et al.
titleChemical abundances of 1111 FGK stars from the HARPS
GTO planet search program. Galactic stellar populations and planets.
journal volume545,
pagesA32 (year2012).

[85]
authorAdibekyan, V. et al.
titleIdentifying the best iron-peak and
-capture elements for chemical tagging: The impact of
the number of lines on measured scatter.
journal volume583,
pagesA94 (year2015).

[86]
authorCastelli, F. & authorKurucz, R. L.
editorPiskunov, N., editorWeiss, W. W. &
editorGray, D. F. (eds) titleNew Grids of
ATLAS9 Model Atmospheres.
(eds editorPiskunov, N., editorWeiss,
W. W. & editorGray, D. F.)
booktitleModelling of Stellar Atmospheres, Vol.
volume210 of seriesIAU Symposium,
pagesA20 (year2003).

[87]
authorAllard, F.
editorBooth, M., editorMatthews, B. C. &
editorGraham, J. R. (eds) titleThe BT-Settl
Model Atmospheres for Stars, Brown Dwarfs and Planets.
(eds editorBooth, M., editorMatthews,
B. C. & editorGraham, J. R.)
booktitleExploring the Formation and Evolution of Planetary
Systems, Vol. volume299, pages271–272

[88]
authorBlackwell, D. E. & authorShallis, M. J.
titleStellar angular diameters from infrared photometry.
Application to Arcturus and other stars; with effective temperatures.
journal volume180,
pages177–191 (year1977).

[89]
authorSchanche, N. et al.
titleWASP-186 and WASP-187: two hot Jupiters discovered
by SuperWASP and SOPHIE with additional observations by TESS.
journal volume499 (1),
pages428–440 (year2020).

[90]
authorWilson, T. G. et al.
titleA pair of sub-Neptunes transiting the bright K-dwarf
TOI-1064 characterized with CHEOPS.
journal volume511 (1),
pages1043–1071 (year2022).

[91]
authorLindegren, L. et al.
titleGaia Early Data Release 3. Parallax bias versus
magnitude, colour, and position.
journal volume649,
pagesA4 (year2021).

[92]
authorBonfanti, A. et al.
titleCHEOPS observations of the HD 108236 planetary
system: a fifth planet, improved ephemerides, and planetary radii.
journal volume646,
pagesA157 (year2021).

[93]
authorBonfanti, A., authorOrtolani, S.,
authorPiotto, G. & authorNascimbeni, V.
titleRevising the ages of planet-hosting stars.
journal volume575,
pagesA18 (year2015).

[94]
authorBonfanti, A., authorOrtolani, S. &
authorNascimbeni, V.
titleAge consistency between exoplanet hosts and field
journal volume585,
pagesA5 (year2016).

[95]
authorMarigo, P. et al.
titleA New Generation of PARSEC-COLIBRI Stellar
Isochrones Including the TP-AGB Phase.
journal volume835,
pages77 (year2017).

[96]
authorScuflaire, R. et al.
titleCLÉS, Code Liégeois d'Évolution
journal volume316,
pages83–91 (year2008).

[97]
authorSalmon, S. J. A. J., authorVan Grootel, V.,
authorBuldgen, G., authorDupret, M. A. &
authorEggenberger, P.
titleReinvestigating  Centauri AB in
light of asteroseismic forward and inverse methods.
journal volume646,
pagesA7 (year2021).

[98]
authorReddy, B. E., authorLambert, D. L. &
authorAllende Prieto, C.
titleElemental abundance survey of the Galactic thick
journal volume367 (4),
pages1329–1366 (year2006).

[99]
authorForeman-Mackey, D. et al.
titledfm/exoplanet: exoplanet v0.2.1

[100]
authorDelrez, L. et al.
titleTransit detection of the long-period volatile-rich
super-Earth $^{2}$ Lupi d with CHEOPS.
journalNat. Astron.
volume5, pages775–787

[101]
authorClaret, A.
titleA new method to compute limb-darkening coefficients
for stellar atmosphere models with spherical symmetry: the space missions
TESS, Kepler, CoRoT, and MOST.
journal volume618,
pagesA20 (year2018).

[102]
authorClaret, A.
titleLimb and Gravity-darkening Coefficients for the
Space Mission CHEOPS.
journalResearch Notes of the American Astronomical
Society volume5 (1), pages13

[103]
authorVan Eylen, V. & authorAlbrecht, S.
titleEccentricity from Transit Photometry: Small Planets
in Kepler Multi-planet Systems Have Low Eccentricities.
journal volume808 (2),
pages126 (year2015).

[104]
authorXie, J.-W. et al.
titleExoplanet orbital eccentricities derived from
LAMOST-Kepler analysis.
journal volume113 (41),
pages11431–11435 (year2016).

[105]
authorHadden, S. & authorLithwick, Y.
titleKepler Planet Masses and Eccentricities from TTV
journal volume154 (1),
pages5 (year2017).

[106]
authorOsborn, H. P.
titleMonoTools: Planets of uncertain periods detector and
howpublishedAstrophysics Source Code Library, record
ascl:2204.020 (year2022).

[107]
authorKipping, D.
titleThe Orbital Period Prior for Single Transits.
journalResearch Notes of the American Astronomical
Society volume2 (4), pages223

[108]
authorVan Eylen, V. et al.
titleThe orbital eccentricity of small planet systems.
journalThe Astronomical Journal
volume157 (2), pages61
(year2019) .

[109]
authorOsborn, H. P. et al.
titleUncovering the true periods of the young
sub-Neptunes orbiting TOI-2076.
journal volume664,
pagesA156 (year2022).

[110]
authorMills, S. M. et al.
titleA resonant chain of four transiting, sub-Neptune
journal volume533 (7604),
pages509–512 (year2016).

[111]
authorSiegel, J. C. & authorFabrycky, D.
titleResonant Chains of Exoplanets: Libration Centers for
Three-body Angles.
journal volume161 (6),
pages290 (year2021).

[112]
authorLopez, T. A. et al.
titleExoplanet characterisation in the longest known
resonant chain: the K2-138 system seen by HARPS.
journal volume631,
pagesA90 (year2019).

[113]
authorRein, H. & authorLiu, S. F.
titleREBOUND: an open-source multi-purpose N-body code
for collisional dynamics.
journal volume537,
pagesA128 (year2012).

[114]
authorDelisle, J.-B.
titlesamsam: Scaled Adaptive Metropolis SAMpler.
howpublishedAstrophysics Source Code Library, record
ascl:2207.011 (year2022).

[115]
authorLeleu, A. et al.
titleRemoving biases on the density of sub-Neptunes
characterised via transit timing variations. Update on the mass-radius
relationship of 34 Kepler planets.
journal volume669,
pagesA117 (year2023).

[116]
authorParviainen, H. & authorAigrain, S.
titleldtk: Limb Darkening Toolkit.
journal volume453 (4),
pages3821–3826 (year2015).

[117]
authorSalvatier, J., authorWiecki, T. V. &
authorFonnesbeck, C.
titleProbabilistic programming in python using pymc3.
journalPeerJ Comput. Sci.
volume2, pagese55 (year2016) .

[118]
authorWatanabe, S. & authorOpper, M.
titleAsymptotic equivalence of bayes cross validation and
widely applicable information criterion in singular learning theory.
journalJ. Mach. Learn. Res.
volume11 (12) (year2010) .

[119]
authorVehtari, A., authorGelman, A. &
authorGabry, J.
titlePractical bayesian model evaluation using
leave-one-out cross-validation and waic.
journalStat. Comput.
volume27, pages1413–1432

[120]
authorArviZ Developers.
titleArviZ: Exploratory analysis of Bayesian models.
howpublishedAstrophysics Source Code Library, record
ascl:2004.012 (year2020).

[121]
authorZechmeister, M. & authorKürster, M.
titleThe generalised Lomb-Scargle periodogram. A new
formalism for the floating-mean and Keplerian periodograms.
journal volume496,
pages577–584 (year2009).

[122]
authorSaar, S. H. & authorDonahue, R. A.
titleActivity-Related Radial Velocity Variation in Cool
journal volume485 (1),
pages319–327 (year1997).

[123]
authorHatzes, A. P.
titleStarspots and exoplanets.
journalAstron. Nachrichten
volume323, pages392–394

[124]
authorMeunier, N., authorDesort, M. &
authorLagrange, A. M.
titleUsing the Sun to estimate Earth-like planets
detection capabilities . II. Impact of plages.
journal volume512,
pagesA39 (year2010).

[125]
authorDumusque, X., authorBoisse, I. &
authorSantos, N. C.
titleSOAP 2.0: A Tool to Estimate the Photometric and
Radial Velocity Variations Induced by Stellar Spots and Plages.
journal volume796 (2),
pages132 (year2014).

[126]
authorQueloz, D. et al.
titleNo planet for HD 166435.
journal volume379,
pages279–287 (year2001).

[127]
authorBoisse, I. et al.
titleStellar activity of planetary host star HD 189
journal volume495 (3),
pages959–966 (year2009).

[128]
authorDumusque, X.
titleRadial velocity fitting challenge. I. Simulating the
data set including realistic stellar radial-velocity signals.
journal volume593,
pagesA5 (year2016).

[129]
authorSimola, U., authorDumusque, X. &
authorCisewski-Kehe, J.
titleMeasuring precise radial velocities and
cross-correlation function line-profile variations using a Skew Normal
journal volume622,
pagesA131 (year2019).

[130]
authorSimola, U. et al.
titleAccounting for stellar activity signals in
radial-velocity data by using change point detection techniques.
journal volume664,
pagesA127 (year2022).

[131]
authorBonfanti, A. et al.
titleTOI-1055 b: Neptunian planet characterised with
HARPS, TESS, and CHEOPS.
journalarXiv e-prints
pagesarXiv:2302.10664 (year2023).

[132]
authorBonfanti, A. & authorGillon, M.
titleMCMCI: A code to fully characterise an exoplanetary
journal volume635,
pagesA6 (year2020).

[133]
authorSchwarz, G.
titleEstimating the Dimension of a Model.
journalAnn. Stat.
volume6 (2), pages461–464
(year1978) .

[134]
authorGelman, A. & authorRubin, D. B.
titleInference from Iterative Simulation Using Multiple
journalStat. Sci.
volume7, pages457–472

[135]
authorRajpaul, V., authorAigrain, S.,
authorOsborne, M. A., authorReece, S. &
authorRoberts, S.
titleA Gaussian process framework for modelling stellar
activity signals in radial velocity data.
journal volume452 (3),
pages2269–2291 (year2015).

[136]
authorBarragán, O., authorAigrain, S.,
authorRajpaul, V. M. & authorZicher, N.
titlePYANETI - II. A multidimensional Gaussian process
approach to analysing spectroscopic time-series.
journal volume509 (1),
pages866–883 (year2022).

[137]
authorBarragán, O. et al.
titleThe young HD 73583 (TOI-560) planetary system: two
10-M$_{{\ensuremath{\oplus}}}$ mini-Neptunes transiting a 500-Myr-old,
bright, and active K dwarf.
journal volume514 (2),
pages1606–1627 (year2022).

[138]
authorZicher, N. et al.
titleOne year of AU Mic with HARPS - I. Measuring the
masses of the two transiting planets.
journal volume512 (2),
pages3060–3078 (year2022).

[139]
authorBarragán, O., authorGandolfi, D. &
authorAntoniciello, G.
titlePYANETI: a fast and powerful software suite for
multiplanet radial velocity and transit fitting.
journal volume482,
pages1017–1030 (year2019).

[140]
authorCale, B. L. et al.
titleDiving Beneath the Sea of Stellar Activity:
Chromatic Radial Velocities of the Young AU Mic Planetary System.
journal volume162 (6),
pages295 (year2021).

[141]
authorBlunt, S. et al.
titleOverfitting Affects the Reliability of Radial
Velocity Mass Estimates of the V1298 Tau Planets.
journalarXiv e-prints
pagesarXiv:2306.08145 (year2023).

[142]
authorDorn, C. et al.
titleCan we constrain the interior structure of rocky
exoplanets from mass and radius measurements?
journalAstron. Astrophys.
volume577, pagesA83 (year2015).

[143]
authorDorn, C. et al.
titleA generalized Bayesian inference method for
constraining the interiors of super Earths and sub-Neptunes.
journalAstron. Astrophys.
volume597, pagesA37 (year2017).

[144]
authorHaldemann, J., authorAlibert, Y.,
authorMordasini, C. & authorBenz, W.
titleAQUA: a collection of H$_{2}$O equations of state
for planetary models.
journalAstron. Astrophys.
volume643, pagesA105 (year2020).

[145]
authorHakim, K. et al.
titleA new ab initio equation of state of hcp-Fe and its
implication on the interior structure and mass-radius relations of rocky
journal volume313,
pages61–78 (year2018).

[146]
authorSotin, C., authorGrasset, O. &
authorMocquet, A.
titleMass radius curve for extrasolar Earth-like planets
and ocean planets.
journal volume191 (1),
pages337–351 (year2007).

[147]
authorLopez, E. D. & authorFortney, J. J.
titleUnderstanding the Mass-Radius Relation for
Sub-neptunes: Radius as a Proxy for Composition.
journalAstrophys. J.
volume792 (1), pages1

[148]
authorThiabaud, A. et al.
titleFrom stellar nebula to planets: The refractory
journalAstron. Astrophys.
volume562, pagesA27 (year2014).

[149]
authorMarboeuf, U., authorThiabaud, A.,
authorAlibert, Y., authorCabral, N. &
authorBenz, W.
titleFrom planetesimals to planets: volatile molecules.
journalAstron. Astrophys.
volume570, pagesA36 (year2014).

[150]
authorVenturini, J., authorGuilera, O. M.,
authorHaldemann, J., authorRonco, M. P. &
authorMordasini, C.
titleThe nature of the radius valley. Hints from
formation and evolution models.
journal volume643,
pagesL1 (year2020).

[151]
authorEmsenhuber, A. et al.
titleThe New Generation Planetary Population Synthesis
(NGPPS). II. Planetary population of solar-like stars and overview of
statistical results.
journal volume656,
pagesA70 (year2021).

[152]
authorIzidoro, A. et al.
titleThe Exoplanet Radius Valley from Gas-driven Planet
Migration and Breaking of Resonant Chains.
journal volume939 (2),
pagesL19 (year2022).

[153]
authorHu, R. et al.
titleUnveiling Shrouded Oceans on Temperate sub-Neptunes
via Transit Signatures of Solubility Equilibria versus Gas Thermochemistry.
journal volume921 (1),
pagesL8 (year2021).

[154]
authorTsai, S.-M. et al.
titleInferring Shallow Surfaces on Sub-Neptune Exoplanets
with JWST.
journal volume922 (2),
pagesL27 (year2021).

[155]
authorHarris, C. R. et al.
titleArray programming with NumPy.
journal volume585 (7825),
pages357–362 (year2020).

[156]
authorHunter, J. D.
titleMatplotlib: A 2D Graphics Environment.
journalComputing in Science and Engineering
volume9 (3), pages90–95

[157]
authorAstropy Collaboration et al.
titleThe Astropy Project: Sustaining and Growing a
Community-oriented Open-source Project and the Latest Major Release (v5.0) of
the Core Package.
journal volume935 (2),
pages167 (year2022).

[158]
authorForeman-Mackey, D., authorHogg, D. W.,
authorLang, D. & authorGoodman, J.
titleemcee: The MCMC Hammer.
journal volume125 (925),
pages306 (year2013).

[159]
authorMacDonald, M. G., authorShakespeare, C. J. &
authorRagozzine, D.
titleA Five-Planet Resonant Chain: Reevaluation of the
Kepler-80 System.
journal volume162 (3),
pages114 (year2021).

[160]
authorCannon, A. J. & authorPickering, E. C.
titleThe Henry Draper catalogue 0h, 1h, 2h, and 3h.
journalAnnals of Harvard College Observatory
volume91, pages1–290 (year1918).

[161]
authorGaia Collaboration et al.
titleGaia Early Data Release 3. Summary of the contents
and survey properties.
journal volume649,
pagesA1 (year2021).

[162]
authorYoss, K. M. & authorGriffin, R. F.
titleRadial Velocities and DDO, BV Photometry of Henry
Draper G5-M Stars Near the North Galactic Pole.
journalJournal of Astrophysics and Astronomy
volume18, pages161 (year1997).

[163]
authorSkrutskie, M. F. et al.
titleThe Two Micron All Sky Survey (2MASS).
journal volume131,
pages1163–1183 (year2006).

[164]
authorDelisle, J. B.
titleAnalytical model of multi-planetary resonant chains
and constraints on migration scenarios.
journal volume605,
pagesA96 (year2017).

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 ground-based 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 [139, 136]. We can share the code used in the data reduction or data analysis on request.

We acknowledge the use of public TESS data from pipelines at the TESS Science Office 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, offices, 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ón of the Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033 and the ERDF “A way of making Europe” through projects PID2019-107061GB-C61, PID2019-107061GB-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ències 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ón 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ánchez 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ómico Nacional on the Sierra de San Pedro Mártir (OAN-SPM), Baja California, México.
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-LP-105. 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 agreement 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ón y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comunicações (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 staff at the F. L. Whipple Observatory for their assistance in 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 Office (BELSPO) in the framework of the PRODEX Program, and by the University of Liège through an ARC grant for Concerted Research Actions financed by the Wallonia-Brussels Federation.
L.De. is an F.R.S.-FNRS Postdoctoral Researcher.
This work was supported by FCT - Fundação para a Ciência e a Tecnologia through national funds and by FEDER through COMPETE2020 - Programa Operacional Competitividade e Internacionalizacão by these grants: UID/FIS/04434/2019, UIDB/04434/2020, UIDP/04434/2020, PTDC/FIS-AST/32113/2017 & POCI-01-0145-FEDER- 032113, PTDC/FIS-AST/28953/2017 & POCI-01-0145-FEDER-028953, 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.