# An unusually low density ultra-short period super-Earth and three mini-Neptunes around the old star TOI-561

G. Lacedelli, L. Malavolta, L. Borsato, G. Piotto, D. Nardiello, A. Mortier, M.

Stalport, A. Collier Cameron, E. Poretti, L. A. Buchhave, et al.

► To cite this version:

G. Lacedelli, L. Malavolta, L. Borsato, G. Piotto, D. Nardiello, et al. An unusually low density ultra-short period super-Earth and three mini-Neptunes around the old star TOI-561. Monthly Notices of the Royal Astronomical Society, 2021, 501, pp.4148-4166. 10.1093/mnras/staa3728. hal-03585883

HAL Id: hal-03585883

https://hal.science/hal-03585883

Submitted on 24 Feb 2022

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

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

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

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

# An unusually low density ultra-short period super-Earth and three mini-Neptunes around the old star TOI-561

G. Lacedelli<sup>1,2,*</sup> L. Malavolta<sup>1,2</sup> L. Borsato<sup>1,2</sup> G. Piotto,<sup>1,2</sup> D. Nardiello<sup>2,3</sup> A. Mortier<sup>1,4,5</sup>

M. Stalport,<sup>6</sup> A. Collier Cameron<sup>1,7</sup> E. Poretti,<sup>8,9</sup> L. A. Buchhave<sup>1,10</sup> M. López-Morales,<sup>11</sup>

V. Nascimbeni<sup>1,2</sup> T. G. Wilson<sup>1,7</sup> S. Udry,<sup>6</sup> D. W. Latham,<sup>11</sup> A. S. Bonomo,<sup>12</sup> M. Damasso<sup>1,12</sup>

X. Dumusque,<sup>6</sup> J. M. Jenkins,<sup>13</sup> C. Lovis,<sup>6</sup> K. Rice<sup>14,15</sup> D. Sasselov,<sup>11</sup> J. N. Winn,<sup>16</sup> G. Andreuzzi,<sup>8,17</sup>

R. Cosentino,<sup>8</sup> D. Charbonneau,<sup>11</sup> L. Di Fabrizio,<sup>8</sup> A. F. Martnez Fiorenzano,<sup>8</sup> A. Ghedina,<sup>8</sup>

A. Harutyunyan,<sup>8</sup> F. Lienhard<sup>1,4</sup> G. Micela,<sup>18</sup> E. Molinari,<sup>19</sup> I. Pagano,<sup>20</sup> F. Pepe,<sup>6</sup> D. F. Phillips,<sup>11</sup>

M. Pinamonti,<sup>12</sup> G. Ricker,<sup>21</sup> G. Scandariato,<sup>20</sup> A. Sozzetti<sup>12</sup> and C. A. Watson<sup>22</sup>

*Affiliations are listed at the end of the paper*

Accepted 2020 November 27. Received 2020 November 26; in original form 2020 September 4

## ABSTRACT

Based on HARPS-N radial velocities (RVs) and TESS photometry, we present a full characterization of the planetary system orbiting the late G dwarf TOI-561. After the identification of three transiting candidates by TESS, we discovered two additional external planets from RV analysis. RVs cannot confirm the outer TESS transiting candidate, which would also make the system dynamically unstable. We demonstrate that the two transits initially associated with this candidate are instead due to single transits of the two planets discovered using RVs. The four planets orbiting TOI-561 include an ultra-short period (USP) super-Earth (TOI-561 b) with period $P_b = 0.45$ d, mass $M_b = 1.59 \pm 0.36 M_⊕$ and radius $R_b = 1.42 \pm 0.07 R_⊕$, and three mini-Neptunes: TOI-561 c, with $P_c = 10.78$ d, $M_c = 5.40 \pm 0.98 M_⊕$, $R_c = 2.88 \pm 0.09 R_⊕$; TOI-561 d, with $P_d = 25.6$ d, $M_d = 11.9 \pm 1.3 M_⊕$, $R_d = 2.53 \pm 0.13 R_⊕$; and TOI-561 e, with $P_e = 77.2$ d, $M_e = 16.0 \pm 2.3 M_⊕$, $R_e = 2.67 \pm 0.11 R_⊕$. Having a density of $3.0 \pm 0.8 \text{ g cm}^{-3}$, TOI-561 b is the lowest density USP planet known to date. Our N-body simulations confirm the stability of the system and predict a strong, anti-correlated, long-term transit time variation signal between planets d and e. The unusual density of the inner super-Earth and the dynamical interactions between the outer planets make TOI-561 an interesting follow-up target.

**Key words:** techniques: photometric – techniques: radial velocities – planets and satellites: composition – planets and satellites: detection – stars: individual: TOI-561 (TIC 377064495, Gaia DR2 3850421005290172416).

## 1 INTRODUCTION

The *Transiting Exoplanet Survey Satellite* (TESS, Ricker et al. 2014) is a NASA all-sky survey designed to search for transiting planets around bright and nearby stars and particularly targeting stars that could reveal planets with radii smaller than Neptune. Since the beginning of its observations in 2018, TESS has discovered more than 66 exoplanets, including about a dozen multiplanet systems (e.g. Dragomir et al. 2019; Dumusque et al. 2019; Günther et al. 2019). Multiplanet systems, orbiting the same star and having formed from the same protoplanetary disc, offer a unique opportunity for comparative planetology. They allow for investigations of the formation and evolution processes, i.e. through studies of relative planet sizes and orbital separations, orbital inclinations relative to the star’s rotation axis, mutual inclination of the orbits, etc. In order to obtain a complete characterization of a system, knowledge of the orbital architecture and the bulk composition of the planets are essential. To obtain such information, transit photometry needs to be

combined with additional techniques that allow for the determination of the planetary masses, i.e. radial velocity (RV) follow-up or transit time variation (TTV) analysis. Up to now, the large majority of known planetary systems have been discovered by the *Kepler* space telescope (Borucki et al. 2010), which has led to an unprecedented knowledge of the ensemble properties of multiple systems (e.g. Latham et al. 2011; Millholland, Wang & Laughlin 2017; Weiss et al. 2018), their occurrence rate (e.g. Fressin et al. 2013), and their dynamical configurations (e.g. Lissauer et al. 2011; Fabrycky et al. 2014). However, many of the *Kepler* targets are too faint for RV follow-up, so most of the planets do not have a mass measurement, preventing a comprehensive understanding of their properties, and of the planetary system. Thanks to the TESS satellite, which targets brighter stars, an increasing number of candidates suitable for spectroscopic follow-up campaigns are being discovered. These new objects will increase the number of well-characterized systems, and will provide a valuable observational counterpart to the theoretical studies on the formation and evolution processes of planetary systems (e.g. Morbidelli et al. 2012; Baruteau et al. 2014, 2016; Davies et al. 2014; Helled et al. 2014; Raymond et al. 2014). In this paper, we combine TESS photometry (Section 2.1)

<footnote>* E-mail: gaia.lacedelli@phd.unipd.it</footnote>
<table><caption>**Table 1.** Astrometric and photometric parameters of TOI-561.</caption><thead><tr><th>Property</th><th>Value</th><th>Source</th></tr></thead><tbody><tr><td colspan="3"><em>Other target identifiers</em></td></tr><tr><td>TIC</td><td>377064495</td><td>A</td></tr><tr><td>Gaia DR2</td><td>3850421005290172416</td><td>B</td></tr><tr><td>2MASS</td><td>J09524454 + 0612589</td><td>C</td></tr><tr><td colspan="3"><em>Astrometric parameters</em></td></tr><tr><td>RA (J2015.5; h:m:s)</td><td>09:52:44.44</td><td>B</td></tr><tr><td>Dec (J2015.5; d:m:s)</td><td>06:12:57.97</td><td>B</td></tr><tr><td>μ<sub>α</sub> (mas yr<sup>-1</sup>)</td><td>-108.432 ± 0.088</td><td>B</td></tr><tr><td>μ<sub>δ</sub> (mas yr<sup>-1</sup>)</td><td>-61.511 ± 0.094</td><td>B</td></tr><tr><td>Systemic velocity (km s<sup>-1</sup>)</td><td>79.54 ± 0.56</td><td>B</td></tr><tr><td>Parallax<sup>a</sup> (mas)</td><td>11.6768 ± 0.0672</td><td>B</td></tr><tr><td>Distance (pc)</td><td>85.80<sup>+0.50</sup><sub>-0.49</sub></td><td>D</td></tr><tr><td colspan="3"><em>Photometric parameters</em></td></tr><tr><td>TESS (mag)</td><td>9.527 ± 0.006</td><td>A</td></tr><tr><td>Gaia (mag)</td><td>10.0128 ± 0.0003</td><td>B</td></tr><tr><td>V (mag)</td><td>10.252 ± 0.006</td><td>A</td></tr><tr><td>B (mag)</td><td>10.965 ± 0.082</td><td>A</td></tr><tr><td>J (mag)</td><td>8.879 ± 0.020</td><td>C</td></tr><tr><td>H (mag)</td><td>8.504 ± 0.055</td><td>C</td></tr><tr><td>K (mag)</td><td>8.394 ± 0.019</td><td>C</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>W3 (mag)</td><td>8.375 ± 0.023</td><td>E</td></tr><tr><td>W4 (mag)</td><td>7.971 ± 0.260</td><td>E</td></tr><tr><td>Spectral type</td><td>G9V</td><td>F</td></tr></tbody></table>

Notes. A) TESS Input Catalogue Version 8 (TICv8, Stassun et al. 2018).

B) Gaia DR2 (Gaia Collaboration et al. 2018).

C) Two Micron All Sky Survey (2MASS, Cutri et al. 2003).

D) Bailer-Jones et al. (2018).

E) Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010).

F) Based on Pecaut & Mamajek (2013), assuming Gaia DR2, Johnson, 2MASS and WISE colour indexes.

<sup>a</sup>Gaia DR2 parallax is corrected by +50 ± 7 μas (with the error added in quadrature) as suggested by Khan et al. (2019).

and high-precision RVs gathered with the HARPS-N spectrograph (Section 2.2) to characterize the multiplanet system orbiting the star TOI-561. The TESS pipeline identified three candidate planetary signals, namely an ultra-short period (USP) candidate (P ~ 0.45 d), and two additional candidates with periods of ~10.8 and ~16.4 d. We determined the stellar properties (Section 3) using three independent methods. Based on our activity analysis, we concluded that TOI-561 is an old, quiet star, and therefore quite appropriate for the study of a complex planetary system. After assessing the planetary nature of the transit-like features (Section 4), we performed a series of analysis – with the tools described in Section 5 – to determine the actual system configuration (Section 6). We further address the robustness of our final solution based on a comparison with other possible models (Section 7). We finally compare the resulting planetary densities with the distribution of known planets in the mass–radius diagram and we predict the expected TTV signal for the planets in the system (Section 8).

# 2 OBSERVATIONS

## 2.1 TESS photometry

TOI-561 was observed by TESS in 2-min cadence mode during observations of sector 8, between 2019 February 2 and 27. The astrometric and photometric parameters of the star are listed in Table 1. Considering the download time, and the loss of 3.26 d of data

due to an interruption in communications between the instrument and the spacecraft that occurred during sector 8,<fnref n="1" /> a total of 20.22 d of science data were collected. The photometric observations for TOI-561 were reduced by the Science Processing Operations Center (SPOC) pipeline (Jenkins et al. 2016; Jenkins 2020), which detected three candidate planetary signals, with periods of 10.8 (TOI-561.01), 0.4 (TOI-561.02), and 16.4 d (TOI-561.03), respectively. The pipeline identified 55 transits of TOI-561.02, 2 transits of TOI-561.01, and 2 transits of TOI-561.03, with depths of 290, 1207, and 923 ppm and signal-to-noise-ratios (SNR) of 10.0, 9.8, and 9.2, respectively. For our photometric analysis, we used the light curve based on the Pre-search Data Conditioning Simple Aperture Photometry (PDCSAP, Smith et al. 2012; Stumpe et al. 2012, 2014). We downloaded the 2-min cadence PDCSAP light curve from the Mikulski Archive for Space Telescopes (MAST),<fnref n="2" /> and removed all the observations encoded as NaN or flagged as bad-quality (DQUALITY > 0) points by the SPOC pipeline.<fnref n="3" /> We performed outliers rejection by doing a cut at 3σ for positive outliers and 5σ (i. e. larger than the deepest transit) for negative outliers. We removed the low-frequency trends in the light curve using the biweight time-windowed slider implemented in the WOTAN package (Hippke et al. 2019), with a window of 1.5 d, and masking the in-transit points to avoid modifications of the transit shape. In order to obtain an independent confirmation of the signals detected in the TESS light curve, we performed an iterative transit search on the detrended light curve using the Transit Least Squares (TLS) algorithm (Hippke & Heller 2019). The first three significant identified signals nicely matched the TESS suggested periods (P<sub>TLS</sub> = 10.78 d, 0.44 d, 16.28 d). In addition, we also extracted the 30-min cadence light curve from the TESS Full-Frame Images (FFIs) using the PATHOS pipeline (Nardiello et al. 2019), in order to obtain an independent confirmation of the detected signals (Section 4).

## 2.2 HARPS-N spectroscopy

