# Investigating the architecture and internal structure of the TOI-561 system planets with CHEOPS, HARPS-N, and TESS

G. Lacedelli, T. G. Wilson, L. Malavolta, M. J. Hooton, A. Collier Cameron, Y. Alibert, A. Mortier, A. Bonfanti, R. D. Haywood, S. Hoyer, et al.

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

## To cite this version:

G. Lacedelli, T. G. Wilson, L. Malavolta, M. J. Hooton, A. Collier Cameron, et al., Investigating the architecture and internal structure of the TOI-561 system planets with CHEOPS, HARPS-N, and TESS. Monthly Notices of the Royal Astronomical Society, 2022, 511, pp.4551-4571. 10.1093/mn-ras/stac199. insu-03647049

HAL Id: insu-03647049

https://insu.hal.science/insu-03647049

Submitted on 20 Apr 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Investigating the architecture and internal structure of the TOI-561 system planets with CHEOPS, HARPS-N, and TESS

G. Lacedelli<sup>1,2*</sup> T. G. Wilson<sup>3</sup> L. Malavolta<sup>1,2</sup> M. J. Hooton<sup>4</sup> A. Collier Cameron<sup>5,3</sup>

Y. Alibert, <sup>4,5</sup> A. Mortier<sup>6,7</sup> A. Bonfanti, <sup>8</sup> R. D. Haywood, <sup>9</sup> S. Hoyer<sup>10</sup> G. Piotto, <sup>1,2</sup> A. Bekkelien, <sup>11</sup>

A. M. Vanderburg, <sup>12,13</sup> W. Benz, <sup>4,5</sup> X. Dumusque, <sup>11</sup> A. Deline, <sup>11</sup> M. López-Morales, <sup>14</sup> L. Borsato<sup>1,2</sup>

K. Rice<sup>15,16</sup> L. Fossati, <sup>8</sup> D. W. Latham, <sup>14</sup> A. Brandeker, <sup>17</sup> E. Poretti, <sup>18,19</sup> S. G. Sousa<sup>1,20</sup>

A. Sozzetti<sup>21</sup> S. Salmon, <sup>11</sup> C. J. Burke, <sup>12</sup> V. Van Grootel, <sup>22</sup> M. M. Fausnaugh, <sup>12</sup> V. Adibekyan, <sup>20</sup>

C. X. Huang, <sup>12,23</sup> H. P. Osborn<sup>1,5,12</sup> A. J. Mustill<sup>24</sup> E. Pallé, <sup>25</sup> V. Bourrier, <sup>11</sup> V. Nascimbeni<sup>1,2</sup>

R. Alonso, <sup>25,26</sup> G. Anglada, <sup>27,28</sup> T. Bárczy, <sup>29</sup> D. Barrado y Navascues, <sup>30</sup> S. C. C. Barros, <sup>20,31</sup>

W. Baumjohann, <sup>8</sup> M. Beck, <sup>11</sup> T. Beck, <sup>4</sup> N. Billot, <sup>11</sup> X. Bonfils, <sup>32</sup> C. Broeg, <sup>4,5</sup> L. A. Buchhave<sup>1,33</sup>

J. Cabrera, <sup>34</sup> S. Charnoz, <sup>35</sup> R. Cosentino, <sup>18</sup> Sz. Csizmadia<sup>34</sup> M. B. Davies, <sup>36</sup> M. Deleuil, <sup>10</sup>

L. Delrez<sup>22,37</sup> O. Demangeon, <sup>20,31</sup> B.-O. Demory, <sup>5</sup> D. Ehrenreich, <sup>11</sup> A. Erikson, <sup>34</sup>

E. Esparza-Borges, <sup>25,26</sup> H.-G. Florén, <sup>17,38</sup> A. Fortier, <sup>4,5</sup> M. Fridlund, <sup>39,40</sup> D. Futyan, <sup>11</sup> D. Gandolfi<sup>1,41</sup>

A. Ghedina, <sup>18</sup> M. Gillon, <sup>37</sup> M. Güdel, <sup>42</sup> P. Guterman, <sup>10,43</sup> A. Harutyunyan, <sup>18</sup> K. Heng, <sup>5,44</sup> K. G. Isaak, <sup>45</sup> J.

M. Jenkins, <sup>46</sup> L. Kiss, <sup>47,48</sup> J. Laskar, <sup>49</sup> A. Lecavelier des Etangs, <sup>50</sup> M. Lendl<sup>11</sup> C. Lovis, <sup>11</sup> D. Magrin, <sup>2</sup>

L. Marafatto, <sup>2</sup> A. F. Martinez Fiorenzano, <sup>18</sup> P. F. L. Maxted<sup>51</sup> M. Mayor, <sup>11</sup> G. Micela, <sup>52</sup> E. Molinari, <sup>53</sup>

F. Murgas, <sup>25</sup> N. Narita, <sup>25,54,55</sup> G. Olofsson, <sup>17</sup> R. Ottensamer, <sup>42</sup> I. Pagano, <sup>56</sup> A. Pasetti, <sup>57</sup> M. Pedani, <sup>18</sup>

F. A. Pepe, <sup>11</sup> G. Peter, <sup>58</sup> D. F. Phillips, <sup>14</sup> D. Pollacco, <sup>44</sup> D. Queloz, <sup>6,11</sup> R. Ragazzoni, <sup>1,2</sup> N. Rando, <sup>59</sup>

F. Ratti, <sup>59</sup> H. Rauer, <sup>34,60,61</sup> I. Ribas, <sup>27,28</sup> N. C. Santos, <sup>20,31</sup> D. Sasselov, <sup>14</sup> G. Scandariato, <sup>56</sup> S. Seager, <sup>12,62,63</sup>

D. Ségransan, <sup>11</sup> L. M. Serrano<sup>41</sup> A. E. Simon, <sup>4</sup> A. M. S. Smith, <sup>34</sup> M. Steinberger, <sup>8</sup> M. Steller, <sup>8</sup>

Gy. Szabó<sup>64,65</sup> N. Thomas, <sup>4</sup> J. D. Twicken, <sup>46,66</sup> S. Udry, <sup>11</sup> N. Walton<sup>67</sup> and J. N. Winn<sup>68</sup>

*Affiliations are listed at the end of the paper*

Accepted 2022 January 18. Received 2022 January 14; in original form 2021 November 25

## ABSTRACT

We present a precise characterization of the TOI-561 planetary system obtained by combining previously published data with TESS and CHEOPS photometry, and a new set of 62 HARPS-N radial velocities (RVs). Our joint analysis confirms the presence of four transiting planets, namely TOI-561 b ($P = 0.45$ d, $R = 1.42$ R<sub>⊕</sub>, $M = 2.0$ M<sub>⊕</sub>), c ($P = 10.78$ d, $R = 2.91$ R<sub>⊕</sub>, $M = 5.4$ M<sub>⊕</sub>), d ($P = 25.7$ d, $R = 2.82$ R<sub>⊕</sub>, $M = 13.2$ M<sub>⊕</sub>), and e ($P = 77$ d, $R = 2.55$ R<sub>⊕</sub>, $M = 12.6$ M<sub>⊕</sub>). Moreover, we identify an additional, long-period signal (>450 d) in the RVs, which could be due to either an external planetary companion or to stellar magnetic activity. The precise masses and radii obtained for the four planets allowed us to conduct interior structure and atmospheric escape modelling. TOI-561 b is confirmed to be the lowest density ($\rho_b = 3.8 \pm 0.5$ g cm<sup>-3</sup>) ultra-short period (USP) planet known to date, and the low metallicity of the host star makes it consistent with the general bulk density-stellar metallicity trend. According to our interior structure modelling, planet b has basically no gas envelope, and it could host a certain amount of water. In contrast, TOI-561 c, d, and e likely retained an H/He envelope, in addition to a possibly large water layer. The inferred planetary compositions suggest different atmospheric evolutionary paths, with planets b and c having experienced significant gas loss, and planets d and e showing an atmospheric content consistent with the original one. The uniqueness of the USP planet, the presence of the long-period planet TOI-561 e, and the complex architecture make this system an appealing target for follow-up studies.

**Key words:** techniques: photometric – techniques: radial velocities – planets and satellites: fundamental parameters – planets and satellites: interiors – stars: individual: TOI-561 (TIC 377064495, Gaia EDR3 3850421005290172416).
# 1 INTRODUCTION

Since the announcement of the first exoplanet orbiting a Sun-like star (Mayor & Queloz 1995), the growing number of discoveries in exoplanetary science have yielded a surprising variety of exoplanets and exoplanetary systems. The field has benefited hugely from dedicated space-based missions, such as CoRoT, Kepler, K2 (Baglin et al. 2006; Borucki et al. 2010; Howell et al. 2014), and recently TESS (Ricker et al. 2014). With more than 170 confirmed planets, and ~4000 planet candidates, the majority of which will likely turn out to be planets, TESS has increased the census of confirmed exoplanets to more than 4500.<fnref n="1" /> Alongside the aforementioned missions, which are designed to discover a large number of exoplanets by searching for transit-like signatures around hundreds of thousands of stars, new characterization missions, with a specific focus on the detailed study of known exoplanets, are now starting to operate. Among them, the Characterising Exoplanet Satellite (CHEOPS, Benz et al. 2021), launched on 18 December 2019, is a 30 cm telescope which is collecting ultra-high precision photometry of known exoplanets, aiming at their precise characterization. CHEOPS met its precision requirements both on bright and faint stars, achieving a noise level of ~15 ppm per 6 h intervals for V ~ 9 mag stars, and 75 ppm per 3 h for V ~ 12 mag stars (Benz et al. 2021). The importance of such a high photometric precision is reflected in CHEOPS' first scientific results, which span a variety of different fields (Lendl et al. 2020; Borsato et al. 2021; Delrez et al. 2021; Leleu et al. 2021; Maxted et al. 2021; Swayne et al. 2021; Szabó et al. 2021; Van Grootel et al. 2021; Morris et al. 2021a; Bonfanti et al. 2021b; Hooton et al. 2022; Barros et al. 2022; Deline et al. 2022; Wilson et al. 2022). As part of its main scientific goals, CHEOPS is refining the radii of known exoplanets to achieve the precision on the bulk density needed for internal structure and atmospheric evolution modelling. To fulfil this aim, radial velocity (RV) follow-ups using high-precision spectrographs are essential to provide the precise planetary masses that can be combined with radii measurements to determine accurate densities. Among the exoplanets having both radius and mass measurements, the ones in well-characterised multiplanetary systems are of particular interest, since they allow for investigation of their formation and evolution processes through comparative planetology, e.g. by comparing their individual inner bulk compositions (e.g. Guenther et al. 2017; Prieto-Arranz et al. 2018), by studying their mutual inclinations and eccentricities (e.g. Fabrycky et al. 2014; Van Eylen et al. 2019; Mills et al. 2019), and by investigating the correlations between their relative sizes, masses and orbital separations (e.g. Lissauer et al. 2011; Ciardi et al. 2013; Millholland, Wang & Laughlin 2017; Weiss et al. 2018; Adams et al. 2020; Jiang, Xie & Zhou 2020).

Within this context, TOI-561, announced simultaneously by Lacedelli et al. (2021) and Weiss et al. (2021) (L21 and W21 hereafter, respectively), is a particularly interesting system, both from the stellar (Section 2.1) and planetary (Section 2.2) perspective. The low stellar metallicity, the presence of an ultra-short period (USP) planet, where USP planets are meant here as planets with periods shorter than one day and radii smaller than 2 R⊕, and the complexity of its planetary configuration make TOI-561 an appealing target for in-depth investigations. In this study, we combine literature data with new TESS observations (Section 3.1), CHEOPS photometry (Section 3.2), and HARPS-N RVs (Section 3.3) to shed

light on the planetary architecture and infer the internal structure of the transiting planets. After assessing the planetary configuration using CHEOPS observations (Section 4.1) and performing a thorough analysis of the global RV data set (Section 4.2), we jointly modelled the photometric and spectroscopic data to obtain the planetary parameters (Section 5). We used our derived stellar and planetary properties to model the internal structures of the transiting planets (Section 6) and their atmospheric evolution (Section 7), before discussing our results and presenting our conclusions (Section 8).

# 2 THE TOI-561 SYSTEM

## 2.1 The host star

TOI-561 is an old, metal-poor, thick disc star (L21, W21), slightly smaller and cooler than the Sun, located ~84 pc away from the Solar System. We report the main astrophysical properties of the star in Table 1.

We adopted the spectroscopic parameters and stellar abundances from L21 (Table 1), which were derived exploiting the high SNR, high-resolution HARPS-N co-added spectrum (L21, Section 3.1) through an accurate analysis using three independent methods, namely the ARES+MOOG equivalent width method (Sousa 2014; Mortier et al. 2014), the Stellar Parameter Classification (Buchhave et al. 2012, 2014) and the CCFPAMS method (Malavolta et al. 2017).

Taking advantage of the updated parameters coming from the Gaia EDR3 release (Gaia Collaboration et al. 2021), we then used the L21 spectral parameters as priors on spectral energy distribution selection to infer the stellar radius ($R_*$) of TOI-561 using the infrared flux method (IRFM, Blackwell & Shallis 1977). The IRFM compares optical and infrared broad-band fluxes and synthetic photometry of stellar atmospheric models, and uses known relationships between stellar angular diameter, $T_{\text{eff}}$ and parallax to derive $R_*$, in a MCMC fashion as detailed in Schanche et al. (2020). For this study, we retrieved from the most recent data releases the Gaia G, $G_{\text{BP}}$, $G_{\text{RP}}$ (Gaia Collaboration et al. 2021), 2MASS J, H, K (Skrutskie et al. 2006), and WISE W1, W2 (Wright et al. 2010) broad-band photometric magnitudes, and we used the stellar atmospheric models from the ATLAS Catalogues (Castelli & Kurucz 2003) and the Gaia EDR3 parallax with the offset of Lindegren et al. (2021) applied, to obtain $R_* = 0.843 \pm 0.005 R_⊙$.

