# Characterization of the HD 108236 system with CHEOPS and TESS Confirmation of a fifth transiting planet

S. Hoyer, A. Bonfanti, A. Leleu, L. Acuña, L. M. Serrano, M., Deleuil, A.

Bekkelien, C. Broeg, H. -G. Florén, D. Queloz, et al.

► To cite this version:

S. Hoyer, A. Bonfanti, A. Leleu, L. Acuña, L. M. Serrano, et al., Characterization of the HD 108236 system with CHEOPS and TESS Confirmation of a fifth transiting planet. Astronomy and Astrophysics - A&A, 2022, 668, 10.1051/0004-6361/202243720. [hal-03935714](https://hal.science/hal-03935714)

HAL Id: hal-03935714

https://hal.science/hal-03935714

Submitted on 19 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Characterization of the HD 108236 system with CHEOPS and TESS

## Confirmation of a fifth transiting planet*

S. Hoyer¹, A. Bonfanti², A. Leleu³,⁴, L. Acuña¹, L. M. Serrano⁵, M. Deleuil¹, A. Bekkelien³, C. Broeg⁴,⁶, H.-G. Florén⁷,⁸, D. Queloz⁹,¹⁰, T. G. Wilson¹¹, S. G. Sousa¹², M. J. Hooton¹⁰,⁴, V. Adibekyan¹², Y. Alibert⁴, R. Alonso¹³,¹⁴, G. Anglada¹⁵,¹⁶, J. Asquier¹⁷, T. Bárczy¹⁸, D. Barrado¹⁹, S. C. C. Barros¹²,²⁰, W. Baumjohann², M. Beck³, T. Beck⁴, W. Benz⁴,⁶, N. Billot³, F. Biondi²¹,²², X. Bonfils²³, A. Brandeker⁷, J. Cabrera²⁴, S. Charnoz²⁵, A. Collier Cameron¹¹, Sz. Csizmadia²⁴, M. B. Davies²⁶, L. Delrez²⁷,²⁸, O. D. S. Demangeon¹²,²⁰, B.-O. Demory⁶, D. Ehrenreich³, A. Erikson²⁴, A. Fortier⁴,⁶, L. Fossati², M. Fridlund²⁹,³⁰, D. Gandolfi⁵, M. Gillon²⁷, M. Güdel³¹, N. Hara³, K. Heng⁶,³², K. G. Isaak³³, J. M. Jenkins³⁴, L. L. Kiss³⁵,³⁶, J. Laskar³⁷, D. W. Latham³⁸, A. Lecavelier des Etangs³⁹, M. Lendl³, C. Lovis³, A. Luntzer⁴⁰, D. Magrin²², P. F. L. Maxted⁴¹, V. Nascimbeni²², G. Olofsson⁷, R. Ottensamer⁴⁰, I. Pagano⁴², E. Pallé¹³, C. M. Persson³⁰, G. Peter⁴³, D. Piazza⁴, G. Piotto²²,²⁴, D. Pollacco³², R. Ragazzoni²²,²⁴, N. Rando¹⁷, H. Rauer²⁴,⁴⁵,⁴⁶, I. Ribas¹⁵,¹⁶, G. R. Ricker⁴⁷, S. Salmon³, N. C. Santos¹²,²⁰, G. Scandariato⁴², S. Seager⁴⁷,⁴⁸,⁴⁹, D. Ségransan³, A. E. Simon⁴, A. M. S. Smith²⁴, M. Steller², Gy. M. Szabó⁵⁰,⁵¹, N. Thomas⁴, J. D. Twicken⁵²,³⁴, S. Udry³, V. Van Grootel²⁸, R. K. Vanderspek⁴⁷, N. A. Walton⁵³, K. Westerdorff⁴³, and J. N. Winn⁵⁴

*(Affiliations can be found after the references)*

Received 5 April 2022 / Accepted 4 September 2022

## ABSTRACT

**Context.** The HD 108236 system was first announced with the detection of four small planets based on TESS data. Shortly after, the transit of an additional planet with a period of 29.54 d was serendipitously detected by CHEOPS. In this way, HD 108236 ($V = 9.2$) became one of the brightest stars known to host five small transiting planets ($R_p < 3 R_⊕$).

**Aims.** We characterize the planetary system by using all the data available from CHEOPS and TESS space missions. We use the flexible pointing capabilities of CHEOPS to follow up the transits of all the planets in the system, including the fifth transiting body.

**Methods.** After updating the host star parameters by using the results from Gaia eDR3, we analyzed 16 and 43 transits observed by CHEOPS and TESS, respectively, to derive the planets' physical and orbital parameters. We carried out a timing analysis of the transits of each of the planets of HD 108236 to search for the presence of transit timing variations.

**Results.** We derived improved values for the radius and mass of the host star ($R_*$ = 0.876 ± 0.007 $R_⊕$ and $M_*$ = 0.867$^{+0.047}_{-0.046} M_⊕$). We confirm the presence of the fifth transiting planet *f* in a 29.54 d orbit. Thus, the HD 108236 system consists of five planets of $R_b$ = 1.587±0.028, $R_c$ = 2.122±0.025, $R_d$ = 2.629 ± 0.031, $R_e$ = 3.008 ± 0.032, and $R_f$ = 1.89 ± 0.04 [$R_⊕$]. We refine the transit ephemeris for each planet and find no significant transit timing variations for planets *c*, *d*, and *e*. For planets *b* and *f*, instead, we measure significant deviations on their transit times (up to 22 and 28 min, respectively) with a non-negligible dispersion of 9.6 and 12.6 min in their time residuals.

**Conclusions.** We confirm the presence of planet *f* and find no significant evidence for a potential transiting planet in a 10.9 d orbital period, as previously suggested. Further monitoring of the transits, particularly for planets *b* and *f*, would confirm the presence of the observed transit time variations. HD 108236 thus becomes a key multi-planetary system for the study of formation and evolution processes. The reported precise results on the planetary radii – together with a profuse RV monitoring – will allow for an accurate characterization of the internal structure of these planets.

**Key words.** planetary systems – planets and satellites: detection – planets and satellites: fundamental parameters – planets and satellites: individual: HD 108236 – planets and satellites: individual: TOI 1233

## 1. Introduction

Multi-planet systems are key to improving our understanding and studying the formation and evolution of planetary systems, with the ultimate goal of finding resemblances with the history of our own Solar System. The dynamics and architecture

of multi-planet systems is crucial for imposing strong constraints on migration processes (e.g., Mills et al. 2016; Delisle 2017; Leleu et al. 2021a). In addition, depending on the system, gravitational interaction between planets can produce detectable transit timing variations (TTVs), especially when planets are located close to mean motion resonances (Miralda-Escudé 2002; Holman & Murray 2005; Agol et al. 2005). In specific cases, TTVs can be used to estimate the masses of the planets without the requirement of complementary radial velocity observations (e.g., Barros et al. 2015; Gillon et al. 2017b; Freudenthal et al. 2018; Agol et al. 2021).

<footnote>* CHEOPS detrended light curves are only available in electronic format the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/ viz-bin/cat/J/A+A/668/A117</footnote>
<table><caption>**Table 1.** Exoplanetary systems orbiting a bright stellar host (V < 10.5), with more than three planets and at least one transiting planet.</caption><thead><tr><th>System name</th><th>V-mag</th><th>Total number of planets</th><th>Number of transiting planets</th><th>P<sub>inner</sub>-P<sub>outer</sub> [d]</th><th>Reference</th></tr></thead><tbody><tr><td>HD 219134</td><td>5.5</td><td>5(6)<sup>(a)</sup></td><td>2</td><td>3.1–46.8</td><td>Motalebi et al. (2015); Gillon et al. (2017a)</td></tr><tr><td>55 Cnc</td><td>5.95</td><td>5</td><td>1</td><td>0.73–4800</td><td>McArthur et al. (2004)</td></tr><tr><td>Kepler-444</td><td>8.9</td><td>5</td><td>5</td><td>3.6–9.7</td><td>Campante et al. (2015)</td></tr><tr><td>HIP 41378</td><td>8.9</td><td>5</td><td>5</td><td>15.6–324</td><td>Vanderburg et al. (2019)</td></tr><tr><td>HD 108236</td><td>9.2</td><td>5</td><td>5</td><td>3.79–29.5</td><td>Daylan et al. (2021); Bonfanti et al. (2021)</td></tr><tr><td>Kepler-37</td><td>9.8</td><td>4</td><td>4</td><td>13.3–39.8</td><td>Barclay et al. (2013)</td></tr><tr><td>V1298 Tau</td><td>10.1</td><td>4</td><td>4</td><td>8.24–60</td><td>David et al. (2019)</td></tr><tr><td>TOI-561</td><td>10.3</td><td>4(5)<sup>(b)</sup></td><td>4</td><td>0.45–77</td><td>Lacedelli et al. (2022)</td></tr></tbody></table>

**Notes.** The orbital periods range of the planets in each system is given in the *P*<sub>inner</sub>-*P*<sub>outer</sub> column. <sup>(a)</sup>The potential presence of a sixth planet is still controversial. <sup>(b)</sup>Hints of an additional planet in a long-period orbit were observed in radial velocities.

The amenability of the systems for in-depth characterizations, including comparative atmospheric studies through transmission spectroscopy, strongly depends on the brightness of the host star due to signal-to-noise considerations (see, e.g., the TSM metric in Kempton et al. 2018). Thanks to the long baseline and high signal-to-noise ratios (S/N) of their observations, to date, most of the transiting multi-planet systems with four or more planets are part of the yield of the *Kepler+K2* missions (Borucki et al. 2010; Howell et al. 2014), for which the magnitude distribution of the stellar hosts peaks towards V ~15.

In Table 1, we show the seven multi-planet systems known with more than three planets, among which there is at least one transiting a stellar host brighter than V = 10.5. Kepler-444 is the brightest detected star (V = 8.9) to host five transiting planets. It is an ultra-compact system formed by very small planets (radii < 0.8 R<sub>⊕</sub>) in orbital periods between 3.6 and 9.7 d. HIP 41378 (V = 8.9, Vanderburg et al. 2019); it hosts two Sub-Neptune planets in orbits of 15.6 and 31.7 d and three additional large planets in long orbital periods (~130–324 d). Therefore, HD 108236 is the second-brightest stellar host of a multiple transiting-planet system in a compact architecture (periods < 30 d). The HD 108236 system, also consisting of five transiting planets, is suitable for a complete and precise orbital and physical characterization thanks to the brightness of its host star.

The HD 108236 system was discovered by the TESS survey (TESS Object of Interest ID TOI 1233; Guerrero et al. 2021), and it was originally reported to consist of four transiting planets with orbital periods between 3.7 and 19.6 d and radii between 1.6 and 3.12 R<sub>⊕</sub> (Daylan et al. 2021). Its magnitude and location in the sky, makes it a perfect target for CHEOPS (Benz et al. 2021). As a follow-up mission, its goal is to provide ultra-high precision photometric observations of known exoplanet host stars for obtaining precise determinations of transiting exoplanet sizes. Bonfanti et al. (2021) used the first three CHEOPS light curves of HD 108236 system to refine the planetary transit ephemeris one year after the discovery observations. In addition, the analysis of these data revealed the presence of a fifth planet with 29.5 d orbital period. Here, we further characterize the HD 108236 system based on the analysis of the full CHEOPS observation dataset consisting of 15 light curves; together with all the available TESS data (i.e., with the addition of Sector 37). In this work, we confirm the presence of the fifth planet, *f*, and we search for the 10.9 d transit signal of a putative sixth planet, which was suggested in the discovery paper.

In Sect. 2, we present the updated stellar host properties based on the updated input from Gaia eDR3. In Sect. 3, we

describe the photometric observations used in this work, reporting 12 new CHEOPS light curves and all the available TESS data. The analysis of the data is described in Sect. 4. We describe the timing analysis of the transits in Sect. 5 and compare different approaches to estimate the planetary masses in Sect. 6. The discussion of the results and the conclusions of the work are presented in Sects. 7 and 8, respectively.

## 2. Host star properties

Bonfanti et al. (2021) already provided a thorough characterization of HD 108236. In this work, we assume the same spectroscopic parameters and elemental abundances, and we recall that they were derived from 13 high-resolution spectra (program 1102.C-0923, PI: Gandolfi) acquired by the High Accuracy Radial velocity Planet Searcher (HARPS, Mayor et al. 2003). In particular, *T*<sub>eff</sub>, log *g*, and [Fe/H] were inferred from the ARES+MOOG tools (Sousa 2014), while the elemental abundances were computed following the procedure detailed in Adibekyan et al. (2012, 2015). These values are reported in Table 2.

In taking advantage of the new stellar parallax from the third Gaia early data release (eDR3, Gaia Collaboration 2021), we were able to update the stellar radius always through the infrared flux method (IRFM Blackwell & Shallis 1977), which has already been summarized in Bonfanti et al. (2021). The slight refinement brought us to *R*<sub>⋆</sub> = 0.876 ± 0.007 R<sub>☉</sub>, which is fully compatible with what reported in Bonfanti et al. (2021). Using *T*<sub>eff</sub>, [Fe/H], and the new *R*<sub>⋆</sub> estimate as input parameters, we then updated also the isochronal mass *M*<sub>⋆</sub> and age *t*<sub>⋆</sub> estimates. A first pair (*M*<sub>⋆,1</sub>, *t*<sub>⋆,1</sub>) was computed employing the isochrone placement technique (Bonfanti et al. 2015, 2016) within PARSEC<sup>1</sup> v1.2S evolutionary models (Marigo et al. 2017), while a second pair (*M*<sub>⋆,2</sub>, *t*<sub>⋆,2</sub>) was outputted by the CLES code (Code Liègeois d'Évolution Stellaire, Scuflaire et al. 2008), following the fitting procedure described in Salmon et al. (2021). We finally combined the two respective estimates after carefully checking their mutual consistency through the χ<sup>2</sup>-based criterion outlined in Bonfanti et al. (2021), ending up with the following results: *M*<sub>⋆</sub> = 0.867<sup>+0.047</sup><sub>-0.046</sub> M<sub>☉</sub> and *t*<sub>⋆</sub> = 6.7<sup>+3.3</sup><sub>-3.4</sub> Gyr. All the relevant stellar parameters are listed in Table 2 and we note that all the updated parameters are consistent with what computed by Bonfanti et al. (2021) well within 1σ.

<footnote n="1">¹ Padova And TRieste Stellar Evolutionary Code: http://stev.oapd.inaf.it/cgi-bin/cmd</footnote>
<table><caption>**Table 2.** Stellar properties of HD 108236.</caption><thead><tr><th colspan="3">HD 108236</th></tr></thead><tbody><tr><td>Alternative<br>names</td><td colspan="2">TOI-1233<br>HIP 60689<br>TIC 260647166<br>Gaia DR2 6125644402384918784</td></tr><tr><th>Parameter</th><th>Value</th><th>Source</th></tr><tr><td>RA [J2000]</td><td>12:26:17.89</td><td>Gaia DR2</td></tr><tr><td>Dec [J2000]</td><td>-51:21:46.21</td><td>Gaia DR2</td></tr><tr><td>V [mag]</td><td>9.24</td><td>Simbad</td></tr><tr><td>G [mag]</td><td>9.0875</td><td>Simbad</td></tr><tr><td>J [mag]</td><td>8.046</td><td>Simbad</td></tr><tr><td>T<sub>eff</sub> [K]</td><td>5660 ± 61</td><td>Spectroscopy</td></tr><tr><td>log g [cgs]</td><td>4.49 ± 0.11</td><td>Spectroscopy</td></tr><tr><td>[Fe/H] [dex]</td><td>-0.28 ± 0.04</td><td>Spectroscopy</td></tr><tr><td>[Mg/H] [dex]</td><td>-0.27 ± 0.03</td><td>Spectroscopy</td></tr><tr><td>[Si/H] [dex]</td><td>-0.29 ± 0.02</td><td>Spectroscopy</td></tr><tr><th colspan="3">Updated values</th></tr><tr><td>d [pc]</td><td>64.66 ± 0.05</td><td>Gaia parallax <sup>(a)</sup></td></tr><tr><td>θ<sup>(b)</sup> [mas]</td><td>0.1261 ± 0.0010</td><td>IRFM</td></tr><tr><td>R<sub>⋆</sub> [R<sub>☉</sub>]</td><td>0.876 ± 0.007</td><td>IRFM</td></tr><tr><td>M<sub>⋆</sub> [M<sub>☉</sub>]</td><td>0.867<sup>+0.047</sup><sub>-0.046</sub></td><td>Isochrones</td></tr><tr><td>t<sub>⋆</sub> [Gyr]</td><td>6.7<sup>+3.3</sup><sub>-3.4</sub></td><td>Isochrones</td></tr><tr><td>L<sub>⋆</sub> [L<sub>☉</sub>]</td><td>0.707 ± 0.032</td><td>from R<sub>⋆</sub> and T<sub>eff</sub></td></tr><tr><td>ρ<sub>⋆</sub> [g cm<sup>-3</sup>]</td><td>1.79 ± 0.11</td><td>from R<sub>⋆</sub> and M<sub>⋆</sub></td></tr></tbody></table>

Notes. (a) Parallax offset from Lindegren et al. (2021) applied. (b) Angular diameter.

## 3. Observations

Here, we describe the photometric data used for the modeling of the system from both CHEOPS and TESS space missions.

### 3.1. CHEOPS

We obtained 15 CHEOPS observation runs (visits) of HD 108236, accumulating ~8.5 d on target, of which ~3.9 d were not covered due to Earth occultations or South Atlantic Anomaly (SAA) crossings by the satellite. With CHEOPS, we observed a total of 16 planetary transits of the system, as reported in the observation log in Table 3. The raw data of each visit were automatically processed by CHEOPS Data Reduction Pipeline (DRP, v13 Hoyer et al. 2020). In short, the DRP performs the instrumental calibration (bias, gain, linearization, and flat-fielding correction) and the environmental correction (cosmic ray hits, background, and smearing correction) before extracting the photometric signal of the target in four different aperture sizes. Three of these apertures have a fixed radius (RINF=22.5 pix, DEFAULT=25 pix, RSUP=30 pix) while a fourth, ROPT, is chosen automatically based on the contamination present in the Field of View (FoV). The contamination is computed by DRP using simulations of the FoV based on the Gaia catalog (Gaia Collaboration 2018). This information, together with other correction times series (e.g., background, smearing), are provided in the final DRP products along with the photometric time series, to be used in the analysis and detrending of the light curves. In this work, we use the light curves obtained with the DEFAULT aperture, as it produces the raw curves with the lowest dispersion when compared to the other photometric apertures.

### 3.2. TESS

Similarly to Daylan et al. (2021) and Bonfanti et al. (2021), our analysis includes the TESS data of HD 108236, consisting of the 2 min light curves obtained during the monitoring of Sec tors 10 and 11. Additionally, we incorporated the 2 min Sector 37 data, collected by TESS camera 2, between 2021-04-02 and 2021-04-28, during its second pass of the Southern Hemisphere. Here, we also use the Presearch Data Conditioning Simple Aperture Photometry (PDC-SAP, Stumpe et al. 2012, 2014; Smith et al. 2012) as provided by the Science Processing Operations Center (SPOC, Jenkins et al. 2016). In the TESS data, we identified and analyzed 18 transits of planet b, 12 transits of planet c, 6 transits of planet d, 4 transits of planet e, and 3 transits of planet f. In order to reduce processing times in the analysis, in this work we used only the portions of the TESS light curves with detected transits (leaving enough out-of-transit data for further detrending), ending up with 34 single light curves which we treated individually; among them, 4 light curves present consecutive transits of two planets and other 5 light curves show overlapped transits (transit events of two planets occurring simultaneously).

## 4. System modeling

### 4.1. Detrending of CHEOPS light curves

CHEOPS rotates around the line of sight while observing. Together with the extended shape of the CHEOPS PSF, this can produce different effects that perturb the aperture photometry (see e.g., Hoyer et al. 2020; Lendl et al. 2020; Bonfanti et al. 2021; Maxted et al. 2022). Some of these effects are the varying contaminating flux from nearby rotating background stars, moving smear trails of bright stars, and the reflections imprinted in the detector produced by different astrophysical sources close to CHEOPS pointing direction. All of these effects vary with the roll angle of the satellite and as a result, the flux of the raw target's light curve is correlated with the roll angle. Therefore, in most cases, a proper detrending against this parameter is fundamental to exploit the exquisite CHEOPS photometry. Following Bonfanti et al. (2021), we modeled the flux versus the roll angle pattern of each CHEOPS light curve through Gaussian processes (GPs, Rasmussen & Williams 2005). In particular, we adopted a Matérn-3/2 kernel and used the GP predictor implemented in the celerite package (Foreman-Mackey et al. 2017), which provides the GP model and its variance. Finally, we obtained the roll-angle detrended light curves, while also adjusting the flux error bars to account for the further variance introduced by the GP modeling. In order to reduce the number of fitted parameters, we used these roll angle detrended CHEOPS light curves in the following transits modeling. The CHEOPS raw and detrended light curves are reported in Appendix A (Figs. A.1–A.3).

### 4.2. Transit modeling

For the analysis of CHEOPS and TESS light curves, we used the juliet package (Espinoza et al. 2019), which implements the modeling of transit light curves from batman (Kreidberg 2015) in a Bayesian framework using the dynasty nested sampling tool (Speagle 2020). We performed the modeling of the full system (star+five planets) using all the CHEOPS (15 roll-angle detrended light curves, as described in Sect. 4.1) and TESS (34) light curves. In this configuration, we are able to retrieve, for example, an estimation of the orbital period (P) and reference time ($T_0$) for each planet. We used as reference the values
<table><caption>**Table 3.** Log of the CHEOPS observations of HD 108236.</caption><thead><tr><th>id<br>#</th><th>Planets</th><th>File Key</th><th>OBSID</th><th>Start date<br>[UTC]</th><th>Duration<br>[h]</th><th>Exp. time<br>[s]</th><th>Efficiency<br>[%]</th></tr></thead><tbody><tr><td>1<sup>(a)</sup></td><td>c, e, (g)</td><td>CH_PR300046_TG000101_V0200</td><td>1015572</td><td>2020-03-10T18:09:15.91</td><td>18.33</td><td>42.0</td><td>61.1</td></tr><tr><td>2</td><td>b</td><td>CH_PR100031_TG015701_V0200</td><td>1084970</td><td>2020-04-19T08:48:11.23</td><td>17.04</td><td>49.0</td><td>59.6</td></tr><tr><td>3<sup>(a)</sup></td><td>d</td><td>CH_PR100031_TG015401_V0200</td><td>1088211</td><td>2020-04-28T07:06:11.05</td><td>18.63</td><td>49.0</td><td>54.8</td></tr><tr><td>4</td><td>c</td><td>CH_PR100031_TG015601_V0200</td><td>1088443</td><td>2020-04-29T16:02:11.28</td><td>11.24</td><td>49.0</td><td>53.7</td></tr><tr><td>5<sup>(a)</sup></td><td>b, f</td><td>CH_PR100031_TG015702_V0200</td><td>1087957</td><td>2020-04-30T17:06:11.04</td><td>16.40</td><td>49.0</td><td>54.5</td></tr><tr><td>6</td><td>c</td><td>CH_PR100031_TG015602_V0200</td><td>1096738</td><td>2020-05-05T21:36:10.33</td><td>10.57</td><td>49.0</td><td>57.1</td></tr><tr><td>7</td><td>c</td><td>CH_PR100031_TG022401_V0200</td><td>1107339</td><td>2020-05-12T02:10:11.27</td><td>07.08</td><td>49.0</td><td>57.5</td></tr><tr><td>8</td><td>b</td><td>CH_PR100031_TG023001_V0200</td><td>1111865</td><td>2020-05-19T23:34:51.37</td><td>08.91</td><td>49.0</td><td>56.6</td></tr><tr><td>9</td><td>d, (g)</td><td>CH_PR100031_TG024301_V0200</td><td>1117140</td><td>2020-05-26T09:18:11.28</td><td>26.37</td><td>49.0</td><td>51.0</td></tr><tr><td>10</td><td>b</td><td>CH_PR100031_TG023002_V0200</td><td>1117572</td><td>2020-05-27T14:31:10.32</td><td>08.35</td><td>49.0</td><td>47.6</td></tr><tr><td>11</td><td>e</td><td>CH_PR100031_TG022701_V0200</td><td>1117376</td><td>2020-05-28T05:21:10.39</td><td>13.14</td><td>49.0</td><td>50.4</td></tr><tr><td>12</td><td>b</td><td>CH_PR100031_TG023003_V0200</td><td>1124008</td><td>2020-05-31T09:24:11.59</td><td>08.89</td><td>49.0</td><td>52.8</td></tr><tr><td>13</td><td>c</td><td>CH_PR100031_TG022901_V0200</td><td>1124302</td><td>2020-06-05T22:42:30.53</td><td>10.38</td><td>49.0</td><td>51.6</td></tr><tr><td>14</td><td>f</td><td>CH_PR110045_TG002601_V0200</td><td>1458282</td><td>2021-04-20T06:58:13.84</td><td>10.95</td><td>49.0</td><td>62.2</td></tr><tr><td>15</td><td>(g)</td><td>CH_PR110045_TG002701_V0200</td><td>1456848</td><td>2021-04-29T18:51:51.55</td><td>17.98</td><td>49.0</td><td>53.9</td></tr></tbody></table>

**Notes.** The second column shows the planetary transits observed during each CHEOPS visit. The last column shows the observation efficiency of each visit, considering the gaps produced by Earth occultations or SAA crossings. We report the unique identifiers of each CHEOPS visit: OBSID and File Key. <sup>(a)</sup>Observations presented in Bonfanti et al. (2021).

published in Bonfanti et al. (2021) to define our priors. In general, we used normal distributions centered in their results, but wide enough to ensure they are weakly informative. For the priors of the transit reference time ($T_0$) of each planet, when possible, we used TESS epochs which do not correspond to overlapped transits to mitigate potential time offsets. We also adopted the same constraints in the eccentricity ($e < 0.1$) resulting from dynamical stability simulations of the system; and the values for the quadratic limb darkening coefficients ($q_1, q_2$) for CHEOPS and TESS passbands (see Bonfanti et al. 2021 for details). In addition, juliet allows us to directly fit the stellar density ($\rho_*$), instead of fitting the scaled semi-major axis for each planet. Therefore, we used the value of stellar density derived from the updated stellar parameters (described in Sect. 2) as a truncated-normal prior in the modeling. We also fit for a photometric offset (mflux) and a “jitter” noise term ($\sigma_w$) for each light curve. As the contamination by external sources is small and already corrected in CHEOPS light curves by the roll angle detrending (described in Sect. 4.1), we fixed the dilution factor of each light curve, $mdilution = 1$. In addition, we fit for a linear and quadratic temporal term ($c_1(t)$ and $c_2(t^2)$) to take care of any residual trend remaining in the CHEOPS transits. Any additional free parameters of our joint fit were the orbital period (P), reference time ($T_0$), planet-to-star radius ratio ($R_p/R_*$), and the impact parameter ($b$) of each planet. The eccentricity ($e$) and the argument of periastron passage ($\omega$) of the orbit were estimated using the parametrization ($\sqrt{e} \sin(\omega), \sqrt{e} \cos(\omega)$). For these two parameters and the impact parameter, $b$, we used Uniform prior distributions. The value of the scaled semi-major axis of the orbit ($a/R_*$) for each planet was derived directly from the fitted $\rho_*$, P, and Kepler's third law. All the priors used in the modeling are listed in Table B.1. The fitted and derived planetary parameters from the posterior distributions of the system are shown in Table 4, while the phased light curves together with the best fit model are presented in Figs. 1 and 2 for CHEOPS and TESS data, respectively. The individual detrended light curves with the best fit model overplotted are shown in Figs. A.2 and A.3.

To confirm our results, we also performed an independent analysis of all the data with the Allesfitter package (Günther & Daylan 2021, 2019), following a similar procedure as described in Bonfanti et al. (2021). We imposed Normal priors on the mean stellar density $\rho_*$ and on the LD coefficients. For all the other system parameters, we adopted wide uniform priors. An upper limit of 0.1 on the orbital eccentricity was also imposed. The temporal detrending benefited of the GP treatment (Matérn-3/2 kernel) for the three TESS sectors, while we used splines for the 15 CHEOPS observations (already pre-detrended against roll-angle following Sect. 4.1) to avoid a dramatic increase of the number of fitted parameters. The resulting parameters are fully consistent, well within $\pm 1\sigma$ level in most of the cases, with the results obtained from the juliet analysis. We notice that orbital eccentricities listed in Table 4 are high when compared to what is expected for multiplanetary systems. We performed a series of tests to explore the possible cause, and finally we attribute the results to the parameterizations used by juliet. In any case, a more detailed RV analysis will provide a final estimate of the eccentricities, something that the photometry presented in this work cannot solve on its own. We describe the performed tests in Appendix C. We also confirm that the rest of the other modeled parameters are robust, no matter the approach used to fit $e$ and $\omega$.

### 4.3. Confirmation of planet *f*

As described in Bonfanti et al. (2021), an additional transit-like feature was serendipitously detected in the CHEOPS light curve when monitoring a transit of planet *b* (see top panel in Fig. 3). They also found signals in Sectors 10 and 11 of the TESS data that are consistent with a planetary candidate *f* with an orbital period of 29.54 d. We detected transits of planet *f* in the TESS Sectors 10 and 11, as well as in Sector 37, which are consistent with the 29.54 d period. The fact that two consecutive transits were identified in the first TESS Sector (as reported in Table E.2) and no other transit consistent with this planet was
<table><caption>**Table 4.** Fitted and derived parameters of the HD 108236 system from CHEOPS and TESS transits.</caption><thead><tr><th colspan="6">HD 108236</th></tr><tr><th>Quadratic L. D. Coeff.</th><th>q<sub>1</sub></th><th>q<sub>2</sub></th><th></th><th>q<sub>1</sub></th><th>q<sub>2</sub></th></tr></thead><tbody><tr><td>CHEOPS</td><td>0.446<sup>+0.016</sup><sub>-0.018</sub></td><td>0.341<sup>+0.017</sup><sub>-0.014</sub></td><td>TESS</td><td>0.277<sup>+0.018</sup><sub>-0.023</sub></td><td>0.230<sup>+0.016</sup><sub>-0.021</sub></td></tr><tr><th>Fitted parameter</th><th>Planet b</th><th>Planet c</th><th>Planet d</th><th>Planet e</th><th>Planet f</th></tr><tr><td>Orbital period, P<sup>(a)</sup> [d]</td><td>3.7958785<sup>+0.0000076</sup><sub>-0.0000070</sub></td><td>6.2036717<sup>+0.0000092</sup><sub>-0.000010</sub></td><td>14.175818<sup>+0.000024</sup><sub>-0.000023</sub></td><td>19.5901277<sup>+0.000030</sup><sub>-0.000030</sub></td><td>29.53935<sup>+0.00012</sup><sub>-0.00014</sub></td></tr><tr><td>T<sub>0</sub><sup>(a)</sup> [BJD-2450000]</td><td>8572.11017<sup>+0.00062</sup><sub>-0.00068</sub></td><td>8572.39445<sup>-0.00067</sup><sub>-0.00065</sub></td><td>8571.33678<sup>+0.00061</sup><sub>-0.00059</sub></td><td>8586.56741<sup>+0.00056</sup><sub>-0.00055</sub></td><td>8616.0403<sup>+0.0014</sup><sub>-0.0014</sub></td></tr><tr><td>R<sub>p</sub>/R<sub>★</sub></td><td>0.01661<sup>+0.00025</sup><sub>-0.00026</sub></td><td>0.02221<sup>+0.00020</sup><sub>-0.00018</sub></td><td>0.02751<sup>+0.00025</sup><sub>-0.00024</sub></td><td>0.03148<sup>+0.00023</sup><sub>-0.00023</sub></td><td>0.01979 ± 0.00043</td></tr><tr><td>√e sin(ω)</td><td>0.221<sup>+0.043</sup><sub>-0.055</sub></td><td>0.154<sup>+0.043</sup><sub>-0.068</sub></td><td>-0.159<sup>+0.051</sup><sub>-0.067</sub></td><td>0.257<sup>+0.025</sup><sub>+0.035</sub></td><td>-0.028<sup>+0.072</sup><sub>-0.076</sub></td></tr><tr><td>√e cos(ω)</td><td>0.120<sup>+0.066</sup><sub>-0.065</sub></td><td>-0.256<sup>+0.047</sup><sub>-0.028</sub></td><td>-0.217<sup>+0.084</sup><sub>-0.051</sub></td><td>-0.048<sup>+0.062</sup><sub>-0.070</sub></td><td>0.129<sup>+0.072</sup><sub>-0.076</sub></td></tr><tr><td>Impact parameter, b [R<sub>★</sub>]</td><td>0.322<sup>+0.060</sup><sub>-0.057</sub></td><td>0.308<sup>+0.055</sup><sub>-0.046</sub></td><td>0.519<sup>+0.029</sup><sub>-0.026</sub></td><td>0.392<sup>+0.032</sup><sub>-0.029</sub></td><td>0.701<sup>+0.020</sup><sub>-0.021</sub></td></tr><tr><th>Derived parameter</th><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Rel. semi-major axis, a/R<sub>★</sub></td><td>11.080<sup>+0.046</sup><sub>-0.043</sub></td><td>15.373<sup>+0.064</sup><sub>-0.059</sub></td><td>26.67<sup>+0.11</sup><sub>-0.10</sub></td><td>33.10<sup>+0.14</sup><sub>-0.13</sub></td><td>43.51<sup>+0.18</sup><sub>-0.17</sub></td></tr><tr><td>Eccentricity, e</td><td>0.067<sup>+0.025</sup><sub>-0.026</sub></td><td>0.091<sup>+0.021</sup><sub>-0.033</sub></td><td>0.075<sup>+0.021</sup><sub>-0.022</sub></td><td>0.073<sup>+0.014</sup><sub>-0.017</sub></td><td>0.024<sup>+0.024</sup><sub>-0.016</sub></td></tr><tr><td>Arg. of periastron, ω [°]</td><td>61<sup>+15</sup><sub>-14</sub></td><td>149<sup>+11</sup><sub>-7</sub></td><td>144<sup>+12</sup><sub>-22</sub></td><td>101<sup>+15</sup><sub>-14</sub></td><td>24<sup>+34</sup><sub>-17</sub></td></tr><tr><td>Orbital inclination, i [°]</td><td>88.23 ± 0.32</td><td>88.8 ± 0.2</td><td>88.8 ± 0.1</td><td>89.3 ± 0.1</td><td>89.07 ± 0.03</td></tr><tr><td>Planet radius [R<sub>⊕</sub>]</td><td>1.587 ± 0.028</td><td>2.122 ± 0.025</td><td>2.629 ± 0.031</td><td>3.008 ± 0.032</td><td>1.89 ± 0.04</td></tr><tr><td>Semi-major axis, a [AU]</td><td>0.0451 ± 0.0004</td><td>0.0626 ± 0.0006</td><td>0.1087 ± 0.0010</td><td>0.1348 ± 0.0012</td><td>0.1773 ± 0.0016</td></tr><tr><td>Equilibrium temp. T<sub>eq</sub><sup>(b)</sup> [K]</td><td>1202 ± 13</td><td>1021 ± 11</td><td>775 ± 8</td><td>696 ± 8</td><td>607 ± 7</td></tr><tr><td>Transit duration, T<sub>I-IV</sub> [h]</td><td>2.45 ± 0.06</td><td>2.93 ± 0.06</td><td>3.52 ± 0.08</td><td>4.16 ± 0.07</td><td>3.8 ± 0.1</td></tr></tbody></table>

**Notes.** The fitted parameters were obtained via the simultaneous modeling of CHEOPS and TESS light curves with juliet package (Sect. 4). For the estimation of the derived parameters, we used the astronomical constants as defined in the python library `astropy constants`. (a) The final and adopted values for P and $T_0$, which are based on the timing analysis of the transits are reported in Table 5. (b) Mean temperature calculated using $A_B = 0$.

detected in the TESS light curves in that time excludes the viability of the 29.54 d period aliases. The TESS phased light curve is presented at the bottom-right panel of Fig. 2. To confirm the presence of this planet, we scheduled a dedicated CHEOPS visit (observation 14 in Table 3) to observe a single transit of this planet. We present in Fig. 3 the light curve obtained from this observation together with the best-fit model resulting from the analysis of CHEOPS and TESS data. In Fig. 3, we also show the light curve reported in Bonfanti et al. (2021) with a transit of both b and f planets. The phased version of the combined two CHEOPS light curves is shown in the bottom-left panel of Fig. 2.

## 4.4. Non-detection of a putative planet g

Daylan et al. (2021) reported a hint of a ~10.9 d periodic transit signal in TESS data with amplitude of 230 ppm (without excluding a possible instrumental origin). Based on the ephemeris reported in that work ($T_0 = 2458570.6781$ BJD and $P = 10.9113$ d), we can see that this planet should have transited during the CHEOPS observations 1 and 9 (see the light curves in upper and middle panels of Fig. 4, respectively), but no clear evidence of this has been detected. Observation 1 has already been published by Bonfanti et al. (2021), who noticed that the transit of this putative planet would have occurred at $0 = 2458919.84$ BJD, being partially overlapped with the transit of planet c. They performed a comparative analysis of a five-planet versus six-planet scenario, finding that there was no strong evidence to reject the five-planet scenario based on the Bayes factors of the two models. According to its ephemeris, the hypothetical planet g would transit at $T_0 = 2458996.22$ BJD during observation 9 (Fig. 4, middle panel). In the time of that CHEOPS

visit, we actually noticed a short flux decrease, but it was shifted 0.13 d from the linear ephemeris at ~2458996.35 BJD. The amplitude of this flux decrease is ~125 ppm, which is around half of what was reported by Daylan et al. (2021) in TESS data. To further probe this periodic signal, we scheduled a ~18 h CHEOPS observation, taking special care to avoid any overlap with the transits of the other five planets in the system (observation 15 in Table 3). In Fig. 4 (bottom panel), we show the detrended light curve of this observation where we found no significant evidence of any transit-like signal. We observed a shallow flux decrease centered at ~2459334.53 BJD (linear ephemeris would predict $T_0 = 2459334.47$ BJD), but also consistent with other features present in the light curve, which likely are residuals of the detrending of the well-known systematics of CHEOPS photometry as a function of its rotation angle. Despite this, we performed additional statistical tests to confirm that the detection of this transit-like feature is not significant (see Appendix D for details). We conclude that with the data we have on hand, we cannot confirm the presence of this additional planet g in the system. In fact, a planet with a 10.9 d period would have a transit duration between 3.7–4 h (assuming a circular orbit and a zero impact parameter). The rms of the light curve of the CHEOPS observation 15 (bottom panel in Fig. 4) at these timescales corresponds to 17–20 ppm, which translates into an upper limit in the size of an undetected planet of 0.42 $R_⊕$ (with $e = b = 0$) at $S/N = 1$.

Taking into account that a transit of this putative planet was not observed in any of the CHEOPS light curves, we can assess the period range around $P = 10.9113$ d. Using only CHEOPS data we can rule out a period in this range assuming that the errors in the ephemeris are given by the orbital period only, that is, fixing
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Phased CHEOPS (**left**) and TESS (**right**) light curves (blue dots) and its 10-min bins (yellow points). The best transit models (solid black curves) for planets b, c and d (from top to bottom) obtained in the global analysis of the data are shown; together with the respective residuals (thin bottom panels). Additional transit events of other planets were removed before plotting the final phased curve.</figcaption></figure>

$T_0$ to the value reported by Daylan et al. (2021). Thus, we can discard periods in the range $P = 10.9113_{-0.018}^{+0.052}$ d using CHEOPS data only. We note that the accumulation of these uncertainties in the orbital period will translate into offsets of more than 20 h at CHEOPS epochs. Moreover, no transits were detected at the expected times in the TESS Sector 37. Additionally, we used the transit least-squares tool (TLS, Hippke & Heller 2019) to search for other transits in this TESS light curve. After masking the transits of the known planets, no significant signals were found with periods between 0.52 and 24.5 d.

**5. Timing analysis**

To perform a detailed timing analysis for each planet, we used the results of the combined analysis in order to determine the transit times for each individual transit event. Therefore, in this case we do not fit for a global $P$ and $T_0$ for each planet but for the central time, $T_c$, of each transit event. Thus, we used as priors the results of the combined fit for all the parameters except for the single transit midtime. With these transit times and the ephemeris (described in Sect. 4 and displayed in Table 4)
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Phased CHEOPS (**left**) and TESS (**right**) light curves (blue dots) and its 10-min bins (yellow points). The best transit models (solid black curves) for planet *e* (**top**) and *f* (**bottom**) obtained in the global analysis of the data. Same details as in Fig. 1.</figcaption></figure>

from the global fit, we generated the observed-minus-calculated (O–C) diagrams of the five planets in the system. A first analysis showed us that the times of the transits derived from overlapping transit events in TESS light curves correspond to large outliers in the O–C diagrams. Thus, for such light curves, we re-fit the two respective planets simultaneously. In this way, we were able to considerably reduce the temporal offsets and the uncertainties of the central times of these transits. The resulting individual transit times are reported in Tables E.1 and E.2.

We checked if we were able to fit a linear function to the timing residuals of each planet, meaning the ephemeris equation can potentially be updated. We evaluate the significance of the *P* or *T*₀ offsets by comparing them with their uncertainties from the global fit. No correction was needed for the orbital periods. We only found significant temporal offsets of 114.5 s (2 σ), -58 s (1.1 σ) and 717 s (5.9 σ) for the *T*₀ of planets *b*, *d* and *f*, respectively. Therefore, we updated the *T*₀ values of the ephemeris equation for these three planets (Table 5). For comparison, in Table 5 we also show the mean error of the transit times, σ̄*τ*<sub>c</sub>, for CHEOPS and TESS transits. The improvement on our estimations of the transit times from CHEOPS light curves in comparison to those from TESS transits is consistent with the results of previous works (e.g., Borsato et al. 2021; Bonfanti et al. 2021). The final O–C diagram of each planet is shown in Fig. 5. There, the overlapped transits are marked with the vertical

dashed lines, and the ±1–2σ uncertainties of the final ephemeris are represented by the gray regions in the diagrams. We noticed non-negligible variations on the transit times of planets *b* and *f*, reflected for example on the rms and mean amplitude values of their O–C residuals. We also noticed hints of an anticorrelation in the transit times of planets *e* and *f*. Thus, in order to investigate further on the results in the transit times, we performed a detailed transit time variations (TTVs) analysis of the system, described in Sect. 6.3.

## 6. Planetary masses

In the absence of radial velocity (RV) measurements, model-dependent estimates of the planetary masses have been previously delivered by Daylan et al. (2021) and Bonfanti et al. (2021). Later, Teske et al. (2021) presented 33 RVs measurements obtained with the Planet Finder Spectrograph (PFS) on the Magellan II telescope at Las Campanas Observatory (Chile). They performed the first Keplerian mass estimation of the HD 108236 system. Among their results for this particular system, it is worth mentioning that they were only able to set an upper limit for the mass of planet *d* and they predicted an eccentricity of ~0.39 for planet *e*. Here, we take advantage of the refined parameters of the system to estimate the masses of the five planets by, first, performing a new fit of the RV data and; in
<table><caption>**Table 5.** Results of the timing analysis of the transits.</caption><thead><tr><th>Planet</th><th>O-C rms<br>[min]</th><th>O-C mean<br>amplitude [min]</th><th>σ<sup>CHEOPS</sup><sub>T<sub>C</sub></sub><br>[min]</th><th>σ<sup>TESS</sup><sub>T<sub>C</sub></sub><br>[min]</th><th>Updated T<sub>0</sub><br>[BJD<sub>TDB</sub> - 2 450 000]</th></tr></thead><tbody><tr><td>b</td><td>9.6</td><td>8.0</td><td>3.6</td><td>6.6</td><td>8572.11149 ± 0.00065</td></tr><tr><td>c</td><td>4.6</td><td>3.4</td><td>2.0</td><td>3.3</td><td>—</td></tr><tr><td>d</td><td>2.4</td><td>2.0</td><td>4.0</td><td>3.1</td><td>8571.33610 ± 0.00060</td></tr><tr><td>e</td><td>2.4</td><td>2.4</td><td>1.0</td><td>2.1</td><td>—</td></tr><tr><td>f</td><td>12.6</td><td>17.8</td><td>4.5</td><td>8.8</td><td>8616.0486 ± 0.0014</td></tr></tbody></table>

**Notes.** Results of the analysis described in Sect. 5. The rms and the mean value of the time residuals of each planet's transit are shown (2nd and 3rd columns). For comparison, we show the mean value of the error bars for the transit times of CHEOPS and TESS transits (4th and 5th columns). We also report the updated values of $T_0$ for planet b, d and f based on the weighted fit of the time residuals (6th column).

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** CHEOPS light curve of HD 108236 during a transit of planet *f* (~2 458 970.5 BJD) and *b* (~2 458 970.7 BJD) (top panel), as well as during a dedicated observation of a planet *f* transit (bottom panel). Detrended photometric points and its 10 min bins are represented by light and dark blue dots, respectively. The best model for each planet is represented by the yellow curve.</figcaption></figure>

addition, by using the M-R relations of Chen & Kipping (2017; CK17 hereafter). Finally, we present our TTV analysis based on the transit timing information of Sect. 5, aimed at exploring whether any constraints can be placed on the planetary masses.

## 6.1. Radial velocity estimations

For the fitting of the RV observations, we also used *juliet*, which implements the RADVEL code (Fulton et al. 2018) linked to the *dynesty* sampler. We set the priors of all the orbital parameters ($P_i$, $T_{0,i}$, $b_i$, $e_i$ and $\omega_i$, i:1..5) as normal distributions defined by the results of the photometric fit. For the five RV

<caption>**Fig. 4.** CHEOPS light curves of HD 108236 during the expected transit time of a potential planet *g* (marked with the yellow triangle). These light curves correspond to observations 1, 9, and 15 in Table 3 (top-to-bottom). Light and dark blue dots represent the detrended and 10 min bins photometric points. No significant statistical evidence of a transit signal with ~10.9 d period is found in these photometric time series.</caption>
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Observed minus calculated diagrams of the transit times of the five planets of the HD 108236 system, obtained by using the ephemeris of the global fit or the updated $T_0$ values reported in Table 5. Transit times were obtained from the individual analysis of the transits (see Sect. 4). Blue and orange symbols correspond, respectively, to TESS and CHEOPS observations. The gray regions represent the ±1–2σ uncertainties of the ephemeris equations. Vertical dashed lines mark overlapped transits in the light curves.</figcaption></figure>

semi-amplitudes, $K_i$, we set a wide uniform prior between 0 and 20 m s⁻¹. We also fit for a jitter and a systemic RV term (σ_PFS2 and µ_PFS2, respectively). The fitted *K*'s and the derived masses for each planet are presented in Table 6, where we also show the masses estimations from the literature. The RV data and the best fitted models are presented in Figs. F.1 and F.2. Due to the paucity and poor sampling of the observations, the RV masses should be taken with caution, and therefore we treated them as only indicative values. It is known that the accurate measurement of masses of small planets in multiple systems requires considerable large amount of data points (e.g., He et al. 2021), which is not the case of the HD 108236 system to date. In fact, the posterior distribution of the RV semi-amplitude for planet *d*, $K_d$, peaks around $K_d \sim 0.41$ m s⁻¹ (see Fig. F.3). As noted in Teske et al. (2021), if we allow a fitting for $K < 0$ values, we then would

obtain a negative-K ($K_d \sim -0.8$ m s⁻¹). This was the reason why Teske et al. (2021) adopted only an upper limit for $K_d$ defined by the root-mean-square (rms) of the RV residuals after removing the best Keplerian model of the other four planets. We confirmed that the results of the other planets were almost unaffected by removing planet *d* from the fit. Therefore, we emphasize that our reported RV result for planet *d* is not reliable.

## 6.2. Statistical estimations

To compare the RV masses with the values derived with a statistical approach, we use the CK17 M-R relations, implemented within the Python probabilistic tool forecaster<fnref n="2" />. We used as

<footnote n="2">² https://github.com/chenjj2/forecaster</footnote>
input the derived planetary radius (and its uncertainties) of the five planets obtained from our photometric fit. The masses predicted by CK17 and their uncertainties are presented in Table 6. We also checked that these values are consistent with the predictions from more recent M-R relations from Otegi et al. (2020), although the uncertainties of the later are 2.3–2.7 times smaller. As these M-R relations are intended to provide a broad estimation of the planetary masses, the values we derived with them should be taken with caution.

We note that for the outer two planets, the RV derived masses are unrealistically high when compared to the models outcomes, particularly for planet *f*. This can likely be explained by the fact that the phase of planet *f* is not well sampled by the RV measurements, and that the fit of planet *e* seems to show a phase offset, likely resulting in the poor mass determination. Since these two planets are located close to a mutual 3:2 mean motion resonance (MMR), it is possible to explore whether the dynamical analysis can provide additional constraints on the mass estimates, as described in Sect. 6.3.

## 6.3. TTV analysis

The period ratio between consecutive pairs of planets ($P_{i+1}/P_i$) in the HD 108236 system is ($P_c/P_b, P_d/P_c, P_e/P_d, P_f/P_e$) = (1.634, 2.285, 1.382, 1.508). Thus, the outer planets, *f* and *e*, are close to the exact 3/2 commensurability ($P_f/P_e \approx 1.508$). The proximity to a two-body mean motion resonances (MMRs), $P_{\text{out}}/P_{\text{in}} \sim (p+q)/p$, with *p* and *q* integers, can generate TTVs that can be observed and modeled to constrain the orbit and masses of the involved planets (e.g., Agol et al. 2005; Lithwick et al. 2012; Nesvorný & Vokrouhlický 2016). We do not currently have enough TTVs measurement to estimate masses. Here, we report the mass posterior of this TTV analysis with the sole purpose to compare them with the CK17 masses from Table 6 (that we use as prior) in order to check whether these masses are consistent with the observed TTVs.

We fit the timings reported in Tables E.1 and E.2 using the code presented in Leleu et al. (2021b): the transit timings are estimated using the TTVfast algorithm (Deck et al. 2014), and the *samsam*<fnref n="3" /> MCMC algorithm (see Delisle et al. 2018) is used to sample the posterior. The mean longitudes, periods, arguments of periastron and eccentricities of the planets have flat priors. The masses and eccentricities posteriors of the TTV analysis are shown in Table 6, and 1000 randomly chosen samples of the fitted model are shown in Fig. 6.

The prior-posterior comparison is made in Fig. 7. For planets *b*, *d*, *e*, and *f*, we see that the mass posterior points toward a lower value than the prior, although the difference lies close to (or within) the $1\sigma$ interval. Planet *c* is the only one for which TTVs mass is larger than the prior, which appears to be due to the tentative chopping in the TTV signal of planet *b* (see Fig. 6).

For the outer two, near-resonant planets *e* and *f*, there are hints of anti-correlated TTVs, which are also apparently not fully consistent with the mass priors. For the derived set of masses, significant TTVs are nonetheless expected to come in the next years for these two planets, as can be seen on the TTV projections in Fig. 6 at dates 9600 [BJD-2450 000] onward. Overall, the TTVs observed for HD 108236 seem to be somewhat at odds with the expected masses and eccentricities, although the current phase coverage is not sufficient to draw strong conclusions and warrant additional observations.

<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Posterior of the transit timings fit presented in Tables E.1 and E.2.</figcaption></figure>

In addition, in order to better constrain the architecture of the system, we used the TTVs posterior to check if the planets *e* and *f* are inside the 3:2 MMR using the Hamiltonian formulation of the second fundamental model for Resonance of Henrard & Lemaitre (1983). We applied the “one degree of freedom” model of first-order resonances presented in Deck et al. (2013), and computed the value of the Hamiltonian parameter $\Gamma'$ over the posterior of the TTV analysis (Eq. (36) in the

<footnote n="3">³ https://gitlab.unige.ch/Jean-Baptiste.Delisle/samsam</footnote>
<table><caption>**Table 6.** Mass estimations of the planets of HD 108236.</caption><thead><tr><th rowspan="3">Planet</th><th colspan="4">This work</th><th colspan="3">Literature</th></tr><tr><th colspan="2">CK17</th><th colspan="2">RVs</th><th rowspan="2">T21 <sup>(a)</sup> (RVs)</th><th rowspan="2">B21 <sup>(b)</sup></th><th rowspan="2">D21 <sup>(c)</sup> (CK17)</th></tr><tr><th>Mass [M<sub>⊕</sub>]</th><th>K [m s<sup>-1</sup>]</th><th>Mass [M<sub>⊕</sub>]</th><th>Density [gr cm<sup>-3</sup>]</th><th>Mass [M<sub>⊕</sub>]</th></tr></thead><tbody><tr><td>b</td><td>3.8<sup>+2.6</sup><sub>-1.3</sub></td><td>1.63<sup>+0.74</sup><sub>-0.72</sub></td><td>3.6<sup>+1.7</sup><sub>-1.6</sub></td><td>5.0 ± 2.0</td><td>4.33 ± 1.54</td><td>4.23 ± 0.40</td><td>5 ± 2</td></tr><tr><td>c</td><td>5.5<sup>+4.0</sup><sub>-2.2</sub></td><td>1.94<sup>+0.79</sup><sub>-0.70</sub></td><td>5.1<sup>+2.1</sup><sub>-1.8</sub></td><td>3.0 ± 1.0</td><td>3.97 ± 1.77</td><td>8.90 ± 0.66</td><td>7 ± 2</td></tr><tr><td>d</td><td>7.3<sup>+5.6</sup><sub>-3.0</sub></td><td>0.46<sup>+0.54</sup><sub>-0.31</sub></td><td>1.6<sup>+1.8</sup><sub>-1.1</sub></td><td>0.5 ± 0.5</td><td>&lt;7.75</td><td>7.75 ± 0.80</td><td>10 ± 2</td></tr><tr><td>e</td><td>9.9<sup>+7.4</sup><sub>-4.3</sub></td><td>3.55<sup>+0.78</sup><sub>-0.76</sub></td><td>13.6<sup>+3.0</sup><sub>-2.9</sub></td><td>3.0 ± 1.0</td><td>19.10 ± 3.91</td><td>8.2<sup>+3.8</sup><sub>-1.2</sub></td><td>13 ± 2</td></tr><tr><td>f</td><td>5.1<sup>+3.8</sup><sub>-2.0</sub></td><td>4.2<sup>+1.8</sup><sub>-1.6</sub></td><td>18.4<sup>+7.9</sup><sub>-7.1</sub></td><td>15.0 ± 6.0</td><td>10.85 ± 5.55</td><td>3.95 ± 0.4</td><td>—</td></tr></tbody></table>

**Notes.** The masses from the radial velocities (RVs) and from the probabilistic method based on the M-R relations of Chen & Kipping (2017; CK17) (Sect. 7.1) are shown. The fitted RV semi-amplitudes, *K*, and the derived planet densities are also shown. We recall in the last 3 columns the values from the literature: <sup>(a)</sup>Teske et al. (2021), <sup>(b)</sup>Bonfanti et al. (2021), and <sup>(c)</sup>Daylan et al. (2021).

<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Mass posterior vs. mass prior of the TTV analysis. The error bars show the $1\sigma$ interval.</figcaption></figure>

aforementioned paper). The resonant state appears for $\Gamma' \ge 1.5$ (Deck et al. 2013), and for this pair of planets (*e* and *f*) we calculated $\Gamma' = -2.76 \pm 0.38$, implying that the planets are outside the 3:2 MMR.

# 7. Discussion

## 7.1. Planetary parameters

As our final planetary parameters, we adopted those resulting from the combined analysis of CHEOPS and TESS data described in Sect. 4 and reported in Table 4. All the parameters are consistent with those presented in Bonfanti et al. (2021) within $1-2\sigma$. Notably, the uncertainties we report, which represent the 68% intervals of the resulting posterior distributions of each parameter, are reduced between 30–80% with respect to previously published uncertainties. A diagram of the orbital configuration, together with a comparative of the sizes of the planets, is presented in Fig. 8.

The uncertainties on the transit ephemeris also improved when compared to Bonfanti et al. (2021) thanks to the extended baseline of the monitored transits. Moreover, we also performed a transit timing analysis (see Sect. 5), which provided us with the updated final $T_0$ reported in Table 4 for planets *b*, *d*, and *f*. For the other planets, we still referred to the values in Table 4.

<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Diagram of the HD 108236 system. We show the scaled semi-major axis, $a/R_*$ (bottom axis) and the respective orbital period of the planets (in units of planet *f*'s period, top axis). We mark the relevant main mean motion resonances with vertical lines. The relative orbital inclination ($90^\circ - i$) of each planet is shown in the left axis. The symbol size scales with the planetary radius in units of Earth radius; and the color with the equilibrium temperature of the planet (with $A_B = 0$).</figcaption></figure>

Despite the fact that none of the three approaches we used to estimate the planetary masses are conclusive, we noticed that (as shown in Table 6) the masses of the three most internal planets are consistent among the different estimations (except for RV mass of planet *d* as discussed in Sect. 6.1). For the two outer planets, the RV masses seem unrealistically high when compared to the other methods, although they are poorly constrained by the data. Below, we use the RV masses to probe the planetary internal structure.

## 7.2. Rocky or gaseous planets

Taking advantage of the precision obtained in the planets' sizes and our RVs estimations of their masses, we placed the five planets of the system in an M-R diagram (Fig. 9). For planet *d*, we used $7.75 M_\oplus$ value as an upper limit based on the estimations from Teske et al. (2021) and Bonfanti et al. (2021). We also plotted the M-R relations from Otegi et al. (2020), red and blue thick curves, representing the “rocky” and “gaseous” planet populations, respectively. Despite the results of our RV fit are
<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Mass-radius diagram for planets with radius and masses below $3.5 R_⊕$ and $30 M_⊕$, respectively. The radius of the HD 108236's planets obtained from the transits modeling are represented by the gray horizontal regions, while the mass values suggested by the RV fit are represented by the black symbols. The upper limit in the mass of planet *d* (based on estimations from Teske et al. 2021 and Bonfanti et al. 2021) is marked with the red arrow. We compared the HD 108236 planets with the current planet population (with estimates of mass and radius better than 30%) represented by the gray symbols and the MR relations from Otegi et al. (2020) (blue and red thick lines). The planet interiors models from Brugger et al. (2017) and Acuña et al. (2021) are represented by the thin curves.</figcaption></figure>

not conclusive, and therefore, the calculated masses and their uncertainties should be treated as indicative values only, we are still able to shed light on the nature of the HD 108236 planets. For this purpose, we overplotted in Fig. 9 the planetary internal models from Brugger et al. (2017) and Acuña et al. (2021). These models consider a planetary interior structure divided in three layers: a Fe-rich core, a silicate mantle and a volatile layer (Brugger et al. 2017). The volatile layer is composed of water, which, depending on the surface pressure and temperature, presents different phases, including liquid, ice, steam, and supercritical. Thus, the models in Fig. 9 represent different ratios of core, mantle and volatile rich envelope with water either in the supercritical (SW) or liquid (LW) form. In the case of the SW mass-radius relationships, we have considered an equilibrium temperature of 1200 K, which is approximately the equilibrium temperature of planets *b* and *c*.

Thus, based on its location on the M-R diagram, planets *d* and *e* likely have a significant volatile envelope surrounding a large mantle. Unfortunately, with the data in hand, for these planets, it is not possible to say anything further about the relative fraction of these two layers. For planet *f*, for which the RVs suggest a rocky composition although the proportion between mantle and core components is not possible to constrain due to the large uncertainties we still have on its mass.

On the contrary, despite their mass uncertainties, it is possible to discriminate between a rocky planet with no volatiles and a volatile-rich planet for planets *b* and *c*. For these planets, we performed a more detailed interior structure analysis within a MCMC Bayesian framework, as described in Acuña et al. (2021). For this analysis, we considered the masses estimated

from the RVs (Table 6) and the radii presented in Table 4. With an equilibrium temperature of about 1200 K, these planets are highly irradiated; therefore, in our interior structure model, we assumed that the volatile layer is not water in a condensed form, but in supercritical and steam phases (Mousis et al. 2020; Acuña et al. 2021). The upper region of the volatile layer consists of an atmosphere in radiative-convective equilibrium that is coupled with the interior at 300 bar. The atmospheric temperature at the bottom of the volatile layer, $T_{300}$, the Bond albedo, and the atmospheric thickness were calculated with a 1D k-correlated model (Pluriel et al. 2019; Marcq et al. 2017). In addition, we also considered the Fe/Si mole ratio as an observable in our MCMC analysis, which is calculated as described in Sotin et al. (2007) and Brugger et al. (2017), using our HD 108236's abundance estimations of Fe, Si, and Mg (Table 2). We obtain a Fe/Si = 0.79 ± 0.08, which is below the solar value Fe/Si$_{⊕}$ = 0.96 (Sotin et al. 2007). The core mass fraction (CMF), the water mass fraction (WMF), and the atmospheric parameters, resulting from our interior structure analysis, are shown in Table 7.

We found that the densities of all planets in the system can be accounted by a core and a mantle with a Fe/Si mole ratio equal to that of the stellar host, together with a volatile layer on top. Under this assumption, in the Fe/Si mole ratio, planet *b* could have a thin atmosphere (<300 bar) or an envelope that could constitute up to 10% of its mass. To explore the possibility of planet *b* being a “dry” super-Earth, that is, composed of a core and a mantle only, we removed the Fe/Si mole ratio as a constraint in our MCMC Bayesian analysis, leaving only the observable constraints imposed by mass and radius. In addition, we fixed the WMF to zero. In this case, planet *b* is compatible with a dry
<table><caption>**Table 7.** Compositional and atmospheric parameters for HD 108236 planets *b* and *c*.</caption><thead><tr><th>Planet</th><th>CMF<sup>(a)</sup></th><th>WMF<sup>(b)</sup></th><th><i>T</i><sub>300</sub> <sup>(c)</sup> [K]</th><th><i>z</i><sub>atm</sub> <sup>(d)</sup> [km]</th><th><i>A</i><sub>B</sub> <sup>(e)</sup></th></tr></thead><tbody><tr><td><i>b</i></td><td>0.24 ± 0.06</td><td>0.05 ± 0.05</td><td>4071 ± 135</td><td>1035 ± 395</td><td>0.22 ± 0.01</td></tr><tr><td><i>c</i></td><td>0.20 ± 0.02</td><td>0.21 ± 0.05</td><td>3804 ± 123</td><td>733 ± 179</td><td>0.22 ± 0.01</td></tr></tbody></table>

**Notes.** We report here the parameters retrieved by our interior-atmospheric MCMC Bayesian analysis: <sup>(a)</sup>core mass fraction; <sup>(b)</sup>water mass fraction; <sup>(c)</sup>atmospheric temperature at 300 bar; <sup>(d)</sup>atmospheric thickness from 300 bar to 20 bar; and <sup>(e)</sup>Bond albedo.

<figure><img src="image_10.png" /><figcaption>**Fig. 10.** Ternary diagram representing the internal structure of HD 108236 *b* and *c*, where we show the estimated core mass fraction (CMF), water mass fraction (WMF) and mantle mass fraction (MMF) defined as $MMF = 1 - CMF - WMF$. The colored contours represent the $1\sigma$ uncertainties derived from the marginal posterior distribution function for each planet. For comparison, the green dot and red square indicate the position of Earth and Mercury in the ternary diagram, respectively.</figcaption></figure>

planet with $CMF = 0.04_{-0.04}^{+0.22}$ and a Fe/Si = $0.67 \pm 0.54$, having a low Fe content compared to the Earth value ($CMF_{⊕} = 0.32$). We observe that planet *c* is more volatile-rich than planet *b* (see Table 7 and Fig. 10). On the other hand, it would be necessary to better constrain the masses of planets *d* and *e* to narrow down the mean value and uncertainties of their volatile mass fractions. If their volatile contents are confirmed to be higher than that of planet *c*, the HD 108236 system could present an increasing volatile mass fraction trend with increasing semi-major axis, which has been noted in other multi-planetary systems with low-mass planets (e.g., Acuña et al. 2022). In addition, if the rocky nature of planet *f* is confirmed, for example via a more complete RV follow-up, it would be very interesting to explore the processes in the formation and/or evolution of the system, such as a Jeans atmospheric escape, required to explain why the more massive and less irradiated planet does not present a gaseous envelope unlike its two precedent planets.

## 8. Conclusions

In this work, we present the most comprehensive characterization of the HD 108236 system based on the photometric data gathered by CHEOPS and TESS space missions. We confirm the existence of a fifth transiting planet with an orbital period of 29.54 d. On the contrary, we found no indication of transit signals of a putative additional planet in a 10.9 d orbit in CHEOPS dedicated observations we performed. A planet larger than $0.42 R_{⊕}$ (with $e = b = 0$) would have been detectable in our dataset at $S/N = 1$. The general picture of the HD 108236 system is summarized in Fig. 8.

With our analysis, we greatly improved the characterization of the planetary properties, in particular, we refined the size estimations of the five planets in the system with uncertainties between 1.5 and 3%. Also, the extensive and detailed timing analysis facilitated an update for the ephemeris of all the planets and to report that TTVs were observed for planets *b* and *f*. Moreover, in the absence of a profuse RV monitoring to date, we performed mass estimations of the planets by different approaches. As illustrated in Table 6, the best agreement between the different applied methods was obtained for the inner planets. On the other hand, a large spread in the mass estimations is seen for the two outer planets. In this case, the constraints set by the timing analysis for that pair of planets in a 3:2 MMR, suggest a low density configuration for these bodies. A TTV analysis that would include additional transits for these two bodies would allow stronger constraints to be placed on their gravitational interaction and, thus, their density as well as the architecture of the system.

Thanks to its brightness, HD 108236 system is fitting for an intense and precise radial velocity follow-up, which together with our precise estimations of the planetary radii will provide definitive constraints on the density of the transiting bodies.

Based on the refined planetary properties, we computed the transmission spectroscopy metric (TSM) of Kempton et al. (2018) to estimate its suitability for a potential atmospheric follow-up with James Webb Space Telescope (JWST, Gardner et al. 2006). By using the RV masses (Table 6) and the upper limit of planet *d* ($M_d = 7.75 \pm 1.0$) from the estimations of Teske et al. (2021) and Bonfanti et al. (2021), we obtained the following TSM values for each planet (*b* to *f*): $54 \pm 26, 77 \pm 30, 73 \pm 10, 56 \pm 13$, and $9 \pm 4$, where the errors are dominated by the uncertainties on the planetary mass, with all of them below the suggested threshold for this range of planet sizes (TSM > 90). The potential for an JWST atmospheric follow-up will be finally defined when the densities of these planets are better determined. So far, our analysis suggests the presence of significant atmospheres in planets *c* to *e*. On the contrary, planet *b* is compatible with the presence of a thin atmosphere, or a Fe-poor dry planet, while planet *f* seems to have a rocky structure.

**Acknowledgements.** 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. Funding for the TESS mission is provided by NASA's Science Mission Directorate. 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. This paper includes data collected by the TESS mission that are publicly available from the Mikulski Archive for Space Telescopes (MAST). SH gratefully acknowledges CNES funding through the grant 837319. The authors acknowledge support from the Swiss NCCR PlanetS
and the Swiss National Science Foundation. LMS gratefully acknowledges financial support from the CRT foundation under Grant No. 2018.2323 'Gaseous or rocky? Unveiling the nature of small worlds'. This project was supported by the CNES. This work was also partially supported by a grant from the Simons Foundation (PI Queloz, grant number 327127). ACC and TW acknowledge support from STFC consolidated grant numbers ST/R000824/1 and ST/V000861/1, and UKSA grant number ST/R003203/1. S.G.S. acknowledges support from FCT through FCT contract nr. CEECIND/00826/2018 and POPH/FSE (EC). YA and MJH acknowledge the support of the Swiss National Fund under grant 200020_172746. We acknowledge support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grants ESP2016-80435-C2-1-R, ESP2016-80435-C2-2-R, PGC2018-098153-B-C33, PGC2018-098153-B-C31, ESP2017-87676-C5-1-R, MDM-2017-0737 Unidad de Excelencia María de Maeztu-Centro de Astrobiología (INTA-CSIC), as well as the support of the Generalitat de Catalunya/CERCA programme. The MOC activities have been supported by the ESA contract No. 4000124370. S.C.C.B. acknowledges support from FCT through FCT contracts nr. IF/01312/2014/CP1215/CT0004. XB, SC, DG, MF and AJ acknowledge their role as ESA-appointed CHEOPS science team members. AJR was supported by the SNSA. ACC acknowledges support from STFC consolidated grant numbers ST/R000824/1 and ST/V000861/1, and UKSA grant number ST/R003203/1. The Belgian participation to CHEOPS has been supported by the Belgian Federal Science Policy Office (BELSPO) in the framework of the PRODEX Program, and by the University of Liège through an ARC grant for Concerted Research Actions financed by the Wallonia-Brussels Federation. L.D. 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 Internacionalizaçã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.D. is supported in the form of work contract (DL 57/2016/CP1364/CT0004) funded by national funds through FCT. B.-O.D. acknowledges support from the Swiss National Science Foundation (PP00P2-190008). This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (project FOUR ACES, grant agreement No 724427). It has also been carried out in the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation (SNSF). D.E. acknowledges financial support from the Swiss National Science Foundation for project 200021_200726. M.F. and C.M.P. gratefully acknowledge the support of the Swedish National Space Agency (DNR 65/19, 174/18). D.G. gratefully acknowledges financial support from the CRT foundation under Grant No. 2018.2323 'Gaseousor rocky? Unveiling the nature of small worlds'. M.G. is an F.R.S.-FNRS Senior Research Associate. KGI is the ESA CHEOPS Project Scientist and is responsible for the ESA CHEOPS Guest Observers Programme. She does not participate in, or contribute to, the definition of the Guaranteed Time Programme of the CHEOPS mission through which observations described in this paper have been taken, nor to any aspect of target selection for the programme. This work was granted access to the HPC resources of MesoPSL financed by the Rgion Île-de-France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche. M.L. acknowledges support of the Swiss National Science Foundation under grant number PCEFP2_194576. P.M. acknowledges support from STFC research grant number ST/M001040/1. G.Sc., G.Pi., I.Pa., L.Bo., V.Na. and R.Ra. acknowledge the funding support from Italian Space Agency (ASI) regulated by "Accordo ASI-NAIF n. 2013-016-R.0 del 9 luglio 2013 e integrazione del 9 luglio 2015 CHEOPS Fasi A/B/C". I.R.I. acknowledges support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grant PGC2018-098153-B-C33, as well as the support of the Generalitat de Catalunya/CERCA programme. S.S. has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 833925, project STAREX). Gy.M.Sz. acknowledges the support of the Hungarian National Research, Development and Innovation Office (NKFIH) grant K-125015, a PRODEX Institute Agreement between the ELTE Eötvös Loránd University and the European Space Agency (ESA-D/SCI-LE-2021-0025), the Lendület LP2018-7/2021 grant of the Hungarian Academy of Science and the support of the city of Szombathely. V.V.G. is an F.R.S-FNRS Research Associate. N.A.W. acknowledges UKSA grant ST/R004838/1. This research has made use of computing facilities operated by CeSAM data center at LAM, Marseille, France. We also thank J.C. Meunier and J.C. Lamber for the support provided in the use of these facilities. Software: We gratefully acknowledge the open-source software which made this work possible: **astropy (Astropy Collaboration 2013, 2018)**, **python (Pérez & Granger 2007)**, **numpy (Harris et al. 2020)**, **scipy (Jones et al. 2001)**, **plotlib (Hunter 2007)**, **jupyter (Kluyver et al. 2016)**.

## References

Acuña, L., Deleuil, Magali, Mousis, Olivier, et al. 2021, A&A, 647, A53

Acuña, L., Lopez, T. A., Morel, T., et al. 2022, A&A, 660, A102

Adibekyan, V. Z., Sousa, S. G., Santos, N. C., et al. 2012, A&A, 545, A32

Adibekyan, V., Figueira, P., Santos, N. C., et al. 2015, A&A, 583, A94

Agol, E., Steffen, J., Sari, R., & Clarkson, W. 2005, MNRAS, 359, 567

Agol, E., Dorn, C., Grimm, S. L., et al. 2021, PSJ, 2, 1

Anderson, D. R., Collier Cameron, A., Hellier, C., et al. 2011, ApJ, 726, L19

Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

Barclay, T., Rowe, J. F., Lissauer, J. J., et al. 2013, Nature, 494, 452

Barros, S. C. C., Almenara, J. M., Demangeon, O., et al. 2015, MNRAS, 454, 4267

Benz, W., Broeg, C., Fortier, A., et al. 2021, Exp. Astron., 51, 109

Blackwell, D. E., & Shallis, M. J. 1977, MNRAS, 180, 177

Bonfanti, A., Ortolani, S., Piotto, G., & Nascimbeni, V. 2015, A&A, 575, A18

Bonfanti, A., Ortolani, S., & Nascimbeni, V. 2016, A&A, 585, A5

Bonfanti, A., Delrez, L., Hooton, M. J., et al. 2021, A&A, 646, A157

Borsato, L., Piotto, G., Gandolfi, D., et al. 2021, MNRAS, 506, 3810

Borucki, W. J., Koch, D., Basri, G., et al. 2010, Science, 327, 977

Brugger, B., Mousis, O., Deleuil, M., & Deschamps, F. 2017, ApJ, 850, 93

Campante, T. L., Barclay, T., Swift, J. J., et al. 2015, ApJ, 799, 170

Chen, J., & Kipping, D. 2017, ApJ, 834, 17

David, T. J., Petigura, E. A., Luger, R., et al. 2019, ApJ, 885, L12

Daylan, T., Pinglé, K., Wright, J., et al. 2021, AJ, 161, 85

Deck, K. M., Payne, M., & Holman, M. J. 2013, ApJ, 774, 129

Deck, K. M., Agol, B. S., & Agol, E. 2013, PASP, 125, 83

Delisle, J. B. 2017, A&A, 605, A96

Delisle, J. B., Ségransan, D., Dumusque, X., et al. 2018, A&A, 614, A133

Eastman, J., Gaudi, B. S., & Agol, E. 2013, PASP, 125, 83

Espinoza, N., Kossakowski, D., & Brahm, R. 2019, MNRAS, 490, 2262

Feroz, F., & Hobson, M. P. 2008, MNRAS, 384, 449

Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2019, Open J. Astrophys., 2, 10

Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220

Freudenthal, J., von Essen, C., Dreizler, S., et al. 2018, A&A, 618, A41

Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. 2018, PASP, 130, 044504

Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1

Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1

Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Sci. Rev., 123, 485

Gillon, M., Demory, B. O., Triaud, A. H. M. J., et al. 2009, A&A, 506, 359

Gillon, M., Triaud, A. H. M. J., Demory, B.-O., Van Grootel, V., et al. 2017a, Nat. Astron., 1, 0056

Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017b, Nature, 542, 456

Guerrero, N. M., Seager, S., Huang, C. X., et al. 2021, ApJS, 254, 39

Günther, M. N., & Daylan, T. 2019, Astrophysics Source Code Library [record ascl:1903.003]

Günther, M. N., & Daylan, T. 2021, ApJS, 254, 13

Hara, N. C., Unger, N., Delisle, J.-B., Díaz, R. F., & Ségransan, D. 2022, A&A, 663, A14

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

He, M. Y., Ford, E. B., & Ragozzine, D. 2021, AJ, 162, 216

Henrard, J., & Lemaitre, A. 1983, Celest. Mech., 30, 197

Hippke, M., & Heller, R. 2019, A&A, 623, A39

Holman, M. J., & Murray, N. W. 2005, Science, 307, 1288

Howell, S. B., Sobeck, C., Haas, M., et al. 2014, PASP, 126, 398

Hoyer, S., Guterman, P., Demangeon, O., et al. 2020, A&A, 635, A24

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, SPIE Conf. Ser., 9913, 99133E

Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open source scientific tools for Python

Kass, R. E., & Raftery, A. E. 1995, J. Am. Stat. Assoc., 90, 773

Kempton, E. M. R., Bean, J. L., Louie, D. R., et al. 2018, PASP, 130, 114401

Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, eds. F. Loizides, & B. Schmidt (Amsterdam: IOS Press), 87

Kreidberg, L. 2015, PASP, 127, 1161

Lacedelli, G., Wilson, T. G., Malavolta, L., et al. 2022, MNRAS, 511, 4551

Leleu, A., Alibert, Y., Hara, N. C., et al. 2021a, A&A, 649, A26

Leleu, A., Chatel, G., Udry, S., et al. 2021b, A&A, 655, A66

Lendl, A., Csizmadia, S., Deline, A., et al. 2020, A&A, 643, A94

Lindegren, L., Bastian, U., Biermann, M., et al. 2021, A&A, 649, A4

Lithwick, Y., Xie, J., & Wu, Y. 2012, ApJ, 761, 122

Marcq, E., Salvador, A., Massol, H., & Davaille, A. 2017, J. Geophys. Res. Planets, 122, 1539

Marigo, P., Girardi, L., Bressan, A., et al. 2017, ApJ, 835, 77
Maxted, P. F. L., Ehrenreich, D., Wilson, T. G., et al. 2022, MNRAS, 514, 77

Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20

McArthur, B. E., Endl, M., Cochran, W. D., et al. 2004, ApJ, 614, L81

Mills, S. M., Fabrycky, D. C., Migaszewski, C., et al. 2016, Nature, 533, 509

Miralda-Escudé, J. 2002, ApJ, 564, 1019

Motalebi, F., Udry, S., Gillon, M., et al. 2015, A&A, 584, A72

Mousis, O., Deleuil, M., Aguichine, A., et al. 2020, ApJ, 896, L22

Nesvorný, D., & Vokrouhlický, D. 2016, ApJ, 823, 72

Otegi, J. F., Bouchy, F., & Helled, R. 2020, A&A, 634, A43

Pérez, F., & Granger, B. E. 2007, Comput. Sci. Eng., 9, 21

Pluriel, W., Marcq, E., & Turbet, M. 2019, Icarus, 317, 583

Rasmussen, C. E., & Williams, C. K. I. 2005, Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning) (Cambridge: The MIT Press)