We collected 82⁴ spectra using HARPS-N at the Telescopio Nazionale Galileo (TNG), in La Palma (Cosentino et al. 2012, 2014), with the goal of precisely determining the masses of the three candidate planets and to search for additional planets. The observations started on 2019 November 17 and ended on 2020 June 13 with an interruption between the end of March and the end of April due to the shut down of the TNG because of Covid-19. In order to precisely characterize the signal of the USP candidate, we collected 6 points per night on 2020 February 4 and 6, thus, covering the whole phase curve of the planet, and two points per night (when weather allowed) during the period of maximum visibility of the target (2020 February–March). The exposure time was set to 1800 s, which resulted in an SNR at 550 nm of 77 ± 20 (median ± standard deviation) and a measurement uncertainty of 1.2 ± 0.6 m s<sup>-1</sup>. We reduced the data using the standard HARPS-N Data Reduction Software (DRS) using a G2 flux template (the closest match to the spectral type of our target) to correct for variations in the flux distribution as a function of the wavelength, and a G2 binary mask to compute the cross-correlation function (CCF, Baranne et al. 1996;

<footnote n="1">¹See TESS Data Release Notes: Sector 8, DR10 (https://archive.stsci.edu/terms/tess_drn.html).</footnote>

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

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

<footnote n="4">⁴62 spectra were collected within the Guaranteed Time Observations (GTO) time (Pepe et al. 2013), while the remaining 20 spectra were collected within the A40_TAC23 program.</footnote>
<table><caption>**Table 2.** HARPS-N RV 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>BIS<br>(m s<sup>-1</sup>)</th><th>FWHM<br>(km s<sup>-1</sup>)</th><th>V<sub>asy</sub></th><th>ΔV<br>(km s<sup>-1</sup>)</th><th>log R'<sub>HK</sub></th><th>H α<br>(dex)</th></tr></thead><tbody><tr><td>2458804.70779</td><td>79700.63</td><td>1.27</td><td>-39.98</td><td>6.379</td><td>0.048</td><td>-0.039</td><td>-5.005</td><td>0.203</td></tr><tr><td>2458805.77551</td><td>79703.74</td><td>0.97</td><td>-36.25</td><td>6.380</td><td>0.049</td><td>-0.036</td><td>-4.984</td><td>0.200</td></tr><tr><td>2458806.76768</td><td>79701.71</td><td>1.05</td><td>-31.81</td><td>6.378</td><td>0.045</td><td>-0.033</td><td>-5.000</td><td>0.200</td></tr><tr><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.

Pepe et al. 2002). All the observations were gathered with the second fibre of HARPS-N illuminated by the Fabry–Perot calibration lamp to correct for the instrumental RV drift, except for the night of 2020 May 31. This observation setting prevented us from using the second fibre to correct for Moon contamination. However, we note that the difference between the systemic velocity of the star and the Moon is always greater than 15 km s⁻¹, therefore preventing any contamination of the stellar CCF (as empirically found by Malavolta et al. 2017a and subsequently demonstrated through simulations by Roy et al. 2020), as the average full width at half-maximum (FWHM) of the CCF for TOI-561 is 6.380 ± 0.004 km s⁻¹. The RV data with their 1σ uncertainties and the associated activity indices (see Section 3.3 for more details) are listed in Table 2. Before proceeding with the analysis, we removed from the total data set five RV measurements, with associated errors greater than 2.5 m s⁻¹ from spectra with SNR < 35, that may affect the accuracy of our results. The detailed procedure performed to identify these points is described in Appendix B1.

# 3 STELLAR PARAMETERS

## 3.1 Photospheric parameters

We derived the photospheric stellar parameters using three different techniques: the curve-of-growth approach, spectral synthesis match, and empirical calibration. The first method minimizes the trend of iron abundances (obtained from the equivalent width, EW, of each line) with respect to excitation potential and reduced EW respectively, to obtain the effective temperature and the micro-turbulent velocity, $\xi_t$. The gravity log g is obtained by imposing the same average abundance from neutral and ionized iron lines. We obtained the EW measurements using ARESv2⁵ (Sousa et al. 2015). We used the local thermodynamic equilibrium (LTE) code MOOG<fnref n="6" /> (Sneden 1973) for the line analysis, together with the ATLAS9 grid of stellar model atmosphere from Castelli & Kurucz (2003). The whole procedure is described in more detail in Sousa (2014). We performed the analysis on a co-added spectrum (SNR > 600), and after applying the gravity correction from Mortier et al. (2014) and adding systematic errors in quadrature (Sousa et al. 2011), we obtained $T_{eff} = 5346 \pm 69$ K, log g = 4.60 ± 0.12, [Fe/H] = -0.40 ± 0.05, and $\xi_t = 0.78 \pm 0.08$ km s⁻¹. The spectral synthesis match was performed using the Stellar Parameters Classification tool (SPC, Buchhave et al. 2012, 2014). It determines effective temperature, surface gravity, metallicity, and line broadening by performing a cross-correlation of the observed spectra with a library of synthetic spectra, and interpolating the correlation peaks to determine the best-matching parameters. For technical reasons, we

<table><caption>**Table 3.** Derived astrophysical stellar parameters.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td><i>T</i><sub>eff</sub><sup>a</sup></td><td>5372 ± 70</td><td>K</td></tr><tr><td>log <i>g</i><sup>a</sup><sub>spec</sub></td><td>4.50 ± 0.12</td><td>–</td></tr><tr><td>[Fe/H]<sup>a</sup><sub>spec</sub></td><td>-0.40 ± 0.05</td><td>–</td></tr><tr><td><i>T</i><sub>eff</sub><sup>b</sup></td><td>5455<sup>+65</sup><sub>-47</sub></td><td>K</td></tr><tr><td>log <i>g</i><sup>b</sup></td><td>4.47 ± 0.01</td><td>–</td></tr><tr><td>[Fe/H]<sup>b</sup></td><td>-0.33<sup>+0.10</sup><sub>-0.05</sub></td><td>–</td></tr><tr><td><i>R</i><sub>*</sub></td><td>0.849 ± 0.007</td><td>R<sub>☉</sub></td></tr><tr><td><i>M</i><sub>*</sub></td><td>0.785 ± 0.018</td><td>M<sub>☉</sub></td></tr><tr><td><i>ρ</i><sub>*</sub></td><td>1.285 ± 0.040</td><td><i>ρ</i><sub>☉</sub></td></tr><tr><td><i>ρ</i><sub>*</sub></td><td>1.809 ± 0.056</td><td>g cm<sup>-3</sup></td></tr><tr><td><i>A</i><sub>V</sub></td><td>0.12<sup>+0.08</sup><sub>-0.06</sub></td><td>mag</td></tr><tr><td><i>v sin i</i></td><td>&lt;2</td><td>km s<sup>-1</sup></td></tr><tr><td><i>age</i><sup>c</sup></td><td>&gt;5</td><td>Gyr</td></tr><tr><td>log R'<sub>HK</sub></td><td>-5.003 ± 0.012</td><td>–</td></tr><tr><td>[Na/H]</td><td>-0.28 ± 0.06</td><td>–</td></tr><tr><td>[Mg/H]</td><td>-0.17 ± 0.05</td><td>–</td></tr><tr><td>[Si/H]</td><td>-0.22 ± 0.05</td><td>–</td></tr><tr><td>[Ca/H]</td><td>-0.27 ± 0.06</td><td>–</td></tr><tr><td>[Ti/H]</td><td>-0.12 ± 0.03</td><td>–</td></tr><tr><td>[Cr/H]</td><td>-0.33 ± 0.08</td><td>–</td></tr><tr><td>[Ni/H]</td><td>-0.37 ± 0.04</td><td>–</td></tr></tbody></table>

Notes. <sup>a</sup> Weighted average of the three spectroscopic methods.

<sup>b</sup> Value inferred from the isochrone fit.

<sup>c</sup> Conservative lower limit.

ran the SPC on the 62 GTO spectra only<fnref n="7" />: the SNR is so high that the spectra are anyway dominated by systematic errors, and including the A40TAC_23 spectra would not change the results. We averaged the values measured for each exposure, and we obtained $T_{eff} = 5389 \pm 50$ K, log g = 4.49 ± 0.10, [M/H] = -0.36 ± 0.08, and $v sin i < 2$ km s⁻¹. We finally used CCFpams,<fnref n="8" /> a method based on the empirical calibration of temperature, metallicity, and gravity on several CCFs obtained with subsets of stellar lines with different sensitivity to temperature (Malavolta et al. 2017b). We obtained $T_{eff} = 5293 \pm 70$ K, log g = 4.50 ± 0.15, and [Fe/H] = -0.40 ± 0.05, after applying the same gravity and systematic corrections as for the EW analysis. We list the final spectroscopic adopted values, i. e. the weighted averages of the three methods, in Table 3. From the co-added HARPS-N spectrum, we also derived the chemical abundances for several refractory elements (Na, Mg, Si, Ca, Ti, Cr, Ni). We used the ARES + MOOG method assuming LTE, as described earlier. The reference for solar values was taken from Asplund et al. (2009), and all values in Table 3 are given relative to the Sun. Details on the method and line lists are described in Adibekyan et al. (2012).

<footnote n="7">⁷SPC runs on a server with access to GTO data only, and the required technical effort to enable the use of A40_TAC23 data, complicated by the global Covid-19 sanitary emergency, was not justified by the negligible scientific gain.</footnote>

<footnote n="8">⁸Available at https://github.com/LucaMalavolta/CCFpams</footnote>

<footnote n="5">⁵Available at http://www.astro.up.pt/~sousasag/ares/</footnote>

<footnote n="6">⁶Available at http://www.as.utexas.edu/~chris/moog.html</footnote>
and Mortier et al. (2013). This analysis shows that this iron-poor star is alpha-enhanced. Using the average abundances of magnesium, silicon, and titanium to represent the alpha-elements and the iron abundance from the ARES + MOOG method (for consistency), we find that [α/Fe] = 0.23.

## 3.2 Mass, radius, and density of the star

For each set of photospheric parameters, we determined the stellar mass and radius using *isochrones* (Morton 2015), with posterior sampling performed by *MultiNest* (Feroz & Hobson 2008; Feroz, Hobson & Bridges 2009; Feroz et al. 2019). We provided as input the parallax of the target from the Gaia DR2 catalogue, after adding an offset of +50 ± 7 µas (with the error added in quadrature to the parallax error) as suggested by Khan et al. (2019), plus the photometry from the TICv8, 2MASS, and *WISE* (Table 1). We used two evolutionary models, [the MESA Isochrones & Stellar Tracks (MIST), Paxton et al. 2011; Choi et al. 2016; Dotter 2016] and the Dartmouth Stellar Evolution Data base (Dotter et al. 2008). For all methods, we assumed $\sigma_{T_{\text{eff}}} = 70$ K, $\sigma_{\text{log g}} = 0.12$, $\sigma_{\text{[Fe/H]}} = 0.05$ (except for SPC, where we kept the original error of 0.08) as a good estimate of the systematic errors regardless of the internal error estimates, to avoid favouring one technique over the others when deriving the stellar mass and radius. We also imposed an upper limit on the age of 13.8 Gyr, i. e. the age of the Universe (Planck Collaboration I 2018). From the mean and standard deviation of all the posterior samplings, we obtained $M_* = 0.785 \pm 0.018 M_{\odot}$ and $R_* = 0.849 \pm 0.007 R_{\odot}$. We derived the stellar density $\rho_* = 1.285 \pm 0.040 \rho_{\odot}$ ($\rho_* = 1.809 \pm 0.056 \text{ g cm}^{-3}$) directly from the posterior distributions of $M_*$ and $R_*$. We summarize the derived astrophysical parameters of the star in Table 3, which also reports temperature, gravity, and metallicity obtained from the posteriors distributions resulting from the *isochrone* fit. A lower limit on the age of ~10 Gyr is obtained considering the 15.86th-percentile of the distribution of the combined posteriors, as for the other parameters. We note, however, that an *isochrone* fit performed through *EXO-FASTv2* (Eastman et al. 2019), assuming the photometric parameters in Table 1 and the spectroscopic parameters in Table 3, using only the MIST evolutionary set, returned a lower limit on the age of 5 Gyr, while all the other parameters were consistent with the results quoted in Table 3. Thus, we decided to assume 5 Gyr as a conservative lower limit for the age of the system. The old stellar age and the sub-solar metallicity suggest that TOI-561 may belong to an old Galactic population, a hypothesis that is also supported by our kinematic analysis. In fact, we derived the Galactic space velocities using the astrometric properties reported in Table 1. For the calculations, we used the *ASTROPY* package, and we assumed the Gaia DR2 RV value of 79.54 km s⁻¹, obtaining the heliocentric velocity components (*U*, *V*, *W*) = (−60.0, −70.8, 16.7) km s⁻¹, in the directions of the Galactic Centre, Galactic rotation, and north Galactic pole, respectively. The derived *UVW* velocities point toward a thick-disc star, as confirmed by the probability membership derived following Bensby, Feltzing & Oey (2014), that implies a ~70 per cent probability that the star belongs to the thick disc, a ~29 per cent probability of being a thin-disc star, and a ~0.0004 per cent probability of belonging to the halo.

## 3.3 Stellar activity

The low value of the log $R'_{\text{HK}}$ index (−5.003 ± 0.012), derived using the calibration by Lovis et al. (2011) and assuming $B-V = 0.71$, low value of the log ($R_{\text{HK}}$) index that TOI-561 is a relatively quiet star. Given its distance of ~86 pc, the lack of

interstellar absorption near the Na D doublet in the HARPS-N co-added spectrum, and the total extinction in the V band from the isochrone fit (0.1 mag), we do not expect any significant effect of the interstellar medium on the log $R'_{\text{HK}}$ index (Fossati et al. 2017). Nevertheless, it is important to check whether the star is showing any sign of activity in all the activity diagnostics at our disposal. In addition to the log $R'_{\text{HK}}$ index, FWHM, and bisector span (BIS) computed by the HARPS-N DRS, we included in our analysis the $V_{\text{asy}}$ (Figueira et al. 2013) and $\Delta V$ (Nardetto et al. 2006) asymmetry indicators, as implemented by Lanza et al. (2018), and the chromospheric activity indicator $H\alpha$ (Gomes da Silva et al. 2011).

The Generalized Lomb-Scargle (GLS, Zechmeister & Kürster 2009) periodograms of the above-mentioned indexes, computed within the frequency range 0.0005–0.5 d⁻¹, i. e. 2–2000 d, are shown in Fig. 1, together with the periodograms of the RVs and *TESS* photometry. For each periodogram, we also report the power threshold corresponding to a False Alarm Probability (FAP) of 1 and 0.1 per cent, computed with a bootstrap approach. The periodogram of the RVs reveals the presence of significant peaks at ≃25, ≃180, ≃10 d (corresponding to one of the transiting planet candidates), and ≃78 d, ordered decreasingly according to their power. None of these peaks has a counterpart in the activity diagnostics here considered, as no signals with a FAP lower than 2.4 per cent can be identified, strongly supporting that the signals in the RVs are not related to stellar activity. We note that the GLS periodogram of the *TESS* light curve identified a periodicity around 3.5 d with an amplitude of 0.13 ppt and a power of 0.014, that is, above the 0.1 per cent FAP threshold. However, it is unlikely that such variability is associated with stellar activity, since a rotational period of just a few days would be extremely atypical for a star older than 1 Gyr (e.g. Douglas et al. 2019), and in contrast with the lack of any signal in all the other above-mentioned activity indicators. Indeed, the rotational period estimated from the log $R'_{\text{HK}}$ using the calibrations of Noyes et al. (1984) and Mamajek & Hillenbrand (2008) supports this assertion, indicating a value around 33 d. We note that this value of the rotational period should be considered as a rough estimate, also because these calibrations are not well tested for old and alpha-enhanced stars like TOI-561. Further evidence against a ~3.5 d rotational period is provided by the low value of the $v \sin i$ (<2 km s⁻¹), that suggests a rotational period >21.5 d, assuming the stellar radius listed in Table 3 and an inclination of 90°. In any case, we verified with a periodogram analysis that our light curve flattening procedure correctly removed the here identified signal at 3.5 d.

In addition, we performed an auto correlation analysis, following the prescription by McQuillan, Aigrain & Mazeh (2013), on the *TESS* light curve (with the transits filtered out), and the ASAS-SN *V* and *g* photometry (Shappee et al. 2014; Kochanek et al. 2017), after applying a 5σ filtering, but no significant periodicity could be identified. A periodogram analysis of the ASAS-SN light curves in each band, either by taking the full data set or by analysing each observing season individually, confirmed these results.

In conclusion, if any activity is present, its signature must be below 0.8 ppt in the short period (rotationally induced activity, <30 d), and 20 ppt in the long term period (magnetic cycles, >100 d), from the RMS of *TESS* and ASAS-SN photometry, respectively. Incidentally, the former is close to the photometric variations of the Sun during the minimum at the end of Solar Cycle 25, when the Sun also reached a log $R'_{\text{HK}}$ very close to the one measured for TOI-561 (Collier Cameron et al. 2019; Milbourne et al. 2019). By comparing our target to the Sun, and in general by taking into account the results of Isaacson & Fischer (2010), it is expected that the contribution to the
<figure><img src="image_3.png" /><figcaption>**Figure 1.** GLS periodogram of the RVs, the *TESS* photometry (PDCSAP), and the spectroscopic activity indexes under analysis. The main peak of each periodogram is highlighted with an orange vertical line. The grey vertical lines represent the signals corresponding to the transit-like signals with periods 10.8 and 16.3 d, and the additional signals identified in the RVs (Section 6) at ≈25, ≈78, and ≈180 d. The dashed and dotted horizontal lines show the 1 and 0.1 per cent FAP levels, respectively. The *TESS* periodogram shows a series of peaks below 10 d, unlikely to be associated with stellar activity given the old age of the star. The FWHM and the log $R'_{HK}$ periodograms have the main peak at 244 and 220 d, respectively, so there is no correspondence with the 180 d signal. Moreover, both of them are below the 1 per cent FAP. The bottom panel shows the window function of the data.</figcaption></figure>

RVs due to the magnetic activity of our star is likely below $1-2 \text{ m s}^{-1}$. Since this value is quite close to the median internal error of our RVs, no hint of the rotational period is provided by either the photometry or the spectroscopic activity diagnostics, and the low activity level is consistent with our derived stellar age (>5 Gyr), we do not include

any activity contributions in the remaining of our analysis, except for an uncorrelated jitter term ($\sigma_{\text{jitter}}$).

## 4 RULING OUT FALSE-POSITIVE SCENARIOS

Previous experience with *Kepler* shows that candidates in multiple systems have a much lower probability of being false positives (Latham et al. 2011; Lissauer et al. 2012). Nevertheless, it is always appropriate to perform a series of checks in order to exclude the possibility of a false positive.

We notice that the star has a good astrometric Gaia DR2 solution (Gaia Collaboration et al. 2018), with zero excess noise and a renormalized unit weight error (RUWE) of 1.1, indicating that the single-star model provides a good fit to the astrometric observations. This likely excludes the presence of a massive companion that could contribute to the star's orbital motion in the Gaia DR2 astrometry, a fact that agrees with the absence of long-term trends in our RVs (see Section 6.1).

Moreover, the overall RV variation below $25 \text{ m s}^{-1}$ and the shape of the CCFs of our HARPS-N spectra exclude the eclipsing binary (EB) scenario, which would be the most likely alternative explanation for the USP planet.

A further confirmation comes from the speckle imaging on the Southern Astrophysical Research (SOAR) telescope that Ziegler et al. (2020) performed on some of the *TESS* planet candidate hosts. According to their analysis (see Tables 3 and 6 therein), no companion is detected around TOI-561 (being the resolution limit for the star 0.041 arcsec, and the maximum detectable $\Delta\text{mag}$ at separation of 1 arcsec 4.76 mag). Still, the 21 arcsec *TESS* pixels and the few-pixels wide point spread function (PSF) can cause the light from neighbours over an arc-minute away to contaminate the target light curve. In the case of neighbouring EBs, eclipses can be diluted and mimic shallow planetary transits. For example, events at ~1 mmag level as in TOI-561.01 and TOI-561.03 can be mimicked by a nearby EB within the *TESS* aperture with a 0.5 per cent eclipse, but no more than 7 mag fainter. This condition is not satisfied in our case, as the only three sources within 100 arcsec from TOI-561 are all fainter than $T = 19.25$ mag and at a distance greater than 59 arcsec, according to the Gaia DR2 catalogue.

An independent confirmation was provided by the analysis of the in-/out-of-transit difference centroids on the *TESS* FFIs (Fig. 2), adopting the procedure described in Nardiello et al. (2020). The analysis of the in-/out-of-transit stacked difference images confirms that, within a box of $10 \times 10 \text{ pixels}^2$ (~200 $\times$ 200 arcsec$^2$) centred on TOI-561, the transit events associated with candidates .01 and .03 occur on our target star, while candidate .02 has too few in-transit points in the 30-min cadence images for this kind of analysis – in any case, its planetary nature will be confirmed by the RV signal of TOI-561 in Section 6.

Finally, in order to exclude the possibility that the transit-like features were caused by instrumental artefacts, we performed some additional checks on the light curve. We visually inspected the FFIs to spot possible causes (including instrumental effects) inducing transit-like features, and we could not find any. We re-extracted the short cadence light curve using the PYTHON package *LIGHTKURVE*<fnref n="9" /> (Lightkurve Collaboration et al. 2018) with different photometric masks and apertures, and we corrected them by using the *TESS* Cotrending Basis Vectors (CBVs); the final results were in agreement with the *TESS*-released PDCSAP light curve. We checked for

<footnote n="9">⁹https://github.com/KeplerGO/lightkurve</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** In-/out-of-transit difference centroid analysis of the transit events associated with the candidates TOI-561.01 (transits 2 and 3) and TOI-561.03 (transits 1 and 4). The star is centred at (0,0), and the grey circles are all the other stars in the Gaia DR2 catalogue, with dimension proportional to their apparent magnitude.</figcaption></figure>

systematics in every light curve pixel, and we found none. Ultimately, we checked for correlations between the flux, the local background, the (X,Y)-position from the PSF-fitting, and the FWHM, with no results. Therefore, we conclude that all the transit-like features in the light curve are real and likely due to planetary transits.

# 5 DATA ANALYSIS TOOLS

We performed the analysis presented in the next sections using *PyORBIT*<fnref n="10" /> (Malavolta et al. 2016, 2018), a convenient wrapper for the analysis of transit light curves and RVs.

In the analysis of the light curve, for each planet we fitted the central time of transit ($T_0$), period ($P$), planetary to stellar radius ratio ($R_p/R_*$), and impact parameter $b$. In order to reduce computational time, we set a narrow, but still uninformative, uniform prior for period and time of transit, as defined by a visual inspection. We fitted a common value for the stellar density $\rho_*$, imposing a Gaussian prior based on the value from Table 3. We included a quadratic limb-darkening law with Gaussian priors on the coefficients $u_1, u_2$, obtained through a bilinear interpolation of limb darkening profiles by Claret (2018).<fnref n="11" /> We initially calculated the standard errors on $u_1, u_2$ using a Monte Carlo approach that takes into account the errors on $T_{\text{eff}}$ and $\log g$ as reported in Table 3, obtaining $u_1 = 0.393 \pm 0.007$ and $u_2 = 0.204 \pm 0.001$. We, however, decided to conservatively increase the error on both coefficients to 0.05. In the fit, we employed the parametrization ($q_1, q_2$) introduced by Kipping (2013). Finally, we included a jitter term to take into account possible *TESS* systematics and short-term stellar activity noise. We assumed uniform, uninformative priors for all the other parameters, although the prior on the stellar density will inevitably affect the other orbital parameters. All the transit models were computed with the *BATMAN*

package (Kreidberg 2015), with an exposure time of 120 s and an oversampling factor of 10 (Kipping 2010).

In the analysis of the RVs, we allowed the periods to span between 2 and 200 d (i. e. the time span of our data set) for the non-transiting planets, while we allowed the semi-amplitude $K$ to vary between 0.01 and $100 \text{ m s}^{-1}$ for all the candidate planets. These two parameters were explored in the logarithmic space. For the transiting candidates, we used the results from the photometric fit (see Appendix A) to impose Gaussian priors on period and time of transit on RV analysis alone, while using the same uninformative priors as for the photometric fit when including the photometric data as well.

For all the signals except the USP candidate, we assumed eccentric orbits with a half-Gaussian zero-mean prior on the eccentricity (with variance 0.098) according to Van Eylen et al. (2019), unless stated otherwise.

We computed the Bayesian evidence using the *MultiNest* nested-sampling algorithm (Feroz & Hobson 2008; Feroz et al. 2009, 2019) with the *PYTHON* wrapper *pyMultiNest* (Buchner et al. 2014). In the specific case of the joint light curve and RV analysis (Section 7), we employed the *dynesty* nested-sampling algorithm (Skilling 2004, 2006; Speagle 2020), which allowed for the computation of the Bayesian evidence in a reasonable amount of time thanks to its easier implementation of the multiprocessing mode. We performed a series of test on a reduced data set, and we verified that the two algorithms provided consistent results with respect to each other. For all the analyses, we assumed 1000 live points and a sampling efficiency of 0.3, including a jitter term for each data set considered in the model.

Global optimization of the parameters was performed using the differential evolution code *PyDE*.<fnref n="12" /> The output parameters were used as a starting point for the Bayesian analysis performed with the *EMCEE* package (Foreman-Mackey et al. 2013), a Markov chain Monte Carlo (MCMC) algorithm with an affine invariant ensemble sampler (Goodman & Weare 2010). We ran the chains with $2n_{\text{dim}}$ walkers, where $n_{\text{dim}}$ is the dimensionality of the model, for a number of steps adapted to each fit, checking the convergence with the Gelman-Rubin statistics (Gelman & Rubin 1992), with a threshold value of $\hat{R} = 1.01$. We also performed an autocorrelation analysis of the chains: if the chains were longer than 100 times the estimated autocorrelation time and this estimate changed by less than 1 per cent, we considered the chains as converged. In each fit, we conservatively set the burn-in value as a number larger than the convergence point as just defined, and we applied a thinning factor of 100.

# 6 UNVEILING THE SYSTEM ARCHITECTURE

## 6.1 Planetary signals in the RV data

Before proceeding with a global analysis, we checked whether we could independently recover the signals identified by the *TESS* pipeline (Section 2.1) in our RV data only. The periodogram analysis of the RVs in Section 3.3 highlighted the presence of several peaks not related to the stellar activity. In particular, an iterative frequency search, performed subtracting at each step the frequency values previously identified, supplied the frequencies $f_1 = 0.039 \text{ d}^{-1}$ ($P_1 \simeq 25.6$ d), $f_2 = 0.006 \text{ d}^{-1}$ or $0.013 \text{ d}^{-1}$ ($P_2 \simeq 170$ d or $\simeq 78$ d) with the two frequencies being related to each other (i. e., removing one of them implies the vanishing of the other one),

<footnote n="10">¹⁰https://github.com/LucaMalavolta/PyORBIT, version 8.1</footnote>

<footnote n="11">¹¹https://vizier.u-strasbg.fr/viz-bin/VizieR?-source = J/A + A/618/A20</footnote>

<footnote n="12">¹²https://github.com/hpparvi/PyDE</footnote>
<table><caption>**Table 4.** Logarithmic Bayesian evidences for the different models under exam. Model 0 corresponds to the model with no additional RVs signal other than the signals from the three transiting candidates i. e. TOI-561.01, .02, .03. Model 1, 2, and 3 correspond to the models with the three transiting candidates plus one, two, and three additional planets, respectively. All the values are expressed with respect to Model 0. We note that the reported errors, as obtained from the nested sampling algorithm, are likely underestimated (Nelson et al. 2020).</caption><thead><tr><th></th><th>Model 0</th><th>Model 1</th><th>Model 2</th><th>Model 3</th></tr></thead><tbody><tr><td>ln Z</td><td>0.0 ± 0.1</td><td>13.4 ± 0.2</td><td>26.1 ± 0.2</td><td>28 ± 0.2</td></tr></tbody></table>

$f_3 = 0.093 \text{ d}^{-1}$ ($P_3 \simeq 10.8 \text{ d}$, corresponding to the TOI-561.01 candidate), and $f_4 = 2.239 \text{ d}^{-1}$ ($P_4 \simeq 0.45 \text{ d}$, corresponding to the TOI-561.02 candidate). After removing these four signals, no other clear dominant frequency emerged in the residuals. Since any attempt to perform a fit of the RVs to characterize the transiting candidates without accounting for additional dominant signals would lead to unreliable results, we decided to test the presence of additional planets in a Bayesian framework. We considered four models, the first one (Model 0) assuming the three transiting candidates only, i. e. TOI-561.01, .02, .03, and then including an additional planet in each of the successive models, i. e. TOI-561.01, .02, .03 plus one (Model 1), two (Model 2), and three (Model 3) additional signals, respectively. We computed the Bayesian evidence for each model using the MultiNest nested-sampling algorithm, following the prescriptions as specified in Section 5. We report the obtained values in Table 4. According to this analysis, we concluded that the model with two additional signals, i. e. Model 2 (with no trend), is strongly favoured over the others, with a difference in the logarithmic Bayes factor $2 \Delta \ln Z > 10$ (Kass & Raftery 1995), both compared to the case with one or no additional signals. In the case of a third additional signal (Model 3), the difference with respect to the two-signal model was less than 2, indicating that there was no strong evidence to favour this more complex model over the simpler model with two additional signals only (Kass & Raftery 1995). We repeated the analysis first including a linear and then a quadratic trend in each of the four models. In all cases, the Bayesian evidence systematically disfavoured the presence of any trend.<fnref n="13" />

The first additional signal was associated with a candidate with $f \simeq 0.04 \text{ d}^{-1}$ ($P \simeq 25.6 \text{ d}$), which corresponds to the strongest peak in the RVs periodogram. Concerning the second additional signal, the MultiNest run highlighted the presence of two clusters of solutions, peaked at about $f = 0.013$ or $0.013 \text{ d}^{-1}$, i. e. $P = 78$ and $180 \text{ d}$, respectively. The frequency analysis confirmed that the signals are aliases of each other, since when we subtract one of them, the other one also disappears. The alias peak is visible in the low-frequency regime of the spectral window (Fig. 1, bottom panel). We should also consider that the longer period is close to the time baseline of our data. In order to disentangle the real frequency from its alias, we computed the Bayesian evidence of the two possible solutions, first allowing the period to vary between 50 and 100 d, and then between 100 and 200 d. The Bayesian evidence slightly favoured the solution with $P \sim 78 \text{ d}$, even if not with strong significance ($\Delta \ln Z \simeq 2$). Since, we could not definitely favour one solution over the other, we decided to perform all the subsequent analyses using both sets of parameters.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Posterior distributions (in the top panels, the blue, red, and green lines, respectively) of the retrieved RV signal of TOI-561.03 according to different injected values for the RV semi-amplitudes of candidates .01 and .03. The black line in the top panels corresponds to the observed posterior of the RV semi-amplitude of candidate .03. Median and $1\sigma$ values are marked with vertical dashed and dotted lines, respectively.</figcaption></figure>

Another important outcome of our frequency search is the absence of a signal with a periodicity of ~16 d, that is, the transiting candidate TOI-561.03. Therefore, in order to test our ability to recover the planetary signals, we performed a series of injection/retrieval simulations, thoroughly explained in Appendix B2. The results of this injection/retrieval test are summarized in Fig. 3. We found that the injected RV amplitude of .01 is not significantly affecting the retrieved value for .03, i. e. the cross-talk between the two signals is negligible. We verified that the same conclusion applies to the other signals as well. More importantly, any attempt to retrieve a null signal at the periodicity of the candidate planet .03 would result in an upper limit of $\approx 0.5 \text{ m s}^{-1}$ as we actually observe with the real data set, when exploring the $K$ parameter in logarithmic space. Any signal equal or higher than $1 \text{ m s}^{-1}$ would have been detected ($> 2\sigma$), even if marginally. A signal with amplitude of $0.5 \text{ m s}^{-1}$ would not lead to the detection of the planet (intended as a $3\sigma$ detection), but the retrieved posterior is expected to differ substantially from the observed one, especially on the lower tail of the distribution. We conclude that the planetary candidate TOI-561.03 is undetected in our RV data set, with an upper limit on the semi-amplitude of $0.5 \text{ m s}^{-1}$ ($M_p < 2.0 M_\oplus$).

## 6.2 Transit attribution

Given the non-detection of the planetary candidate TOI-561.03 in the RV data, we investigated more closely the transit-like features associated with this candidate in the TESS light curve, at $T_0^{14} \simeq 1521.9 \text{ d}$ and $T_0 \simeq 1538.2 \text{ d}$, referred from now on as transits 1 and 4, respectively, given their sequence in the TESS light curve (when excluding the transits of the USP candidate). From our preliminary

<footnote n="13">¹³For the model with three additional signals and a quadratic trend, the calculation of the Bayesian evidence did not converge.</footnote>

<footnote n="14">¹⁴All the $T_0$s in this section are expressed in BJD-2457000.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Transits 1 ($T_0 \simeq 1521.9$ d) and 4 ($T_0 \simeq 1538.2$ d) in the TESS detrended light curve associated with the candidate TOI-561.03. The best-fitting transit model from the three-planet model photometric fit is overplotted (black solid line). The black dots are the data points binned over 15 min. With respect to transit 1, the duration of transit 4 looks underestimated by the global model, with a systematic offset in the residuals, especially in the pre-transit phase.</figcaption></figure>

three-planet photometric fit (Fig. A1), we noted that, with respect to the other candidates, TOI-561.03 appears to have a longer transit duration compared to the model, and the residuals show some deviations in the ingress/egress phases. To better understand the cause of these deviations, we checked how the model fits each transit. As Fig. 4 shows, the global model appears to better reproduce the first transit associated with TOI-561.03 (transit 1) than the second transit (transit 4), that has a duration that looks underestimated by the model. Moreover, a two-sample Kolmogorov–Smirnov statistical test<sup>15</sup> (Hodges 1958) on the residuals of transits 1 and 4 suggests that the two residual samples are not drawn from the same distribution (threshold level $\alpha = 0.05$, statistics KS = 0.178, p-value $\ll 0.01$).

Therefore, we hypothesized that the two transit-like features may be unrelated, i. e. they correspond to the transits of two distinct planets. Since two additional planets are actually detected in the RV data, and their periods are longer than the TESS light-curve interval (i.e. that TESS can detect, at most, only one transit for each

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparison between period (left-hand panels) and $T_0$ (right-hand panels) obtained from the RV fit and from the fit of each single transit. Top and bottom panels refer to transits 1 and 4, respectively. Each panel shows the posterior distribution of the analysed parameter, and the shaded area indicates the region within the 68.27th percentile from the mode of the distribution. The vertical solid lines indicate the inferred best-fitting value of the parameter, with thickness proportional to the associated error.</figcaption></figure>

of them), we tested the possibility that the two transits previously associated with TOI-561.03 could indeed be due to the two additional planets inferred from the RV analysis. To check our hypothesis, we first analysed the RV data set with a model encompassing four planets, of which only .01 and .02 have period and time of transit constrained by TESS. In other words, we performed the same RV analysis as described in Appendix B2, but without including TOI-561.03 in the model. We repeated the analysis twice in order to disentangle the periodicity at 78 d from its alias at 180 d, and vice versa. We used the posteriors of the fit to compute the expected time of transit of the outer planets. We then performed two independent fits of transits 1 and 4 with PYORBIT, following the prescriptions as specified in Section 5. We imposed a lower boundary on the period of 22 d, in order to exclude the periods that would imply a second transit of the same planet in the TESS light curve, and an upper limit of 200 d. As a counter-measure against the degeneracy between eccentricity and impact parameter in a single-transit fit, we kept the Van Eylen et al. (2019) eccentricity prior knowing that high eccentricities for such a compact, old system are quite unlikely (Van Eylen et al. 2019). Finally, we compared the posteriors of period and time of transit from the photometric fit with those from RVs, knowing that the former will provide extremely precise transit times, but a broad distribution in period, while RVs give us precise periods, but little information on the transit times. The results are summarized in Fig. 5: the 25.7 ± 0.3 d signal detected in the RVs is located in the vicinity of the main peak of transit 1 period distribution, while the 78.6<sup>+1.8</sup><sub>-2.5</sub> d signal is close to the main peak in transit 4 period distribution. Moreover, Fig. 5 definitely confirms that both the conjunction times inferred from the RV fit corresponding to the ~25 and ~78 d signals, respectively $T_0 = 1520_{-6}^{+3}$ d and $T_0 = 1532_{-9}^{+12}$ d, are consistent with the (much more precise) $T_0$s inferred from the

<footnote n="15"><sup>15</sup>We used the PYTHON version implemented in scipy.stats.ks_2samp.</footnote>
<table><caption>**Table 5.** Final parameters of the TOI-561 system.</caption><thead><tr><th>Parameter</th><th>TOI-561b</th><th>TOI-561c</th><th>TOI-561d</th><th>TOI-561e</th></tr></thead><tbody><tr><td>P (d)</td><td>0.446578 ± 0.000017</td><td>10.779 ± 0.004</td><td>25.62 ± 0.04</td><td>77.23 ± 0.39</td></tr><tr><td>T<sub>0</sub><sup>a</sup> (d)</td><td>1517.498 ± 0.001</td><td>1527.060 ± 0.004</td><td>1521.882 ± 0.004</td><td>1538.181 ± 0.004</td></tr><tr><td>a/R<sub>*</sub></td><td>2.646 ± 0.031</td><td>22.10 ± 0.26</td><td>39.35 ± 0.46</td><td>82.13 ± 0.99</td></tr><tr><td>a (AU)</td><td>0.01055 ± 0.00008</td><td>0.08809 ± 0.0007</td><td>0.1569 ± 0.0012</td><td>0.3274<sup>+0.0028</sup><sub>-0.0027</sub></td></tr><tr><td>R<sub>p</sub>/R<sub>*</sub></td><td>0.0152 ± 0.0007</td><td>0.0308 ± 0.0009</td><td>0.0271 ± 0.0014</td><td>0.0286 ± 0.0011</td></tr><tr><td>R<sub>p</sub> (R<sub>⊕</sub>)</td><td>1.423 ± 0.066</td><td>2.878 ± 0.096</td><td>2.53 ± 0.13</td><td>2.67 ± 0.11</td></tr><tr><td>b</td><td>0.14<sup>+0.13</sup><sub>-0.10</sub></td><td>0.18<sup>+0.16</sup><sub>-0.12</sub></td><td>0.32<sup>+0.17</sup><sub>-0.19</sub></td><td>0.34<sup>+0.13</sup><sub>-0.20</sub></td></tr><tr><td>i (deg)</td><td>87.0<sup>+2.1</sup><sub>-2.8</sub></td><td>89.53<sup>+0.32</sup><sub>-0.39</sub></td><td>89.54<sup>+0.28</sup><sub>-0.21</sub></td><td>89.75<sup>+0.14</sup><sub>-0.08</sub></td></tr><tr><td>T<sub>14</sub> (h)</td><td>1.327<sup>+0.021</sup><sub>-0.030</sub></td><td>3.77<sup>+0.07</sup><sub>-0.15</sub></td><td>4.85<sup>+0.20</sup><sub>-0.35</sub></td><td>6.96<sup>+0.34</sup><sub>-0.38</sub></td></tr><tr><td>e</td><td>0 (fixed)</td><td>0.060<sup>+0.067</sup><sub>-0.042</sub></td><td>0.051<sup>+0.064</sup><sub>-0.036</sub></td><td>0.061<sup>+0.051</sup><sub>-0.042</sub></td></tr><tr><td>ω (deg)</td><td>90 (fixed)</td><td>200<sup>+55</sup><sub>-49</sub></td><td>246<sup>+67</sup><sub>-124</sub></td><td>155 ± 83</td></tr><tr><td>K<sup>b</sup> (m s<sup>-1</sup>)</td><td>1.56 ± 0.35</td><td>1.84 ± 0.33</td><td>3.06 ± 0.33</td><td>2.84 ± 0.41</td></tr><tr><td>M<sub>p</sub><sup>b</sup> (M<sub>⊕</sub>)</td><td>1.59 ± 0.36</td><td>5.40 ± 0.98</td><td>11.95 ± 1.28</td><td>16.0 ± 2.3</td></tr><tr><td>ρ<sub>p</sub> (ρ<sub>⊕</sub>)</td><td>0.55 ± 0.14</td><td>0.23 ± 0.05</td><td>0.74 ± 0.14</td><td>0.84 ± 0.16</td></tr><tr><td>ρ<sub>p</sub> (g cm<sup>-3</sup>)</td><td>3.0 ± 0.8</td><td>1.3 ± 0.3</td><td>4.1 ± 0.8</td><td>4.6 ± 0.9</td></tr><tr><td colspan="5" style="text-align: center;">Common parameter</td></tr><tr><td>ρ<sub>*</sub> (ρ<sub>⊕</sub>)</td><td>1.248 ± 0.043</td><td></td><td></td><td></td></tr><tr><td>u<sub>1</sub></td><td>0.401 ± 0.048</td><td></td><td></td><td></td></tr><tr><td>u<sub>2</sub></td><td>0.208 ± 0.049</td><td></td><td></td><td></td></tr><tr><td>σ<sup>c</sup><sub>jitter,ph</sub></td><td>0.000024<sup>+0.000018</sup><sub>-0.000011</sub></td><td></td><td></td><td></td></tr><tr><td>σ<sup>d</sup><sub>jitter</sub> (m s<sup>-1</sup>)</td><td>1.29 ± 0.23</td><td></td><td></td><td></td></tr><tr><td>γ<sup>e</sup> (m s<sup>-1</sup>)</td><td>79702.58 ± 0.29</td><td></td><td></td><td></td></tr></tbody></table>

Notes. a BJD<sub>TDB</sub>-2457000.

b The here reported values of planet b correspond to the weighted mean between the values inferred from the floating chunk offset method ($K_b = 1.80 \pm 0.38 \text{ m s}^{-1}$, $M_b = 1.83 \pm 0.39 \text{ M}_\oplus$) and from the joint photometric and RV fit ($K_b = 1.39 \pm 0.32 \text{ m s}^{-1}$, $M_b = 1.42 \pm 0.33 \text{ M}_\oplus$).

c Photometric jitter term. d Uncorrelated RV jitter term. e RV offset.

individual fit of transits 1 ($T_0 = 1521.885 \pm 0.004 \text{ d}$) and 4 ($T_0 = 1538.178 \pm 0.006 \text{ d}$), respectively. Regarding the alias at 182 ± 7 d, while the RV period is consistent with the corresponding posterior from the transit fit, the conjunction time $T_0 = 1628 \pm 13 \text{ d}$ that is derived from our analysis is not compatible with any of the transits in the TESS light curve. We also note that the proportion of the orbital period covered by the TESS photometry is ~2.3 times larger for the candidate with 78 d period, thus increasing the chance of getting a transit of it. In conclusion, taking into account both photometric and RV observations, the most plausible solution for the TOI-561 system is a four-planet configuration in which transits 1 and 4 are associated with the planets that have periods of ~25 and ~78 d detected in the RV data, and the 180 d signal is considered an alias of the 78 d signal.

Given this final configuration, hereafter, we will refer to the planets with period ~0.45, ~10.8, ~25, and ~78 d as planets b, c, d, and e, respectively.

## 6.3 The system architecture

Given the presence of two single-transit planets in our data, a joint photometric and RV modelling is necessary in order to characterize the orbital parameters of all members of the TOI-561 system in the best possible way. We considered a four-planet model, with a circular orbit for the USP planet and allowing non-zero-eccentricity orbits for the others. We performed the PyORBIT fit as specified in Section 5, running the chains for 150 000 steps, and discarding the first 50 000 as burn-in. We summarize the results of our best-fitting model in Table 5, and show the transit models, the phase folded RVs, and the global RV model in Figs 6, 7, and 8, respectively. We

obtained a robust detection of the USP planet (planet b) RV semi-amplitude ($K_b = 1.39 \pm 0.32 \text{ m s}^{-1}$), that corresponds to a mass of $M_b = 1.42 \pm 0.33 \text{ M}_\oplus$, while for the 10.8 d period planet (planet c), we obtained $K_c = 1.84 \pm 0.33 \text{ m s}^{-1}$, corresponding to $M_c = 5.40 \pm 0.98 \text{ M}_\oplus$. We point out that the here reported value of $K_b$ and $M_b$ is obtained from the joint photometric and RV fit. However, the final value of $K_b$ and $M_b$ that we decided to adopt (see Section 6.4 for more details) is the weighted mean between the values obtained from the joint fit reported in this section and from the floating chunk offset method described in the next section. In addition, we inferred the presence of two additional planets, with periods of 25.62 ± 0.04 d (planet d) and 77.23 ± 0.39 d (planet e), and robustly determined semi-amplitudes of $K_d = 3.06 \pm 0.33 \text{ m s}^{-1}$ ($M_d = 11.95 \pm 1.28 \text{ M}_\oplus$) and $K_e = 2.84 \pm 0.41 \text{ m s}^{-1}$ ($M_e = 16.0 \pm 2.3 \text{ M}_\oplus$). Both planets show a single transit in the TESS light curve, previously attributed to a transiting planet with period ~16 d, whose presence has, however, been ruled out by our analysis. This allowed us to infer a planetary radius of $R_d = 2.53 \pm 0.13 \text{ R}_\oplus$ and $R_e = 2.67 \pm 0.11 \text{ R}_\oplus$ for planet d and e, respectively.

We performed the stability analysis of our determined solution, computing the orbits for 100 Kyr with the `whfast` integrator (with fixed time-step of 0.1 d) implemented within the `REBOUND` package (Rein & Liu 2012; Rein & Tamayo 2015). During the integration, we checked the dynamical stability of the solution with the Mean Exponential Growth factor of Nearby Orbits (MEGNO or $\langle Y \rangle$) indicator developed by Cincotta & Simó (2000) and implemented within `REBOUND` by Rein & Tamayo (2016). We ran 10 simulations with initial parameters drawn from a Gaussian distribution centred on the best-fitting parameters and standard deviation derived in this
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Top panel: 2-min cadence flattened light curve of TOI-561. The transits of planet b (P ~ 0.45 d), c (P ~ 10.8 d), d (P ~ 25.6 d), e (P ~ 77.2 d) are highlighted with blue, orange, red, and purple triangles, respectively. Bottom panels: TOI-561 phase-folded 2-min light curves over the best-fitting models (solid lines) for the four planets. The light-curve residuals are shown in the bottom panels.</figcaption></figure>

section. All the 10 runs resulted in a MEGNO value of 2, indicating that the family of solutions is stable.

Finally, we checked the presence of any additional signal in the RVs residuals after removing the four-planet model contribution. The GLS periodogram showed a non-significant peak at ~2.5 d, with a normalized power of 0.20, that is, below the 1 per cent FAP threshold (0.26). As a supplemental confirmation, we ran a *PyORBIT* fit of the RVs, assuming first a four-planet model plus an additional signal, and then a four-planet model adding a Gaussian Process (GP) regression. For the latter approach, we employed the quasi-periodic kernel as formulated by Grunblatt, Howard & Haywood (2015), with no priors on the GP hyper-parameters, since we could not identify any activity-related signal in the ancillary data sets (see Section 3.3).<fnref n="16" /> In both cases, the (hyper-)parameters of the additional signal did not reach convergence, while the results for the four transiting planets were consistent with those reported above.

Considering these results, we adopt the parameters and configuration determined in this section as the representative ones for the TOI-561 system, with the only exception of the mass and semi-amplitude of TOI-561 b, that we discuss in the next section.

## 6.4 Alternative characterization of the USP planet

If the separation between the period of the planet and all the other periodic signals is large enough, and the RV signal has a similar or larger semi-amplitude, it is possible to determine the RV semi-amplitude for a USP planet without any assumptions about the number of planets in the system or the activity of the host star.

Under such conditions, during a single night, the influence of any other signal is much smaller than the measurement error and, thus, it can be neglected. If two or more observations are gathered during the same night and they span a large fraction of the orbital phase, the RV semi-amplitude of the USP planet can be precisely measured by just applying nightly offsets to remove all the other signals (e.g. Hatzes et al. 2010; Howard et al. 2013; Pepe et al. 2013; Frustagli et al. 2020 for a recent example). Such an approach, also known as floating chunk offset method (FCO; Hatzes 2014), has proven extremely reliable even in the presence of complex activity signals, as shown by Malavolta et al. (2018). In our case, the shortest, next periodic signal (i. e. TOI-561 c at 10.78 d) is ≈24 times the period of TOI-561 b (i. e. the USP planet at 0.45 d), with similar predicted RV semi-amplitude, making this target suitable for the FCO approach. Thanks to our observational strategy (see Section 2.2), we could use 10 different nights for this analysis. Most notably, during two nights, we managed to gather six observations spanning nearly 5 h, i. e. more than 40 per cent of the orbital period of TOI-561 b, at opposite orbital phases, thus, providing a good coverage in phase of the RV curve. We did not include RV measurements with an associated error greater than 2.5 m s⁻¹ (see Appendix B1). We performed the analysis with *PyORBIT* as specified in Section 5, assuming a circular orbit for the USP planet and including a RV jitter as a free parameter to take into account possible short-term stellar variability and any underestimation of the error bars. From our analysis, we obtained a RV semi-amplitude of $K_p = 1.80 \pm 0.38$ m s⁻¹, corresponding to a mass of $M_p = 1.83 \pm 0.39$ M⊕. The resulting RV jitter is $j < 0.9$ m s⁻¹ (84.13th percentile of the posterior). We show the phase folded RVs of the USP planet in Fig. 9. Since the greater reliability of this method over a full fit of the RV data set is counter-balanced by the smaller number of RVs, we decided not to privilege one over the other. Therefore, we assumed as final semi-amplitude and mass of TOI-561 b the weighted mean of the values obtained from the two

<footnote n="16">¹⁶We are well aware that this is a sub-optimal use of GP regression, and that this approach may be justified in this specific case only as an attempt to identify additional signals.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Phase-folded RV fit with residuals from the joint four-planet photometric and RV analysis. Planets b, c, d, and e are shown in blue, orange, red, and purple, respectively. The reported error bars include the jitter term, added in quadrature.</figcaption></figure>

methods (FCO approach and joint photometric and RV fit), i. e. $K_b = 1.56 \pm 0.35$ m s⁻¹, corresponding to a mass of $M_b = 1.59 \pm 0.36$ M⊕. Table 5 lists the above-mentioned values for TOI-561 b.

# 7 COMPARISON WITH OTHER MODELS

Our final configuration is quite different from the initial one suggested by the *TESS* automatic pipeline. However, the analyses

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Four-planet model from the joint photometric and RV analysis. The grey curve is the the best-fitting model, and the blue points are the HARPS-N data. The residuals are shown in the bottom panel. The reported error bars include the jitter term, added in quadrature.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Phase folded RVs of the 10 nights used to model the RV semi-amplitude of the USP planet using the FCO approach.</figcaption></figure>

performed on the currently available data clearly disfavour the scenario with a ~16 d period candidate. In fact, in addition to the previous analyses, we also performed a joint photometric and RV fit assuming a five-planet model including the 16 d period candidate, and assuming that the two additional signals seen in the RVs were caused by two non-transiting planets, the inner one with period of ~25 d and the outer one both in the case of ~78 and ~180 d period. Such a model, including the TOI-561.01, .02, .03 candidates plus two additional signals, corresponds to the favoured model (Model 2) identified in Section 6.1, and is therefore representative of the best-fitting solution when assuming the *TESS* candidate attribution. In fact, Table 4 suggests that in this case two additional signals need to be added to the three transiting candidates to best reproduce the RV data set, and therefore the five-planet model should be considered also in the joint photometric and RV modelling.

According to the Bayesian evidence (Table 6), computed with the *dynesty* algorithm as specified in Section 5, the four-planet model is strongly favoured with respect to the five-planet model in both cases, with a difference in the logarithmic Bayes factor $2 \Delta \ln \mathcal{Z} \gg 10$ (Kass & Raftery 1995).
<table><caption>**Table 6.** Logarithmic Bayesian evidences for the models considered in Section 7. Model 0 corresponds to the four-planet model, that includes TOI-561.01, .02, and the two additional planets identified in the RVs, showing a single transit each. Models 1 and 2 correspond the five-planet model, i. e. including TOI-561.01, .02, .03, and the two additional RV planets (assumed in this case not to transit), in the case of an outer planet at ~78 and ~180 d period, respectively (see Section 6.1). All the values are expressed with respect to Model 0. We note that the reported errors, as obtained from the nested sampling algorithm, are likely underestimated (Nelson et al. 2020).</caption><thead><tr><th></th><th>Model 0</th><th>Model 1</th><th>Model 2</th></tr></thead><tbody><tr><th>ln Z</th><td>0.0 ± 0.9</td><td>-77.8 ± 1.0</td><td>-76.9 ± 1.0</td></tr></tbody></table>

Moreover, we checked the stability of the five-planet model solutions as described in Section 6.3, with the external planet both on an orbit of 78 and 180 d. For all the planetary parameters, including the mass of the 16 d period planet,<fnref n="17" /> we used the values and standard deviations derived from the joint photometric and RV fit, except for the inclination of the two external planets, that we fixed to 90°.

All of 10 runs yielded unstable solutions, with a close encounter or an ejection occurring within the integration time. In order to assess the origin of the instability of the system, we tested a four-planet configuration following the same procedure as above, removing one planet each time. We found that the orbital configuration of the system could be stable only if we remove the candidate with period of ~16 d. Therefore, the stability analysis additionally confirms our determined four-planet configuration, ruling out the presence of a ~16 d period planet.

## 8 DISCUSSION AND CONCLUSIONS

According to our analysis, TOI-561 hosts four transiting planets, including a USP planet, a ~10.8 d period planet, and two external planets with periods of ~25.6 and ~77.2 d. The latter were initially detected in the RVs data only, but based on our subsequent analyses, we were able to identify a single transit of each planet in the *TESS* light curve; those transits were initially associated with a candidate planet with period of ~16 d, whose presence we ruled out. As a ‘lesson learned’, we would suggest that caution should be taken when candidate planets, detected by photometric pipelines, are based on just two transits. In such cases, one should not hesitate to consider alternative scenarios.

TOI-561 joins the sample of 88 confirmed systems with four or more planets,<fnref n="18" /> and it is one of the few multiplanet systems with both a mass and radius estimate for all the planets. Our global photometric and RV model allowed us to determine the masses and densities of all the planets with high precision, with a significance of ~4.4σ for planet b and >5σ for planets c, d, and e. In Fig. 10, we show the position of TOI-561 b, c, d, and e in the mass–radius diagram of exoplanets with masses and radii measured with a precision better than 30 per cent. The comparison with the theoretical mass–radius curves excludes an Earth-like composition (~33 per cent iron and 67 per cent silicates) for all planets in the system, whose internal structure we further analyse in the following sections.

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Mass–radius diagram for known exoplanets with mass and radius measurements more precise than 30 per cent, colour-coded according to their incidental flux in Earth units. The TOI-561 planets are labelled and represented with coloured diamonds. The USP planets are highlighted with black thick contours. The solid coloured lines represent the theoretical mass-radius curves for various chemical compositions according to Zeng et al. (2019). The shaded grey region marks the maximum value of iron content predicted by collisional stripping (Marcus et al. 2010). The planetary data are taken from the The Extrasolar Planets Encyclopaedia catalogue (http://exoplanet.eu/catalog/) updated to 2020 August 17.</figcaption></figure>

## 8.1 TOI-561 b

The density ($\rho_b = 3.0 \pm 0.8 \text{ g cm}^{-3}$) of the USP planet is consistent with a 50 per cent (or even more) water composition. Such a composition may be compatible with a water-world scenario, where ‘water worlds’ are planets with massive water envelopes, in the form of high pressure H$_2$O ice, comprising > 5 per cent of the total mass. Even assuming the higher mass value inferred with the FCO method ($M_b = 1.83 \pm 0.39 \text{ M}_\oplus$, implying a density of $\rho_b = 3.5 \pm 0.9 \text{ g cm}^{-3}$), TOI-561 b would be located close to the 25 per cent water composition theoretical curve in the mass–radius diagram, and it would be consistent with a rocky composition only at a confidence level greater than 2σ in both radius and mass. Given its proximity to the host star (incident flux $F_p \simeq 5100 F_\oplus$), the presence of any thick H–He envelope has to be excluded due the photo-evaporation processes that such old close-in planets are expected to suffer (e.g. Lopez 2017). Nevertheless, the possibility of a water-world scenario is an intriguing one. An H$_2$O-dominated composition would imply that the planet formed beyond the snow line, accreted a considerable amount of condensed water, and finally migrated inwards (Zeng et al. 2019). While the determination of the precise interior composition of TOI-561 b is beyond the scope of this work, if such an interpretation is proven trustworthy by future observational campaigns, TOI-561 b would support the hypothesis that the formation of super-Earths with a significant amount of water is indeed possible. However, an important caveat should be considered while investigating this scenario. If TOI-561 b was a water world, being more irradiated than the runaway greenhouse irradiation limit, the planet would present a massive and very extended steam atmosphere. Such an atmosphere would substantially increase the measured radius compared to a condensed water world (Turbet et al. 2020). Therefore, a comparison with the condensed water-world theoretical curves should be used

<footnote n="17">¹⁷The mass of the 16 d period planet obtained from the fit was 0.62 ± 1.03 M⊕ and 1.19 ± 1.27 M⊕ for the ~78 and ~180 d external planet period, respectively. Obviously, when selecting the 10 samples, the mass was constrained to positive values.</footnote>

<footnote n="18">¹⁸According to the https://exoplanetarchive.ipac.caltech.edu/</footnote>
with caution, since in this case it could lead to an overestimation of
the bulk water content (Turbet et al. 2020).

Finally, we note that the USP planet is located on the opposite side
of the radius valley, i. e. the gap in the distribution of planetary radii
at ~1.7–2 R<sub>⊕</sub> (Fulton et al. 2017), with respect to all the other planets
in the system. The origin of the so-called radius valley is likely due
to a transition between rocky and non-rocky planets with extended
H–He envelopes, with several physical mechanisms proposed as
explanation, i.e. photo-evaporation (Chen & Rogers 2016; Owen &
Wu 2017; Jin & Mordasini 2018; Lopez & Rice 2018), core-powered
mass-loss (Ginzburg, Schlichting & Sari 2018; Gupta & Schlichting
2019), or superposition of rocky and non-rocky planet populations
(Lee & Chiang 2016; Lopez & Rice 2018). In the TOI-561 system,
planet c is located above the radius valley and it indeed appears to
require a thick H–He envelope (see next section). In the same way,
the compositions of planet d and e are consistent with the presence of
a gaseous envelope. However, the density of TOI-561 b is lower
than expected for a planet located below the radius valley, where we
mainly expect rocky compositions. Moreover, TOI-561 b is the first
USP planet with such a low measured density (see Fig. 10). We note
that also the USP planets WASP-47 e and 55 Cnc e are less dense
than an Earth-like rocky planet, even if both of them have higher
densities than TOI-561 b, i. e. ρ<sub>W47c</sub> = 6.4 ± 0.6 g cm<sup>-3</sup> (Vanderburg
et al. 2017), and ρ<sub>55Cnc<sub>e</sub></sub> = 6.3 ± 0.8 g cm<sup>-3</sup> (Demory et al. 2016),
respectively. Vanderburg et al. (2017) proposed the presence of water
envelopes as a possible explanation for the low densities of these two
planets, even though the inferred amount of water was smaller than
the one required to explain TOI-561 b location in the mass–radius
diagram. It should also be considered that both planets are more
massive than TOI-561 b, i. e. M<sub>W47c</sub> = 6.83 ± 0.66 M<sub>⊕</sub> (Vanderburg
et al. 2017) and M<sub>55Cnc<sub>e</sub></sub> = 8.08 ± 0.31 M<sub>⊕</sub> (Demory et al. 2016),
thus increasing their chances of retaining a small envelope of high-
metallicity volatile materials (or water steam) that could explain their
low densities (Vanderburg et al. 2017). Given its smaller mass, this
scenario is less probable for TOI-561 than for WASP-47 e and 55
Cnc e, making the object even more peculiar. With its particular
properties, this planet could be an intriguing case to test also other
extreme planetary composition models. For example, given the
metal-poor alpha-enriched host star, the planet is likely to have a
lighter core composition.

## 8.2 TOI-561 c, d, and e

TOI-561 c, with a density of ρ_c ~ 1.3 g cm⁻³, is located above the threshold of a 100 per cent water composition, and given its position in the mass–radius diagram, we suppose the presence of a significant gaseous envelope surrounding an Earth-like iron core and a silicate mantle, and possibly a significant water layer (high-pressure ice). If the inner USP planet is water-rich, there is no simple planet formation scenario in which the outer three planets are water-poor. It is simpler to assume that all four planets were formed with similar volatile abundances, and that the inner USP planet lost all of its H–He layer, plus much of its water content, while the outer planets could keep them. Following Lopez & Fortney (2014), assuming a rocky Earth- like core and a solar composition H–He envelope, we estimate that an H–He envelope comprising ~ 4.9 per cent of the planet mass could explain the density of TOI-561 c, using our derived stellar and planetary parameters.

Planets TOI-561 d and e are consistent with a >50 per cent water
composition, a feature that may place them among the water worlds.
However, such densities are also consistent with the presence of a
rocky core plus water mantel surrounded by a gaseous envelope.

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Predicted TTV signal of TOI-561 d and e assuming our best-fitting model (see Table 5). The planets show a strong, anticorrelated signal. The signals of the USP planet (<1 s) and of planet c (<1 min) are not reported.</figcaption></figure>

We estimate that a H–He envelope of ~1.8 and ~2.3 per cent of the
planet mass could explain the observed planetary properties.

## 8.3 Dynamical insights

Our analysis shows that the orbital inclinations of planets c, d and e are all consistent within 1σ (see Table 5), and that the difference with the inclination of the USP planet is of the order of Δi ~ 2.5. According to the analysis of Dai, Masuda & Winn (2018), when the innermost planet has a/R<sub>*</sub> < 5, the minimum mutual inclination with other planets in the system often reaches values up to 5°–10°, with larger period ratios (P<sub>c</sub>/P<sub>b</sub> > 5–6) implying a higher mutual inclination. Considering the large period ratio of TOI-561 (P<sub>c</sub>/P<sub>b</sub> ~ 24) and the value of a<sub>b</sub>/R<sub>*</sub> = 2.6, the measured Δi ~ 2.5 in this case is much lower that the expected inclination dispersion of 6.7 ± 0.7 that Dai et al. (2018) inferred for systems with similar orbital configurations, indicating that the TOI-561 system probably evolved through a mechanism that did not excite the inclination of the innermost planet.

We also performed a dynamical N-body simulation to check if
significant TTVs are expected in the TOI-561 system with our
determined configuration. In fact, the period ratio of TOI-561 d and
e indicates that the planets are close to a 3:1 commensurability, hint
of a second-order mean motion resonance (MMR), that may suggest
the presence of a strong dynamical interaction between these planets.
Starting from the initial configuration (as reported in Table 5), we
numerically integrated the orbits using the N-body integrator *ias15*
within the *REBOUND* package (Rein & Liu 2012). We assumed as
reference time the T₀ of the USP planet (see Table 5), that roughly
corresponds to the beginning of the *TESS* observations of TOI-561.
During the integration, we computed the transit times of each planet
following the procedure described in Borsato et al. (2019), and we
compared the inferred transit times with the linear ephemeris in
order to obtain the TTV signal, reported as an observed-calculated
diagram (*O* − *C*, Agol & Fabrycky 2018) in Fig. 11. According to our
simulation, TOI-561 d and e display an anticorrelated TTV signal,
with a very long TTV period of ~4850 days (~13 yr), and TTV
amplitudes of ~62 min (planet d) and ~84 min (planet e), calculated
computing the GLS periodogram of the simulated TTVs. The anticorrelated signal demonstrates that the two planets are expected to dynamically interact (Agol & Fabrycky 2018). In contrast, the predicted TTV amplitude of planet c is extremely low (~0.9 min), being the planet far from any period commensurability, as well as the USP planet, which has a negligible TTV signal (<1 s). With the solution for the planetary system we propose in this paper, TOI-561 is a good target for a TTV follow-up, that will, however, require a very long time baseline in order to tackle the long-period TTV pattern. To better sample such a long-period TTV signal, it could be worth specifically re-observing the target when the deviations from the linear ephemeris are higher, i. e. during the periods corresponding to the *O* − *C* peaks (or dips) in Fig. 11. According to our simulation, the first peak (dip) corresponds to the period between 2020 March–December, while the second one will be between 2026 January–October, i. e. corresponding to the time-spans between ~400–700 and ~2500–3000 d of integration in Fig. 11, respectively. We remark that this calculation is performed assuming the *T*<sub>0</sub>s inferred from single transit observations, thus implying a significant uncertainty in the TTV phase determination. Therefore, additional photometric observations are necessary to refine the linear ephemeris of the planets, and consequently also the prediction of the TTV phase.

**8.4 Prospects for atmospheric characterization**

Given the interesting composition of the planets in the system, we checked if the TOI-561 planets would be accessible targets for atmospheric characterization through transmission spectroscopy, e.g. with the *James Webb Space Telescope (JWST)*. For all the planets in the system, we calculated the Transmission Spectroscopy Metric (TSM, Kempton et al. 2018), which predicts the expected transmission spectroscopy SNR of a 10-hr observing campaign with *JWST/Near Infrared Imager* and *Slitless Spectrograph (NIRISS)* under the assumptions of cloud-free atmospheres, the same atmospheric composition for all planets of a given type, and a fixed mass–radius relation. We obtained TSM values of 19, 107, 24, and 14 for planets b, c, d, and e, respectively. According to Kempton et al. (2018),<fnref n="19" /> this classifies TOI-561 b and c as high-quality atmospheric characterization targets among the *TESS* planetary candidates. However, it should be noted that the TSM metric assumes rocky composition for planets with radius <1.5 R<sub>⊕</sub> and according to our analysis TOI-561 b is not compatible with such a composition. The same caveat holds for planet c, for which the assumptions under which the TSM is calculated may not be totally valid (e.g. the mass obtained from our analysis is not the same as if calculated with the Chen & Kipping (2017) mass–radius relation, that is the relation assumed in Kempton et al. (2018), and that would imply a mass of M<sub>c</sub> ≃ 8.7 M<sub>⊕</sub>). Therefore, this estimate of the atmospheric characterization feasibility should be used with caution, especially as the TSM metric has been conceived to prioritize targets for follow-up, and not to precisely determine the atmospheric transmission properties.

8.5 Summary and conclusions

According to our analysis, TOI-561 hosts a nearly co-planar four-
planet system, with an unusually low density USP super-Earth
(planet b), a mini-Neptune (planet c) with a significant amount of

volatiles surrounding a rocky core, and two mini-Neptunes, which
are both consistent with a water-world scenario or with a rocky core
surrounded by a gaseous envelope, and that are expected to show a
strong, long-term TTV signal. The multiplanetary nature of TOI-561
offers a unique opportunity for comparative exoplanetology. TOI-561
planets may be compared with the known population of multiplanet
systems to understand their underlying distribution and occurrences,
and to give insights on the formation and evolution processes of
close-in planets, especially considering the intriguing architecture of
the system, with the presence of a uncommonly low-density USP
super-Earth and three mini-Neptunes on the opposite side of the
radius valley.

Considering the few available data (i. e. 2 transits for planet c,
1 transit for planets d, e), additional observations are needed to
unequivocally confirm our solution. Further, high-precision pho-
tometric (i.e. with TESS, that will re-observe TOI-561 in sector
35 – 2021 February/March, or with the CHEOPS satellite) and
RVs observations will help improving the precision on the planets
parameters, both allowing for the detection of eventual TTVs and
increasing the time-span of the RV data set, that could also unveil
possible additional long-period companions.

ACKNOWLEDGEMENTS

We thank the anonymous referee for the constructive comments and
recommendations which helped improving this paper. 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
Programme. Resources supporting this work were provided by the
NASA High-End Computing (HEC) Programme through the NASA
Advanced Supercomputing (NAS) Division at Ames Research Cen-
ter for the production of the SPOC data products. Based on obser-
vations 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 programme, and
A40TAC-23 programme 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 Scot-
tish 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. Parts of
this work have been supported by the National Aeronautics and
Space Administration under grant no. NNX17AB59G issued through
the Exoplanets Research Program. 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 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. This work is
made possible by a grant from the John Templeton Foundation. The

<footnote n="19">¹⁹The authors suggest to select planets with TSM > 12 for R<sub>p</sub> < 1.5 M<sub>⊕</sub>, TSM > 92 for 1.5 R<sub>⊕</sub> < R<sub>p</sub> < 2.75 R<sub>⊕</sub>, and TSM > 84 for 2.75 R<sub>⊕</sub> < R<sub>p</sub> < 4 R<sub>⊕</sub>.</footnote>
opinions expressed in this publication are those of the authors and do not necessarily reflect the views of the John Templeton Foundation. GL acknowledges support by CARIPARO Foundation, according to the agreement CARIPARO-Università degli Studi di Padova (Pratica n. 2018/0098), and scholarship support by the 'Soroptimist International d'Italia' association (Cortina d'Ampezzo Club). GLa, LBo, GPi, VN, GS, and IPa 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'. DNa acknowledges the support from the French Centre National d'Etudes Spatiales (CNES). AM acknowledges support from the senior Kavli Institute Fellowships. ACC acknowledges support from STFC consolidated grant ST/R000824/1 and UK Space Agency grant ST/R003203/1. ASB and MPi acknowledge financial contribution from the ASI-INAF agreement n.2018-16-HH.0. XD is grateful to the Branco-Weiss Fellowship for continuous support. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No. 851555). JNW thanks the Heising-Simons Foundation for support.