We combined two different sets of stellar evolutionary tracks and isochrones, PARSEC<fnref n="2" /> (PAdova & TRieste Stellar Evolutionary Code, v1.2S; Marigo et al. 2017) and CLES (Code Liègeois d'Évolution Stellaire, Scuflaire et al. 2008), to derive the stellar mass ($M_*$) and age ($t_*$) of TOI-561. As the star is significantly alpha-enhanced, we avoided using [Fe/H] as a proxy for the stellar metallicity; instead, we inserted both [Fe/H] and [α/Fe] in relation (3) provided by Yi et al. (2001), obtaining an overall scaling of metal abundances [M/H] = −0.23 ± 0.06. Besides [M/H], the main input parameters for computing $M_*$ and $t_*$ were $T_{\text{eff}}$ and $R_*$. In addition, we used as inputs $\log R'_{\text{HK}}$ and the upper limit on $v \sin i$ from L21, and the yttrium over magnesium abundance [Y/Mg] = −0.22 ± 0.07, as computed from [Mg/H] and [Y/H] reported by W21. These indices improve the model convergence by discarding unlikely young isochrones, as broadly discussed in section 2.2.3 of Bonfanti & Gillon (2020), and references therein. The PARSEC results were obtained using the isochrone placement algorithm of

<footnote n="1">¹From NASA Exoplanet Archive (Akeson et al. 2013), https://exoplanetarc.hive.ipac.caltech.edu/.</footnote>

<footnote n="2">²http://stev.oapd.inaf.it/cgi-bin/cmd</footnote>
<table><caption>**Table 1.** Stellar properties of TOI-561.</caption><thead><tr><th></th><th>TOI-561</th><th></th></tr></thead><tbody><tr><td>TIC</td><td>377064495</td><td></td></tr><tr><td>Gaia EDR3</td><td>3850421005290172416</td><td></td></tr><tr><td>2MASS</td><td>J09524454+0612589</td><td></td></tr><tr><td>Parameter</td><td>Value</td><td>Source</td></tr><tr><td>RA (J2016; hh:mm:ss.ss)</td><td>09:52:44.43</td><td>A</td></tr><tr><td>Dec (J2016; dd:mm:ss.ss)</td><td>+06:12:57.94</td><td>A</td></tr><tr><td>μ<sub>α</sub> (mas yr<sup>-1</sup>)</td><td>-108.504 ± 0.022</td><td>A</td></tr><tr><td>μ<sub>δ</sub> (mas yr<sup>-1</sup>)</td><td>-61.279 ± 0.019</td><td>A</td></tr><tr><td>γ (km s<sup>-1</sup>)</td><td>79.54 ± 0.56</td><td>A</td></tr><tr><td>Parallax (mas)</td><td>11.8342 ± 0.0208</td><td>A</td></tr><tr><td>Distance (pc)</td><td>84.25 ± 0.12</td><td>B</td></tr><tr><td>TESS (mag)</td><td>9.527 ± 0.006</td><td>C</td></tr><tr><td>G (mag)</td><td>10.0181 ± 0.0028</td><td>A</td></tr><tr><td>G<sub>BP</sub> (mag)</td><td>10.3945 ± 0.0028</td><td>A</td></tr><tr><td>G<sub>RP</sub> (mag)</td><td>9.4692 ± 0.0038</td><td>A</td></tr><tr><td>V (mag)</td><td>10.252 ± 0.006</td><td>C</td></tr><tr><td>B (mag)</td><td>10.965 ± 0.082</td><td>C</td></tr><tr><td>J (mag)</td><td>8.879 ± 0.020</td><td>D</td></tr><tr><td>H (mag)</td><td>8.504 ± 0.055</td><td>D</td></tr><tr><td>K (mag)</td><td>8.394 ± 0.019</td><td>D</td></tr><tr><td>W1 (mag)</td><td>8.337 ± 0.023</td><td>E</td></tr><tr><td>W2 (mag)</td><td>8.396 ± 0.020</td><td>E</td></tr><tr><td>T<sub>eff</sub> (K)</td><td>5372 ± 70</td><td>F</td></tr><tr><td>log g (cgs)</td><td>4.50 ± 0.12</td><td>F</td></tr><tr><td>[Fe/H] (dex)</td><td>-0.40 ± 0.05</td><td>F</td></tr><tr><td>[Mg/H] (dex)</td><td>-0.17 ± 0.05</td><td>F</td></tr><tr><td>[Si/H] (dex)</td><td>-0.22 ± 0.05</td><td>F</td></tr><tr><td>[Ti/H] (dex)</td><td>-0.12 ± 0.03</td><td>F</td></tr><tr><td>[α/Fe] (dex)</td><td>0.23 ± 0.04</td><td>F</td></tr><tr><td>[M/H] (dex)</td><td>-0.23 ± 0.06</td><td>G</td></tr><tr><td>[Y/Mg] (dex)</td><td>-0.22 ± 0.07</td><td>G<sup>a</sup></td></tr><tr><td>log R<sub>HK</sub></td><td>-5.003 ± 0.012</td><td>F</td></tr><tr><td>v sin i (km s<sup>-1</sup>)</td><td>&lt;2</td><td>F</td></tr><tr><td>R<sub>*</sub> (R<sub>☉</sub>)</td><td>0.843 ± 0.005</td><td>G, IRFM</td></tr><tr><td>M<sub>*</sub> (M<sub>☉</sub>)</td><td>0.806 ± 0.036</td><td>G, isochrones</td></tr><tr><td>t<sub>*</sub> (Gyr)</td><td>11.0<sup>+2.8</sup><sub>-3.5</sub></td><td>G, isochrones</td></tr><tr><td>ρ<sub>*</sub> (ρ<sub>☉</sub>)</td><td>1.34 ± 0.06</td><td>G, from R<sub>*</sub> and M<sub>*</sub></td></tr><tr><td>ρ<sub>*</sub> (g cm<sup>-3</sup>)</td><td>1.89 ± 0.09</td><td>G, from R<sub>*</sub> and M<sub>*</sub></td></tr><tr><td>L<sub>*</sub> (L<sub>☉</sub>)</td><td>0.533 ± 0.029</td><td>G, from R<sub>*</sub> and T<sub>eff</sub></td></tr><tr><td>Spectral type</td><td>G9V</td><td>F</td></tr></tbody></table>

Notes. (A) Gaia EDR3 (Gaia Collaboration et al. 2021). (B) Bailer-Jones et al. (2021). (C) TESS Input Catalogue Version 8 (TICv8, Stassun et al. 2018). (D) Two Micron All Sky Survey (2MASS; Cutri et al. 2003). (E) Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010). (F) L21. (G) This work.

Bonfanti et al. (2015), Bonfanti, Ortolani & Nascimbeni (2016), which retrieves the best-fit parameters by interpolating within a pre-computed grid of models, while the CLES algorithm models directly the star through a Levenberg-Marquardt minimization (Salmon et al. 2021). The final adopted values ($M_{*} = 0.806 \pm 0.036 M_{⊙}$, $t_{*} = 11.0^{+2.8}_{-3.5}$ Gyr) are a combination of the outputs from both sets of models, as described in detail in Bonfanti et al. (2021b). The derived mass and radius, listed in Table 1, are consistent within $1\sigma$ with the values reported in L21 ($R_{*} = 0.849 \pm 0.007 R_{⊙}$, $M_{*} = 0.785 \pm 0.018 M_{⊙}$).

## 2.2 The planetary system

The discovery of a multiplanetary system orbiting TOI-561 was announced simultaneously by L21 and W21 in two independent papers. The main planetary parameters from both studies are reported in Table 2.

<table><caption>**Table 2.** Literature parameters of the proposed planets orbiting TOI-561.</caption><thead><tr><th>TOI-561 b</th><th>Lacedelli et al. (2021)</th><th>Weiss et al. (2021)</th></tr></thead><tbody><tr><td>P (d)</td><td>0.446578 ± 0.000017</td><td>0.446573<sup>+0.000032</sup><sub>-0.000021</sub></td></tr><tr><td>T<sub>0</sub> (TBJD)</td><td>1517.498 ± 0.001</td><td>1517.4973 ± 0.0018</td></tr><tr><td>R<sub>p</sub> (R<sub>⊕</sub>)</td><td>1.423 ± 0.066</td><td>1.45 ± 0.11</td></tr><tr><td>K (m s<sup>-1</sup>)</td><td>1.56 ± 0.35</td><td>3.1 ± 0.8</td></tr><tr><td>M<sub>p</sub> (M<sub>⊕</sub>)</td><td>1.59 ± 0.36</td><td>3.2 ± 0.8</td></tr><tr><td>TOI-561 c</td><td></td><td></td></tr><tr><td>P (d)</td><td>10.779 ± 0.004</td><td>10.77892 ± 0.00015</td></tr><tr><td>T<sub>0</sub> (TBJD)</td><td>1527.060 ± 0.004</td><td>1527.05825 ± 0.00053</td></tr><tr><td>R<sub>p</sub> (R<sub>⊕</sub>)</td><td>2.878 ± 0.096</td><td>2.90 ± 0.13</td></tr><tr><td>K (m s<sup>-1</sup>)</td><td>1.84 ± 0.33</td><td>2.4 ± 0.8</td></tr><tr><td>M<sub>p</sub> (M<sub>⊕</sub>)</td><td>5.40 ± 0.98</td><td>7.0 ± 2.3</td></tr><tr><td>TOI-561 d</td><td></td><td></td></tr><tr><td>P (d)</td><td>25.62 ± 0.04</td><td>—</td></tr><tr><td>T<sub>0</sub> (TBJD)</td><td>1521.882 ± 0.004</td><td>—</td></tr><tr><td>R<sub>p</sub> (R<sub>⊕</sub>)</td><td>2.53 ± 0.13</td><td>—</td></tr><tr><td>K (m s<sup>-1</sup>)</td><td>3.06 ± 0.33</td><td>—</td></tr><tr><td>M<sub>p</sub> (M<sub>⊕</sub>)</td><td>11.95 ± 1.28</td><td>—</td></tr><tr><td>TOI-561 e</td><td></td><td></td></tr><tr><td>P (d)</td><td>77.23 ± 0.39</td><td>—</td></tr><tr><td>T<sub>0</sub> (TBJD)</td><td>8538.181 ± 0.004</td><td>—</td></tr><tr><td>R<sub>p</sub> (R<sub>⊕</sub>)</td><td>2.67 ± 0.11</td><td>—</td></tr><tr><td>K (m s<sup>-1</sup>)</td><td>2.84 ± 0.41</td><td>—</td></tr><tr><td>M<sub>p</sub> (M<sub>⊕</sub>)</td><td>16.0 ± 2.3</td><td>—</td></tr><tr><td>TOI-561 f<sup>a</sup></td><td></td><td></td></tr><tr><td>P (d)</td><td>—</td><td>16.287 ± 0.005</td></tr><tr><td>T<sub>0</sub> (TBJD)</td><td>—</td><td>1521.8828 ± 0.0035</td></tr><tr><td>R<sub>p</sub> (R<sub>⊕</sub>)</td><td>—</td><td>2.32 ± 0.16</td></tr><tr><td>K (m s<sup>-1</sup>)</td><td>—</td><td>0.9 ± 0.6</td></tr><tr><td>M<sub>p</sub> (M<sub>⊕</sub>)</td><td>—</td><td>3.0<sup>+2.4</sup><sub>-1.9</sub></td></tr><tr><td>M<sub>*</sub> (M<sub>☉</sub>)</td><td>0.785 ± 0.018</td><td>0.805 ± 0.030</td></tr></tbody></table>

The two papers presented different RV data sets, collected with HARPS-N and HIRES, respectively, to confirm the planetary nature of three candidates identified by TESS in sector 8, the only available sector at the time of the publications. The TESS-identified signals had periods of ~0.45, ~10.8, and ~16 d. The two inner candidates were confirmed by both L21 and W21, with the names of TOI-561 b (an USP super-Earth, with period $P_{b} \sim 0.4465$ d, and radius $R_{b} \sim 1.4 R_{\oplus}$) and TOI-561 c (a warm mini-Neptune, with $P_{c} \sim 10.779$ d, and $R_{c} \sim 2.9 R_{\oplus}$). However, two different interpretations for the third TESS signal were proposed by the authors. In the scenario presented in L21, the two transits related to the third TESS signal were interpreted as single transits of two distinct planets, TOI-561 d ($P_{d} \sim 25.6$ d, $R_{d} \sim 2.5 R_{\oplus}$), and TOI-561 e ($P_{e} \sim 77$ d, $R_{e} \sim 2.7 R_{\oplus}$). The periods of these two planets were inferred from the RV analysis, which played an essential role in determining the final planetary architecture. In fact, the ephemeris match between the RV and photometric fits (see fig. 5 of L21) and the non-detection of the 16 d signal in the HARPS-N data set, combined with the different durations of the two TESS transits and results from the long-term stability analysis led the authors to converge on a four-planet configuration, presenting robust mass and radius detection for all the four planets in the system (L21, table 5). In contrast, W21 proposed the presence of a single planet at the period suggested by TESS (TOI-561 f, $P_{f} \sim 16.29$ d, $R_{f} \sim 2.3 R_{\oplus}$), based on the analysis of the two available transits. W21 pointed out that the 8.1 d alias of planet f's orbital period is also consistent with the TESS data, with

<footnote><sup>a</sup> Based on W21 abundances.</footnote>

<footnote><sup>a</sup>Referred as TOI-561 d in W21.</footnote>
the even transit falling into the *TESS* download gap, even though in this case the transit duration would be too long compared to what is expected for a 8 d period planet (Section 4.9, W21). However, the authors could not obtain an accurate mass determination for this planet, with the 60 HIRES RVs being consistent with a non-detection (Section 7.2, W21). An additional discrepancy between the two studies is the mass of the USP planet, differing by almost a factor of 2. According to the W21 analysis, TOI-561 b has a mass of 3.2 ± 0.8 M<sub>⊕</sub>, making it consistent with a rocky composition and placing it among the population of typical small (<2 R<sub>⊕</sub>), extremely irradiated USP planets (Sanchis-Ojeda et al. 2015; Dai et al. 2021). Instead, assuming the low mass (M<sub>b</sub> = 1.59 ± 0.36 M<sub>⊕</sub>) inferred from L21 analysis, TOI-561 b is not consistent with a pure rocky composition, and it is the lowest density USP super-Earth known to date, calling for a more complex interpretation (e.g. lighter core composition, deep water reservoirs, presence of a high-metallicity, volatile materials or water steam envelope, etc.).

The complexity of this system and the differences between the two studies demanded further investigations. We therefore decided to collect additional, precise photometric and RV data (Section 3) to shed light on the planetary configuration and on the internal composition of the TOI-561 planets.

# 3 OBSERVATIONS

## 3.1 TESS photometry

During its two-year primary mission (Ricker et al. 2014), *TESS* observed TOI-561 in two-minute cadence mode between 2019 February 2 and February 27 (sector 8). After entering its extended mission, *TESS* re-observed the star in 2-min cadence mode during sector 35, between 2021 February 9 and March 6. At the beginning of the second orbit, the spacecraft dropped out of Fine Pointing mode for 3.44 d, entering Coarse Pointing mode.<fnref n="3" /> Data collected during Coarse Pointing mode were flagged and removed from the Pre-search Data Conditioning Simple Aperture Photometry (PDCSAP; Smith et al. 2012; Stumpe et al. 2012, 2014) light curves, leading to a total of 19.86 d of science data. The photometric observations of TOI-561 were reduced by the Science Processing Operations Center (SPOC) pipeline and searched for evidence of transiting planets (Jenkins et al. 2016; Jenkins 2020). For our photometric analysis, we used the light curves based on the PDCSAP, downloading the 2-min cadence data from the Mikulski Archive for Space Telescopes (MAST),<fnref n="4" /> and removing all the observations encoded as *NaN* or flagged as bad-quality (DQUALITY>0) points by the SPOC pipeline.<fnref n="5" /> We performed outlier rejection by doing a cut at 3σ for positive outliers and 5σ (i.e. larger than the deepest transit) for negative outliers. The resulting *TESS* light curves of sectors 8 and 35 are shown in Fig. 1, and Table 3 summarizes the total number of transits observed by *TESS* for each planet.

To refine the ephemeris of planet d in time for the scheduling of the *CHEOPS* observations (Section 4.1), we also extracted the 10-min cadence light curve of sector 35 using the quick-look *TESS* Full

Frame Images (FFIs) calibrated using the *TESS* Image CAlibrator<fnref n="6" /> package (TICA, Fausnaugh et al. 2020).

## 3.2 CHEOPS photometry

