key: cord-0174507-siyf0sil
authors: Lacedelli, G.; Malavolta, L.; Borsato, L.; Piotto, G.; Nardiello, D.; Mortier, A.; Stalport, M.; Cameron, A. Collier; Poretti, E.; Buchhave, L. A.; L'opez-Morales, M.; Nascimbeni, V.; Wilson, T. G.; Udry, S.; Latham, D. W.; Bonomo, A. S.; Damasso, M.; Dumusque, X.; Jenkins, J. M.; Lovis, C.; Rice, K.; Sasselov, D.; Winn, J. N.; Andreuzzi, G.; Cosentino, R.; Charbonneau, D.; Fabrizio, L. Di; Fiorenzano, A. F. Martinez; Ghedina, A.; Harutyunyan, A.; Lienhard, F.; Micela, G.; Molinari, E.; Pagano, I.; Pepe, F.; Phillips, D. F.; Pinamonti, M.; Ricker, G.; Scandariato, G.; Sozzetti, A.; Watson, C. A.
title: An unusually low density ultra-short period super-Earth and three mini-Neptunes around the old star TOI-561
date: 2020-09-04
journal: nan
DOI: nan
sha: eade7256d3f51ca2a4efa79c436f3296aef40696
doc_id: 174507
cord_uid: siyf0sil

Based on HARPS-N radial velocities (RVs) and TESS photometry, we present a full characterisation of the planetary system orbiting the early K 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_{rm b} = 0.45$ d, mass $M_{rm b} =1.59 pm 0.36$ M$_oplus$ and radius $R_{rm b}=1.42 pm 0.07$ R$_oplus$, and three mini-Neptunes: TOI-561 c, with $P_{rm c} = 10.78$ d, $M_{rm c} = 5.40 pm 0.98$ M$_oplus$, $R_{rm c}= 2.88 pm 0.09$ R$_oplus$; TOI-561 d, with $P_{rm d} = 25.6$ d, $M_{rm d} = 11.9 pm 1.3$ M$_oplus$, $R_{rm d} = 2.53 pm 0.13$ R$_oplus$; and TOI-561 e, with $P_{rm e} = 77.2$ d, $M_{rm e} = 16.0 pm 2.3$ M$_oplus$, $R_{rm e} = 2.67 pm 0.11$ R$_oplus$. Having a density of $3.0 pm 0.8$ 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.

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 multi-planet 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 E-mail: gaia.lacedelli@phd.unipd.it 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 characterisation 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 proper-ties of multiple systems (e.g. Latham et al. 2011; Millholland et al. 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 characterised 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; Raymond et al. 2014; Helled et al. 2014 ; Baruteau et al. 2014 Baruteau et al. , 2016 Davies et al. 2014) . In this paper, we combine TESS photometry (Section 2.1) and high precision RVs gathered with the HARPS-N spectrograph (Section 2.2) to characterise the multi-planet system orbiting the star TOI-561. The TESS pipeline identified three candidate planetary signals, namely an ultra-short period (USP) candidate (P ∼ 0.45 days), and two additional candidates with periods of ∼ 10.8 and ∼ 16.4 days. 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. We first analysed separately the photometric (Section 4) and RVs data (Section 5), also performing some injection/retrieval simulations to test our ability to recover the planetary signals given our dataset sampling and precision. We then performed a series of analyses to determine the actual system configuration (in Sections 6 and 7), and we present our best-fitting solution from the joint photometric and RVs modelling in Section 8. 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 9).