DATA AVAILABILITY

HARPS-N observations and data products are available through the Data & Analysis Center for Exoplanets (DACE) at https://dace.u nige.ch/. TESS data products can be accessed through the official NASA website https://heasarc.gsfc.nasa.gov/docs/tess/data-access. html. All underlying data are available either in the appendix/online supporting material or will be available via VizieR at CDS.

REFERENCES

Adibekyan V. Z., Sousa S. G., Santos N. C., Delgado Mena E., González Hernández J. I., Israelian G., Mayor M., Khachatryan G., 2012, A&A, 545, A32

Agol E., Fabrycky D. C., 2018, in Deeg H., Belmonte J., eds, Handbook of Exoplanets, Springer, Cham, p. 797

Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481

Bailer-Jones C. A. L., Rybizki J., Fouesneau M., Mantelet G., Andrae R., 2018, AJ, 156, 58

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

Baruteau C. et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI, Univ. Arizona Press, Tucson, p. 667

Baruteau C., Bai X., Mordasini C., Mollière P., 2016, Space Sci. Rev., 205, 77

Bensby T., Feltzing S., Oey M. S., 2014, A&A, 562, A71

Borsato L. et al., 2019, MNRAS, 484, 3233

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

Buchner J. et al., 2014, A&A, 564, A125

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

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