Salmon, S. J. A. J., Van Grootel, V., Buldgen, G., Dupret, M. A., & Eggenberger, P. 2021, A&A, 646, A7

Scuflaire, R., Théado, S., Montalbán, J., et al. 2008, Ap&SS, 316, 83

Smith, J. C., Stumpe, M. C., Van Cleve, J. E., et al. 2012, PASP, 124, 1000

Sotin, C., Grasset, O., & Mocquet, A. 2007, Icarus, 191, 337

Sousa, S. G. 2014, Determination of Atmospheric Parameters of B-, A-, F- and G-Type Stars. Series: GeoPlanet: Earth and Planetary Sciences, eds. E. Niemczura, B. Smalley, & W. Pych (Springer International Publishing (Cham)), 297

Speagle, J. S. 2020, MNRAS, 493, 3132

Stumpe, M. C., Smith, J. C., Van Cleve, J. E., et al. 2012, PASP, 124, 985

Stumpe, M. C., Smith, J. C., Catanzarite, J. H., et al. 2014, PASP, 126, 100

Teske, J., Wang, S. X., Wolfgang, A., et al. 2021, ApJS, 256, 33

Vanderburg, A., Huang, C. X., Rodriguez, J. E., et al. 2019, ApJ, 881, L19

Wilcoxon, F. 1945, Biometrics Bull., 1, 80