2.1 TESS photometry TOI-561 was observed by TESS in two-minute cadence mode during observations of sector 8, between 2 February and 27 February 2019. The astrometric and photometric parameters of the star are listed in Table 1 . Considering the download time, and the loss of 3.26 days of data due to an interruption in communications between the instrument and the spacecraft that occurred during sector 8 1 , a total of 20.22 days 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 days (TOI-561.01), 0.4 days (TOI-561.02), and 16.4 days (TOI-561.03), respectively. The pipeline identified 55 transits of TOI-561.02, two transits of TOI-561.01, and two transits of TOI-561.03, with depths of 290, 1207, and 923 1 See TESS Data Release Notes: Sector 8, DR10 (https:// archive.stsci.edu/tess/tess_drn.html). E) Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010 ). a Gaia DR2 parallax is corrected by +50 ± 7 µas (with the error added in quadrature) as suggested by Khan et al. (2019) .

ppm and signal-to-noise-ratios (S/N) of 10.0, 9.8 and 9.2, respectively. For our photometric analysis in Section 4, we used the light curve based on the Pre-search Data Conditioning Simple Aperture Photometry 2 (PDCSAP, Smith et al. 2012; Stumpe et al. 2012 Stumpe et al. , 2014 . We also extracted the 30-minutes 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).

We collected 82 3 spectra using HARPS-N at the Telescopio Nazionale Galileo (TNG), in La Palma (Cosentino et al. 2012 (Cosentino et al. , 2014 , with the goal of precisely determining the masses of the three candidate planets and to search for additional planets. The observations started on November 17, 2019 and ended on June 13, 2020, 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 characterise the signal of the USP candidate, we collected 6 points per night on February 4 and February 6, 2020, 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 (February-March 2020). The exposure time was set to 1800 seconds, which resulted in a S/N at 550 nm of 77 ± 20 (median ± standard deviation) and a measurement uncertainty of 1.2 ± 0.6 m s −1 . 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 crosscorrelation function (CCF, Baranne et al. 1996; 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 May 31, 2020. 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 −1 , 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 −1 .

The RV data with their 1σ uncertainties and the associated activity indices (see Section 3.3 for more details) are listed in Table 2 .

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 microturbulent velocity, ξ t . The gravity log g is obtained by imposing the same average abundance from neutral and ionised iron lines. We obtained the EW measurements using ARESv2 4 (Sousa et al. 2015) . We used the local thermodynamic equilibrium (LTE) code MOOG 5 (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 (S/N> 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 ± 69 K, log g = 4.60 ± 0.12, [Fe/H] = −0.40 ± 0.05 and ξ t = 0.78 ± 0.08 km s −1 .

The spectral synthesis match was performed using the Stellar Parameters Classification tool (SPC, Buchhave et al. 2012 Buchhave et al. , 2014 . It determines effective temperature, surface gravity, metallicity and line broadening by performing a crosscorrelation 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 ran the SPC on the 62 GTO spectra only 6 : the S/N 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 ± 50 K, log g = 4.49 ± 0.10, [M/H] = −0.36 ± 0.08 and v sin i < 2 km s −1 . We finally used CCFpams 7 , 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±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. 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) 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.

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 et al. 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, Dotter 2016; Choi et al. 2016; Paxton et al. 2011 ) and the Dartmouth Stellar Evolution Database (Dotter et al. 2008) . For all methods, we assumed σ T eff = 70 K, σ log g = 0.12, σ [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 . From the mean and standard deviation of all the posterior samplings we obtained M = 0.785 ± 0.018 M and R = 0.849 ± 0.007 R . We derived the stellar density ρ = 1.285 ± 0.040 ρ (ρ = 1.809 ± 0.056 g cm −3 ) directly from the posterior distributions of M and R . We summarise 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.86-th percentile of the distribution of the combined posteriors, as for the other parameters. We note however that an isochrone fit performed through EXOFASTv2 (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 low value of the log R HK index (−5.003 ± 0.012), derived using the calibration by Lovis et al. (2011) and assuming B − V = 0.71, indicates 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 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 HK index, FWHM, and bisector span (BIS) computed by the HARPS-N DRS, we included in our analysis the V asy (Figueira et al. 2013) and ∆V (Nardetto et al. 2006 ) asymmetry indicators, as implemented by Lanza et al. (2018) , and the chromospheric activity indicator Hα (Gomes da .

The Generalized Lomb-Scargle (GLS, Zechmeister & Kürster 2009) periodograms of the above-mentioned indexes are shown in Figure 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%, computed with a bootstrap approach. The periodogram of the RVs 8 reveals the presence of significant peaks at 25 days, 180 days, 10 days (corresponding to one of the transiting candidates), and 78 days, 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% 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 days with an amplitude of 0.13 ppt and a power of 0.014, that is, above the 0.1% 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. Moreover, we performed an auto correlation analysis, following the prescription by McQuillan et al. (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 dataset 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 days), and 20 ppt in the long term period (magnetic cycles, > 100 days), 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 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 RVs due to the magnetic activity of our star is likely below 1-2 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 (σ jitter ).

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. First, we notice that the star has a good astrometric Gaia DR2 solution (Gaia , with zero excess noise and a re-normalised 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 5.2). Moreover, the overall RV variation below 25 m s −1 and the shape of the CCFs of our HARPS-N spectra exclude the eclipsing binary 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 ∆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 eclipsing binaries (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 eclipsing binary within the TESS aperture with a 0.5% eclipse, but no more than 7 magnitudes 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. 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 days, and the additional signals identified in the RVs analysis (Sections 5 and 8) at 25, 78 and 180 days. The dashed and dotted horizontal lines show the 1% and 0.1% FAP levels, respectively. The TESS periodogram shows a series of peaks below 10 days, 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 days, respectively, so there is no correspondence with the 180 days signal. Moreover, both of them are below the 1% FAP.

An independent confirmation was provided by the analysis of the in-/out-of-transit difference centroids on the TESS FFIs (Figure 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 × 10 pixels 2 (∼ 200 × 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-minute 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 5. 

We downloaded the two-minute cadence PDCSAP light curve, and removed all the observations encoded as NaN or flagged as bad-quality (DQUALITY>0) points by the SPOC pipeline 9 . 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 , with a window of 1.5 days, and masking the in-transit points to avoid modifications of the transit shape. We verified with a periodogram analysis that our flattening procedure correctly removed the signal at 3.5 days identified in Section 3.3. The final light curve is shown in Figure 3 . 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 ). The first three significant identified signals nicely matched the TESS suggested periods (P TLS = 10.78 d, 0.44 d, 16.28 d).

In this stage, we fit the transits assuming circular orbits for all the candidate planets, given the uncertainty associated with the eccentricity from the analysis of TESS data alone (Winn 2010) . 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. We fitted a common value for the stellar density ρ , imposing a Gaussian prior based on the value from Table 3 . We included a quadratic limb-darkening law with 9 https://archive.stsci.edu/missions/tess/doc/ EXP-TESS-ARC-ICD-TM-0014.pdf Gaussian priors on the coefficients u 1 , u 2 , obtained through a bilinear interpolation of limb darkening profiles by Claret (2018) 10 . We initially calculated the standard errors on u 1 , u 2 using a Monte Carlo approach that takes into account the errors on T eff and log g as reported in Table 3 , obtaining u 1 = 0.393 ± 0.007 and u 2 = 0.204 ± 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 shortterm 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. We performed the light curve fitting using PyOR-BIT 11 (Malavolta et al. 2016 (Malavolta et al. , 2018 , a convenient wrapper for the analysis of transit light curves and radial velocities. All the transit models were computed with the batman package (Kreidberg 2015) , with an exposure time of 120 seconds and an oversampling factor of 10 (Kipping 2010). Global optimisation of the parameters was performed using the differential evolution code PyDE 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 decided to fit 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 d, 0.45 d, and 16.3 d, respectively. We ran the chains with 40 walkers for 100 000 steps, checking the convergence with the Gelman-Rubin statistic (Gelman & Rubin 1992) , with a threshold value ofR = 1.01. We also performed an auto-correlation analysis of the chains: if the chains were longer than 100 times the estimated auto-correlation time and this estimate changed by less that 1%, we considered the chains as converged. As a conservative choice, we discarded the first 20 000 steps as burn-in, i. e., a number larger than the convergence point as just defined, and then we applied a thinning factor of 100, obtaining a total of 32 000 independent samples for each parameter. We list the obtained parameters in Table 4 and we show the best-fitting transit models in Figure 3 . In order to test whether our light curve flattening affected the inferred parameters of the planetary candidates, we also ran the Py-ORBIT 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 note 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 (Figure 3 ). To better understand the cause of these deviations, we checked how the model fits each transit. As Figure 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 current model. Moreover, a two-sample Kolmogorov-Smirnov statistical test 13 (Hodges 1958 ) on the residuals of transit 1 and 4 suggests that the two residual samples are not drawn from the same distribution (threshold level α = 0.05, statistics KS = 0.178, p−value 0.01). Given this discrepancy, we performed some additional checks on the light curve to exclude the possibility that the transit-like features were caused by instrumental artefacts. 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 14 (Lightkurve ) 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 systematics in every light curve pixel, and we found none. Finally, we checked for correlations be- 13 We used the Python version implemented in scipy.stats.ks_2samp. 14 https://github.com/KeplerGO/lightkurve tween 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. We further investigate and propose a solution for the discrepancy of TOI-561.03 associated transits in the following sections.

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 an 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 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.01 at 10.78 days) is 24 times the period of TOI-561.02 (i. e., the USP planet at 0.45 days), with similar predicted RV semi-amplitude, making this target suitable for the nightly offset approach. Thanks to our observational strategy (see Section 2.2) we could use ten different nights for this analysis. Most notably, during two nights we managed to gather six observations spanning nearly 5 hours, i. e., more than 40% of the orbital period of TOI-561.02, 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 −1 (see Section 5.3). We performed the analysis with PyORBIT, 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 errorbars. We assumed Gaussian priors on the orbital period and the central time of transit for TOI-561.02 following the results obtained in Section 4.2. The MCMC run, convergence checks, and analysis of results were performed as described in that same section. From our analysis, we obtained a RV semi-amplitude of K p = 1.80 ± 0.38 m s −1 , corresponding to a mass of M p = 1.83 ± 0.39 M ⊕ . The resulting RV jitter is j < 0.9 m s −1 (84.13-th percentile of the posterior). We show the phase folded RVs of the USP planet in Figure 5 .

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 d −1 (P 1 25.6 d), f 2 = 0.006 d −1 or 0.013 d −1 (P 2 170 d or 78 d) with the two frequencies being related to each other (i. e., removing one of them implies the vanishing of the other one), f 3 = 0.093 d −1 (P 3 10.8 d, corresponding to the TOI-561.01 candidate), and f 4 = 2.239 d −1 (P 4 0.45 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 characterise 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 assuming the transiting candidates only, and then including an additional planet in each of the successive models. We computed the Bayesian evidence for each model using the MultiNest nested-sampling algorithm (Feroz & Hobson 2008; Feroz et al. 2009 Feroz et al. , 2019 with the Python wrapper py-MultiNest (Buchner, J. et al. 2014) . We assumed 1000 live points and a sampling efficiency of 0.3, including the jitter term in the model. For the transiting candidates, we used the results from Section 4 to impose priors on period and time of transit, while for the additional signals we allowed the periods to span between 2 and 200 days (i. e., the time span of our dataset). 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) , while we allowed the semi-amplitude K to vary between 0.01 and 100 m s −1 for all the candidate planets. We report the obtained Bayesian evidences of the four models in Table 5 . According to this analysis, we concluded that the model with two additional signals (with no trend) is strongly favoured over the others, with a difference in the logarithmic Bayes factor 2 ∆ 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, 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 Model 0 corresponds to the model with no additional RVs signal other than the signal from the transiting candidates, model 1, 2 and 3 to the models with 1, 2 and 3 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 each of the four models. In all cases, the Bayesian evidence systematically disfavoured the presence of any trend 15 . The first additional signal was associated with a candidate with ∼ 25 d period, 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 78 and 180 days respectively, without any preference of one over the other. A periodogram analysis confirmed that the signals are aliases of each other. In fact, when subtracting one of the two signals from the periodogram, the other one also disappears. In order to disentangle the real period from its alias, we computed the Bayesian evidence of the two possible solutions, first allowing the period to vary between 50 and 100 days, and then between 100 and 200 days. The Bayesian evidence slightly favoured the solution with P ∼ 78 d, even if not with strong significance (∆ ln Z 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.

Before proceeding with a more detailed analysis with the selected model, 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 same methodology as described in the next section, and to compare the resulting RV semi-amplitudes with those derived using the full dataset. To reduce computational time, we decided to remove from the dataset 5 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 semiamplitudes in order to identify the anomalous RV measurement. With this approach, we found out that a total of 5 RV measurements, with associated errors greater than 2.5 m s −1 and S/N< 35 were systematically producing a decrease in the semi-amplitude of candidates .01 and .02 by ≈ 0.1−0.2 m s −1 , and we therefore removed these points from our dataset 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 S/N 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. The analysis presented in 5.1 has already been performed without including these points.

We performed a PyDE+emcee fit with PyORBIT, following the same methodology as described in Section 4.2, assuming the model suggested by the Bayesian evidence, i. e. a model with the three transiting candidates plus two additional ones. We used the same priors as specified in the previous section. We performed two independent fits, constraining the period of the outer signal to be shorter or longer than 100 days, in order to disentangle the 78 periodicity from its alias at 180 respectively. The results of this analysis are reported in Tables A1 and A2 in Appendix A. 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 0.5 m s −1 corresponding to a rather nonphysical mass of 2 M ⊕ (at 1σ) for a planet with R p 2.7 R ⊕ . 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 datasets were simulated by assuming the orbital parameters as previously determined for planets .01, .02, and the non-transiting planets, while the RV semi-amplitude of the candidate planet at 16 d was varied between 0.0 m s −1 and 1.5 m s −1 in steps of 0.5 m s −1 . For computational reasons, we performed this analysis only with the 78− day solution for the outer planet. We projected the model onto 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 −1 added in quadrature, while preserving the original value in the analysis. We built 50 different noise realisations and analysed each of them with the same methodology as before, i. e., PyDE+emcee through PyORBIT, but for a shorter chain length 16 to reduce computing time. The posteriors of each parameter were then obtained by putting together the individual posterior distributions from each noise realisation. We finally repeated the same analysis but varying the RV semi-amplitude of 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 −1 with respect to the value of 1.7 m s −1 used in the previous analysis. The results of this injection/retrieval test are summarised in Figure 6 . As a first result, we can see 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 16 10 000 steps after convergence, reached at approximately 15 000 steps. of the candidate planet .03 would result in an upper limit of ≈ 0.5 m s −1 as we actually observe with the real dataset, when exploring the K parameter in logarithmic space. Any signal equal or higher than 1 m s −1 would have been definitely detected (> 2σ), even if marginally. A signal with amplitude of 0.5 m s −1 would not lead to the detection of the planet (intended as a 3-σ detection), but the retrieved posterior is expected to differ substantially from the observed one, especially on the lower tail of the distribution. We can thus affirm that the planet is undetected, with an upper limit on the semi-amplitude of 0.5 m s −1 .

We performed the stability analysis of the solutions of the favoured models from Section 4 and 5, that is three candidate transiting planets with two additional RV signals. We computed 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 Y ) indicator developed by Cincotta & Simó (2000) and implemented within rebound by Rein & Tamayo (2016) . For each configuration, with the external planet on an orbit of 78 d and 180 d, we ran 10 simulations with initial parameters drawn from a Gaussian distribution centred on the best-fitting parameters and standard deviation from previous sections 17 . All of them 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.

Considering these results, we questioned the existence of the candidate at ∼ 16 d period, TOI-561.03. Our photometric analysis (Section 4, in particular Figure 4 ) already highlighted the discrepancy between the two transit-like features attributed to the ∼ 16 d period candidate. Moreover, despite the high precision of our RV dataset, the signal of TOI-561.03 remained undetected. Our injection/retrieval simulations (Section 5.4) demonstrated that only a candidate with a semi-amplitude K < 0.5 m s −1 (M p < 2.0 M ⊕ ) would be undetected using our current dataset. Even assuming the most optimistic case of 2.0 M ⊕ , such a candidate would have an extremely low density (ρ < 0.10 ρ ⊕ ) when considering our inferred radius from the light curve fit (R p = 2.70 ± 0.09 R ⊕ ). The instability of the system caused by the 16 d candidate planet finally convinced us that all these observational facts were not just coincidences. Therefore, we decided to reject the hypothesis of a 16 d period candidate (TOI-561.03), and we investigated a different scenario.

After disclaiming the presence of the transiting candidate TOI-561.03, we hypothesised that the two transit-like features at T 0 18 1521.9 d and T 0 1538.2 d (transit 1 and 4, indicated with green triangles in Figure 3 , respectively) 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 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 dataset 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 repeated the analysis of Section 5 without including TOI-561.03 in the model. Once again, 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 transit 1 and 4 with PyORBIT, following again the same prescriptions as Section 4. We imposed a lower boundary on the period of 22 days, 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 days. 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 radial velocities, 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 summarised in Figure 7 : 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 +1.8 −2.5 d signal is close to the main peak in transit 4 period distribution. Moreover, Figure 7 definitely confirms that both the conjunction times (T c s) inferred from the RV fit corresponding to the ∼ 25 and ∼ 78 days signals, respectively T c = 1520 +3 −6 d and T c = 1532 +12 −9 d, are consistent with the (much more precise) T 0 s inferred from the individual fit of transit 1 (T 0 = 1521.885 ± 0.004 d) and 4 (T 0 = 1538.178 ± 0.006 d) respectively. Regarding the alias at 182 ± 7 days, while the RV period is consistent with the corresponding posterior from the transit fit, the conjunction time T c of 1628 ± 13 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, and the stability analysis, the most plausible solution for the TOI-561 system is a fourplanet configuration in which transits 1 and 4 are associated with the planets that have periods of ∼ 25 d 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 days as planets b, c, d and e, respectively.

Given the presence of two single-transit planets in our data, a joint photometric and RV modelling is necessary in order to characterise the orbital parameters of all members of the TOI-561 system in the best possible way. We assumed a fourplanet model, with a circular orbit for the USP planet and Keplerian orbits for the others, and uniform priors for all the parameters, except for the limb darkening coefficients and the stellar density, for which we assumed the same priors as in Section 4, and for the eccentricities of the three external planets, for which we assumed the priors specified in Van Eylen et al. (2019) for compact multi-planet systems. Once again we used PyORBIT, with a global optimisation PyDE run followed by an emcee analysis. We ran 64 walkers for 150 000 steps, discarding the first 50 000 as burn-in and assuming a thinning factor of 100. The choice of the burn-in value was driven by the results of the auto-correlation analysis, as described in 4.2. We summarise the results of our best-fitting model in Table 6 , and show the transit models, the phase folded RVs, and the global RV model in Figures 8, 9 , and 10 respectively. We obtained a robust detection of the USP planet (planet b) RV semi-amplitude (K b = 1.56±0.35 m s −1 ), that corresponds to a mass of M b = 1.59 ± 0.36 M ⊕ 19 , while for the 10.8 d period planet (planet c) we obtained K c = 1.84 ± 0.33 m s −1 , corresponding to M c = 5.40 ± 0.98 M ⊕ . Moreover, we inferred the presence of two additional planets, with periods of 25.62 ± 0.04 days (planet d) and 77.23 ± 0.39 days (planet e), and robustly determined semi-amplitudes of K d = 3.06 ± 0.33 m s −1 (M d = 11.95 ± 1.28 M ⊕ ) and K e = 2.84 ± 0.41 m s −1 (M e = 16.0 ± 2.3 M ⊕ ). 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 ± 0.13 R ⊕ and R e = 2.67 ± 0.11 R ⊕ for planet d and e respectively. We checked the stability of the system with the newly determined configuration as described in Section 6, and all the 10 runs resulted in a MEGNO value of 2, indicating that the family of solutions is stable. As a final test, we performed a joint photometric and RV fit assuming a five-planet model including the 16 d period planet, and assuming that the two additional signals seen in the RVs were caused by two non-transiting planets, both in the case of a ∼ 78 d and ∼ 180 d external period planet. According to the Bayesian Information Criterion (BIC), the four-planet model was favoured with respect to both the five-planet models, having a ∆BIC 10 (Kass & Raftery 1995) in both cases (∆BIC 78d = 77, ∆BIC 180d = 84). 19 The final adopted value of K b and M b is the the weighed mean between the values obtained from the nightly offset method (Section 5.1) and from the joint photometric and RV fit (see Table 6 ).

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 days, with a normalised power of 0.20, that is, below the 1% 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 et al. (2015) , with no priors on the GP hyper-parameters, since we could not identify any activity-related signal in the ancillary datasets (see Section 3.3) 20 . 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 all the previous analyses and results, we adopt the parameters and configuration determined in this section as the representative ones for the TOI-561 system.

According to our analysis, TOI-561 hosts four transiting planets, including an USP planet, a ∼ 10.8 d period planet and two external planets with periods of ∼ 25.6 and ∼ 77.2 days. 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 4 or more planets 21 , and it is one of the few multi-planet 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 Figure 11 we show the position of TOI-561 b, c, d and e in the mass-radius diagram of exoplanets with masses and radiii measured with a precision better than 30%. The comparison with the theoretical massradius curves excludes an Earth-like composition (∼ 33% iron and 67% silicates) for all planets in the system. The density (ρ b = 3.0 ± 0.8 g cm −3 ) of the USP planet is consistent with a 50% (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% of the total mass. Even assuming the higher mass value inferred with the nightly offset method 22 (M b = 1.83 ± 0.39 M ⊕ , implying a density of ρ b = 3.5 ± 0.9 g cm −3 ), TOI-561 b would be located close to the 25% 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 5100 F ⊕ ), 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 with caution, since in this case it could lead to an overestimation of the bulk water content (Turbet et al. 2020) . TOI-561 c, with a density of ρ c ∼ 1.3 g cm −3 , is located above the threshold of a 100% 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 waterrich, 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% 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% 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. We estimate that a H-He envelope of ∼ 1.8% and ∼ 2.3% of the planet mass could explain the observed planetary properties. 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 ⊕ (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. photoevaporation (Chen & 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. 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 Figure 11 ). 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., ρ W47 e = 6.4 ± 0.6 g cm −3 (Vanderburg et al. 2017) and ρ 55Cnc e = 6.3 ± 0.8 g cm −3 (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 W47 e = 6.83 ± 0.66 M ⊕ (Vanderburg et al. 2017 ) and M 55Cn e = 8.08±0.31 M ⊕ (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.

As an additional remark, we note that the orbital inclinations of planets c, d and e are all consistent within 1σ, and that the difference with the inclination of the USP planet is of the order of ∆I ∼ 2.5 deg. According to the analysis of Dai et al. (2018) , when the innermost planet has a/R < 5, the minimum mutual inclination with other planets in the system often reaches values up to 5-10 deg, with larger period ratios (P c /P b > 5-6) implying an higher mutual inclination. Considering the large period ratio of TOI-561 (P c /P b ∼ 24) and the value of a b /R = 2.6, the measured ∆I ∼ 2.5 deg in this case is much lower that the expected inclination dispersion of 6.7 ± 0.7 deg 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.

Therefore, 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. Given the interesting composition of the planets in the system, we checked if the TOI-561 planets Figure 9 . Phase-folded RV fit with residuals from the joint fourplanet photometric and RV analysis. Planets b, c, d, and e are shown in blue, orange, red and purple, respectively. The reported errorbars include the jitter term, added in quadrature.

would be accessible targets for atmospheric characterisation 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-hour 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) 23 , this classifies TOI-561 b and c as high-quality atmospheric characterisation 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 ⊕ 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 c 8.7 M ⊕ ). Therefore, this estimate of the atmospheric characterisation feasibility should be used with caution, especially as the TSM metric has been conceived to prioritise targets for follow-up, and not to precisely determine the atmospheric transmission properties.

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 highprecision photometric (i.e. with TESS, that will re-observe TOI-561 in sector 35 -February/March 2021, 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 dataset, that could also unveil eventual additional long-period companions. In this framework, we 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 6 ), 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 0 of the USP planet (see Table 6 ), 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 Figure 12 . According to our simulation, TOI-561 d and e display an anti-correlated TTV signal, with a very long TTV period of ∼ 4850 days (∼ 13 yr), and TTV amplitudes 24 of ∼ 62 minutes (planet d) and ∼ 84 minutes (planet e). The anti-correlated 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 sec). 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 Figure 12 . According to our simulation, the first peak (dip) corresponds to the period between March-December 2020, while the second one will be between January-October 2026, i. e., corresponding to the time-spans between ∼ 400-700 and ∼ 2500-3000 days of integration in Figure 12 respectively. We remark that this calculation is performed assuming the T 0 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.

The multi-planetary nature of TOI-561 offers a unique opportunity for comparative exoplanetology. TOI-561 planets may be compared with the known population of multi-planet 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. 24 We calculated the TTV period and amplitudes computing the GLS periodogram of the simulated TTVs. Figure 11 . Mass-radius diagram for known exoplanets with mass and radius measurements more precise than 30%, 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 August 17, 2020. TOI-561d TOI-561e Figure 12 . Predicted TTV signal of TOI-561 d and e assuming our best-fitting model (see Table 6 ). The planets show a strong, anti-correlated signal. The signals of the USP planet (< 1 sec) and of planet c (< 1 min) are not reported. 

Handbook of Exoplanets

Protostars and Planets VI. p

IAU Symposium

Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series

Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series. p. 8

VizieR Online Data Catalog

Protostars and Planets VI. p

EXOFASTv2: A public, generalized, publication-quality exoplanet modeling code

Protostars and Planets VI. p

Kepler Data Processing Handbook

Software and Cyberinfrastructure for

Lightkurve: Kepler and TESS time series analysis in Python

isochrones: Stellar model grid package

Protostars and Planets VI. p

ARES + MOOG: A Practical Overview of an Equivalent Width (EW) Method to Derive Stellar Parameters

Proceedings of the National Academy of Science

Università degli Studi di Padova, Vicolo dell'Osservatorio 3, I-35122, Padova, Italy 2 INAF -Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, I-35122

Cambridge CB3 0HE, UK 5 Kavli Institute for Cosmology

DK-2800 Kgs. Lyngby, Denmark 11 Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA, 02138, USA 12 INAF -Osservatorio Astrofisico di Torino, Via Osservatorio 20, I-10025

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

HARPS-N observations and data products are available through the Data & Analysis Center for Exoplanets (DACE) at https://dace.unige.ch/. TESS data products can be accessed through the official NASA website https://heasarc. gsfc.nasa.gov/docs/tess/data-access.html. All underlying data are available either in the appendix/online supporting material or will be available via VizieR at CDS. This paper has been typeset from a T E X/L A T E X file prepared by the author.