Chen H., Rogers L. A., 2016, ApJ, 831, 180

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

Cincotta P. M., Simó C., 2000, A&AS, 147, 205

Claret A., 2018, A&A, 618, A20

Cloutier R. et al., 2019, A&A, 621, A49

Collier Cameron A. et al., 2019, MNRAS, 487, 1082

Cosentino R. et al., 2012, Proc. SPIE Conf. Ser. Vol. 8446, Harps-N: the new planet hunter at TNG, SPIE, Bellingham, p. 1

Cosentino R. et al., 2014, HARPS-N @ TNG, two year harvesting data: performances and results, Proc. SPIE Conf. Ser. Vol. 9147, SPIE, Bellingham, p. 8

Cutri R. M. et al., 2003, 2MASS All-Sky Catalog of Point Sources, VizieR Online Data Catalog, II/246

Dai F., Masuda K., Winn J. N., 2018, ApJ, 864, L38

Davies M. B., Adams F. C., Armitage P., Chambers J., Ford E., Morbidelli A., Raymond S. N., Veras D., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI, Univ. Arizona Press, Tucson, p. 787

Demory B.-O., Gillon M., Madhusudhan N., Queloz D., 2016, MNRAS, 455, 2018

Dotter A., 2016, ApJS, 222, 8