Wilson, T. G., Goffo, E., Alibert, Y., et al. 2022, MNRAS, 511, 1043

¹ Aix-Marseille Univ., CNRS, CNES, LAM, 38 rue Frédéric Joliot-Curie, 13388 Marseille, France
e-mail: sergio.hoyer@lam.fr

² Space Research Institute, Austrian Academy of Sciences, Schmiedl-Strasse 6, 8042 Graz, Austria

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

⁴ Physikalisches Institut, University of Bern, Gesellsschaftstrasse 6, 3012 Bern, Switzerland

⁵ Dipartimento di Fisica, Universita degli Studi di Torino, via Pietro Giuria 1, 10125 Torino, Italy

⁶ Center for Space and Habitability, Gesellsschaftstrasse 6, 3012 Bern, Switzerland

⁷ Department of Astronomy, Stockholm University, AlbaNova University Center, 10691 Stockholm, Sweden

⁸ Department of Astronomy, Stockholm University, 10691 Stockholm, Sweden

⁹ ETH Zurich, Department of Physics, Wolfgang-Pauli-Strasse 2, 8093 Zurich, Switzerland

¹⁰ Cavendish Laboratory, JJ Thomson Avenue, Cambridge CB3 0HE, UK

¹¹ Centre for Exoplanet Science, SUPA School of Physics and Astronomy, University of St Andrews, North Haugh, St Andrews KY16 9SS, UK