To confirm the planetary architecture and improve the planetary parameters, we obtained three visits of TOI-561 with *CHEOPS*, the ESA small class mission dedicated to the characterization of known exoplanets (Benz et al. 2021). The observations, collected within the Guaranteed Time Observing (GTO) programme, were carried out between 2021 January 23 and April 15, for a total of 73.85 h on target. During the three visits, we observed a total of eight transits of TOI-561 b, two transits of TOI-561 c, and one transit of TOI-561 d. The three *CHEOPS* light curves have an observing efficiency, i.e. the actual time spent observing the target with respect to the total visit duration, of 64 per cent, 75 per cent, and 61 per cent, respectively. The observing efficiency is linked to data gaps, which are intrinsically present in all *CHEOPS* light curves (see e.g. Bonfanti et al. 2021b; Delrez et al. 2021; Leleu et al. 2021), and are related to *CHEOPS*'s low-Earth orbit. In fact, during (1) South Atlantic Anomaly (SAA) crossing, (2) target occultation by the Earth, and (3) too high stray light contamination, no data are downlinked. This results in data gaps, whose number and extension depend on the target sky position (Benz et al. 2021). For all the visits, we adopted an exposure time of 60 s. The summary log of the *CHEOPS* observations is reported in Table 4.

Data were reduced using the latest version of the *CHEOPS* automatic Data Reduction Pipeline (DRP v13; Hoyer et al. 2020), which performs aperture photometry of the target after calibrating the raw images (event flagging, bias, gain, non-linearity, dark current, and flat-field) and correcting them for instrumental and environmental effects (smearing trails, cosmic rays, de-pointing, stray light, and background). The target flux is obtained for a set of three fixed-radius apertures, namely *R* = 22.5 arcsec (RINF), 25.0 arcsec (DEFAULT), 30.0 arcsec (RSUP), plus an additional one specifically computed to optimize the radius based on the instrumental noise and contamination level of each target (OPTIMAL). Moreover, the DRP estimates the contamination in the photometric aperture due to nearby targets using the sources listed in the Gaia DR2 catalogue (Gaia Collaboration et al. 2018) to simulate the *CHEOPS* Field-of-View (FoV) of the target, as described in detail in Hoyer et al. (2020). No strong contaminants are present in the TOI-561 FoV, and the main contribution to the contamination is due to the smearing trails of a *G* = 10.20 mag star at a projected sky distance of ~117.9 arcsec, which rotates around the target inside the CCD window because of the *CHEOPS* field rotation (Benz et al. 2021). During the third visit three telegraphic pixels (pixels with a non-stable and abnormal behaviour during the visit) appeared within the *CHEOPS* aperture, one of them inside the *CHEOPS* PSF (Fig. 2). A careful treatment, described in detail in Appendix A, was applied to correct for their effect. In the subsequent analysis, we adopted for all the visits the RINF photometry (see Fig. A1 in Appendix A), which minimized the light-curve root mean square (RMS) dispersion, and we removed the outliers by applying a 4σ clipping.

Finally, a variety of non-astrophysical sources, such as varying background, nearby contaminants or others, can produce short-term photometric trends in the *CHEOPS* light curves on the timescale of one orbit, due to the rotation of the *CHEOPS* FoV around the target and due to the nature of the spacecraft orbit. To correct for these effects, we detrended the light curves using the basis vectors