Dotter A., Chaboyer B., Jevremović D., Kostov V., Baron E., Ferguson J. W., 2008, ApJS, 178, 89

Douglas S. T., Curtis J. L., Agüeros M. A., Cargile P. A., Brewer J. M., Meibom S., Jansen T., 2019, ApJ, 879, 100

Dragomir D. et al., 2019, ApJ, 875, L7

Dumusque X. et al., 2019, A&A, 627, A43

Eastman J. D., EXOFASTv2: Generalized publication-quality exoplanet modeling code, Astrophysics Source Code Library, ascl:1710.003

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

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

Feroz F., Hobson M. P., Bridges M., 2009, MNRAS, 398, 1601

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

Figueira P., Santos N. C., Pepe F., Lovis C., Nardetto N., 2013, A&A, 557, A93

Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306

Fossati L. et al., 2017, A&A, 601, A104

Fressin F. et al., 2013, ApJ, 766, 81

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

Fulton B. J. et al., 2017, AJ, 154, 109

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

Gelman A., Rubin D. B., 1992, Stat. Sci., 7, 16

Ginzburg S., Schlichting H. E., Sari R., 2018, MNRAS, 476, 759

Gomes da Silva J., Santos N. C., Bonfils X., Delfosse X., Forveille T., Udry S., 2011, A&A, 534, A30