¹² Instituto de Astrofísica e Ciencias do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150-762 Porto, Portugal

¹³ Instituto de Astrofísica de Canarias, 38200 La Laguna, Tenerife, Spain

¹⁴ Departamento de Astrofísica, Universidad de La Laguna, 38206 La Laguna, Tenerife, Spain

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

¹⁶ Institut d'Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain

¹⁷ ESTEC, European Space Agency, 2201AZ Noordwijk, NL, The Netherlands

¹⁸ Admatis, K. and Kálmán Street, 3534 Miskolc, Hungary

¹⁹ Depto. de Astrofísica, Centro de Astrobiología (CSIC-INTA), ESAC campus, 28692 Villanueva de la Cañada (Madrid), Spain

²⁰ Departamento de Física e Astronomia, Faculdade de Ciencias, Universidade do Porto, Rua do Campo Alegre, 4169-007 Porto, Portugal

²¹ Max Planck Institute for Extraterrestrial Physics, Gießenbachstraße 1, 85748 Garching bei München, Germany

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

²³ Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

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

²⁵ Université de Paris, Institut de physique du globe de Paris, CNRS, 75005 Paris, France

²⁶ Centre for Mathematical Sciences, Lund University, Box 118, 221 00 Lund, Sweden

²⁷ Astrobiology Research Unit, Université de Liège, Allée du 6 Août 19C, 4000 Liège, Belgium

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