<footnote n="3">³See *TESS* Data Release Notes: Sector 35, DR51 (https://archive.stsci.edu/tess/tess_drn.html).</footnote>

<footnote n="4">⁴https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html</footnote>

<footnote n="5">⁵https://archive.stsci.edu/missions/tess/doc/EXP-TESS-ARC-ICD-TM-0014-Rev-F.pdf</footnote>

<footnote n="6">⁶https://archive.stsci.edu/hlsp/tica</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** TESS sector 8 (left) and 35 (right) PDCSAP light curves of TOI-561. In the top panel, the dark red solid line shows the best-fitting transit and Matérn-3/2 kernel Gaussian Process (GP) model, as detailed in Section 5. The central panel shows the flattened light curve after the removal of the GP component, with the best-fitting transit model superimposed (dark red solid line). The transits of planets c, d, and e are labelled and highlighted with orange, green, and red vertical lines, respectively. The expected locations of the transits of planets c and d occurring during the data gaps of sector 35 are marked with pale, dashed orange and green lines, respectively. Planet e is not expected to transit in sector 35. The transits of the USP planet are too shallow to be individually visible and are not indicated. Light-curve residuals are shown in the bottom panel.</figcaption></figure>

<table><caption>**Table 3.** Number of TOI-561 transits observed by TESS.</caption><thead><tr><th></th><th>TOI-561 b</th><th>TOI-561 c</th><th>TOI-561 d</th><th>TOI-561 e</th></tr></thead><tbody><tr><td>Sector 8</td><td>41</td><td>2</td><td>1</td><td>1</td></tr><tr><td>Sector 35</td><td>43</td><td>1</td><td>-</td><td>-</td></tr></tbody></table>

provided by the DRP, as detailed in Section 5. The resulting detrended light curves are shown in Fig. 3.

## 3.3 HARPS-N spectroscopy

In addition to the 82 RVs published in L21, we collected 62 high-resolution spectra using HARPS-N at the Telescopio Nazionale Galileo (TNG) in La Palma (Cosentino et al. 2012, 2014). These were used to refine the planetary masses and confirm the system configuration. The new observations were collected between 2020 November 15 and 2021 June 1. Following the same strategy of the previous season (L21), in addition to 30 single observations, we collected six points per night on 2021 February 8 and 10, and two points per night on ten additional nights, specifically targeting the USP planet. The exposure time for all the observations was set to 1800 s, resulting in a signal-to-noise ratio (SNR) at 550 nm of 83 ± 20 (median ± standard deviation) and a RV measurement uncertainty of 1.0 ± 0.4 m s⁻¹. All the observations were gathered with the second

HARPS-N fibre illuminated by the Fabry-Perot calibration lamp to correct for the instrumental RV drift.

We reduced the global HARPS-N data set (144 RVs in total) using the new version of the HARPS-N Data Reduction Software based on the ESPRESSO pipeline (DRS, version 2.3.1; see Dumusque et al. 2021 for more details). We used a G2 flux template to correct for variations in the flux distribution as a function of wavelength, and a G2 binary mask to compute the cross-correlation function (CCF; Baranne et al. 1996; Pepe et al. 2002). We report the RVs and the associated activity indices (see Section 4.2) with their 1σ uncertainties in Table 5. As in L21, we removed from the first season data set five RVs with associated errors >2.5 m s⁻¹ from spectra with SNR <35 (see Appendix B in L21). All the RV uncertainties of the second season data set were below 2.5 m s⁻¹, so no points were removed.

## 3.4 HIRES spectroscopy

We included in our analysis 60 high-resolution spectra collected with the W.M. Keck Observatory HIRES instrument on Mauna Kea, Hawaii between 2019 May and 2020 October. The data set was published in W21, and we refer to that paper for details regarding the observing and data reduction procedures. The HIRES data set has an RMS of 5 m s⁻¹, and a median individual RV uncertainty of 1.4 m s⁻¹ (W21).
<table><caption>**Table 4.** Log of TOI-561 CHEOPS observations.</caption><thead><tr><th>Visit (No.)</th><th>File key</th><th>Starting date (UTC)</th><th>Duration (h)</th><th>Data points (No.)</th><th>Efficiency (Per cent)</th><th>Exposure time (s)</th><th>Planets</th></tr></thead><tbody><tr><td>1</td><td>CH_PR100031_TG037001_V0200</td><td>2021-01-23T15:29:07</td><td>15.67</td><td>604</td><td>64</td><td>60</td><td>b,c</td></tr><tr><td>2</td><td>CH_PR100008_TG000811_V0200</td><td>2021-03-29T10:19:08</td><td>4.42</td><td>207</td><td>75</td><td>60</td><td>b,c</td></tr><tr><td>3</td><td>CH_PR100031_TG039301_V0200</td><td>2021-04-12T23:52:28</td><td>53.76</td><td>1978</td><td>61</td><td>60</td><td>b,d</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Extraction of 60 × 60 arcsec of the CHEOPS FoV during the third visit centred on TOI-561. The dashed black circle represents the RINF photometric aperture surrounding the CHEOPS PSF, whose centroid is marked by the black cross. The positions of the three identified telegraphic pixels, including the one located within the CHEOPS PSF (see Appendix A), are highlighted by the red, circled crosses.</figcaption></figure>

# 4 PROBING THE SYSTEM ARCHITECTURE

## 4.1 CHEOPS confirmation of TOI-561 d

To solve the discrepancy among the planetary architectures proposed by L21 and W21 (Section 2.2), we initially looked for the transits of TOI-561 d (~25 d) and TOI-561 f (~16 d) in the TESS sector 35 light curve, whereas TOI-561 e (~77 d) was not expected to transit during those TESS observations. However, as shown in the top panel of Fig. 4, the transits of planet d and f occurred during the light curve gap (Section 3.1), and so we could not use the new TESS data to conclusively discriminate between the two planetary configurations. None the less, these observations ruled out the planet f alias at ~8.1 d mentioned in W21, since no transit events were detected at its predicted transit times.

We therefore decided to probe the L21 scenario collecting a transit of TOI-561 d using CHEOPS. We opted for the scheduling of the last seasonal observing window, in April, in order to take advantage of the most updated ephemeris to optimize the scheduling. For this reason, we performed a global fit adding to the literature data a partial set of the new HARPS-N RVs, as of 2021 March 16, and including the TESS sector 35 light curve extracted from the second data release of the TICA FFIs in 2021 March. Even if no transit was detected, the new TESS sector helped to reduce the time window in which to search. In fact, the TESS data partially covered the 3σ-uncertainty transit window, enabling us to exclude some time-spans in the computation of the CHEOPS visit. Thanks to the ephemeris update, the CHEOPS 3σ observing window shrank from ~7.4 to ~2.2 d, demonstrating the importance of the early TESS data releases

in the scheduling of follow-up observations. The bottom panel of Fig. 4 shows the CHEOPS visit scheduled to observe TOI-561 d, whose transit occurred almost exactly at the predicted time, so confirming the planetary period and giving further credence to the four-planet scenario proposed by L21.

Even updating the ephemeris using the partial new HARPS-N data set, the last possible CHEOPS observing window of TOI-561 e in the 2021 season was still longer than seven days because of ephemeris uncertainties. Even including the full set of RVs would have not helped as the target was no longer observable with CHEOPS when the HARPS-N campaign finished. Given the high pressure on the CHEOPS schedule, we therefore plan the TOI-561 e observations for the 2022 observing season. The ephemeris for the 2022 CHEOPS observations will be updated using the TESS Sectors 45 and 46 observations in November–December 2021, and the results will be presented in a future publication.

## 4.2 Additional signals in the RV data

Before proceeding with the global modelling, we analysed the RV data sets in order to confirm the robustness of the L21 scenario and search for potential new signals. The $l_1$-periodogram<fnref n="7" /> (Hara et al. 2017) of the combined HARPS-N and HIRES RVs (Fig. 5) shows four significant peaks corresponding to the planetary periods reported in L21, plus hints of a possible longer period signal with a broad peak around 400–600 d. We investigated the presence of this additional signal in a Bayesian framework using PYORBIT<fnref n="8" /> (Malavolta et al. 2016, 2018), a package for light curve and RVs analysis. We employed the DYNESTY nested-sampling algorithm (Skilling 2004, 2006; Speagle 2020), assuming 1000 live points, and including offset and jitter terms for each data set. We first performed a four-planet fit of the combined data sets, using the L21 values to impose Gaussian priors on periods and transit times,<fnref n="9" /> and assuming eccentric orbits with a half-Gaussian zero-mean prior on the eccentricity (with variance 0.098; Van Eylen et al. 2019), except for the circular orbit of the USP planet. We let the semi-amplitude K vary between 0.01 and 100 m s⁻¹. As can be seen in Fig. 6, the RV residuals show an anomalous positive variation at ~9000 BJD-2450000, and the Generalized Lomb-Scargle (GLS, Zechmeister & Kürster 2009) periodogram of the RV residuals revealed the presence of a significant, broad peak at low frequencies. Moreover, the HARPS-N jitter was 1.84 m s⁻¹, which is unusually high when compared to the value reported in L21 ($\sigma_{HARPS-N} = 1.29 \pm 0.23$ m s⁻¹). We therefore performed a second fit including a fifth Keplerian signal, allowing the period to span between 2 and 900 d. According to the Bayesian Evidence, this model is strongly favoured with respect to the four-planet model, with a difference in

<footnote n="7">⁷https://github.com/nathanchara/l1periodogram.</footnote>

<footnote n="8">⁸https://github.com/LucaMalavolta/PyORBIT, V8.1.</footnote>

<footnote n="9">⁹We note that we obtained the same results when using uniform, uninformative priors, also for the five- and six-Keplerian fits.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** CHEOPS detrended light curves of TOI-561. Visits 1, 2, and 3 are shown in the top left, top right, and bottom panel, respectively. The best-fitting model is over-plotted as a red solid line, and residuals are shown for each visit. The transits of planets b, c, and d are highlighted with blue, orange, and green triangles, respectively.</figcaption></figure>

<table><caption>**Table 5.** HARPS-N RV and activity indices measurements.</caption><thead><tr><th>BJD<sub>TDB</sub><br>(d)</th><th>RV<br>(m s<sup>-1</sup>)</th><th>σ<sub>RV</sub><br>(m s<sup>-1</sup>)</th><th>FWHM<br>(km s<sup>-1</sup>)</th><th>σ<sub>FWHM</sub><br>(km s<sup>-1</sup>)</th><th>BIS<br>(m s<sup>-1</sup>)</th><th>σ<sub>BIS</sub><br>(m s<sup>-1</sup>)</th><th>Contrast</th><th>σ<sub>contr</sub></th><th>S-index</th><th>σ<sub>S</sub></th><th>Hα<br>(dex)</th><th>σ<sub>Hα</sub><br>(dex)</th></tr></thead><tbody><tr><td>2458804.70780</td><td>79695.97</td><td>1.13</td><td>6.415</td><td>0.002</td><td>-86.82</td><td>2.26</td><td>59.879</td><td>0.021</td><td>0.1643</td><td>0.0005</td><td>0.2101</td><td>0.0002</td></tr><tr><td>2458805.77552</td><td>79699.66</td><td>0.85</td><td>6.419</td><td>0.002</td><td>-85.13</td><td>1.71</td><td>59.810</td><td>0.016</td><td>0.1702</td><td>0.0003</td><td>0.2124</td><td>0.0001</td></tr><tr><td>2458806.76769</td><td>79697.50</td><td>0.91</td><td>6.415</td><td>0.002</td><td>-83.66</td><td>1.82</td><td>59.861</td><td>0.017</td><td>0.1689</td><td>0.0003</td><td>0.2082</td><td>0.0001</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

Note. This table is available in its entirety in machine-readable form.

the logarithmic evidences $\Delta \ln Z = 19.0$ (Kass & Raftery 1995).<fnref n="10" /> Moreover, the HARPS-N jitter decreased to ~1.37 m s⁻¹. After this fit, the periodogram of the residuals did not show evidence of additional significant peaks (Fig. 6). This is confirmed also by the comparison with a six-Keplerian model that we tested, with the period of the sixth Keplerian free to span between 2 and 900 d, whose Bayesian Evidence differed by less than two from the five-Keplerian model one, indicating that there was no strong evidence to favour a more complex model (Kass & Raftery 1995).

The fitted period of the fifth Keplerian was ~480 d. Such a long-term signal could be induced either by stellar activity, considering that stellar magnetic fields related to magnetic cycles can show variability on timescales of the order of 1–3 yr (e.g. Collier Cameron 2018; Hatzes 2019; Crass et al. 2021), or by an additional long-period

planet. We refer here to an eventual long-period planet because, given the inferred semi-amplitude of ~2 m s⁻¹ (Table 6), an external companion with mass equal to 13 M<sub>j</sub> (assuming this value to be the threshold between planetary and sub-stellar objects) would have an inclination of ~0.03 deg. Such an inclination would imply an almost perpendicular orbit with respect to the orbital plane of the four inner planets, hinting at a very unlikely configuration. Therefore, in the hypothesis of the presence of an external companion, it would most likely be a planetary-mass object.

On one hand, all the five signals, including the long-term one, are recovered in an independent analysis that we performed with the CCF-based SCALPELS algorithm (Collier Cameron et al. 2021). Concisely, SCALPELS projects the RV time-series on to the highest variance principal components of the time series of autocorrelation functions of the CCF, with the aim of distinguishing RV variations caused by orbiting planets from activity-induced distortions on each CCF. The absence of the signal in the SCALPELS shape-driven velocities indicates that the long-term periodicity is not due to shape

<footnote n="10">¹⁰According to Kass & Raftery (1995), a difference $\Delta \ln Z > 5$ sets a strong evidence against the null hypothesis, which in our case corresponds to the four-planet model.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Top: 2-min cadence detrended TESS light curve of sector 35. The predicted transit times of TOI-561 c and d (according to L21 ephemeris), and TOI-561 f (according to W21 ephemeris), are highlighted with orange, green and black vertical solid lines, respectively. The black dashed lines indicate the predicted position of planet f alias at ~8.1 d. The only transit present in the light curve is the one of TOI-561 c at ~9260 BJD−2450000, while no transit events occurred at the predicted times of planet f alias. The transits of planet d and f fall into the time-series gap. Bottom: CHEOPS visit scheduled to observe TOI-561 d. The green vertical solid line indicates the predicted transit time used to compute the CHEOPS observing window after the ephemeris update (Section 4.1). The transit occurred within the 68 per cent highest probability density interval, highlighted by the pale green region. We note that this transit is not consistent with the ephemeris propagation of planet f, which would have transited at 9319.94 BJD−2450000, so almost one day after the observed CHEOPS transit.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** $\ell_1$-periodogram of the combined HARPS-N and HIRES data sets, computed on a grid of frequencies from 0 to 2.5 cycles per day. The total time-span of the observations is 768 d. The code automatically accounts for the offset between HARPS-N and HIRES data by using the mean value of each data set.</figcaption></figure>

changes in the line profiles, supporting the idea of a planetary origin. Moreover, TOI-561 is not expected to be a particularly active star given its old age and low $\log R_{\text{HK}}^*$, as assessed in the L21 and W21 activity analyses. As can be seen in Fig. 7, the GLS periodogram of the majority of the activity indicators extracted with the HARPS-N DRS, i.e. full width at half-maximum (FWHM), bisector span (BIS), contrast and $H_\alpha$, do not show significant peaks, with none of them exceeding the 0.1 False Alarm Probability (FAP) threshold,

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Time-series and GLS periodogram of the RV residuals after the four-planet and five-Keplerian fits as described in Section 4.2. In the residuals plot, the HARPS-N and HIRES RVs are plotted with dark blue diamonds and light blue circles, respectively. In the periodogram plots, the dashed and dotted horizontal lines show the 1 and 0.1 per cent False Alarm Probability (FAP) level, respectively. The red vertical line indicates the main peak of each periodogram. The long-period peak around frequencies 0.0017–0.0025 $d^{-1}$ ($P = 400–600$ d) in the four-planet residuals periodogram is modelled by the fifth Keplerian, and no more significant peaks are identified in the five-Keplerian residuals periodogram. Moreover, the positive variation at ~9000 BJD-2450000 in the four-planet fit residuals disappears in the five-Keplerian fit residuals.</figcaption></figure>

which we computed using a bootstrap approach, at the frequency of interest. On the other hand, the periodogram of the S-index, which is particularly sensitive to magnetically induced activity, shows a significant, broad peak at low frequencies, potentially suggesting that the previously identified long-term variability is related to stellar activity. Considering this, we performed an additional DYNESTY fit assuming a four-planet model and including a Gaussian Process (GP) regression with a quasi-periodic kernel, as formulated in Grunblatt, Howard & Haywood (2015), to account for the long-term signal. We modelled simultaneously the RVs and the S-index time series in order to better inform the GP (Langellier et al. 2021; Osborn et al. 2021), using two independent covariance matrices for each data set with common GP hyper-parameters except for the amplitude of the covariance matrix, assuming uniform, non-informative priors on all of them. The fit suggests a periodicity longer than ~570 d, but the GP model is too flexible to derive a precise period value, considering also that the global RV baseline (~768 d) is comparable with the periodicity of the long-term signal. The inferred semi-amplitudes of the four known planets differed by less than 0.07$\sigma$.
<table><caption>**Table 6.** Parameters of the TOI-561 system, including the fifth Keplerian signal, as determined from the joint photometric and spectroscopic fit described in Section 5.</caption><thead><tr><th></th><th>TOI-561 b</th><th>TOI-561 c</th><th>TOI-561 d</th><th>TOI-561 e</th><th>5<sup>th</sup> Keplerian</th></tr></thead><tbody><tr><td colspan="6" style="text-align: center;">Planetary parameters</td></tr><tr><td>P (d)</td><td>0.4465688<sup>+0.0000007</sup><sub>-0.0000008</sub></td><td>10.778831<sup>+0.000034</sup><sub>-0.000036</sub></td><td>25.7124<sup>+0.0001</sup><sub>-0.0002</sub></td><td>77.03<sup>+0.25</sup><sub>-0.24</sub></td><td>473<sup>+36</sup><sub>-25</sub></td></tr><tr><td>T<sub>0</sub> (TBJD)<sup>a</sup></td><td>2317.7498 ± 0.0005</td><td>2238.4629<sup>+0.0008</sup><sub>-0.0009</sub></td><td>2318.966<sup>+0.003</sup><sub>-0.004</sub></td><td>1538.180<sup>+0.004</sup><sub>-0.005</sub></td><td>1664<sup>+28</sup><sub>-33</sub></td></tr><tr><td>a/R<sub>*</sub></td><td>2.685<sup>+0.024</sup><sub>-0.025</sub></td><td>22.43<sup>+0.20</sup><sub>-0.21</sub></td><td>40.04<sup>+0.36</sup><sub>-0.37</sub></td><td>83.22<sup>+0.77</sup><sub>-0.79</sub></td><td>279<sup>+14</sup><sub>-10</sub></td></tr><tr><td>a (AU)</td><td>0.0106 ± 0.0001</td><td>0.0884 ± 0.0009</td><td>0.158 ± 0.002</td><td>0.328 ± 0.003</td><td>1.1<sup>+0.6</sup><sub>-0.4</sub></td></tr><tr><td>R<sub>p</sub>/R<sub>*</sub></td><td>0.0155 ± 0.0004</td><td>0.0316 ± 0.0004</td><td>0.0306 ± 0.0008</td><td>0.0278<sup>+0.0016</sup><sub>-0.0014</sub></td><td>—</td></tr><tr><td>R<sub>p</sub> (R<sub>⊕</sub>)</td><td>1.425 ± 0.037</td><td>2.91 ± 0.04</td><td>2.82 ± 0.07</td><td>2.55<sup>+0.12</sup><sub>-0.13</sub></td><td>—</td></tr><tr><td>b</td><td>0.13<sup>+0.10</sup><sub>-0.09</sub></td><td>0.12<sup>+0.13</sup><sub>-0.08</sub></td><td>0.45<sup>+0.11</sup><sub>-0.17</sub></td><td>0.28<sup>+0.15</sup><sub>-0.18</sub></td><td>—</td></tr><tr><td>i (deg)</td><td>87.2<sup>+1.9</sup><sub>-2.1</sub></td><td>89.69<sup>+0.21</sup><sub>-0.31</sub></td><td>89.40<sup>+0.21</sup><sub>-0.11</sub></td><td>89.80<sup>+0.13</sup><sub>-0.10</sub></td><td>—</td></tr><tr><td>T<sub>14</sub> (h)</td><td>1.31 ± 0.02</td><td>3.75<sup>+0.05</sup><sub>-0.08</sub></td><td>4.54<sup>+0.32</sup><sub>-0.29</sub></td><td>6.98<sup>+0.24</sup><sub>-0.40</sub></td><td>—</td></tr><tr><td>e</td><td>0 (fixed)</td><td>0.030<sup>+0.035</sup><sub>-0.021</sub></td><td>0.122<sup>+0.054</sup><sub>-0.048</sub></td><td>0.079<sup>+0.058</sup><sub>-0.050</sub></td><td>0.085<sup>+0.083</sup><sub>-0.059</sub></td></tr><tr><td>ω (deg)</td><td>90 (fixed)</td><td>291<sup>+55</sup><sub>-84</sub></td><td>235<sup>+14</sup><sub>-26</sub></td><td>143<sup>+42</sup><sub>-44</sub></td><td>348<sup>+198</sup><sub>-53</sub></td></tr><tr><td>K (m s<sup>-1</sup>)</td><td>1.93 ± 0.21</td><td>1.81<sup>+0.23</sup><sub>-0.22</sub></td><td>3.34<sup>+0.23</sup><sub>-0.22</sub></td><td>2.19 ± 0.23</td><td>1.94 ± 0.27</td></tr><tr><td>M<sub>p</sub> (M<sub>⊕</sub>)</td><td>2.00 ± 0.23</td><td>5.39<sup>+0.69</sup><sub>-0.68</sub></td><td>13.2<sup>+1.0</sup><sub>-0.9</sub></td><td>12.6 ± 1.4</td><td>20 ± 3<sup>b</sup></td></tr><tr><td>ρ<sub>p</sub> (ρ<sub>⊕</sub>)</td><td>0.69 ± 0.10</td><td>0.22 ± 0.03</td><td>0.59 ± 0.06</td><td>0.76 ± 0.14</td><td>—</td></tr><tr><td>ρ<sub>p</sub> (g cm<sup>-3</sup>)</td><td>3.8 ± 0.5</td><td>1.2 ± 0.2</td><td>3.2 ± 0.3</td><td>4.2 ± 0.8</td><td>—</td></tr><tr><td>S<sub>p</sub> (S<sub>⊕</sub>)</td><td>4745 ± 269</td><td>68.2 ± 3.9</td><td>21.4 ± 1.3</td><td>4.96 ± 0.28</td><td>—</td></tr><tr><td>T<sub>eq</sub> (K)</td><td>2310 ± 33</td><td>800 ± 11</td><td>598 ± 9</td><td>415 ± 6</td><td>—</td></tr><tr><td>g<sub>p</sub><sup>d</sup> (m s<sup>-2</sup>)</td><td>9.7 ± 1.2</td><td>6.2 ± 0.8</td><td>16.3 ± 1.5</td><td>19.0 ± 2.9</td><td>—</td></tr><tr><td colspan="6" style="text-align: center;">Common parameters</td></tr><tr><td>R<sub>*</sub><sup>e</sup> (R<sub>☉</sub>)</td><td></td><td>0.843 ± 0.005</td><td></td><td></td><td></td></tr><tr><td>M<sub>*</sub><sup>e</sup> (M<sub>☉</sub>)</td><td></td><td>0.806 ± 0.036</td><td></td><td></td><td></td></tr><tr><td>ρ<sub>*</sub> (ρ<sub>☉</sub>)</td><td></td><td>1.31 ± 0.05</td><td></td><td></td><td></td></tr><tr><td>u<sub>1</sub>, TESS</td><td></td><td>0.33 ± 0.08</td><td></td><td></td><td></td></tr><tr><td>u<sub>2</sub>, TESS</td><td></td><td>0.23 ± 0.09</td><td></td><td></td><td></td></tr><tr><td>u<sub>1</sub>, CHEOPS</td><td></td><td>0.46 ± 0.07</td><td></td><td></td><td></td></tr><tr><td>u<sub>2</sub>, CHEOPS</td><td></td><td>0.22 ± 0.09</td><td></td><td></td><td></td></tr><tr><td>σ<sup>i</sup><sub>HARPS-N</sub> (m s<sup>-1</sup>)</td><td></td><td>1.40<sup>+0.15</sup><sub>-0.14</sub></td><td></td><td></td><td></td></tr><tr><td>σ<sup>f</sup><sub>HIRES</sub> (m s<sup>-1</sup>)</td><td></td><td>2.77<sup>+0.36</sup><sub>-0.31</sub></td><td></td><td></td><td></td></tr><tr><td>γ<sup>g</sup><sub>HARPS-N</sub> (m s<sup>-1</sup>)</td><td></td><td>79700.41 ± 0.26</td><td></td><td></td><td></td></tr><tr><td>γ<sup>g</sup><sub>HIRES</sub> (m s<sup>-1</sup>)</td><td></td><td>-1.20 ± 0.42</td><td></td><td></td><td></td></tr></tbody></table>

<sup>a</sup>TESS Barycentric Julian Date (BJD-2457000). <sup>b</sup>Minimum mass in the hypothesis of a planetary origin. <sup>c</sup>Computed as $T_{eq} = T_{\star} (R_{\star}/2a)^{1/2}$ [$f(1 - A_B)^{1/4}$, assuming $f=1$ and a null Bond albedo ($A_B = 0$). <sup>d</sup> Planetary surface gravity. <sup>e</sup> As determined from the stellar analysis in Section 2.1. <sup>f</sup> RV jitter term. <sup>g</sup> RV offset.

from the five-Keplerian model ones, indicating that the different modelling of the long-term signal is not influencing the results for the known, transiting planets. Finally, as in the case of the 5-Keplerian fit, the HARPS-N jitter is significantly improved ($\sigma_{HARPS-N} \sim 1.30 \text{ m s}^{-1}$) when including the GP model. Therefore, since our Bayesian analyses showed that the modelling of the long-term signal is necessary to obtain the best picture of the system, we decided to perform the global fit assuming a five-Keplerian model, but without drawing conclusions on the origin of the fifth signal. We stress that the five-Keplerian fit does not provide absolute evidence of the presence of a fifth planet, since also poorly sampled stellar activity could be well modelled using a Keplerian (Pepe et al. 2013; Mortier & Collier Cameron 2017; Affer et al. 2016), especially in our case where the RV baseline is of the order of the signal periodicity. Since it is not possible to distinguish a true planetary signal from an activity signal that has not been observed long enough to exhibit a loss of coherence in its phase or amplitude, only a follow-up campaign over several years can allow one to better understand the nature of this long-term signal.

## 5 JOINT PHOTOMETRIC AND RV ANALYSIS

To infer the properties of the TOI-561 planets, we jointly modelled all photometric and spectroscopic data with PYORBIT, using PYDE<fnref n="11" /> + EMCEE (Foreman-Mackey et al. 2013) as described in section 5 of L21, and adopting the same convergence criteria. We ran 96 chains (twice the number of the model parameters) for 250 000 steps, discarding the first 50 000 as burn-in.

Based on the analysis presented in the previous section, we assumed a five-Keplerian model, including four planets plus a fifth Keplerian with period free to span between 2 and 900 d. We fitted a common value for the stellar density, using the value reported in Table 1 as Gaussian prior. We adopted the quadratic limb-darkening law as parametrized by Kipping (2013), putting Gaussian priors on the $u_1$, $u_2$ coefficients, obtained for the TESS and CHEOPS passband through a bilinear interpolation of limb darkening profiles by Claret (2017) and Claret (2021), respectively, and assuming a $1\sigma$

<footnote n="11">¹¹https://github.com/hpparvi/PyDE.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** GLS periodogram of the HARPS-N activity indices. The main peak of each periodogram is highlighted with a red vertical line. The dashed and dotted horizontal lines indicate the 1 and 0.1 per cent FAP levels, respectively. The only peak above the 0.1 FAP level is the low-frequency peak in the S-index periodogram, as discussed in Section 4.2.</figcaption></figure>

uncertainty of 0.1 for each coefficient. We imposed a half-Gaussian zero-mean prior (Van Eylen et al. 2019) on the planet eccentricities, except for the USP planet, whose eccentricity was fixed to zero. We assumed uniform priors for all the other parameters.

To model the long-term correlated noise in the *TESS* light curve, we included in the fit a GP regression with a Matérn-3/2 kernel against time, as shown in Fig. 1, and we added a jitter term to account for possible extra white noise. We pre-decorrelated the *CHEOPS* light curves with the *PYCHEOPS*<fnref n="12" /> package (Maxted et al. 2021), selecting the detrending parameters according to the Bayes factor to obtain the best correlated noise model for each visit. For all the three *CHEOPS* visits, a decorrelation for the first three harmonics of the roll angle was necessary, plus first-order polynomials in time, x-y centroid position, and smearing. We then used the detrended light curves (Fig. 3) for the global PYORBIT fit. In order to check if the detrending was affecting our results for the planetary parameters, we also performed an independent global analysis with the JULIET package (Espinoza, Kossakowski & Brahm 2019), including in the global modelling the basis vectors selected with PYCHEOPS to detrend the data simultaneously. All the results were consistent within 1σ, indicating that the pre-detrending did not significantly alter our inferred results. Finally, for both the HARPS-N and HIRES data sets we included jitter and offset terms as free parameters.

We summarize our best-fitting model results in Table 6, and we show the transit model, phase-folded RVs and global RV model in Figs 8, 9, and 10, respectively. We inferred precise masses and

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Phase-folded *TESS* (left) and *CHEOPS* (right) light curves of TOI-561 b, c, and d. Planet e shows a single transit in the *TESS* light curve, and it has no *CHEOPS* observations. For each planet, the coloured line indicates the best-fitting model, and residuals are shown in the bottom panels. Data points binned over 20 min (planet b) and 30 min (planets c, d and e) are shown with coloured dots.</figcaption></figure>

radii for all the four planets in the system, whose positions in the mass-radius diagram are shown in Fig. 11. With a radius of $R_b = 1.425 \pm 0.037 R_⊕$ and a mass of $M_b = 2.00 \pm 0.23 M_⊕$ (from $K_b = 1.93 \pm 0.21 \text{ m s}^{-1}$), TOI-561 b is located in a region of the mass-radius diagram which is not consistent with a pure rocky composition, as will be also shown in Section 6 by our internal structure modelling. Our analysis confirms TOI-561 b to be the lowest density ($\rho_b = 3.8 \pm 0.5 \text{ g cm}^{-3}$) USP planet known to date (Fig. 12). In order to further confirm the planetary density, we also performed a specific RV analysis of TOI-561 b using the Floating Chunk Offset method (FCO; Hatzes 2014). The FCO analysis, detailed in Appendix B,

<footnote n="12">¹²https://github.com/pmaxted/pycheops.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Phase-folded HARPS-N and HIRES RVs with residuals of TOI-561 b, c, d and e, as resulting from the joint photometric and spectroscopic fit. The error bars include the jitter term added in quadrature.</figcaption></figure>

confirms the low mass inferred for TOI-561 b, and consequently its low density. Thanks to the CHEOPS observations, we also improved significantly the radius of TOI-561 c, for which we obtained a value of $R_c = 2.91 \pm 0.04 R_⊕$. From the semi-amplitude $K_c = 1.81^{+0.23}_{-0.22}$ m s⁻¹, we inferred a mass of $M_c = 5.39^{+0.69}_{-0.68} M_⊕$, implying a density of $\rho_c = 1.2 \pm 0.2 \text{ g cm}^{-3}$. From the combined fit of one TESS and one

<figure><img src="image_13.png" /><figcaption>**Figure 10.** Global model (grey line) with residuals of HARPS-N and HIRES RVs according to the 5-Keplerian photometric and spectroscopic fit. The error bars include the jitter term added in quadrature.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Mass-radius diagram for exoplanets with radii and masses measured with a precision better than 30 per cent, colour coded according to their incidental flux. Data are taken from the Extrasolar Planets Encyclopaedia catalogue (http://exoplanet.eu/catalog/, Schneider et al. 2011) as of 18 October 2021. The TOI-561 planets are labelled, and highlighted with coloured diamonds. The USP planets are emphasized with thick, black-contoured circles. The theoretical mass-radius curves for various chemical compositions (Zeng et al. 2019) are represented by solid coloured lines, while the dashed lines indicate the curves for an Earth-like core surrounded by a H₂ envelope (2 per cent mass fraction) at varying equilibrium temperatures. The forbidden region predicted by collisional stripping (Marcus et al. 2010) is marked by the shaded grey region.</figcaption></figure>

CHEOPS transit, we inferred a radius of $2.82 \pm 0.07 R_⊕$ for planet d, which has a mass of $M_c = 13.2^{+1.0}_{-0.9} M_⊕$ (from $K_d = 3.34^{+0.23}_{-0.22}$ m s⁻¹) and a resulting density of $\rho_d = 3.2 \pm 0.3 \text{ g cm}^{-3}$. Finally, for TOI-561 e, which shows a single transit in TESS sector 8, we derived a radius of $R_e = 2.55^{+0.12}_{-0.13} R_⊕$, a mass of $M_e = 12.6 \pm 1.4 M_⊕$, and an average density of $\rho_d = 4.2 \pm 0.8 \text{ g cm}^{-3}$. Lastly, the period inferred for the fifth Keplerian in the model was $473^{+36}_{-25}$ d, with a
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Mass–radius diagram of confirmed USP planets (*P* < 1 d, *R*<sub>p</sub> < 2 *R*<sub>⊕</sub>) as taken from the Extrasolar Planets Encyclopaedia catalogue in date 18 October 2021. Iso-density lines are plotted in grey. TOI-561 b stands out as the lowest density USP planet known to date (ρ<sub>b</sub> = 3.8 ± 0.5 g cm<sup>-3</sup>).</figcaption></figure>

7.2σ detected semi-amplitude of 1.94 ± 0.27 m s⁻¹. As discussed in Section 4.2, additional data spanning a longer baseline are needed to definitively confirm the planetary nature of this long-term signal.

# 6 INTERNAL STRUCTURE MODELLING

We modelled the internal planetary structure in a Bayesian framework, following the procedure detailed in Leleu et al. (2021). Our model assumes fully-differentiated planets composed of four layers, comprising an iron and sulfur central core, a silicate mantle which includes Si, Mg, and Fe, a water layer, and a pure H/He gas layer. The inner core is modelled assuming the Hakim et al. (2018) equation of state (EOS), the silicate mantle uses the Sotin, Grasset & Mocquet (2007) EOS, and the water layer uses the Haldemann et al. (2020) EOS. The core, mantle, and water layer compose the ‘solid’ part of the planet. The thickness of the gas envelope is computed as a function of stellar age and irradiation, and mass and radius of the solid part, according to the model presented in Lopez & Fortney (2014). We assumed no compression effects of the gas envelope on the solid part, a hypothesis which is justified a *posteriori* given the low mass fraction of gas obtained for each planet (see below).

Our Bayesian model fits the planetary system as a whole, rather than performing an independent fit for each planet, in order to account for the correlations between the absolute planetary masses and radii, which depend on the stellar properties. The model fits the stellar (mass, radius, effective temperature, age, chemical abundances of Fe, Mg, Si) and planetary properties (RV semi-amplitudes, transit depths, orbital periods) to derive the posterior distributions of the internal structure parameters. The internal structure parameters modelled for each planet are the mass fractions of the core, mantle and water layer, the mass of the gas envelope, the iron molar fraction in the core, the silicon and magnesium molar fraction in the mantle, the equilibrium temperature, and the age of the planet (equal to the age of the star). For a more extensive discussion on the relation among input data and derived parameters we refer to Leleu et al. (2021). We assumed the mass fraction of the inner core, mantle, and water layer to be uniform on the simplex (the surface on which they add up to one), with the water mass fraction having an upper boundary of 0.5 (Thiabaud et al. 2014; Marboeuf et al. 2014). For the mass of the gas envelope, we

assumed a uniform prior in logarithmic space. Finally, we assumed the Si/Mg/Fe molar ratios of each planet to be equal to the stellar atmospheric values (even though Adibekyan et al. (2021a) recently showed that the stellar and planetary abundances may not be always correlated in a one-to-one relation). We emphasize the fact that, as in many Bayesian analyses, the results presented below in terms of planet internal structure depend to some extent on the selection of the priors, which we chose following i.e. Dorn et al. (2017), Dorn et al. (2018), and Leleu et al. (2021). Analysing the same data with very different priors (e.g. non-uniform core/mantle/water mass fraction or gas fraction uniform in linear scale) would lead to different conclusions.

We show the results of the internal structure modelling for the four planets in Fig. 13. As expected from its closeness to the host star, planet b has basically no H/He envelope, while the other three planets show a variable amount of gas mass. Planet c hosts a relatively massive gaseous envelope, with a gas mass of (5 and 95 per cent quantiles) $M_{gas,c} = 0.07_{-0.02}^{+0.04}$ M<sub>⊕</sub> (1.3<sub>-0.4</sub><sup>+0.8</sup> weight per cent wt per cent). Planet d hosts the most massive envelope ($M_{gas,d} = 0.10_{-0.07}^{+0.13}$ M<sub>⊕</sub>), which, considering the total mass of the planet, correspond to a smaller relative mass fraction of 0.8<sub>-0.5</sub><sup>+1.0</sup> wt per cent, while TOI-561 e's envelope spans a range between -10.7 < log $M_{gas,e}$ < -1.0, implying an upper limit on the gas mass of 0.11 M<sub>⊕</sub> (<0.9 wt per cent). As expected from its low density, TOI-561 b could host a significant amount of water, having a water mass of $M_{H_2O,b} = 0.62_{-0.44}^{+0.32}$ M<sub>⊕</sub> (31<sub>-22</sub><sup>+16</sup> wt per cent). We stress that this result is highly dependent on the caveat of including only a solid water layer in the model. In fact, a massive water layer, if present on a planet with such a high equilibrium temperature, would imply the presence of a massive steam atmosphere (Turbet et al. 2020). This would in turn considerably change the inferred water mass fraction with respect to a model that includes only a solid water layer. Due to the presence of the gas envelope, the amount of water in both planet c and d is almost unconstrained ($M_{H_2O,c} = 1.29_{-1.14}^{+1.24}$ M<sub>⊕</sub>, i.e. 24<sub>-21</sub><sup>+23</sup> wt per cent; $M_{H_2O,d} = 3.56_{-3.18}^{+2.78}$ M<sub>⊕</sub>, i.e. 27<sub>-24</sub><sup>+21</sup> wt per cent), while TOI-561 e modelling points toward a massive water layer, with $M_{H_2O,e} = 4.50_{-3.65}^{+1.69}$ M<sub>⊕</sub> (36<sub>-29</sub><sup>+13</sup> wt per cent).

# 7 ATMOSPHERIC EVOLUTION

We employed the system parameters derived in this work to constrain the evolution of the stellar rotation period, which we use as a proxy for the evolution of the stellar high-energy emission affecting atmospheric escape, and the predicted initial atmospheric mass fraction of the detected transiting planets $f_{\text{atm}}^{\text{start}}$, that is the mass of the planetary atmosphere at the time of the dispersal of the protoplanetary disc. To this end, we used the planetary atmospheric evolution code PASTA described by Bonfanti et al. (2021a), which is an updated version of the original code presented by Kubyshkina et al. (2019b), Kubyshkina et al. (2019a). The code models the evolution of the planetary atmospheres combining a model predicting planetary atmospheric escape rates based on hydrodynamic simulations (this has the advantage over other commonly used analytical estimates to account for both XUV-driven and core-powered mass loss; Kubyshkina et al. 2018), a model of the stellar high-energy (X-ray plus extreme ultraviolet; XUV) flux evolution (Bonfanti et al. 2021a), a model relating planetary parameters and atmospheric mass (Johnstone et al. 2015b), and stellar evolutionary tracks (Choi et al. 2016). The main assumptions of the framework are that planet migration did not occur after the dispersal of the protoplanetary disk, and that the planets hosted at some point in the past or still host a hydrogen-dominated atmosphere.
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Posterior distributions of the main parameters describing the internal structure of TOI-561 b (**top left**), c (**top right**), d (**bottom left**), and e (**bottom right**). Each corner plot shows the mass fraction of the inner core and of the water layer, the molar fractions of silicon and magnesium in the mantle, the iron molar fraction in the inner core, and the mass of gas in logarithmic scale. On top of each column are printed the mean and the 5 per cent and 95 per cent quantiles values. For each planet, and we show an illustration of the radius fractions of the inner core+mantle (dark grey), water layer (dark blue), and gas envelope (light blue), corresponding to the medians of the posterior distributions. The coloured rectangles indicate the uncertainty on the corresponding layer thickness, while the black dashed outer rings represent the uncertainty on the total radius. Equilibrium temperature and planetary surface gravity are reported for each planet.</figcaption></figure>

For each planet, the evolution calculations begin at an age of 5 Myr, which is the age assumed in the code for the dispersal of the protoplanetary disc. At each time step, the framework derives the mass-loss rate from the atmospheric escape model employing the stellar flux and the system parameters, and uses it to update the atmospheric mass fraction. This procedure is then repeated until the age of the system is reached or the planetary atmosphere has completely escaped. The free parameters of the algorithm are the initial atmospheric mass fraction at the time of the dispersal of the

protoplanetary disc, and the indexes of the power law controlling the stellar rotation period (see Bonfanti et al. 2021a, for a detailed description of the mathematical formulation of the power law), that we use as proxy for the stellar XUV emission.

The free parameters are constrained by implementing the atmospheric evolution algorithm in a Bayesian framework employing the MCMC tool presented by Cubillos et al. (2017). The framework uses the system parameters with their uncertainties as input priors. It then computes millions of forward planetary evolutionary tracks, varying
<figure><img src="image_17.png" /><figcaption>**Figure 14.** *From left to right:* Posterior distributions (dark blue lines) of the stellar rotation period at an age of 150 Myr and of the initial atmospheric mass fractions of TOI-561 b, c, d, and e. In each panel, the purple region represents the 68 per cent highest probability density intervals. In the left-hand panel, the black thin line shows the rotation period distribution of stars member of open clusters with ages around 150 Myr. Data are taken from Johnstone et al. (2015a), who report the rotation period of ~ 2000 stars belonging to the Pleiades, M50, M35, and NGC 2516, whose ages are between 125 and 150 Myr. To generate the black histogram we selected a sub-sample of 578 stars, which have masses that deviate from $M_*$ less than 0.1 $M_\odot$. In the other panels, the horizontal orange lines mark the uniform prior used in the fit, scaled to the highest peak of each posterior distribution for better visualization. The light blue lines indicate the current atmospheric mass fraction of each planet determined as described in Section 6.</figcaption></figure>

the input parameters according to the shape of the prior distributions, and varying the free parameters within pre-defined ranges, fitting the current planetary atmospheric mass fractions obtained as described in Section 6. The fit is done at the same time for all planets, thus simultaneously constraining the rotational period, and the results are posterior distributions of the free parameters. In particular, we opted for fitting for the planetary atmospheric mass fractions instead of the planetary radii. This enables the code to be more accurate by avoiding the continuous conversion of the atmospheric mass fraction into planetary radius, given the other system parameters (see also Delrez et al. 2021).

Fig. 14 shows the results of the planetary atmospheric evolution simulations. As a proxy for the evolution of the stellar rotation period, in Fig. 14, we show the posterior distribution of the stellar rotation period at an age of 150 Myr, further comparing it to the distribution of stellar rotation periods observed in stars member of young clusters of comparable age and with masses that deviate from $M_*$ less than 0.1 $M_\odot$ (from Johnstone et al. 2015a). The inferred posterior distribution for the rotation period is consistent with membership of the slowly rotating period-colours sequence in clusters of this age. However, this comparison should be taken with some caution, since there are no comprehensive studies on the rotation-colour distributions of 150 Myr-old clusters with the same metallicity as TOI-561. The initial atmospheric mass fractions of planets b and c are rather broad and peak at about one planetary mass. This is because both planets are close enough to the host star and have a small enough mass to have been subject to significant atmospheric escape. Therefore, to enable the presence of a thin hydrogen atmosphere, as predicted by the internal structure model, both planets had to host a significant hydrogen envelope after the formation and atmospheric accretion processes. Instead, planets d and e are far from the host star and massive enough not to have been subject to significant atmospheric escape, which is why we obtain an initial atmospheric mass fraction that resembles the current one. We also find that the posterior distributions of all input parameters match well the inserted priors (not shown here). As a whole, the results indicate that the currently observed system parameters are compatible with a scenario in which migration happened (if at all) exclusively inside the protoplanetary disk. Otherwise the code would have led to mismatches between the prior and posterior of the input parameters (particularly for what concerns the planetary masses and/or the stellar mass and age), in addition to showing incoherent results in the posterior distribution of the output parameters. This is

for example the case of the TOI-1064 system, which is composed by two transiting planets with comparable masses and irradiation levels, but significantly different radii (Wilson et al. 2022). In our framework in which planets do not migrate after the dispersal of the protoplanetary nebula, reproducing the physical parameters of the planets composing the TOI-1064 system requires different evolutions of the stellar rotation rate, which is not possible, thus calling for a post-nebula migration.

## 8 DISCUSSION AND CONCLUSIONS

In this study, we confirm the presence of four transiting planets around TOI-561, with orbital periods of approximately 0.44, 10.8, 25.7, and 77 d (Table 6). Our analysis disproves the presence of the previously suggested planet TOI-561 f ($P \sim 16.3$ d; W21). TOI-561 is one of the few four-planet systems having precise radius and mass measurements for all the planets. Thanks to our global photometric and RV analysis, we refined all masses and radii with respect to the L21 values, and we precisely determined the planetary bulk densities, with uncertainties of 14.4 per cent, 13.6 per cent, 10.2 per cent, and 18.4 per cent for planets b, c, d, and e, respectively. The higher uncertainty on planet e reflects the lower precision in the radius determination (5 per cent uncertainty), which is based on the analysis of a single TESS transit, and highlights the importance of the high-precision CHEOPS photometry. In fact, with a single CHEOPS transit we managed to decrease the uncertainty on the radius of planet d from 5.1 per cent (L21, based on one TESS transit) to 2.5 per cent. Including also the improvement on the mass, this implied a decrease on the density uncertainty from 18.9 per cent to 10.2 per cent. We expect a similar improvement for planet e with future CHEOPS observations scheduled for 2022. The improvement in the radius of TOI-561 e is particularly important, since the planet is an interesting target for the study of the internal structure of cold sub-Neptunes. Its long period ($P_d = 77.03_{-0.24}^{+0.25}$ d) implies an insolation flux of $S_e = 4.96 \pm 0.28$ $S_⊕$ and a relatively cool zero Bond albedo equilibrium temperature of $T_{eq. e} = 415 \pm 6$ K. As shown in Fig. 15, TOI-561 e is one of the few cool, long-period planets orbiting a star bright enough for precise RV characterization, and it is therefore an optimal test-case to refine tools and models that will be useful to characterize targets of future long-staring missions like PLATO.

TOI-561 hosts one of the most intriguing USP planets discovered to date. As initially suggested by L21, our analysis confirms that TOI-561 b is the lowest density ($\rho_b = 3.8 \pm 0.5$ g cm$^{-3}$) USP super-Earth
<figure><img src="image_18.png" /><figcaption>**Figure 15.** V magnitude versus planetary periods for confirmed transiting exoplanets as reported in the Extrasolar Planets Encyclopaedia catalogue (http://exoplanet.eu/catalog/) in date 2021 October 18. The dashed vertical line marks V = 12 mag. TOI-561 e is one of the few long-period planets orbiting a star bright enough for precise RV characterization.</figcaption></figure>

that we know of (see Fig. 12), and it paves the way for in-depth studies of interior composition, and formation and evolution processes of USP planets. Even though now the mass values are consistent within $1\sigma$, contrary to what proposed by W21 (see Section 2.2) TOI-561 b is not consistent with a pure rocky composition, and to explain the planetary density our internal structure modelling (Section 6) predicts basically no H/He envelope, and a massive water layer. In this regard, an important point to consider is that, with an insolation flux of $S_b \simeq 4745 S_⊕$, the planet receives more irradiation from the star than the theoretical runaway greenhouse limit (Kasting, Whitmire & Reynolds 1993; Goldblatt & Watson 2012; Kopparapu et al. 2013). In this case, a large water content would imply the presence of an extended steam atmosphere, which in turn would increase the measured radius with respect to a purely condensed water world, leading in our model to an overestimation of the bulk water content (Turbet et al. 2020). The presence of a water steam envelope could eventually be tested with the *James Webb Space Telescope* (JWST). In fact, with an Emission Spectroscopy Metric (ESM; Kempton et al. 2018) value of 8.2, TOI-561 b is a promising target for secondary eclipse and phase curve observations. More complex models, including a lighter core compositions (i.e. a Ca/Al enriched core), the modelling of water steam envelopes, or wet-melt solid interiors related to deep water reservoirs (Dorn & Lichtenberg 2021), could be an interesting step forward in the understanding of the planet structure and composition. The low density of TOI-561 b could also be related to the fact that the host star is a metal-poor, thick-disk star. Adibekyan et al. (2021a) showed that the composition of the rocky planets reflects the chemical abundances of the host star (even though not in a one-to-one relation), so implying a lighter composition for TOI-561 b with respect to other USP planets that orbit more metal-rich stars.<fnref n="13" /> According to Adibekyan et al. (2021a), the low density of TOI-561 b is consistent with the general

$\rho/\rho_{\text{Earth-like}} - f_{\text{iron}}^{\text{star}}$ trend and dispersion inferred from the sample of rocky planets analysed by the authors (see Figs. 2, 3 therein), where $\rho/\rho_{\text{Earth-like}}$ is the planetary density normalised to that expected for an Earth-like composition, and $f_{\text{iron}}^{\text{star}}$ is the iron-to-silicate mass fraction of the protoplanetary disk as inferred from the stellar properties. An additional interesting remark concerns the Galactic kinematics of the host star. According to our analysis, performed as described in Mustill, Lambrechts & Davies (2021), TOI-561 is located in a low-density region of the 6-dimensional Galactic phase space (see Winter et al. 2020, Mustill et al. 2021, and Kruijssen et al. 2021 for definition and discussion), which is not surprising given that TOI-561 is a thick disk star (Mustill et al. 2021). Kruijssen, Longmore & Chevance (2020) showed that stars in low-density regions seem to host no super-Earths, but only sub-Neptunes, i.e. planets having a significant H/He envelope and therefore located above the radius gap. In this context, TOI-561 b is an interesting object that runs counter to this finding. We point out that this result should be taken with some caution, since the Kruijssen et al. (2020) sample does not include planets with periods shorter than one day, and it excludes stars with ages >4.5 Gyr.<fnref n="14" />

All the four planets seem to host a large water layer (Section 6), although with high uncertainties, especially for planet c and d, due to the degeneracy related to the possible presence of a gas envelope. Also in this case, the presence of a considerable amount of water could be linked with the stellar properties. In fact, Santos et al. (2017) showed that metal-poor, thick disk stars are expected to form planetary building blocks with a higher water mass fraction (~76 per cent) compared to metal-rich, thin disk stars (~58 per cent). Therefore, we would expect these stars to produce water-rich planets, a result that is in agreement with our findings on the TOI-561 system.

Except for TOI-561 b, all the other planets are suggested to host a non-negligible H/He envelope. In particular, the gas content of planet c (~1.3 wt per cent, the highest mass fraction among the four planets) implies a much lower density with respect to the density of planet d, even though the two planets have a similar size. This is reflected in the different positions of the planets in the mass-radius diagram (Fig. 11). The two planets show hints of a different evolution for what concerns their gas content. In fact, our atmospheric evolution model (Section 7) suggests that planet c underwent a strong envelope loss after the atmospheric accretion and the dispersal of the protoplanetary nebula, while planet d (as well as planet e) did not experience strong atmospheric escape, with a current gas content that is comparable to the original one. The surprising difference in gas mass fraction between planets c, d and e, not only at present time but also at the end of their formation phase, takes probably its origin in the conditions that prevailed during the protoplanetary disk phase. Planet c is indeed likely sub-critical because of its low mass, where sub-critical planets are those with masses below the critical value required to initiate runaway gas accretion (see Helled et al. 2014 for a recent review on the core accretion model), whereas planets d and e never accreted large amounts of gas as demonstrated in Section 7, and so they also remained always below the critical mass. The interpretation of the different gas mass fractions could therefore result from the structure of sub-critical planets. In this case, the gas mass fraction depends on the core mass, the thermodynamical properties in the disc, and more importantly the accretion rate of

<footnote n="14">¹⁴We note however that the stellar ages used in Kruijssen et al. (2020) are quite inhomogeneous, coming directly from the NASA Exoplanet Archive, and can therefore show a large scatter with respect to a homogeneous determination (Adibekyan et al. 2021b).</footnote>

<footnote n="13">¹³All the USP planets shown in Fig. 11 have [Fe/H] > −0.14.</footnote>
solids (lower accretion rate translating in larger gas mass fraction). Interpreting the internal structure of the four planets of the system in a global planetary system formation model could therefore constrain these parameters.

With its derived properties, TOI-561 c has a Transmission Spectroscopy Metric (TSM, Kempton et al. 2018) of 110.4, and is therefore a suitable target for atmospheric characterization with JWST.<fnref n="15" /> Instead, planets d and e have lower TSM values of 30.7 and 16.2, respectively. As the TSM is proportional to the equilibrium temperature, it is not surprising to obtain lower values for the two planets, given their longer periods.

In addition to the characterization of the four planets, we also identified a significant long-term signal ($P \sim 473$ d) in the RVs. On the basis of our current data set, we are not able to distinguish between a stellar (magnetically-induced) or planetary origin. Long-term monitoring using both spectroscopic ground-based facilities and future long-staring missions like the PLATO spacecraft will allow us to shed light on the nature of this additional signal and to potentially find new outer companions. It is worth noting that, if the above-mentioned signal proves in future to be of planetary origin, there is a non-zero chance that, under the assumption of co-planarity, such a planet would transit. In fact, assuming the same inclination of planet e and using the semimajor axis $a/R_\odot = 279_{-10}^{+14}$ derived from our global fit, we infer an impact parameter of $0.97_{-0.63}^{+0.49}$. Moreover, the planet would orbit in TOI-561's empirical habitable zone ($175 \le P \le 652$ d), as originally defined by Kasting et al. (1993) using a 1D climate model, and later updated in Kopparapu et al. (2013), Ramirez & Kaltenegger (2016) for main-sequence stars with $2600 < T_{\text{eff}} < 10\,000$ K.

This work bears witness to the fruitful results that can be obtained by the timely combination of data coming from different instruments. It adds to the works (Bonfanti et al. 2021b; Leleu et al. 2021; Delrez et al. 2021) that prove the potential of CHEOPS in precisely characterizing TESS-discovered exoplanets, as well as demonstrating the key role of high-precision spectrographs such as HARPS-N when working in synergy with space-based facilities.

## ACKNOWLEDGEMENTS

We thank the referee for the useful comments that helped improving the quality of the manuscript. We gratefully acknowledge the MuSCAT2 team for the availability in collaborating with the CHEOPS Consortium on target monitoring, and the NGTS team for providing ground-based photometry that helped the scheduling of our observations. 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 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 the mission. This work is based on observations made with the Italian Telescopio Nazionale Galileo (TNG) operated on the island of La Palma by the Fundación Galileo Galilei of the INAF at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias (GTO program, and A40TAC_23 program from INAF-TAC). The HARPS-N project was funded by the Prodex Program of the Swiss Space Office

(SSO), the Harvard University Origin of Life Initiative (HUOLI), the Scottish Universities Physics Alliance (SUPA), the University of Geneva, the Smithsonian Astrophysical Observatory (SAO), and the Italian National Astrophysical Institute (INAF), University of St. Andrews, Queen's University Belfast and University of Edinburgh. This paper includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST). Funding for the TESS mission is provided by the NASA Explorer Program. 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 research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This research has made use of data obtained from the portal http://www.exoplanet.eu/ of The Extrasolar Planets Encyclopaedia. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. GL acknowledges support by CARIPARO Foundation, according to the agreement CARIPARO-Università degli Studi di Padova (Pratica n. 2018/0098). TW and ACC acknowledge support from STFC consolidated grant numbers ST/R000824/1 and ST/V000861/1, and UKSA grant number ST/R003203/1. YA, MJH, B-OD, and ML acknowledge the support of the Swiss National Fundation under grant numbers 200020_172746, PP00P2-190080, and PCEFP2_194576. SH gratefully acknowledges CNES funding through the grant number 837319. GPi, VNa, GSs, IPa, LBo, and RRa acknowledge the funding support from Italian Space Agency (ASI) regulated by 'Accordo ASI-INAF n. 2013-016-R.0 del 9 luglio 2013 e integrazione del 9 luglio 2015 CHEOPS Fasi A/B/C'. ADe acknowledges support from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (project FOUR ACES, grant agreement No. 724427), and from the National Centre for Competence in Research 'PlanetS' supported by the Swiss National Science Foundation (SNSF). KR is grateful for support from the UK STFC via grant number ST/V000594/1. This work has been supported by the National Aeronautics and Space Administration under grant number NNX17AB59G, issued through the Exoplanets Research Program. S.S. has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 833925, project STAREX). MG is an F.R.S.-FNRS Senior Research Associate. VVG is an F.R.S-FNRS Research Associate. LD is an F.R.S.-FNRS Postdoctoral Researcher. This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. AMu and MF acknowledge support from the Swedish National Space Agency (career grant numbers 120/19C, DNR 65/19, 174/18). ABr was supported by the SNSA. We acknowledge support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through

<footnote n="15">¹⁵Kempton et al. (2018) suggest to select planets with TSM >92 for 1.5 $R_⊕$ < $R_p$ <2.75 $R_⊕$, and TSM >84 for 2.75 $R_⊕$ < $R_p$ <4 $R_⊕$.</footnote>
grant numbers 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 Maria 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. SGS, SCCB, and VA acknowledge support from FCT through FCT contract nr. CEECIND/00826/2018, POPH/FSE (EC), nr. IF/01312/2014/CP1215/CT0004, and IF/00650/2015/CP1273/CT0001, respectively. ODSD is supported in the form of work contract (DL 57/2016/CP1364/CT0004) funded by national funds through FCT. XB, SC, DG, MF, and JL acknowledge their role as ESA-appointed CHEOPS science team members. 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. 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. 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). DG and LMS gratefully acknowledge financial support from the CRT foundation under grant number 2018.2323 'Gaseous rocky? Unveiling the nature of small worlds'. 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 Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche. PM acknowledges support from STFC research grant number ST/M001040/1. This work was also partially supported by a grant from the Simons Foundation (PI: Queloz, grant number 327127). GyMSz 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. This work is partly supported by JSPS KAKENHI grant number JP18H05439, JST CREST grant number JPMJCR1761, the Astrobiology Center of National Institutes of Natural Sciences (NINS) (grant number AB031010). EE-B acknowledges financial support from the European Union and the State Agency of Investigation of the Spanish Ministry of Science and Innovation (MICINN) under the grant number PRE2020-093107 of the Pre-Doc Program for the Training of Doctors (FPI-SO) through FSE funds.

DATA AVAILABILITY

HARPS-N observations and data products are available through the Data & Analysis Center for Exoplanets (DACE) at

https://dace.unige.ch/. TESS data products can be accessed through the official NASA website https://heasarc.gsfc.nasa.gov/docs/tess/da ta-access.html. All underlying data are available either in the online supporting material or will be available via VizieR at CDS, where we will provide both the raw and detrended CHEOPS light curves.

REFERENCES

Adams F. C., Batygin K., Bloch A. M., Laughlin G., 2020, MNRAS, 493, 5520

Adibekyan V. et al., 2021a, Science, 374, 330

Adibekyan V. et al., 2021b, A&A, 649, A111

Affer L. et al., 2016, A&A, 593, A117

Akeson R. L. et al., 2013, PASP, 125, 989

Baglin A. et al., 2006, in 36th COSPAR Scientific Assembly, 36, 3749

Bailer-Jones C. A. L., Rybizki J., Fouesneau M., Demleitner M., Andrae R., 2021, VizieR Online Data Catalog, I/352

Baranne A. et al., 1996, A&AS, 119, 373

Barros S. C. C. et al., 2022, A&A, 657, A52

Benz W. et al., 2021, Exp. Astron., 51, 109

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

Bonfanti A., Fossati L., Kubyshkina D., Cubillos P. E., 2021a, A&A, 656, A157

Bonfanti A., Gillon M., 2020, A&A, 635, A6

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

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

Bonfanti A. et al., 2021b, A&A, 646, A157

Borsato L. et al., 2021, MNRAS, 506, 3810

Borucki W. J. et al., 2010, Science, 327, 977

Buchhave L. A. et al., 2012, Nature, 486, 375

Buchhave L. A. et al., 2014, Nature, 509, 593

Castelli F., Kurucz R. L., 2003, in Piskunov N., Weiss W. W., Gray D. F., eds, IAU Symp. Vol. 210, Modelling of Stellar Atmospheres, Astron. Soc. Pac., San Francisco (USA), p. A20

Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ, 823, 102

Ciardi D. R., Fabrycky D. C., Ford E. B., Gautier T. N. I., Howell S. B., Lissauer J. J., Ragozzine D., Rowe J. F., 2013, ApJ, 763, 41

Caret A., 2017, A&A, 600, A30

Caret A., 2021, Res. Notes Am. Astron. Soc., 5, 13

Collier Cameron A. et al., 2021, MNRAS, 505, 1699

Collier Cameron A., 2018, The Impact of Stellar Activity on the Detection and Characterization of Exoplanets. Springer International Publishing, Cham, p. 1791 doi:10.1007/978-3-319-55333-7_23, https://doi.org/10.1 007/978-3-319-55333-7_23

Cosentino R. et al., 2012, in McLean I.S., Ramsay S.K., Takami H. (eds), Proc. SPIE Conf. Ser. Vol. 8446, Ground-based and Airborne Instrumentation for Astronomy IV, SPIE, Bellingham

Cosentino R. et al., 2014, in Ramsay S. K., McLean I. S., Hideki T., eds, Proc SPIE Conf. Ser. Vol. 9147, Ground-based and Airborne Instrumentation for Astronomy V, SPIE, Bellingham

Crass J. et al., 2021, Extreme Precision Radial Velocity Working Group Final Report. preprint (arXiv:2107.14291)

Cubillos P., Harrington J., Loredo T. J., Lust N. B., Blecic J., Stemm M., 2017, AJ, 153, 3

Cutri R. M. et al., 2003, VizieR Online Data Catalog, II/246

Dai F. et al., 2021, AJ, 162, 62

Deline A. et al., 2022, preprint (arXiv:2201.04518)

Delrez L. et al., 2021, Nature Astron., 5, 775

Dorn C. et al., 2017, A&A, 597, A37

Dorn C., Lichtenberg T., 2021, ApJL, 922, L4

Dorn C., Mosegaard K., Grimm S. L., Alibert Y., 2018, ApJ, 865, 20

Dumusque X. et al., 2021, A&A, 648, A103

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

Fabrycky D. C. et al., 2014, ApJ, 790, 146

Fausnaugh M. M., Burke C. J., Ricker G. R., Vanderspek R., 2020, Res. Notes Am. Astron. Soc., 4, 251
Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

Frustagli G. et al., 2020, A&A, 633, A133

Gaia Collaboration et al., 2018, A&A, 616, A1

Gaia Collaboration et al., 2021, A&A, 649, A1

Goldblatt C., Watson A. J., 2012, Phil. Trans. R. Soc. London Ser. A, 370, 4197

Grunblatt S. K., Howard A. W., Haywood R. D., 2015, ApJ, 808, 127

Guenther E. W. et al., 2017, A&A, 608, A93

Hakim K., Rivoldini A., Van Hoolst T., Cottenier S., Jaeken J., Chust T., Steinle-Neumann G., 2018, Icarus, 313, 61

Haldemann J., Alibert Y., Mordasini C., Benz W., 2020, A&A, 643, A105

Hara N. C., Boué G., Laskar J., Correia A. C. M., 2017, MNRAS, 464, 1220

Hatzes A. P., 2014, A&A, 568, A84

Hatzes A. P., 2019, The Doppler Method for the Detection of Exoplanets. IOP Publishing, Bristol (UK), p. 2514

Helled R. et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI, The University of Arizona Press, Tucson (USA), p. 643

Hooton M. J. et al., 2022, A&A, 658, A75

Howard A. W. et al., 2013, Nature, 503, 381

Howell S. B. et al., 2014, PASP, 126, 398

Hoyer S., Guterman P., Demangeon O., Sousa S. G., Deleuil M., Meunier J. C., Benz W., 2020, A&A, 635, A24

Jenkins J. M. et al., 2016, in Chiozzi G., Guzman J. C., eds, Proc. SPIE Conf. Ser. Vol. 9913, Software and Cyberinfrastructure for Astronomy IV. SPIE, Bellingham, p.1232

Jenkins J. M., 2020, Kepler Data Processing Handbook, Kepler Science Document KSCI-19081-003, Kepler Project Office, Moffett Field, California, https://archive.stsci.edu/kepler/documents.html

Jiang C.-F., Xie J.-W., Zhou J.-L., 2020, AJ, 160, 180

Johnstone C. P. et al., 2015b, ApJ, 815, L12

Johnstone C. P., Güdel M., Brott I., Lüftinger T., 2015a, A&A, 577, A28

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

Kasting J. F., Whitmire D. P., Reynolds R. T., 1993, Icarus, 101, 108

Kempton E. M. R. et al., 2018, PASP, 130, 114401

Kipping D. M., 2013, MNRAS, 435, 2152

Kopparapu R. K. et al., 2013, ApJ, 765, 131

Kruijssen J. M. D., Longmore S. N., Chevance M., 2020, ApJ, 905, L18

Kruijssen J. M. D., Longmore S. N., Chevance M., Laporte C. F. P., Motylinski M., Keller B. W., Henshaw J. D., 2021, preprint (arXiv:2109.06182)

Kubyshkina D. et al., 2018, A&A, 619, A151

Kubyshkina D. et al., 2019a, A&A, 632, A65

Kubyshkina D. et al., 2019b, ApJ, 879, 26

Lacedelli G. et al., 2021, MNRAS, 501, 4148 (L21)

Langellier N. et al., 2021, AJ, 161, 287

Leleu et al., 2021, A&A, 649, A26

Lendl M. et al., 2020, A&A, 643, A94

Lindegren L. et al., 2021, A&A, 649, A4

Lissauer J. J. et al., 2011, ApJS, 197, 8

Lopez E. D., Fortney J. J., 2014, ApJ, 792, 1

Malvolta L. et al., 2016, A&A, 588, A118

Malvolta L. et al., 2018, AJ, 155, 107

Malvolta L., Lovis C., Pepe F., Sneden C., Udry S., 2017, MNRAS, 469, 3965

Marboeuf U., Thiabaud A., Alibert Y., Cabral N., Benz W., 2014, A&A, 570, A36

Marcus R. A., Sasselov D., Stewart S. T., Hernquist L., 2010, ApJ, 719, L45
Marigo et al., 2017, ApJ, 835, 77

Maxted P. F. L. et al., 2021, preprint (arXiv:2111.08828)

Mayor M., Queloz D., 1995, Nature, 378, 355

Millholland S., Wang S., Laughlin G., 2017, ApJ, 849, L33

Mills S. M., Howard A. W., Petigura E. A., Fulton B. J., Isaacson H., Weiss L. M., 2019, AJ, 157, 198

Morris B. M. et al., 2021a, A&A, 653, A173

Morris D., Heng K., Brandeker A., Swan A., Lendl M., 2021b, A&A, 651, L12

Mortier A., Collier Cameron A., 2017, A&A, 601, A110

Mortier A., Sousa S. G., Adibekyan V. Z., Brandão I. M., Santos N. C., 2014, A&A, 572, A95

Mustill A. J., Lambrechts M., Davies M. B., 2021, preprint (arXiv:2103.15823)

Osborn H. P. et al., 2021, MNRAS, 502, 4842

Pepe F. et al., 2013, Nature, 503, 377

Pepe F., Mayor M., Galland F., Naef D., Queloz D., Santos N. C., Udry S., Burnet M., 2002, A&A, 388, 632

Prieto-Arranz J. et al., 2018, A&A, 618, A116

Ramirez R. M., Kaltenegger L., 2016, ApJ, 823, 6

Ricker G. R. et al., 2014, J. Astron. Telesc. Instr. Syst., 1, 1

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

Sanchis-Ojeda R. et al., 2015, ApJ, 812, L11

Santos N. C. et al., 2017, A&A, 608, A94

Schanche N. et al., 2020, MNRAS, 499, 428

Schneider J. et al., 2011, A&A, 532, A79

Scuflaire R., Théado S., Montalbán J., Miglio A., Bourge P. O., Godart M., Thoul A., Noels A., 2008, Ap&SS, 316, 83

Skilling J., 2004, in Fischer R., Preuss R., Toussaint U. V., eds, AIP Conf. Ser. Vol. 735, Bayesian Inference and Maximum Entropy Methods in Science and Engineering: 24th International Workshop on Bayesian Inference and Maximum Entropy Methods in Science and Engineering. p. 395 doi:10.1063/1.1835238

Skilling J., 2006, Bayesian Anal., 1, 833

Skrutskie M. F. et al., 2006, AJ, 131, 1163

Smith J. C. et al., 2012, PASP, 124, 1000

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

Sousa S. G., 2014, in Niemczura E., Smalley B., Pych W., eds, Determination of Atmospheric Parameters of B, A, F and G Type Stars, Springer International Publishing, Cham (SWI)

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

Stassun K. G. et al., 2018, AJ, 156, 102

Stumpe M. C. et al., 2012, PASP, 124, 985

Stumpe M. C., Smith J. C., Catanzarite J. H., Van Cleve J. E., Jenkins J. M., Twicken J. D., Girouard F. R., 2014, PASP, 126, 100

Swayne M. I. et al., 2021, MNRAS, 506, 306

Szabó G. M. et al., 2021, A&A, 654, A159

Thiabaud A., Marboeuf U., Alibert Y., Cabral N., Leya I., Mezger K., 2014, A&A, 562, A27

Turbet M., Bolmont E., Ehrenreich D., Gratier P., Leconte J., Selsis F., Hara N., Lovis C., 2020, A&A, 638, A41

Van Eylen V. et al., 2019, AJ, 157, 61

Van Grootel V. et al., 2021, A&A, 650, A205

Weiss L. M. et al., 2018, AJ, 156, 254

Weiss L. M. et al., 2021, AJ, 161, 56 (W21)

Wilson T. G. et al., 2022, preprint (arXiv:2201.03570)

Winter A. J., Kruijssen J. M. D., Longmore S. N., Chevance M., 2020, Nature, 586, 528

Wright E. L. et al., 2010, AJ, 140, 1868

Yi S., Demarque P., Kim Y.-C., Lee Y.-W., Ree C. H., Lejeune T., Barnes S., 2001, ApJS, 136, 417

Zechmeister M., Kürster M., 2009, A&A, 496, 577

Zeng L. et al., 2019, Proc. Natl. Acad. Sci., 116, 9723

SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

**Table 5. HARPS-N RV and activity indices measurements.**

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.
<figure><img src="image_19.png" /><figcaption>**Figure A1.** CHEOPS RINF light curves of TOI-561 as extracted from the DRP, with $4\sigma$-clipping for outliers removal. Visits 1, 2, and 3 are shown from top to bottom.</figcaption></figure>

## APPENDIX A: CHEOPS LIGHT CURVES AND TELEGRAPHIC PIXEL TREATMENT

As described in Section 3.2, the three CHEOPS visits of TOI-561 were reduced via the standard DRP processing. The light curves presented in this study, obtained using the RINF aperture size (RINF = 0.9 × DEFAULT, where DEFAULT = 25 px; see also Section 3.2), are shown in Fig. A1. While for the two initial visits the automatic DRP processes was performed, the appearance of some telegraphic pixels during the third visit required a more in-depth analysis.

In addition to the large number of known hot pixels present in the CHEOPS CCD (some of them visible in Fig. 2), some normal pixels can change their behaviour during the duration of a visit, for example becoming ‘hot’ after a SAA crossing. These pixels, called ‘telegraphic’ for their abnormal behaviour, can affect the photometry if located within the photometric aperture (see e.g. Leleu et al. 2021). During the third CHEOPS visit, we identified an unusual flux bump before the ingress of TOI-561 d transit, at BJD ~2459318.75 (top panel, Fig. A2). After analysing the statistics of each pixel light curve within the photometric aperture, we detected a telegraphic pixel with a large flux variation (see second panel, Fig. A2) located within the CHEOPS PSF. The exact position of this pixel on the CHEOPS CCD is shown in Fig. 2. We masked the pixel flux and repeated the photometric extraction of the visit using the RINF aperture, so removing the flux jump in the light curve (bottom panel, Fig. A2). During this analysis, we detected two additional telegraphic pixels within the photometric aperture, inducing smaller, but still significant variations in the light curve flux (third panel, Fig. A2). We corrected for the effect of these pixels as described above.

While investigating the nature of the flux bump happening during the third visit, we also extracted the light curve using a PSF-photometry approach exploiting the PIPE (PSF Imagette Photometric Extraction) software.<fnref n="16" /> PIPE is a photometric extraction package

<figure><img src="image_20.png" /><figcaption>**Figure A2.** Top panel: TOI-561 RINF original light curve of the third visit (light blue dots) after the removal of $4\sigma$ outliers, with overplotted the 15-min binned light curve (dark blue dots). The start of the flux jump due to the telegraphic pixel is marked with the dashed vertical line. Second panel: Light curve of the telegraphic pixel located within the CHEOPS PSF. Third panel: Light curve of the two additional telegraphic pixels located within the RINF aperture. Bottom panel: Corrected light curve after masking the three telegraphic pixels.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure A3.** Top panel: Comparison between DRP and PIPE-extracted light curve of TOI-561 third visit, before the telegraphic pixel correction and with $4\sigma$ outliers removal. The DRP has an MAD of 371 ppm over the whole visit, while PIPE of 325 ppm. Bottom panel: PIPE light curve after the telegraphic pixel correction. The light curve gets slightly noisier (MAD = 331 ppm) because one less pixel is considered in the reduction, but more reliable thanks to the exclusion of the telegraphic pixel flux. In both panels, the beginning of the flux jump is highlighted with a vertical dashed line.</figcaption></figure>

<footnote n="16">¹⁶https://pipe-cheops.readthedocs.io/en/latest/index.html</footnote>
specifically developed to extract CHEOPS light curves by applying PSF photometry on the 60-pixel imagettes, complementing the official DRP extraction. The use of PSF photometry makes usually easier to filter out the impact of hot pixels and cosmic rays, by either giving them a lower weight or masking them entirely in the fitting process. However, in this case the telegraphic pixel was located inside the CHEOPS PSF, requiring a careful manual masking. As for the DRP light curve, the flux bump in the PIPE photometry is reduced after masking the telegraphic pixel (bottom panel, Fig. A3). The PIPE-extracted light curve resulted in a slightly lower mean absolute deviation (MAD) with respect to the DRP photometry (top panel, Fig. A3), mainly due to the lower number of outliers present in the PSF photometry. For a more detailed comparison between PIPE and DRP photometries, see Morris et al. (2021b). We performed the same global analysis described in Section 5 using the PIPE light curve instead of the DRP one, obtaining consistent results and comparable uncertainties on the transit parameters of both planets b and d. We therefore decided to use the light curve obtained with the official DRP extraction in our final analysis.

## APPENDIX B: FLOATING CHUNK OFFSET METHOD ON TOI-561 B

In order to investigate the literature discrepancy on the mass of TOI-561 b (Section 2.2), we adopted a specific observing strategy with HARPS-N targeting the USP planet (Section 3.3), obtaining multiple observations during the same night for 22 nights. Multiple nightly observations can be used to precisely infer the mass of USP planets using the Floating Chunk Offset method (FCO; Hatzes 2014), which consists in applying a nightly offset to remove all the other signals present in the system, both of planetary and stellar origin (i.e. Howard et al. 2013; Pepe et al. 2013; Malavolta et al. 2018; Frustagli et al. 2020). The FCO method is only applicable when the separation between the USP period and the period of all the other signals is large enough, and the RV semi-amplitude has a similar or larger value with respect to the other signals. As demonstrated in L21, these conditions apply to TOI-561 b, for which the authors derived an FCO semi-amplitude of $K_b, FCO = 1.80 \pm 0.38 \text{ m s}^{-1}$ ($M_b, FCO = 1.83 \pm 0.39 \text{ M}_\oplus$) exploiting multiple observations collected over ten nights.

Here, we applied the FCO method to TOI-561 b on a total of 22 HARPS-N nights, adding 12 novel nights to the 10 nights already presented in L21. Out of the total set, four nights have six multiple observations extending over more than 40 per cent of the orbital period of the planet, and span opposite orbital phases to provide an

<figure><img src="image_22.png" /><figcaption>**Figure B1.** Phase-folded RVs of the 22 HARPS-N nights used to model TOI-561 b with the FCO method. The error bars include the jitter term added in quadrature.</figcaption></figure>

optimal phase coverage. We performed a PYDE + EMCEE fit with PYORBIT, assuming a fixed zero eccentricity and Gaussian priors on period and $T_0$ coming from the global fit, and we added a jitter term to account for possible additional white noise. We derived a semi-amplitude of $K_b = 1.81 \pm 0.31 \text{ m s}^{-1}$, corresponding to a mass of $M_b = 1.86 \pm 0.33 \text{ M}_\oplus$, with a jitter of $0.96_{-0.23}^{+0.25} \text{ m s}^{-1}$. Fig. B1 shows the resulting phase-folded RVs. The derived mass and semi-amplitude are nicely in agreement with the L21 values, and they support the values inferred from our joint photometric and RV modelling (Section 5), being consistent within $1\sigma$. Given the higher number of RVs included in the joint fit, which led to smaller uncertainties on the derived parameters, we decided to adopt as final values for TOI-561 b the ones obtained from the global modelling, i.e. $K_b = 1.93 \pm 0.21 \text{ m s}^{-1}$ and $M_b = 1.99 \pm 0.22 \text{ M}_\oplus$.

<footnote n="1">¹Dipartimento di Fisica e Astronomia 'Galileo Galilei', Università degli Studi di Padova, Vicolo dell'Osservatorio 3, I-35122 Padova, Italy</footnote>

<footnote n="2">²INAF, Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, I-35122 Padova, Italy</footnote>

<footnote n="3">³Centre for Exoplanet Science, SUPA School of Physics and Astronomy, University of St Andrews, North Haugh, St Andrews KY16 9SS, UK</footnote>

<footnote n="4">⁴Physikalisches Institut, University of Bern, Gesellsschaftstrasse 6, CH-3012 Bern, Switzerland</footnote>

<footnote n="5">⁵Center for Space and Habitability, Gesellsschaftstrasse 6, CH-3012 Bern, Switzerland</footnote>

<footnote n="6">⁶Astrophysics Group, Cavendish Laboratory, University of Cambridge, J.J. Thomson Avenue, Cambridge CB3 0HE, UK</footnote>

<footnote n="7">⁷Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="8">⁸Space Research Institute, Austrian Academy of Sciences, Schmiedlstrasse 6, A-8042 Graz, Austria</footnote>

<footnote n="9">⁹Astrophysics Group, University of Exeter, Exeter EX4 2QL, UK</footnote>

<footnote n="10">¹⁰Aix Marseille Univ, CNRS, CNES, LAM, 38 rue Frédéric Joliot-Curie, F-13388 Marseille, France</footnote>

<footnote n="11">¹¹Department of Astronomy, University of Geneva, Chemin Pegasi 51, Versoix, CH-1290, Switzerland</footnote>

<footnote n="12">¹²Department of Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA</footnote>

<footnote n="13">¹³Department of Astronomy, University of Wisconsin-Madison, Madison, WI 53706, USA</footnote>

<footnote n="14">¹⁴Center for Astrophysics | Harvard and Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="15">¹⁵SUPA, Institute for Astronomy, University of Edinburgh, Blackford Hill, Edinburgh EH9 3HJ, Scotland, UK</footnote>

<footnote n="16">¹⁶Centre for Exoplanet Science, University of Edinburgh, Edinburgh EH93FD, UK</footnote>

<footnote n="17">¹⁷Department of Astronomy, Stockholm University, AlbaNova University Center, SE-10691 Stockholm, Sweden</footnote>

<footnote n="18">¹⁸Fundación Galileo Galilei - INAF, Rambla J. A. F. Perez, 7, E-38712 S.C. Tenerife, Spain</footnote>

<footnote n="19">¹⁹INAF – Osservatorio Astronomico di Brera, via E. Bianchi 46, I-23807 Merate (LC), Italy</footnote>

<footnote n="20">²⁰Instituto de Astrofísica e Ciencias do Espaco, Universidade do Porto, CAUP, Rua das Estrelas, P-4150-762 Porto, Portugal</footnote>

<footnote n="21">²¹INAF – Osservatorio Astrofísico di Torino, Via Osservatorio 20, I-10025 Pino Torinese, Italy</footnote>

<footnote n="22">²²Space Sciences, Technologies and Astrophysics Research (STAR) Institute, Université de Liège, Allée du 6 Août 19C, B-4000 Liège, Belgium</footnote>

<footnote n="23">²³University of Southern Queensland, Centre for Astrophysics, West Street, Toowoomba, QLD 4350 Australia</footnote>

<footnote n="24">²⁴Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University, Box 43, SE-22100 Lund, Sweden</footnote>

<footnote n="25">²⁵Instituto de Astrofísica de Canarias, E-38200 La Laguna, Tenerife, Spain</footnote>

<footnote n="26">²⁶Departamento de Astrofísica, Universidad de La Laguna, E-38206 La Laguna, Tenerife, Spain</footnote>
<footnote n="27"><sup>27</sup>Institut de Ciencies de l'Espai (ICE, CSIC), Campus UAB, Can Magrans s/n, E-08193 Bellaterra, Spain</footnote>

<footnote n="28"><sup>28</sup>Institut d'Estudis Espacials de Catalunya (IEEC), E-08034 Barcelona, Spain</footnote>

<footnote n="29"><sup>29</sup>Admatis, 5. Kandó Kálmán Street, 3534 Miskolc, Hungary</footnote>

<footnote n="30"><sup>30</sup>Depto. de Astrofísica, Centro de Astrobiología (CSIC-INTA), ESAC campus, E-28692 Villanueva de la Cañada (Madrid), Spain</footnote>

<footnote n="31"><sup>31</sup>Departamento de Física e Astronomía, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, P-4169-007 Porto, Portugal</footnote>

<footnote n="32"><sup>32</sup>Université Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="33"><sup>33</sup>DTU Space, National Space Institute, Technical University of Denmark, Elektrovej 328, DK-2800 Kgs. Lyngby, Denmark</footnote>

<footnote n="34"><sup>34</sup>Institute of Planetary Research, German Aerospace Center (DLR), Rutherfordstrasse 2, P-12489 Berlin, Germany</footnote>

<footnote n="35"><sup>35</sup>Université de Paris, Institut de physique du globe de Paris, CNRS, F-75005 Paris, France</footnote>

<footnote n="36"><sup>36</sup>Centre for Mathematical Sciences, Lund University, SE-22100 Lund, Sweden</footnote>

<footnote n="37"><sup>37</sup>Astrobiology Research Unit, Université de Liège, Allée du 6 Août 19C, B-4000 Liège, Belgium</footnote>

<footnote n="38"><sup>38</sup>Department of Astronomy, Stockholm University, SE-106 91 Stockholm, Sweden</footnote>

<footnote n="39"><sup>39</sup>Leiden Observatory, University of Leiden, PO Box 9513, NL-2300 RA Leiden, the Netherlands</footnote>

<footnote n="40"><sup>40</sup>Department of Space, Earth and Environment, Chalmers University of Technology, Onsala Space Observatory, SE-43992 Onsala, Sweden</footnote>

<footnote n="41"><sup>41</sup>Dipartimento di Fisica, Università degli Studi di Torino, via Pietro Giuria 1, I-10125, Torino, Italy</footnote>

<footnote n="42"><sup>42</sup>Department of Astrophysics, University of Vienna, Tuerkenschanzstrasse 17, A-1180 Vienna, Austria</footnote>

<footnote n="43"><sup>43</sup>Division Technique INSU, CS20330, F-83507 La Seyne sur Mer cedex, France</footnote>

<footnote n="44"><sup>44</sup>Department of Physics, University of Warwick, Gibbet Hill Road, Coventry CV4 7AL, UK</footnote>

<footnote n="45"><sup>45</sup>Science and Operations Department – Science Division (SCI-SC), Directorate of Science, European Space Agency (ESA), European Space Research and Technology Centre (ESTEC), Keplerlaan 1, NL-2201 AZ Noordwijk, the Netherlands</footnote>

<footnote n="46"><sup>46</sup>NASA Ames Research Center, Moffett Field, CA 94035, USA</footnote>

<footnote n="47"><sup>47</sup>Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, 1121 Budapest, Konkoly Thege Miklós út 15-17, Hungary</footnote>

<footnote n="48"><sup>48</sup>ELTE Eötvös Loránd University, Institute of Physics, Pázmány Péter sétány I/A, 1117 Budapest, Hungary</footnote>

<footnote n="49"><sup>49</sup>IMCCE, UMR8028 CNRS, Observatoire de Paris, PSL Univ., Sorbonne Univ., 77 av. Denfert-Rochereau, F-75014 Paris, France</footnote>

<footnote n="50"><sup>50</sup>Institut d'astrophysique de Paris, UMR7095 CNRS, Université Pierre & Marie Curie, 98bis blvd. Arago, F-75014 Paris, France</footnote>

<footnote n="51"><sup>51</sup>Astrophysics Group, Keele University, Staffordshire ST5 5BG, UK</footnote>

<footnote n="52"><sup>52</sup>INAF – Osservatorio Astronomico di Palermo, Piazza del Parlamento 1, I-90134 Palermo, Italy</footnote>

<footnote n="53"><sup>53</sup>INAF – Osservatorio Astronomico di Cagliari, Via della Scienza 5, I-09047 Selargius, Italy</footnote>

<footnote n="54"><sup>54</sup>Komaba Institute for Science, The University of Tokyo, 3-8-1 Komaba, Meguro, Tokyo 153-8902, Japan</footnote>

<footnote n="55"><sup>55</sup>Astrobiology Center, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan</footnote>

<footnote n="56"><sup>56</sup>INAF, Osservatorio Astrofisico di Catania, Via S. Sofia 78, I-95123 Catania, Italy</footnote>

<footnote n="57"><sup>57</sup>P&P Software GmbH, High Tech Center, Taegerwilen, CH-8274, Switzerland</footnote>

<footnote n="58"><sup>58</sup>Institute of Optical Sensor Systems, German Aerospace Center (DLR), Rutherfordstrasse 2, D-12489 Berlin, Germany</footnote>

<footnote n="59"><sup>59</sup>ESTEC, European Space Agency, NL-2201 AZ Noordwijk, the Netherlands</footnote>

<footnote n="60"><sup>60</sup>Center for Astronomy and Astrophysics, Technical University Berlin, Hardenbergstrasse 36, D-10623 Berlin, Germany</footnote>

<footnote n="61"><sup>61</sup>Institut für Geologische Wissenschaften, Freie Universität Berlin, D-12249 Berlin, Germany</footnote>

<footnote n="62"><sup>62</sup>Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, USA</footnote>

<footnote n="63"><sup>63</sup>Department of Aeronautics and Astronautics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA</footnote>

<footnote n="64"><sup>64</sup>ELTE Eötvös Loránd University, Gothard Astrophysical Observatory, 9700 Szombathely, Szent Imre h. u. 112, Hungary</footnote>

<footnote n="65"><sup>65</sup>MTA-ELTE Exoplanet Research Group, 9700 Szombathely, Szent Imre h. u. 112, Hungary</footnote>

<footnote n="66"><sup>66</sup>SETI Institute, Mountain View, CA 94043, USA</footnote>

<footnote n="67"><sup>67</sup>Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="68"><sup>68</sup>Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA</footnote>

<footnote>This paper has been typeset from a TEX/LATEX file prepared by the author.</footnote>