Goodman J., Weare J., 2010, Commun. Appl. Math. Comput. Sci., 5, 65

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

Günther M. N., Pozuelos F. J., Dittmann J. A. et al., 2019, Nature Astron., 3, 1099

Gupta A., Schlichting H. E., 2019, MNRAS, 487, 24

Hatzes A. P. et al., 2010, A&A, 520, A93

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

Helled R. et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI, Univ. Arizona Press, Tucson, p. 643

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

Hippke M., David T. J., Mulders G. D., Heller R., 2019, AJ, 158, 143

Hodges J. L., 1958, Ark. Mat., 3, 469

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

Isaacson H., Fischer D., 2010, ApJ, 725, 875

Jenkins J. M. et al., 2016, in Chiozzi G., Guzman J. C., eds, 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, Tech. rep., NASA Ames Research Center, Mountain View

Jin S., Mordasini C., 2018, ApJ, 853, 163

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

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

Khan S. et al., 2019, A&A, 628, A35

Kipping D. M., 2010, MNRAS, 408, 1758

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

Kochanek C. S. et al., 2017, PASP, 129, 104502

Kreidberg L., 2015, PASP, 127, 1161

Lanza A. F. et al., 2018, A&A, 616, A155
Latham D. W. et al., 2011, ApJ, 732, L24