²⁹ Leiden Observatory, University of Leiden, PO Box 9513, 2300 RA Leiden, The Netherlands

³⁰ Department of Space, Earth and Environment, Chalmers University of Technology, Onsala Space Observatory, 43992 Onsala, Sweden

³¹ University of Vienna, Department of Astrophysics, Türkenschanzstrasse 17, 1180 Vienna, Austria

³² Department of Physics, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK

³³ Science and Operations Department – Science Division (SCI-SC), Directorate of Science, European Space Agency (ESA), European Space Research and Technology Centre (ESTEC), Keplerlaan 1, 2201-AZ Noordwijk, The Netherlands

³⁴ NASA Ames Research Center, Moffett Field, CA, 94035, USA

³⁵ Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, 1121 Budapest, Konkoly Thege Miklós út 15-17, Hungary

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

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

³⁸ Center for Astrophysics, Harvard & Smithsonian, 60 Garden Street, Cambridge, MA, 02138, USA

³⁹ Institut d'astrophysique de Paris, CNRS, UMR 7095, Sorbonne Université, 98 bis bd Arago, 75014 Paris, France

⁴⁰ Department of Astrophysics, University of Vienna, Tuerkenschanzstrasse 17, 1180 Vienna, Austria

⁴¹ Astrophysics Group, Keele University, Staffordshire, ST5 5BG, UK