Lee E. J., Chiang E., 2016, ApJ, 817, 90

Lightkurve Collaboration et al., 2018, Astrophysics Source Code Library, record ascl:1812.013

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

Lissauer J. J. et al., 2012, ApJ, 750, 112

Lopez E. D., 2017, MNRAS, 472, 245

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

Lopez E. D., Rice K., 2018, MNRAS, 479, 5303

Lovis C. et al., 2011, preprint (arXiv:1107.5325)

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

Malavolta L. et al., 2017a, AJ, 153, 224

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

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

Mamajek E. E., Hillenbrand L. A., 2008, ApJ, 687, 1264

Marcus R. A., Sasselov D., Stewart S. T., Hernquist L., 2010, ApJ, 719, L45

McQuillan A., Aigrain S., Mazeh T., 2013, MNRAS, 432, 1203

Milbourne T. W. et al., 2019, ApJ, 874, 107

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

Morbidelli A., Lunine J., O'Brien D., Raymond S., Walsh K., 2012, Annu. Rev. Earth Planet. Sci., 40, 251

Mortier A., Santos N. C., Sousa S. G., Fernandes J. M., Adibekyan V. Z., Delgado Mena E., Montalto M., Israelian G., 2013, A&A, 558, A106

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

Morton T. D., 2015, Astrophysics Source Code Library, record ascl:1812.013

Nardetto N., Mourard D., Kervella P., Mathias P., Mérand A., Bersier D., 2006, A&A, 453, 309

Nardiello D. et al., 2019, MNRAS, 490, 3806

Nardiello D. et al., 2020, MNRAS, 495, 4924

Nelson B. E. et al., 2020, AJ, 159, 73

Noyes R. W., Hartmann L. W., Baliunas S. L., Duncan D. K., Vaughan A. H., 1984, ApJ, 279, 763

Owen J. E., Wu Y., 2017, ApJ, 847, 29

Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011, ApJS, 192, 3

Pecaut M. J., Mamajek E. E., 2013, ApJS, 208, 9

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

Planck Collaboration, 2020, A&A, 641, A1

Raymond S. N., Kokubo E., Morbidelli A., Morishima R., Walsh K. J., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI, Univ. of Arizona Press, Tucson, p. 595

Rein H., Liu S. F., 2012, A&A, 537, A128

Rein H., Tamayo D., 2015, MNRAS, 452, 376

Rein H., Tamayo D., 2016, MNRAS, 459, 2275

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

Roy A. et al., 2020, AJ, 159, 161

Seager S., Mallen-Ornelas G., 2003, ApJ, 585, 1038

Shappee B. J. et al., 2014, ApJ, 788, 48

Skilling J., 2004, in Fischer R., Preuss R., Toussaint U. V., eds, AIP Conf. Proc. 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, Am. Inst. Phys., New York, p. 395

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

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

Sneden C., 1973, ApJ, 184, 839

Sousa S. G., 2014, ARES + MOOG: A Practical Overview of an Equivalent Width (EW) Method to Derive Stellar Parameters, Springer, Cham, p. 297

Sousa S. G., Santos N. C., Israelian G., Lovis C., Mayor M., Silva P. B., Udry S., 2011, A&A, 526, A99

Sousa S. G., Santos N. C., Adibekyan V., Delgado-Mena E., Israelian G., 2015, A&A, 577, A67

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

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

Vanderburg A. et al., 2017, AJ, 154, 237

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

Winn J. N., 2010, in Seager S., eds, Exoplanets, University of Arizona Press, Tucson, p. 55

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

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

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

Ziegler C., Tokovinin A., Briceño C., Mang J., Law N., Mann A. W., 2020, AJ, 159, 19

## SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

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.

## APPENDIX A: PHOTOMETRIC ANALYSIS

We performed a preliminary light-curve fit of the three candidate planets found by the SPOC pipeline and our independent TLS analysis, that is TOI-561.01, .02, and .03 with periods of about 10.8, 0.45, and 16.3 d, respectively. We fit the transits using *PyORBIT*, as specified in Section 5, but assuming circular orbits for all the candidate planets, given the uncertainty associated with the eccentricity from the analysis of *TESS* data alone (Winn 2010). We ran the chains for 100 000 steps, discarding the first 20 000 as burn-in. We list the obtained parameters in Table A1 and we show the best-fitting transit models in Fig. A1. In order to test whether our light-curve flattening affected the inferred parameters of the planetary candidates, we also ran the *PyORBIT* fit on the original PDCSAP light curve. For all the candidates, the difference between the parameters of the two runs was lower than the error on the parameters themselves, indicating that the flattening did not significantly alter the results.

We stress that, at last, our global analysis disclaimed the presence of the planetary candidate TOI-561.03, linking the transits here associated with this candidate to single transits of two additional planets discovered in the system (see Section 6).
<figure><img src="image_14.png" /></figure>

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

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

<figure><img src="image_17.png" /><figcaption>**Figure A1.** Top panel: 2-min cadence flattened light curve of TOI-561. The transits of candidates TOI-561.02 (P ~ 0.45 d), .01 (P ~ 10.8 d), and .03 (P ~ 16.3 d) are highlighted with blue, orange and, green triangles, respectively. Bottom panels: TOI-561 phase-folded light curves over the best-fitting models (solid lines) for the three planets. The grey points are the TESS 2-minute data, the coloured dots are the data points binned over 15 min. The light-curve residuals are shown in the bottom panel. Note the deviations from zero of the residuals in the ingress/egress phase for TOI-561.03.</figcaption></figure>