⁴² INAF, Osservatorio Astrofisico di Catania, Via S. Sofia 78, 95123 Catania, Italy

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

⁴⁴ Dipartimento di Fisica e Astronomia “Galileo Galilei”, Universita degli Studi di Padova, Vicolo dell'Osservatorio 3, 35122 Padova, Italy

⁴⁵ Zentrum für Astronomie und Astrophysik, Technische Universität Berlin, Hardenbergstr. 36, 10623 Berlin, Germany

⁴⁶ Institut für Geologische Wissenschaften, Freie Universität Berlin, 12249 Berlin, Germany

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

⁴⁸ Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA, 02139, USA

⁴⁹ Department of Aeronautics and Astronautics, MIT, 77 Massachusetts Avenue, Cambridge, MA, 02139, USA

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

⁵¹ MTA-ELTE Exoplanet Research Group, 9700 Szombathely, Szent Imre h. u. 112, Hungary

⁵² SETI Institute, Mountain View, CA 94043, USA

⁵³ Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK

⁵⁴ Department of Astrophysical Sciences, Princeton University, 4 Ivy Ln., Princeton NJ 08544, USA
**Appendix A: Raw and detrended CHEOPS light curves**

In Fig. A.1 we show the raw light curves as output by the CHEOPS Data Reduction Pipeline (DRP-v13, Hoyer et al. 2020). These light curves are usually affected by systematics as a function of the rotation angle of the Field of View along the pointing direction. In Figs. A.2 and A.3 we show the detrended version of the CHEOPS light curves (see Sect. 4.1 for details), together with the best models obtained from the global modeling of the system (Sect. 4).
<figure><img src="image_11.png" /><figcaption>**Fig. A.1.** Normalized CHEOPS raw light curves of HD 108236 used in this work. Large outliers have been clipped out to improve visualization. The CHEOPS label number on top of each panel correspond to the ID given in Table 3.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. A.2.** Normalized CHEOPS detrended light curves of HD 108236 used in this work and its 10 min bins are shown, respectively, with light and dark blue symbols. The best transit models obtained from the global analysis of the system are represented by the yellow curve. The residuals are shown in the bottom panels.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. A.3.** Normalized CHEOPS detrended light curves of HD 108236 used in this work and its 10 min bins are shown, respectively, with light and dark blue symbols. The best transit models obtained from the global analysis of the system are represented by the yellow curve. The residuals are shown in the bottom panels.</figcaption></figure>
## Appendix B: Initial priors of the LC analysis