<table><caption>**Table A1.** Planetary parameters of the three transiting candidates from the initial light-curve fitting.</caption><thead><tr><th>Parameter</th><th>TOI-561.02</th><th>TOI-561.01</th><th>TOI-561.03</th></tr></thead><tbody><tr><td>P (d)</td><td>0.44656 ± 0.00007</td><td>10.780 ± 0.005</td><td>16.309<sup>+0.010</sup><sub>-0.008</sub></td></tr><tr><td>T<sub>0</sub><sup>d</sup> (d)</td><td>1517.4988 ± 0.0019</td><td>1527.060 ± 0.004</td><td>1521.884<sup>+0.003</sup><sub>-0.006</sub></td></tr><tr><td>a/R<sub>*</sub></td><td>2.611 ± 0.030</td><td>21.81 ± 0.25</td><td>28.75 ± 0.33</td></tr><tr><td>a (AU)</td><td>0.01055 ± 0.00008</td><td>0.0881 ± 0.0007</td><td>0.1161 ± 0.0009</td></tr><tr><td>R<sub>p</sub>/R<sub>*</sub></td><td>0.01544 ± 0.0007</td><td>0.0308 ± 0.0009</td><td>0.0285 ± 0.0008</td></tr><tr><td>R<sub>p</sub> (R<sub>⊕</sub>)</td><td>1.46 ± 0.06</td><td>2.91 ± 0.10</td><td>2.70 ± 0.09</td></tr><tr><td>b</td><td>0.16<sup>+0.14</sup><sub>-0.11</sub></td><td>0.17 ± 0.12</td><td>0.07<sup>+0.07</sup><sub>-0.05</sub></td></tr><tr><td>i (deg)</td><td>86.5<sup>+2.7</sup><sub>-3.0</sub></td><td>89.54<sup>+0.30</sup><sub>-0.33</sub></td><td>89.86<sup>+0.10</sup><sub>-0.15</sub></td></tr><tr><td>T<sub>14</sub><sup>b</sup> (hr)</td><td>1.343<sup>+0.022</sup><sub>-0.034</sub></td><td>3.82<sup>+0.06</sup><sub>-0.10</sub></td><td>4.44 ± 0.06</td></tr><tr><td colspan="4">Common parameter</td></tr><tr><td>ρ<sub>*</sub> (ρ<sub>☉</sub>)</td><td>1.200 ± 0.041</td><td></td><td></td></tr><tr><td>u<sub>1</sub></td><td>0.381 ± 0.047</td><td></td><td></td></tr><tr><td>u<sub>2</sub></td><td>0.192 ± 0.050</td><td></td><td></td></tr></tbody></table>

Notes. a BJD<sub>TDB</sub>-2457000.

b Transit duration is derived from the posterior distributions using the formulas in Seager & Mallen-Ornelas (2003).

## APPENDIX B: RV ANALYSIS

### B1 Removal of anomalous points

Before proceeding with a detailed analysis, we verified if any anomalous RV measurement was affecting our analysis. We followed a similar approach to that of Cloutier et al. (2019), but slightly more

sophisticated due to the presence of (possibly up to) five planetary signals. Instead of analysing the power variation of the periodogram's peaks associated with the candidate planets while removing one point at the time, we decided to perform a full RV fit with the methodology described in Section 5, and to compare the resulting RV semi-amplitudes with those derived using the full data set. To reduce
<table><caption>**Table B1.** Best-fitting parameters from the five-planet RV fit, assuming period boundaries of 2–100 d for the outermost planet.</caption><thead><tr><th>Parameter</th><th>TOI-561.02</th><th>TOI-561.01</th><th>TOI-561.03</th><th>TOI-561.04</th><th>TOI-561.05</th></tr></thead><tbody><tr><td>P (d)</td><td>0.44658 ± 0.00001</td><td>10.778 ± 0.004</td><td>16.294 ± 0.008</td><td>25.64<sup>+0.21</sup><sub>-0.18</sub></td><td>77.9 ± 1.9</td></tr><tr><td>T<sub>0</sub><sup>a</sup> (d)</td><td>1517.4983 ± 0.0008</td><td>1527.061 ± 0.003</td><td>1521.883 ± 0.004</td><td>1521<sup>+3</sup><sub>-5</sub></td><td>1535<sup>+9</sup><sub>-10</sub></td></tr><tr><td>e</td><td>0 (fixed)</td><td>0.069<sup>+0.068</sup><sub>-0.048</sub></td><td>0.069<sup>+0.074</sup><sub>-0.048</sub></td><td>0.073<sup>+0.078</sup><sub>-0.051</sub></td><td>0.061<sup>+0.068</sup><sub>-0.043</sub></td></tr><tr><td>ω (deg)</td><td>90 (fixed)</td><td>178 ± 75</td><td>235<sup>+135</sup><sub>-100</sub></td><td>275<sup>+60</sup><sub>-80</sub></td><td>100<sup>+93</sup><sub>-113</sub></td></tr><tr><td>K (m s<sup>-1</sup>)</td><td>1.41 ± 0.33</td><td>1.73 ± 0.36</td><td>&lt;0.37</td><td>3.12 ± 0.36</td><td>2.78 ± 0.44</td></tr><tr><td>M<sub>p</sub> (M<sub>⊕</sub>)</td><td>1.43 ± 0.33</td><td>5.1 ± 1.0</td><td>&lt;1.27</td><td>12.2 ± 1.4</td><td>15.7 ± 2.5</td></tr><tr><td colspan="6">Common parameter</td></tr><tr><td>σ<sub>jitter</sub><sup>b</sup> (m s<sup>-1</sup>)</td><td>1.32 ± 0.23</td><td></td><td></td><td></td><td></td></tr><tr><td>γ<sup>c</sup> (m s<sup>-1</sup>)</td><td>79702.58 ± 0.30</td><td></td><td></td><td></td><td></td></tr></tbody></table>

Notes. <sup>a</sup> BJD<sub>TDB</sub>-2457000.

<sup>b</sup> Uncorrelated jitter term.

<sup>c</sup> RV offset.

<table><caption>**Table B2.** Best-fitting parameters from the five-planet RV fit, assuming period boundaries of 100–200 d for the outermost planet.</caption><thead><tr><th>Parameter</th><th>TOI-561.02</th><th>TOI-561.01</th><th>TOI-561.03</th><th>TOI-561.04</th><th>TOI-561.05</th></tr></thead><tbody><tr><td>P (d)</td><td>0.44658 ± 0.00001</td><td>10.779 ± 0.004</td><td>16.294 ± 0.007</td><td>25.82 ± 0.19</td><td>179.5<sup>+8.3</sup><sub>-7.4</sub></td></tr><tr><td>T<sub>0</sub><sup>a</sup> (d)</td><td>1517.4983 ± 0.0009</td><td>1527.061 ± 0.003</td><td>1521.883 ± 0.004</td><td>1518 ± 3</td><td>1632<sup>+13</sup><sub>-15</sub></td></tr><tr><td>e</td><td>0 (fixed)</td><td>0.067<sup>+0.072</sup><sub>-0.047</sub></td><td>0.064<sup>+0.070</sup><sub>-0.045</sub></td><td>0.072<sup>+0.071</sup><sub>-0.051</sub></td><td>0.058<sup>+0.064</sup><sub>-0.041</sub></td></tr><tr><td>ω (deg)</td><td>90 (fixed)</td><td>148<sup>+118</sup><sub>-107</sub></td><td>189<sup>+118</sup><sub>-127</sub></td><td>287<sup>+67</sup><sub>-73</sub></td><td>128<sup>+98</sup><sub>-113</sub></td></tr><tr><td>K (m s<sup>-1</sup>)</td><td>1.57 ± 0.32</td><td>0.69<sup>+0.41</sup><sub>-0.46</sub></td><td>&lt;0.54</td><td>3.10 ± 0.36</td><td>3.17 ± 0.49</td></tr><tr><td>M<sub>p</sub> (M<sub>⊕</sub>)</td><td>1.59 ± 0.33</td><td>2.01<sup>+1.20</sup><sub>-1.35</sub></td><td>&lt;1.91</td><td>12.1 ± 1.4</td><td>23.7 ± 3.7</td></tr><tr><td colspan="6">Common parameter</td></tr><tr><td>σ<sub>jitter</sub><sup>b</sup> (m s<sup>-1</sup>)</td><td>1.34 ± 0.23</td><td></td><td></td><td></td><td></td></tr><tr><td>γ<sup>c</sup> (m s<sup>-1</sup>)</td><td>79703.86 ± 0.25</td><td></td><td></td><td></td><td></td></tr></tbody></table>

Notes. <sup>a</sup> BJD<sub>TDB</sub>-2457000.

<sup>b</sup> Uncorrelated jitter term.

<sup>c</sup> RV offset.

computational time, we decided to remove from the data set five consecutive observations at once (i. e. performing 17 iterations rather than 82), and then performed the leave-one-out cross-validation on those subsets showing deviating RV semi-amplitudes in order to identify the anomalous RV measurement. With this approach, we found out that a total of five RV measurements, with associated errors greater than 2.5 m s⁻¹ and SNR < 35 were systematically producing a decrease in the semi-amplitude of candidates .01 and .02 by ≈0.1–0.2 m s⁻¹, and we therefore removed these points from our data set in order to improve the accuracy of our results, even if the total variation in RV semi-amplitude was within the error bars. We note that these observations are clearly outliers at more than 2σ in both the SNR of the spectra and the RV error distributions (see Section 2.2), which is simply the consequence of having been gathered in sub-optimal weather conditions. A much simpler sigma-clipping selection would have led to the exclusion of the same data points. The complex approach we employed in this work can thus be avoided in future analysis involving HARPS-N data.

## B2 RV modelling and injection/retrieval tests

Given the results of the frequency analysis in Section 6.1, we performed a PyDE + emcee RV fit with PyORBIT, following the methodology as described in Section 5, and assuming the model suggested by the Bayesian evidence computed in Section 6.1 (see Table 4), i. e. a model with the three transiting candidates plus two additional ones. We performed two independent fits, constraining the

period of the outer signal to be shorter or longer than 100 d, in order to disentangle the 78 periodicity from its alias at 180, respectively. We ran the chains for 150 000 steps, discarding the first 50 000 as burn-in. The results of this analysis are reported in Tables B1 and B2.

In all our RV fits, regardless of the assumed period of the outermost planet, TOI-561.03 (i. e. the candidate with period of ~16.3 d) remains undetected with an upper limit of $K \lesssim 0.5$ m s⁻¹, corresponding to a rather non-physical mass of $\lesssim 2$ M<sub>⊕</sub> (at $1\sigma$) for a planet with $R_p \simeq 2.7$ R<sub>⊕</sub>. We, thus, performed a series of injection/retrieval simulations in order to assess the influence of the observational sampling and of the precision in the mass measurements of the other planets. In a first run, the synthetic data sets were simulated by assuming the orbital parameters as previously determined in the RV fits for the candidate planets .01, .02, and the non-transiting candidates, while the RV semi-amplitude of the candidate planet at 16 d was varied between 0.0 and 1.5 m s⁻¹ in steps of 0.5 m s⁻¹. For computational reasons, we performed this analysis only with the 78-d solution for the outer planet. We projected the model on to the real epochs of observation and then, we added a Gaussian noise corresponding to the measured error plus an RV jitter of 1.0 m s⁻¹ added in quadrature, while preserving the original value in the analysis. We built 50 different noise realizations and analysed each of them with the same methodology as before, i. e. PyDE + emcee through PyORBIT, but for a shorter chain length<fnref n="20" /> to reduce computing time. The posteriors of each parameter were then obtained by putting together

<footnote n="20">²⁰10 000 steps after convergence, reached at approximately 15 000 steps.</footnote>
the individual posterior distributions from each noise realization. We
finally repeated the same analysis but varying the RV semi-amplitude
of the candidate planet .01, i. e. the closest signal in frequency space
and the one with the most uncertain RV semi-amplitude measurement
other than the USP candidate, by ±0.5 m s⁻<fnref n="1" /> with respect to the value
of 1.7 m s⁻<fnref n="1" /> used in the previous analysis.

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

<footnote n="1">¹Department of Physics and Astronomy, 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">³Aix-Marseille Université, CNRS, CNES, LAM, F-13013 Marseille, France</footnote>

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

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

<footnote n="6">⁶Observatoire Astronomique de l'Université de Genève, 51 Chemin des
Maillettes, CH-1290 Versoix, Switzerland</footnote>

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

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

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

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

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

<footnote n="12">¹²INAF – Osservatorio Astrofisico di Torino, Via Osservatorio 20, I-10025
Pino Torinese, Italy</footnote>

<footnote n="13">¹³NASA Ames Research Center, Moffett Field, CA 94035, USA</footnote>

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

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

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

<footnote n="17">¹⁷INAF – Osservatorio Astronomico di Roma, Via Frascati 33, I-00078 Monte
Porzio Catone, Italy</footnote>

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

<footnote n="19">¹⁹INAF – Osservatorio Astronomico di Cagliari, via della Scienza 5, I-09047
Selargius, Italy</footnote>

<footnote n="20">²⁰INAF – Osservatorio Astrofisico di Catania, Via S. Sofia 78, I-95123
Catania, Italy</footnote>

<footnote n="21">²¹Department of Earth, Atmospheric and Planetary Sciences, and Kavli
Institute for Astrophysics and Space Research, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA</footnote>

<footnote n="22">²²Astrophysics Research Centre, School of Mathematics and Physics, Queen's
University Belfast, Belfast BT7 1NN, UK</footnote>