In Table B.1, we list the priors distributions used for the global fit of all the photometric data described in Sect. 4, consisting of 34 TESS and 15 CHEOPS light curves.

<table><caption>**Table B.1.** Priors for each of the parameters fitted in the global photometric analysis of CHEOPS and TESS data.</caption><thead><tr><th>Parameter</th><th>Prior</th></tr></thead><tbody><tr><td></td><td><strong>Star</strong></td></tr><tr><td>ρ<sub>★</sub> [kg m<sup>-3</sup>]</td><td>TN(1790, 110, 1680, 1900)</td></tr><tr><td>q<sub>1,TESS</sub></td><td>N(0.34, 0.05)</td></tr><tr><td>q<sub>2,TESS</sub></td><td>N(0.27, 0.05)</td></tr><tr><td>q<sub>1,CHEOPS</sub></td><td>N(0.46, 0.05)</td></tr><tr><td>q<sub>2,CHEOPS</sub></td><td>N(0.32, 0.05)</td></tr><tr><td></td><td><strong>Planet b</strong></td></tr><tr><td>R<sub>p</sub>/R<sub>★</sub></td><td>N(0.01688, 0.01)</td></tr><tr><td>T<sub>0-2</sub> 450 000. [BJD]</td><td>N(8572.1128, 0.01)</td></tr><tr><td>P [d]</td><td>N(3.796, 0.001)</td></tr><tr><td></td><td><strong>Planet c</strong></td></tr><tr><td>R<sub>p</sub>/R<sub>★</sub></td><td>N(0.0217, 0.01)</td></tr><tr><td>T<sub>0-2</sub> 450 000. [BJD]</td><td>N(8572.3949, 0.01)</td></tr><tr><td>P [d]</td><td>N(6.203, 0.001)</td></tr><tr><td></td><td><strong>Planet d</strong></td></tr><tr><td>R<sub>p</sub>/R<sub>★</sub></td><td>N(0.0266, 0.01)</td></tr><tr><td>T<sub>0-2</sub> 450 000. [BJD]</td><td>N(8571.3368, 0.01)</td></tr><tr><td>P [d]</td><td>N(14.1758, 0.001)</td></tr><tr><td></td><td><strong>Planet e</strong></td></tr><tr><td>R<sub>p</sub>/R<sub>★</sub></td><td>N(0.0322, 0.01)</td></tr><tr><td>T<sub>0-2</sub> 450 000. [BJD]</td><td>N(8586.5677, 0.01)</td></tr><tr><td>P [d]</td><td>N(19.59, 0.001)</td></tr><tr><td></td><td><strong>Planet f</strong></td></tr><tr><td>R<sub>p</sub>/R<sub>★</sub></td><td>N(0.0211, 0.01)</td></tr><tr><td>T<sub>0-2</sub> 450 000. [BJD]</td><td>N(8616.039, 0.01)</td></tr><tr><td>P [d]</td><td>N(29.5398, 0.001)</td></tr><tr><td></td><td><strong>All planets</strong></td></tr><tr><td>b [R<sub>★</sub>]</td><td>U(0.001, 0.99)</td></tr><tr><td>√esin(ω)</td><td>U(-0.3, 0.3)</td></tr><tr><td>√ecos(ω)</td><td>U(-0.3, 0.3)</td></tr><tr><td></td><td><strong>Instrumental</strong></td></tr><tr><td>mflux_TESS<sub>i</sub></td><td>N(0.0, 0.1)</td></tr><tr><td>σ<sub>w_TESS<sub>i</sub></sub> [ppm]</td><td>logU(0.1, 1000)</td></tr><tr><td>mflux_CHEOPS<sub>i</sub></td><td>N(0.0, 0.1)</td></tr><tr><td>σ<sub>w_CHEOPS<sub>i</sub></sub> [ppm]</td><td>logU(0.1, 1000)</td></tr><tr><td>c<sub>1</sub>(t) CHEOPS<sub>i</sub></td><td>U(-1.0, 1.0)</td></tr><tr><td>c<sub>2</sub>(t<sup>2</sup>) CHEOPS<sub>i</sub></td><td>U(-1.0, 1.0)</td></tr></tbody></table>

**Notes.** The following distributions were used: Normal (N(μ₀, σ₀)), Truncated Normal (TN(μ₀, σ₀, limit_low, limit_upper)), Uniform (U(limit_low, limit_upper)), log-Uniform (logU(limit_low, limit_upper)). The instrumental parameters are defined for each fitted light curve.

## Appendix C: Role of the parametrization in the fitted parameters: e, ω and b.

We explored different configurations of the system modeling with juliet to compare the resulting values of the orbits' eccentricities and any possible effect it might have in other planetary parameters. As described in Sect. 4.2, we used the parametrization $f_s = \sqrt{e} \sin(\omega)$ and $f_c = \sqrt{e} \cos(\omega)$ to fit for the eccentricity (e) and the argument of the periastron ($\omega$). In addition, we per-

formed two extra modelings by either fitting e and ω directly or assuming e=0 and ω = 90°. The different models were also compared with the fit performed using Allesfitter. We found that, in general, all the fitted parameters were consistent well within 1σ. The largest differences are found in the eccentricities as illustrated in Fig. D.2. We notice that, for planets b, c and d the eccentricities derived using ($f_s, f_c$) parametrization are systematically large in comparison with the other fits. We found that in this case, eccentricities are correlated with the impact parameters, which were also discrepant with the values obtained with the other fits (see right panels in Fig. D.2). These differences likely raise because juliet does not fit for $b_{circ}$ (the impact parameter of the circular orbit, see e.g., Gillon et al. 2009), but for the eccentric $b'$, hence producing correlations between $b'$ and the resulting e as:

$$b' = b_{\text{circ}}(1 + e \sin(\omega))/(1 - e^2). \quad (C.1)$$

Allesfitter, instead, assumes ($R_p+R_s)/a$ and $cos(i)$ as priors, whose combination essentially gives $b_{circ}$. Thus, with Allesfitter, we avoided $b'$ versus e correlations; we also obtained non-zero eccentricities, but only $1.5\sigma$ away from the zero (although Allesfitter's uncertainties are systematically larger). It is also known that fitting directly e and ω biases towards low eccentricities (e.g Anderson et al. 2011; Eastman et al. 2013), behavior that we also confirm. Therefore, we attribute the larger eccentricities reported in Table 4 to the parametrization used by juliet. We recall the difficulty in obtaining a robust estimate of the orbital eccentricities with photometry of primary transits alone and thus a detailed RV analysis is required for achieving a final conclusion in this regards.

## Appendix D: Ruling out the putative planet g

Confirming the five-planet architecture of the system or preferring the six-planet scenario is an example of model selection problem. Before acquiring the last CHEOPS light curve of the putative planet g (Fig. 4, bottom panel), we compared the two scenarios by performing several pairs of comprehensive analyses, which involved all the other CHEOPS light curves (where planet g is expected to transit twice) and the two TESS sectors, 10 and 11, (where planet g is expected to transit five times). To address the model selection problem, the analyses were performed through the dynamic nested sampling technique (see e.g., Feroz & Hobson 2008; Feroz et al. 2019) implemented in the dynasty package (Speagle 2020) and incorporated in Allesfitter. In fact, dynasty outputs the Bayesian evidence $Z$ of the adopted model, so that it is then straightforward to decide whether to prefer a more complicated model (e.g., $M_1$) against a simpler one ($M_0$) by computing the Bayes factor $B_{10}$ as defined in Kass & Raftery (1995):

$$B_{10} = \frac{Z_1}{Z_0}, \quad (D.1)$$

where $Z_i$ ($i = 0, 1$) is the Bayesian evidence of the $i^{th}$ model. As reported by Kass & Raftery (1995), we recall that a very strong evidence against the null hypothesis (i.e. the more complicated $M_1$ model is accepted) occurs when $\ln B_{10} > 5$.

Following a first analysis, where we included only planets from b to f (model $M_5$), we added the putative planet g to our input file, adopting uniform priors in agreement with the guess values from Daylan et al. (2021) (model $M_6$). The code was actually able to model also the sixth planet, but we obtained a transit depth $\delta = 127 \pm 21$ ppm, which is half the
value suggested by Daylan et al. (2021) and ~6σ away. Moreover, the posterior distributions of P and T₀ have heavy tails, with the uncertainty on T₀ that is one order of magnitude larger than what we found for the other planets. This indicates a poor convergence, probably due to the code attempts of modeling transit-like wiggles. We also looked at the rms of the residuals coming from the five-planet and six-planet analyses (rms₅ and rms₆, respectively), considering whether rms₆ values were statistically lower than rms₅. To this end we applied the Wilcoxon test (Wilcoxon 1945), whose null hypothesis is that the two samples comes from a single population (i.e., they are statistically equal) obtaining *p*-values > 0.79, which provides no evidence for preferring the six-planet scenario. On the other hand, we computed the Bayes factor *B*₆₅ comparing *M*₆ vs. *M*₅ obtaining ln *B*₆₅ = +8.9, which would strongly favor the existence of planet *g*.

We then decided to test whether the Bayes factor on its own is actually a reliable indicator in predicting the existence of a shallow transit-like signals. To this end, we repeated the analysis of all the light curves mentioned above and replaced the sixth planet candidate from Daylan et al. (2021) with a fake planet having the same characteristics, except for its T₀. The fake transits were placed at those timings where we could see a residual noise compatible with weak transit-like features. For the first stress test, we set T₀ = 2458959.08 BJD simulating a transit during CHEOPS observation 2 before the transit of planet *b* (Fig. A.2, right-hand side panel of the first row): model *M*<sub>f2</sub>. After the Allesfitter run, the Bayes factor comparing *M*<sub>f2</sub> with *M*<sub>5</sub> gave ln *B*<sub>f25</sub> = -0.3, which states that the existence of this sixth fake planet is not supported. We then performed another stress test, setting T₀ = 2458968.0 BJD so to mimic a transit during CHEOPS observation 3 (Fig. A.2, left-hand side panel of the second row) before the transit of planet *d*: model *M*<sub>f3</sub>. This time the comparison with *M*<sub>5</sub> yielded to ln *B*<sub>f35</sub> = +4.5, hence favoring the six-planet scenario containing a fake planet though. This is quite surprising, as it means that the transit-like wiggles of the order of 100 ppm occurring at the transit timings of this fake planet drive the Bayes factor in supporting *M*<sub>f3</sub>. As a double check, we repeated the latter analysis, but looking for a deeper transit (δ = 500 ppm) and the new proposed model (*M*<sub>fd3</sub>) was finally rejected by the Bayes factor criterion as expected (ln *B*<sub>fd35</sub> = -0.3).

As the answer from the Bayes factor is not definitive, we decided to schedule one last CHEOPS observation of the putative planet *g* (observation 15), after ensuring its transit would not overlap with any of the other planets. The detrended light curve is shown in Fig. 4 (bottom planet) and our global fitting did not detect any relevant transit feature.

To further assess the likelihood of a transiting body within this CHEOPS observation we conduct a more statistically rigorous analysis of this visit. Firstly, we noticed that this dataset suffers from the so-called ramp effect at the beginning of the observations that likely comes from changes in telescope temperature (Maxted et al. 2022) that manifests as point-spread function (PSF) shape changes. Moreover, there appears to be roll angle trends that may result as PSF shape changes due to the rotating field of view. To combat these effects, we utilize a novel PSF-based detrending method recently reported in Wilson et al. (2022) to remove these effects in CHEOPS observations of TOI-1064. This tool performs a principal component analysis (PCA) on the auto-correlation function of the CHEOPS subarray images in order to assess subtle changes in the PSF shape. The significant components to be used for decorrelation are then selected via a leave-one-out-cross-validation method. The result of the detrending using these vectors is shown in Fig. D.1 with

<figure><img src="image_14.png" /><figcaption>**Fig. D.1.** CHEOPS light curve of the final visit of HD 108236; (top panel) the DEFAULT DRP fluxes in blue, (second panel) the raw data in blue plus the linear model built from the PCA PSF method in black, (third panel) detrended fluxes in blue and the best fit transit model in black, and (bottom panel) the residuals to the fit in blue.</figcaption></figure>

the constructed linear model presented in the second panel. Satisfied with the correction of the systematic noise sources, we subsequently computed the true and false inclusion probabilities (TIP and FIP; Hara et al. 2022) of the presence of a transit in this dataset by fitting 0 and 1 planet models combined with the linear model determined via the PSF PCA method mentioned above and calculating the TIP and FIP values via comparison of the Bayesian evidence and posterior distributions of the fits. For all transit T₀ values within the observing window of the light curve, we find FIP~1 that statistically confirms that there is no transit within this window and rules out the presence of a planet *g* on this period. A definitive response may likely come from future RV measurements; in any case, the available data ultimately do not support the existence of planet *g*.
<figure><img src="image_15.png" /><figcaption>**Fig. D.2.** Comparison of fitted values using different parameterizations/priors. We show the fitted eccentricity, *e* (left panels), argument of periastron *ω* (middle panels), and impact parameter *b* (right panels) of each planet. In the x-axis we denote the used method: the (*f*<sub>s</sub>,*f*<sub>c</sub>) parametrization, fixing the eccentricity to zero (*e* = 0) and by jumping directly on *e* and *ω* (*e*,*ω*). We also show the resulting values of the modeling with Allesfiter code (A).</figcaption></figure>

**Appendix E: Transit times**

The results of the timing analysis of the transits as described in
Sect. 5 are shown in Tables E.1 and E.2. These transit times are
the result of the light curve fits for each planet separately, using
as priors the results of the global photometric fit described in
Sect. 4.
<table><caption>**Table E.1.** Central time of each of the transits of planet *b* and *c* analyzed in this work.</caption><thead><tr><th rowspan="2">Epoch</th><th colspan="3">Planet b</th><th rowspan="2">Source</th></tr><tr><th>T<sub>c</sub></th><th>+σ</th><th>-σ</th></tr></thead><tbody><tr><td>0</td><td>8572.1080</td><td>0.0038</td><td>0.0032</td><td>TESS</td></tr><tr><td>1</td><td>8575.9000</td><td>0.0055</td><td>0.0060</td><td>TESS</td></tr><tr><td>2</td><td>8579.7025</td><td>0.0069</td><td>0.0082</td><td>TESS</td></tr><tr><td>4</td><td>8587.2939</td><td>0.0025</td><td>0.0039</td><td>TESS</td></tr><tr><td>5</td><td>8591.0970</td><td>0.0078</td><td>0.0043</td><td>TESS</td></tr><tr><td>6</td><td>8594.8917</td><td>0.0039</td><td>0.0038</td><td>TESS</td></tr><tr><td>7</td><td>8598.6734</td><td>0.0039</td><td>0.0028</td><td>TESS</td></tr><tr><td>8</td><td>8602.4686</td><td>0.0052</td><td>0.0040</td><td>TESS</td></tr><tr><td>9</td><td>8606.2730</td><td>0.0036</td><td>0.0050</td><td>TESS</td></tr><tr><td>11</td><td>8613.8544</td><td>0.0068</td><td>0.0044</td><td>TESS</td></tr><tr><td>12</td><td>8617.6598</td><td>0.0021</td><td>0.0021</td><td>TESS</td></tr><tr><td>13</td><td>8621.4498</td><td>0.0027</td><td>0.0042</td><td>TESS</td></tr><tr><td>102</td><td>8959.3020</td><td>0.0020</td><td>0.0019</td><td>CHEOPS</td></tr><tr><td>105</td><td>8970.6821</td><td>0.0028</td><td>0.0027</td><td>CHEOPS</td></tr><tr><td>110</td><td>8989.6601</td><td>0.0028</td><td>0.0024</td><td>CHEOPS</td></tr><tr><td>112</td><td>8997.2494</td><td>0.0011</td><td>0.0013</td><td>CHEOPS</td></tr><tr><td>113</td><td>9001.0425</td><td>0.0026</td><td>0.0052</td><td>CHEOPS</td></tr><tr><td>194</td><td>9308.5086</td><td>0.0044</td><td>0.0037</td><td>TESS</td></tr><tr><td>195</td><td>9312.3224</td><td>0.0042</td><td>0.0043</td><td>TESS</td></tr><tr><td>196</td><td>9316.0992</td><td>0.0052</td><td>0.0035</td><td>TESS</td></tr><tr><td>198</td><td>9323.6931</td><td>0.0026</td><td>0.0028</td><td>TESS</td></tr><tr><td>199</td><td>9327.4805</td><td>0.0040</td><td>0.0078</td><td>TESS</td></tr><tr><td>200</td><td>9331.2915</td><td>0.0059</td><td>0.0079</td><td>TESS</td></tr></tbody></table><table><thead><tr><th rowspan="2">Epoch</th><th colspan="3">Planet c</th><th rowspan="2">Source</th></tr><tr><th>T<sub>c</sub></th><th>+σ</th><th>-σ</th></tr></thead><tbody><tr><td>0</td><td>8572.3921</td><td>0.0015</td><td>0.0013</td><td>TESS</td></tr><tr><td>1</td><td>8578.5998</td><td>0.0018</td><td>0.0013</td><td>TESS</td></tr><tr><td>2</td><td>8584.7979</td><td>0.0038</td><td>0.0038</td><td>TESS</td></tr><tr><td>3</td><td>8591.0124</td><td>0.0034</td><td>0.0036</td><td>TESS</td></tr><tr><td>5</td><td>8603.4095</td><td>0.0025</td><td>0.0019</td><td>TESS</td></tr><tr><td>6</td><td>8609.6151</td><td>0.0022</td><td>0.0012</td><td>TESS</td></tr><tr><td>7</td><td>8615.8166</td><td>0.0041</td><td>0.0047</td><td>TESS</td></tr><tr><td>8</td><td>8622.0309</td><td>0.0013</td><td>0.0015</td><td>TESS</td></tr><tr><td>56</td><td>8919.7989</td><td>0.0013</td><td>0.0016</td><td>CHEOPS</td></tr><tr><td>64</td><td>8969.4285</td><td>0.0020</td><td>0.0022</td><td>CHEOPS</td></tr><tr><td>65</td><td>8975.6328</td><td>0.0005</td><td>0.0005</td><td>CHEOPS</td></tr><tr><td>66</td><td>8981.8337</td><td>0.0021</td><td>0.0019</td><td>CHEOPS</td></tr><tr><td>70</td><td>9006.6546</td><td>0.0008</td><td>0.0008</td><td>CHEOPS</td></tr><tr><td>119</td><td>9310.6308</td><td>0.0024</td><td>0.0038</td><td>TESS</td></tr><tr><td>120</td><td>9316.8357</td><td>0.0014</td><td>0.0013</td><td>TESS</td></tr><tr><td>121</td><td>9323.0385</td><td>0.0015</td><td>0.0014</td><td>TESS</td></tr><tr><td>122</td><td>9329.2407</td><td>0.0014</td><td>0.0013</td><td>TESS</td></tr></tbody></table>

**Notes.** The transit times are expressed in BJD<sub>TDB</sub>-2 450 000.

<table><caption>**Table E.2.** Central time of each of the transits of planet *d*, *e* and *f* analyzed in this work.</caption><thead><tr><th rowspan="2">Epoch</th><th colspan="3">Planet d</th><th rowspan="2">Source</th></tr><tr><th>T<sub>c</sub></th><th>+σ</th><th>-σ</th></tr></thead><tbody><tr><td>0</td><td>8571.3355</td><td>0.0016</td><td>0.0014</td><td>TESS</td></tr><tr><td>1</td><td>8585.5137</td><td>0.0019</td><td>0.0014</td><td>TESS</td></tr><tr><td>2</td><td>8599.6872</td><td>0.0015</td><td>0.0011</td><td>TESS</td></tr><tr><td>3</td><td>8613.8654</td><td>0.0034</td><td>0.0036</td><td>TESS</td></tr><tr><td>28</td><td>8968.2577</td><td>0.0011</td><td>0.0009</td><td>CHEOPS</td></tr><tr><td>30</td><td>8996.6104</td><td>0.0049</td><td>0.0043</td><td>CHEOPS</td></tr><tr><td>52</td><td>9308.4794</td><td>0.0028</td><td>0.0025</td><td>TESS</td></tr><tr><td>53</td><td>9322.6590</td><td>0.0022</td><td>0.0023</td><td>TESS</td></tr><tr><th rowspan="2">Epoch</th><th colspan="3">Planet e</th><th rowspan="2">Source</th></tr><tr><th>T<sub>c</sub></th><th>+σ</th><th>-σ</th></tr><tr><td>0</td><td>8586.5685</td><td>0.0014</td><td>0.0013</td><td>TESS</td></tr><tr><td>1</td><td>8606.1591</td><td>0.0012</td><td>0.0011</td><td>TESS</td></tr><tr><td>17</td><td>8919.5979</td><td>0.0006</td><td>0.0006</td><td>CHEOPS</td></tr><tr><td>21</td><td>8997.9600</td><td>0.0008</td><td>0.0008</td><td>CHEOPS</td></tr><tr><td>37</td><td>9311.4060</td><td>0.0011</td><td>0.0013</td><td>TESS</td></tr><tr><td>38</td><td>9330.9942</td><td>0.0023</td><td>0.0020</td><td>TESS</td></tr><tr><th rowspan="2">Epoch</th><th colspan="3">Planet f</th><th rowspan="2">Source</th></tr><tr><th>T<sub>c</sub></th><th>+σ</th><th>-σ</th></tr><tr><td>-1</td><td>8586.4885</td><td>0.0044</td><td>0.0038</td><td>TESS</td></tr><tr><td>0</td><td>8616.0337</td><td>0.0059</td><td>0.0050</td><td>TESS</td></tr><tr><td>12</td><td>8970.5272</td><td>0.0015</td><td>0.0018</td><td>CHEOPS</td></tr><tr><td>24</td><td>9324.9808</td><td>0.0054</td><td>0.0037</td><td>CHEOPS</td></tr><tr><td>24</td><td>9324.9857</td><td>0.0088</td><td>0.0087</td><td>TESS</td></tr></tbody></table>

**Notes.** The transit times are expressed in BJD<sub>TDB</sub>-2 450 000.
Appendix F: Radial velocities fit

The radial velocity measurements reported by Teske et al. (2021) were fitted using juliet package. The priors of the orbital parameters were defined as normal distributions using the fitted values and uncertainties of the photometric fit (Sect. 4) reported in Table 4. The full RV model as a function of time, and the phased model for each planet are shown in Figures F.1 and F.2, respectively. The corner plot of the posterior distributions of the fitted parameters are shown in Fig. F.3.

<figure><img src="image_16.png" /><figcaption>**Fig. F.1.** Radial velocity measurements from Teske et al. (2021) (blue symbols) with the best fitted model (black solid curve) and its uncertain-
ties (red regions). The bottom panels are closer views of the top panel
for the relevant time coverage of the observations. The RV error bars
include the estimated jitter.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. F.2.** Phased radial velocities with the full model removed except for the component of the plotted planet (solid curve).</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. F.3.** Corner plot of the posterior distributions of the parameters estimated by the modeling of the RV measurements of HD 108236.</figcaption></figure>