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, M. López-Morales, V. Nascimbeni, T. G. Wilson, S. Udry, D. W. Latham, A. S. Bonomo, M. Damasso, X. Dumusque, J. M. Jenkins, C. Lovis, K. Rice, D. Sasselov, J. N. Winn, G. Andreuzzi, R. Cosentino, D. Charbonneau, L. Di Fabrizio, A. F. Martinez Fiorenzano, A. Ghedina, A. Harutyunyan, F. Lienhard, G. Micela, E. Molinari, I. Pagano, F. Pepe, D. F. Phillips, M. Pinamonti, G. Ricker, G. Scandariato, A. Sozzetti, C. A. Watson
MMNRAS , 1–20 (2020) Preprint 30 November 2020 Compiled using MNRAS L A TEX style file v3.0
An unusually low density ultra-short period super-Earth and threemini-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 , M. López-Morales , V. Nascimbeni , T. G. Wilson ,S. Udry , D. W. Latham , A. S. Bonomo , M. Damasso , X. Dumusque , J. M. Jenkins , C. Lovis ,K. Rice , , D. Sasselov , J. N. Winn , G. Andreuzzi , , R. Cosentino , D. Charbonneau ,L. Di Fabrizio , A. F. Martinez Fiorenzano , A. Ghedina , A. Harutyunyan , F. Lienhard , G. Micela ,E. Molinari , I. Pagano , F. Pepe , D. F. Phillips , M. Pinamonti , G. Ricker , G. Scandariato ,A. Sozzetti , C. A. Watson 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 characterisation of the planetary systemorbiting the late G dwarf TOI-561. After the identification of three transiting candidates by
TESS , we discovered two additionalexternal planets from RV analysis. RVs cannot confirm the outer
TESS transiting candidate, which would also make the systemdynamically unstable. We demonstrate that the two transits initially associated with this candidate are instead due to single transitsof 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 𝑃 b = .
45 d, mass 𝑀 b = . ± .
36 M ⊕ and radius 𝑅 b = . ± .
07 R ⊕ , and three mini-Neptunes:TOI-561 c, with 𝑃 c = .
78 d, 𝑀 c = . ± .
98 M ⊕ , 𝑅 c = . ± .
09 R ⊕ ; TOI-561 d, with 𝑃 d = . 𝑀 d = . ± . ⊕ , 𝑅 d = . ± .
13 R ⊕ ; and TOI-561 e, with 𝑃 e = . 𝑀 e = . ± . ⊕ , 𝑅 e = . ± .
11 R ⊕ . Having a density of3 . ± . − , TOI-561 b is the lowest density USP planet known to date. Our N-body simulations confirm the stability of thesystem and predict a strong, anti-correlated, long-term transit time variation signal between planets d and e. The unusual densityof the inner super-Earth and the dynamical interactions between the outer planets make TOI-561 an interesting follow-up target. Key words: planets and satellites: detection – planets and satellites: composition – star: individual: TOI-561 (TIC 377064495,
Gaia
DR2 3850421005290172416) – techniques: photometric – techniques: radial velocities
The
Transiting Exoplanet Survey Satellite ( TESS , Ricker et al. 2014)is a NASA all-sky survey designed to search for transiting planetsaround bright and nearby stars, and particularly targeting stars thatcould reveal planets with radii smaller than Neptune. Since the be-ginning of its observations in 2018,
TESS has discovered more than66 exoplanets, including about a dozen multi-planet systems (e.g.Dragomir et al. 2019; Dumusque et al. 2019; Günther et al. 2019).Multi-planet systems, orbiting the same star and having formed fromthe same protoplanetary disc, offer a unique opportunity for compar-ative planetology. They allow for investigations of the formation andevolution processes, i.e. through studies of relative planet sizes andorbital separations, orbital inclinations relative to the star’s rotation ★ E-mail: [email protected] axis, mutual inclination of the orbits, etc. In order to obtain a completecharacterisation of a system, knowledge of the orbital architecture andthe bulk composition of the planets are essential. To obtain such in-formation, transit photometry needs to be combined with additionaltechniques 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 systemshave been discovered by the
Kepler space telescope (Borucki et al.2010), which has led to an unprecedented knowledge of the ensembleproperties of multiple systems (e.g. Latham et al. 2011; Millhollandet al. 2017; Weiss et al. 2018), their occurrence rate (e.g. Fressinet al. 2013), and their dynamical configurations (e.g. Lissauer et al.2011; Fabrycky et al. 2014). However, many of the
Kepler targets aretoo faint for RV follow-up, so most of the planets do not have a massmeasurement, preventing a comprehensive understanding of theirproperties, and of the planetary system. Thanks to the
TESS satellite, © a r X i v : . [ a s t r o - ph . E P ] N ov G. Lacedelli et al. which targets brighter stars, an increasing number of candidates suit-able for spectroscopic follow-up campaigns are being discovered.These new objects will increase the number of well characterisedsystems, and will provide a valuable observational counterpart tothe theoretical studies on the formation and evolution processes ofplanetary systems (e.g. Morbidelli et al. 2012; Raymond et al. 2014;Helled et al. 2014; Baruteau et al. 2014, 2016; Davies et al. 2014).In this paper, we combine
TESS photometry (Section 2.1) andhigh precision RVs gathered with the HARPS-N spectrograph (Sec-tion 2.2) to characterise the multi-planet system orbiting the starTOI-561. The
TESS pipeline identified three candidate planetary sig-nals, namely an ultra-short period (USP) candidate ( 𝑃 ∼ .
45 days),and two additional candidates with periods of ∼ . ∼ . TESS photometry
TOI-561 was observed by
TESS in two-minute cadence mode duringobservations of sector 8, between 2 February and 27 February 2019.The astrometric and photometric parameters of the star are listed inTable 1. Considering the download time, and the loss of 3 .
26 daysof data due to an interruption in communications between the in-strument and the spacecraft that occurred during sector 8 , a total of20 .
22 days of science data were collected. The photometric observa-tions for TOI-561 were reduced by the Science Processing OperationsCenter (SPOC) pipeline (Jenkins et al. 2016; Jenkins 2020), whichdetected three candidate planetary signals, with periods of 10 . . . . . .
2, respectively. For our photometric analysis, we used thelight curve based on the Pre-search Data Conditioning Simple Aper-ture Photometry (PDCSAP, Smith et al. 2012; Stumpe et al. 2012,2014). We downloaded the two-minute cadence PDCSAP light curvefrom the Mikulski Archive for Space Telescopes (MAST) , and re-moved all the observations encoded as NaN or flagged as bad-quality(
DQUALITY>0 ) points by the SPOC pipeline . We performed outliersrejection by doing a cut at 3 𝜎 for positive outliers and 5 𝜎 (i. e. largerthan the deepest transit) for negative outliers. We removed the lowfrequency trends in the light curve using the biweight time-windowed See
TESS
Data Release Notes: Sector 8, DR10 ( https://archive.stsci.edu/tess/tess_drn.html ). https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html https://archive.stsci.edu/missions/tess/doc/EXP-TESS-ARC-ICD-TM-0014.pdf Table 1.
Astrometric and photometric parameters of TOI-561Property Value Source
Other target identifiers
TIC 377064495 A
Gaia
DR2 3850421005290172416 B2MASS J09524454+0612589 C
Astrometric parameters
RA (J2015.5; h:m:s) 09:52:44.44 BDec (J2015.5; d:m:s) 06:12:57.97 B 𝜇 𝛼 (mas yr − ) − . ± .
088 B 𝜇 𝛿 (mas yr − ) − . ± .
094 BSystemic velocity (km s − ) 79 . ± .
56 BParallax 𝑎 (mas) 11 . ± . . + . − . D Photometric parametersTESS (mag) 9 . ± .
006 A
Gaia (mag) 10 . ± . V (mag) 10 . ± .
006 A B (mag) 10 . ± .
082 A J (mag) 8 . ± .
020 C H (mag) 8 . ± .
055 C K (mag) 8 . ± .
019 C W1 (mag) 8 . ± .
023 E W2 (mag) 8 . ± .
020 E W3 (mag) 8 . ± .
023 E W4 (mag) 7 . ± .
260 ESpectral type G9V FA)
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 color indexes. 𝑎 Gaia
DR2 parallax is corrected by + ± 𝜇 as (with the error addedin quadrature) as suggested by Khan et al. (2019). slider implemented in the wotan package (Hippke et al. 2019), witha window of 1 . TESS light curve, we per-formed an iterative transit search on the detrended light curve usingthe Transit Least Squares (TLS) algorithm (Hippke & Heller 2019).The first three significant identified signals nicely matched the
TESS suggested periods ( 𝑃 TLS = .
78 d, 0 .
44 d, 16 .
28 d).In addition, we also extracted the 30-minutes cadence light curvefrom the
TESS
Full-Frame Images (FFIs) using the
PATHOS pipeline(Nardiello et al. 2019), in order to obtain an independent confirmationof the detected signals (Section 4).
MNRAS000
MNRAS000 , 1–20 (2020)
OI-561: an USP super-Earth and three mini-Neptunes We collected 82 spectra using HARPS-N at the TelescopioNazionale Galileo (TNG), in La Palma (Cosentino et al. 2012, 2014),with the goal of precisely determining the masses of the three can-didate planets and to search for additional planets. The observationsstarted on November 17, 2019 and ended on June 13, 2020, with aninterruption between the end of March and the end of April due tothe shut down of the TNG because of Covid-19. In order to preciselycharacterise the signal of the USP candidate, we collected 6 pointsper night on February 4 and February 6, 2020, thus covering thewhole phase curve of the planet, and two points per night (whenweather allowed) during the period of maximum visibility of thetarget (February-March 2020). The exposure time was set to 1800seconds, which resulted in a S/N at 550 nm of 77 ±
20 (median ± stan-dard deviation) and a measurement uncertainty of 1 . ± . − .We reduced the data using the standard HARPS-N Data ReductionSoftware (DRS) using a G2 flux template (the closest match to thespectral type of our target) to correct for variations in the flux dis-tribution as a function of the wavelength, and a G2 binary mask tocompute the cross-correlation function (CCF, Baranne et al. 1996;Pepe et al. 2002). All the observations were gathered with the secondfibre of HARPS-N illuminated by the Fabry-Perot calibration lamp tocorrect for the instrumental RV drift, except for the night of May 31,2020. This observation setting prevented us from using the secondfibre to correct for Moon contamination. However, we note that thedifference between the systemic velocity of the star and the Moonis always greater than 15 km s − , therefore preventing any contam-ination of the stellar CCF (as empirically found by Malavolta et al.2017a and subsequently demonstrated through simulations by Royet al. 2020), as the average full width at half maximum (FWHM) ofthe CCF for TOI-561 is 6 . ± .
004 km s − .The RV data with their 1 𝜎 uncertainties and the associated activityindices (see Section 3.3 for more details) are listed in Table 2. Be-fore proceeding with the analysis, we removed from the total dataset5 RV measurements, with associated errors greater than 2 . − from spectra with S/N <
35, that may affect the accuracy of ourresults. The detailed procedure performed to identify these points isdescribed in Appendix B1.
We derived the photospheric stellar parameters using three differenttechniques: the curve-of-growth approach, spectral synthesis match,and empirical calibration.The first method minimizes the trend of iron abundances (obtainedfrom the equivalent width, EW, of each line) with respect to excita-tion potential and reduced EW respectively, to obtain the effectivetemperature and the microturbulent velocity, 𝜉 t . The gravity log g isobtained by imposing the same average abundance from neutral andionised iron lines. We obtained the EW measurements using ARESv2 (Sousa et al. 2015). We used the local thermodynamic equilibrium(LTE) code MOOG (Sneden 1973) for the line analysis, together with
62 spectra were collected within the Guaranteed Time Observations (GTO)time (Pepe et al. 2013), while the remaining 20 spectra were collected withinthe A40_TAC23 program. Available at Available at the
ATLAS9 grid of stellar model atmosphere from Castelli & Ku-rucz (2003). The whole procedure is described in more detail inSousa (2014). We performed the analysis on a co-added spectrum(S/N > 𝑇 eff = ±
69 K, log g = . ± . [ Fe / H ] = − . ± .
05 and 𝜉 t = . ± .
08 km s − .The spectral synthesis match was performed using the Stellar Pa-rameters Classification tool ( SPC , Buchhave et al. 2012, 2014). It de-termines effective temperature, surface gravity, metallicity and linebroadening by performing a cross-correlation of the observed spectrawith a library of synthetic spectra, and interpolating the correlationpeaks to determine the best-matching parameters. For technical rea-sons, we ran the
SPC on the 62 GTO spectra only : the S/N is sohigh that the spectra are anyway dominated by systematic errors, andincluding the A40TAC_23 spectra would not change the results. Weaveraged the values measured for each exposure, and we obtained 𝑇 eff = ±
50 K, log g = . ± . [ M / H ] = − . ± .
08 and 𝑣 sin 𝑖 < − .We finally used CCFpams , a method based on the empirical cal-ibration of temperature, metallicity and gravity on several CCFsobtained with subsets of stellar lines with different sensitivity to tem-perature (Malavolta et al. 2017b). We obtained 𝑇 eff = ±
70 K,log g = . ± .
15 and [ Fe / H ] = − . ± .
05, after applying thesame gravity and systematic corrections as for the EW analysis.We list the final spectroscopic adopted values, i. e., the weightedaverages of the three methods, in Table 3.From the co-added HARPS-N spectrum, we also derived the chem-ical abundances for several refractory elements (Na, Mg, Si, Ca, Ti,Cr, Ni). We used the ARES+MOOG method assuming LTE, as de-scribed earlier. The reference for solar values was taken from Asplundet 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 mag-nesium, silicon, and titanium to represent the alpha-elements and theiron abundance from the ARES+MOOG method (for consistency),we find that [ 𝛼 / Fe ] = . For each set of photospheric parameters, we determined the stellarmass and radius using isochrones (Morton 2015), with posteriorsampling performed by
MultiNest (Feroz & Hobson 2008; Ferozet al. 2009, 2019). We provided as input the parallax of the targetfrom the
Gaia
DR2 catalogue, after adding an offset of + ± 𝜇 as(with the error added in quadrature to the parallax error) as suggestedby Khan et al. (2019), plus the photometry from the TICv8, 2MASSand WISE (Table 1). We used two evolutionary models, the MESAIsochrones & 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 𝜎 𝑇 eff =
70 K, 𝜎 log g = . 𝜎 [ Fe / H ] = .
05 (except for
SPC , where we kept theoriginal error of 0 .
08) as a good estimate of the systematic errorsregardless of the internal error estimates, to avoid favouring one SPC runs on a server with access to GTO data only, and the required tech-nical effort to enable the use of A40_TAC23 data, complicated by the globalCovid-19 sanitary emergency, was not justified by the negligible scientificgain. Available at https://github.com/LucaMalavolta/CCFpams
MNRAS , 1–20 (2020)
G. Lacedelli et al.
Table 2.
HARPS-N Radial Velocity Measurements.BJD
TDB RV 𝜎 R 𝑉 BIS FWHM 𝑉 asy Δ 𝑉 log R (cid:48) HK H 𝛼 ( d ) (m s − ) (m s − ) (m s − ) (km s − ) (km s − ) (dex)2458804.70779 79700.63 1.27 -39.98 6.379 0.048 -0.039 -5.005 0.2032458805.77551 79703.74 0.97 -36.25 6.380 0.049 -0.036 -4.984 0.2002458806.76768 79701.71 1.05 -31.81 6.378 0.045 -0.033 -5.000 0.200... ... ... ... ... ... ... ... ...This table is available in its entirety in machine-readable form. technique over the others when deriving the stellar mass and radius.We also imposed an upper limit on the age of 13 . 𝑀 ★ = . ± .
018 M (cid:12) and 𝑅 ★ = . ± .
007 R (cid:12) . We derivedthe stellar density 𝜌 ★ = . ± . 𝜌 (cid:12) ( 𝜌 ★ = . ± .
056 g cm − )directly from the posterior distributions of 𝑀 ★ and 𝑅 ★ .We summarise the derived astrophysical parameters of the starin Table 3, which also reports temperature, gravity and metallic-ity obtained from the posteriors distributions resulting from the isochrone fit. A lower limit on the age of ∼
10 Gyr is obtainedconsidering the 15 . EXOFASTv2 (Eastman et al. 2019),assuming the photometric parameters in Table 1 and the spectro-scopic parameters in Table 3, using only the MIST evolutionary set,returned a lower limit on the age of 5 Gyr, while all the other pa-rameters were consistent with the results quoted in Table 3. Thus,we decided to assume 5 Gyr as a conservative lower limit for theage of the system. The old stellar age and the sub-solar metallic-ity suggest that TOI-561 may belong to an old Galactic population,an hypothesis that is also supported by our kinematic analysis. Infact, we derived the Galactic space velocities using the astromet-ric properties reported in Table 1. For the calculations we used the astropy package, and we assumed the
Gaia
DR2 radial velocityvalue of 79 .
54 km s − , obtaining the heliocentric velocity compo-nents ( 𝑈, 𝑉, 𝑊 ) = (− . , − . , . ) km s − , in the directions ofthe Galactic center, Galactic rotation, and north Galactic pole, re-spectively. The derived 𝑈𝑉𝑊 velocities point toward a thick-diskstar, as confirmed by the probability membership derived followingBensby et al. (2014), that implies a ∼
70% probability that the starbelongs to the thick disc, a ∼
29% probability of being a thin-discstar and a ∼ . The low value of the log R (cid:48) HK index ( − . ± . 𝐵 − 𝑉 = . (cid:39)
86 pc, the lack of interstellar absorption near the Na D doubletin the HARPS-N co-added spectrum, and the total extinction in the V band from the isochrone fit (0 . (cid:48) HK index(Fossati et al. 2017). Nevertheless, it is important to check whetherthe star is showing any sign of activity in all the activity diagnostics atour disposal. In addition to the log R (cid:48) HK index, FWHM, and bisectorspan (BIS) computed by the HARPS-N DRS, we included in ouranalysis the 𝑉 asy (Figueira et al. 2013) and Δ 𝑉 (Nardetto et al. 2006)asymmetry indicators, as implemented by Lanza et al. (2018), and thechromospheric activity indicator H 𝛼 (Gomes da Silva et al. 2011). Table 3.
Derived astrophysical stellar parameters.Parameter Value Unit 𝑇 eff 𝑎 spec ±
70 Klog g 𝑎 spec . ± .
12 - [ Fe / H ] 𝑎 spec − . ± .
05 - 𝑇 eff 𝑏 + − Klog g 𝑏 . ± .
01 - [ Fe / H ] 𝑏 − . + . − . - 𝑅 ★ . ± .
007 R (cid:12) 𝑀 ★ . ± .
018 M (cid:12) 𝜌 ★ . ± . 𝜌 (cid:12) 𝜌 ★ . ± .
056 g cm − 𝐴 𝑉 . + . − . mag 𝑣 sin 𝑖 < − age 𝑐 > (cid:48) HK − . ± .
012 - [ Na/H ] − . ± .
06 - [ Mg/H ] − . ± .
05 - [ Si/H ] − . ± .
05 - [ Ca/H ] − . ± .
06 - [ Ti/H ] − . ± .
03 - [ Cr/H ] − . ± .
08 - [ Ni/H ] − . ± .
04 - 𝑎 Weighted average of the three spectroscopicmethods. 𝑏 Value inferred from the isochrone fit. 𝑐 Conservative lower limit.
The Generalized Lomb-Scargle (GLS, Zechmeister & Kürster2009) periodograms of the above-mentioned indexes, computedwithin the frequency range 0 . . − , i. e., 2–2000 days, areshown in Figure 1, together with the periodograms of the RVs and TESS photometry. For each periodogram, we also report the powerthreshold corresponding to a False Alarm Probability (FAP) of 1%and 0 . (cid:39)
25 days, (cid:39)
180 days, (cid:39)
10 days (corresponding to one of the transiting planetcandidates), and (cid:39)
78 days, ordered decreasingly according to theirpower. None of these peaks has a counterpart in the activity diag-nostics here considered, as no signals with a FAP lower than 2 . TESS light curve identified a periodicity around 3 . .
13 ppt and a power of 0 . . MNRAS000
13 ppt and a power of 0 . . MNRAS000 , 1–20 (2020)
OI-561: an USP super-Earth and three mini-Neptunes would be extremely atypical for a star older than 1 Gyr (e.g. Douglaset al. 2019), and in contrast with the lack of any signal in all theother above-mentioned activity indicators. Indeed, the rotational pe-riod estimated from the log R (cid:48) 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 rota-tional period should be considered as a rough estimate, also becausethese calibrations are not well tested for old and alpha-enhanced starslike TOI-561. Further evidence against a ∼ . 𝑣 sin 𝑖 ( < − ), that suggestsa rotational period > . ◦ . In any case, we verified with a pe-riodogram analysis that our light curve flattening procedure correctlyremoved the here identified signal at 3 . TESS light curve(with the transits filtered out), and the ASAS-SN V and g photom-etry (Shappee et al. 2014; Kochanek et al. 2017), after applying a5- 𝜎 filtering, but no significant periodicity could be identified. Aperiodogram analysis of the ASAS-SN light curves in each band, ei-ther by taking the full dataset or by analysing each observing seasonindividually, confirmed these results.In conclusion, if any activity is present, its signature must be below0 . <
30 days),and 20 ppt in the long term period (magnetic cycles, >
100 days),from the RMS of
TESS and ASAS-SN photometry respectively. In-cidentally, the former is close to the photometric variations of theSun during the minimum at the end of Solar Cycle 25, when the Sunalso reached a log R (cid:48) HK very close to the one measured for TOI-561(Collier Cameron et al. 2019; Milbourne et al. 2019). By comparingour target to the Sun, and in general by taking into account the resultsof Isaacson & Fischer (2010), it is expected that the contributionto the RVs due to the magnetic activity of our star is likely below1-2 m s − . Since this value is quite close to the median internal errorof our RVs, no hint of the rotational period is provided by eitherthe photometry or the spectroscopic activity diagnostics, and the lowactivity level is consistent with our derived stellar age ( > 𝜎 jitter ). Previous experience with
Kepler shows that candidates in multi-ple systems have a much lower probability of being false positives(Latham et al. 2011; Lissauer et al. 2012). Nevertheless, it is alwaysappropriate to perform a series of checks in order to exclude thepossibility 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 re-normalised unit weight error (RUWE) of 1 .
1, indicating that thesingle-star model provides a good fit to the astrometric observations.This likely excludes the presence of a massive companion that couldcontribute to the star’s orbital motion in the
Gaia
DR2 astrometry, afact that agrees with the absence of long-term trends in our RVs (seeSection 6.1).Moreover, the overall RV variation below 25 m s − and the shapeof the CCFs of our HARPS-N spectra exclude the eclipsing binaryscenario, which would be the most likely alternative explanation forthe USP planet.A further confirmation comes from the speckle imaging on theSouthern Astrophysical Research (SOAR) telescope that Ziegler et al. . . . . . G L Sp o w e r RV Period (d) . . . . . . G L Sp o w e r TESS . . . . . . G L Sp o w e r FWHM . . . . . . G L Sp o w e r BIS . . . . . . G L Sp o w e r V asy . . . . . . G L Sp o w e r ∆ V . . . . . . G L Sp o w e r log R HK . . . . . . G L Sp o w e r H α . . . . . Frequency (d − ) . . . . . . G L Sp o w e r Window function
Figure 1.
GLS periodogram of the RVs, the
TESS photometry (PDCSAP)and the spectroscopic activity indexes under analysis. The main peak of eachperiodogram is highlighted with an orange vertical line. The grey vertical linesrepresent the signals corresponding to the transit-like signals with periods10 . . (cid:39) (cid:39)
78 and (cid:39)
180 days. The dashed and dotted horizontal lines showthe 1% and 0 .
1% FAP levels, respectively. The
TESS periodogram shows aseries of peaks below 10 days, unlikely to be associated with stellar activitygiven the old age of the star. The FWHM and the log R (cid:48) HK periodograms havethe main peak at 244 and 220 days, respectively, so there is no correspondencewith the 180 days signal. Moreover, both of them are below the 1% FAP. Thebottom panel shows the window function of the data.MNRAS , 1–20 (2020) G. Lacedelli et al. (2020) performed on some of the
TESS planet candidate hosts. Ac-cording to their analysis (see Tables 3 and 6 therein), no companionis detected around TOI-561 (being the resolution limit for the star0 .
041 arcsec, and the maximum detectable Δ mag at separation of1 arcsec 4 .
76 mag). Still, the 21 arcsec
TESS pixels and the few-pixels wide point spread function (PSF) can cause the light fromneighbours over an arc-minute away to contaminate the target lightcurve. In the case of neighbouring eclipsing binaries (EBs), eclipsescan be diluted and mimic shallow planetary transits. For example,events at ∼ TESS aperture witha 0 .
5% eclipse, but no more than 7 magnitudes fainter. This conditionis not satisfied in our case, as the only three sources within 100 arcsecfrom TOI-561 are all fainter than 𝑇 = .
25 mag and at a distancegreater than 59 arcsec, according to the
Gaia
DR2 catalogue.An independent confirmation was provided by the analysis of thein-/out-of-transit difference centroids on the
TESS
FFIs (Figure 2),adopting the procedure described in Nardiello et al. (2020). Theanalysis of the in-/out-of transit stacked difference images confirmsthat, within a box of 10 ×
10 pixels ( ∼ ×
200 arcsec ) centredon TOI-561, the transit events associated with candidates .01 and .03occur on our target star, while candidate .02 has too few in-transitpoints in the 30-minute cadence images for this kind of analysis — inany case, its planetary nature will be confirmed by the RV signal ofTOI-561 in Section 6.Finally, in order to exclude the possibility that the transit-like fea-tures were caused by instrumental artefacts, we performed some addi-tional checks on the light curve. We visually inspected the FFIs to spotpossible causes (including instrumental effects) inducing transit-likefeatures, and we could not find any. We re-extracted the short cadencelight curve using the python package lightkurve (Lightkurve Col-laboration et al. 2018) with different photometric masks and aper-tures, and we corrected them by using the TESS
Cotrending BasisVectors (CBVs); the final results were in agreement with the
TESS -released PDCSAP light curve. We checked for systematics in everylight curve pixel, and we found none. Ultimately, we checked for cor-relations between the flux, the local background, the (X,Y)-positionfrom the PSF-fitting, and the FWHM, with no results. Therefore, weconclude that all the transit-like features in the light curve are realand likely due to planetary transits.
We performed the analysis presented in the next sections using
Py-ORBIT (Malavolta et al. 2016, 2018), a convenient wrapper for theanalysis of transit light curves and radial velocities.In the analysis of the light curve, for each planet we fitted thecentral time of transit ( 𝑇 ), period ( 𝑃 ), planetary to stellar radiusratio ( 𝑅 p / 𝑅 ★ ), and impact parameter 𝑏 . In order to reduce compu-tational time, we set a narrow, but still uninformative, uniform priorfor period and time of transit, as defined by a visual inspection. Wefitted a common value for the stellar density 𝜌 ★ , imposing a Gaus-sian prior based on the value from Table 3. We included a quadraticlimb-darkening law with Gaussian priors on the coefficients 𝑢 , 𝑢 ,obtained through a bilinear interpolation of limb darkening profilesby Claret (2018) . We initially calculated the standard errors on 𝑢 , https://github.com/KeplerGO/lightkurve https://github.com/LucaMalavolta/PyORBIT , version 8.1 https://vizier.u-strasbg.fr/viz-bin/VizieR?-source=J/A+A/618/A20 Figure 2.
In-/out-of-transit difference centroid analysis of the transit eventsassociated with the candidates TOI-561.01 (transit 2 and 3) and TOI-561.03(transit 1 and 4). The star is centred at (0,0), and the grey circles are all theother stars in the
Gaia
DR2 catalogue, with dimension proportional to theirapparent magnitude. 𝑢 using a Monte Carlo approach that takes into account the errorson 𝑇 eff and log g as reported in Table 3, obtaining 𝑢 = . ± . 𝑢 = . ± . .
05. In the fit we employedthe parametrization ( 𝑞 , 𝑞 ) introduced by Kipping (2013). Finally,we included a jitter term to take into account possible TESS sys-tematics and short-term stellar activity noise. We assumed uniform,uninformative priors for all the other parameters, although the prioron the stellar density will inevitably affect the other orbital parame-ters. All the transit models were computed with the batman package(Kreidberg 2015), with an exposure time of 120 seconds and anoversampling factor of 10 (Kipping 2010).In the analysis of the radial velocities, we allowed the periods tospan between 2 and 200 days (i. e., the time span of our dataset) forthe non-transiting planets, while we allowed the semi-amplitude 𝐾 to vary between 0 .
01 and 100 m s − for all the candidate planets.These two parameters were explored in the logarithmic space. Forthe transiting candidates, we used the results from the photometricfit (see Appendix A) to impose Gaussian priors on period and timeof transit on RV analysis alone, while using the same uninformativepriors as for the photometric fit when including the photometric dataas well.For all the signals except the USP candidate, we assumed eccentricorbits with a half-Gaussian zero-mean prior on the eccentricity (withvariance 0 . MultiNest nested-sampling algorithm (Feroz & Hobson 2008; Feroz et al. 2009,2019) with the Python wrapper pyMultiNest (Buchner, J. 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; Skilling 2006; Speagle 2020), which allowed forthe computation of the Bayesian evidence in a reasonable amount
MNRAS000
MNRAS000 , 1–20 (2020)
OI-561: an USP super-Earth and three mini-Neptunes of time thanks to its easier implementation of the multi-processingmode. We performed a series of test on a reduced dataset, and we ver-ified that the two algorithms provided consistent results with respectto each other. For all the analyses, we assumed 1000 live points anda sampling efficiency of 0 .
3, including a jitter term for each datasetconsidered in the model.Global optimisation of the parameters was performed using thedifferential evolution code
PyDE . The output parameters were usedas a starting point for the Bayesian analysis performed with the em-cee package (Foreman-Mackey et al. 2013), a Markov chain MonteCarlo (MCMC) algorithm with an affine invariant ensemble sampler(Goodman & Weare 2010). We ran the chains with 2 𝑛 dim walkers,where 𝑛 dim is the dimensionality of the model, for a number of stepsadapted to each fit, checking the convergence with the Gelman-Rubinstatistics (Gelman & Rubin 1992), with a threshold value of ˆ 𝑅 = . Before proceeding with a global analysis, we checked whether wecould independently recover the signals identified by the
TESS pipeline (Section 2.1) in our RV data only. The periodogram anal-ysis of the RVs in Section 3.3 highlighted the presence of severalpeaks not related to the stellar activity. In particular, an iterativefrequency search, performed subtracting at each step the frequencyvalues previously identified, supplied the frequencies 𝑓 = .
039 d − ( 𝑃 (cid:39) . 𝑓 = .
006 d − or 0 .
013 d − ( 𝑃 (cid:39)
170 d or (cid:39)
78 d)with the two frequencies being related to each other (i. e., removingone of them implies the vanishing of the other one), 𝑓 = .
093 d − ( 𝑃 (cid:39) . 𝑓 = .
239 d − ( 𝑃 (cid:39) .
45 d, corresponding to the TOI-561.02 can-didate). After removing these four signals, no other clear dominantfrequency emerged in the residuals. Since any attempt to perform a fitof the RVs to characterise the transiting candidates without account-ing for additional dominant signals would lead to unreliable results,we decided to test the presence of additional planets in a Bayesianframework. We considered four models, the first one (Model 0) as-suming the three transiting candidates only,i. e., TOI-561.01, .02,.03, and then including an additional planet in each of the successivemodels,i. e., TOI-561.01, .02, .03 plus one (Model 1), two (Model 2)and three (Model 3) additional signals, respectively. We computedthe Bayesian evidence for each model using the
MultiNest nested-sampling algorithm, following the prescriptions as specified in Sec-tion 5. We report the obtained values in Table 4. According to thisanalysis, 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 Δ ln Z >
10 (Kass& Raftery 1995), both compared to the case with one or no addi-tional signals. In the case of a third additional signal (Model 3), thedifference with respect to the two-signal model was less than 2, indi-cating that there was no strong evidence to favour this more complexmodel over the simpler model with two additional signals only (Kass https://github.com/hpparvi/PyDE Table 4.
Logarithmic Bayesian evidences for the different models underexam. Model 0 corresponds to the model with no additional RVs signal otherthan 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 transitingcandidates plus one, two and three additional planets, respectively. All thevalues 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).Model 0 Model 1 Model 2 Model 3ln Z . ± . . ± . . ± . ± . & Raftery 1995). We repeated the analysis first including a linearand then a quadratic trend in each of the four models. In all cases,the Bayesian evidence systematically disfavoured the presence of anytrend .The first additional signal was associated with a candidate with 𝑓 (cid:39) .
04 d − ( 𝑃 (cid:39) . MultiNest run highlighted the presence of two clusters of solutions,peaked at about 𝑓 = .
013 d − or 0 .
013 d − , i. e., 𝑃 =
78 and 180days respectively. The frequency analysis confirmed that the signalsare aliases of each other, since when we subtract one of them, the otherone also disappears. The alias peak is visible in the low-frequencyregime of the spectral window (Figure 1, bottom panel). We shouldalso consider that the longer period is close to the time baseline ofour 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 days, and thenbetween 100 and 200 days. The Bayesian evidence slightly favouredthe solution with 𝑃 ∼
78 d, even if not with strong significance( Δ ln Z (cid:39) ∼
16 days, that is, the transiting can-didate TOI-561.03. Therefore, in order to test our ability to recoverthe planetary signals, we performed a series of injection/retrievalsimulations, thoroughly explained in Appendix B2. The results ofthis injection/retrieval test are summarised in Figure 3. We foundthat the injected RV amplitude of .01 is not significantly affectingthe retrieved value for .03, i. e. the cross-talk between the two sig-nals is negligible. We verified that the same conclusion applies tothe other signals as well. More importantly, any attempt to retrieve anull signal at the periodicity of the candidate planet .03 would resultin an upper limit of ≈ . − as we actually observe with thereal dataset, when exploring the 𝐾 parameter in logarithmic space.Any signal equal or higher than 1 m s − would have been detected( > 𝜎 ), even if marginally. A signal with amplitude of 0 . − would not lead to the detection of the planet (intended as a 3- 𝜎 de-tection), but the retrieved posterior is expected to differ substantiallyfrom the observed one, especially on the lower tail of the distribution.We conclude that the planetary candidate TOI-561.03 is undetectedin our RV dataset, with an upper limit on the semi-amplitude of0 . − ( 𝑀 p < . ⊕ ). For the model with three additional signals and a quadratic trend, thecalculation of the Bayesian evidence did not converge.MNRAS , 1–20 (2020)
G. Lacedelli et al. K . po s t e r i o r s K . i n j = . m s - K . i n j = . m s - K . i n j = . m s - K .03inj =0.0 ms -1 K .03inj =0.5 ms -1 K .03inj =1.0 ms -1 K .03inj =1.5 ms -1 K . [ m s - ] K .03 [ms -1 ] Figure 3.
Posterior distributions (in the top panels, the blue, red and greenlines respectively) of the retrieved RV signal of TOI-561.03 according todifferent injected values for the RV semi-amplitudes of candidates .01 and.03. The black line in the top panels corresponds to the observed posterior ofthe RV semi-amplitude of candidate .03. Median and 1- 𝜎 values are markedwith vertical dashed and dotted lines respectively. Given the non-detection of the planetary candidate TOI-561.03 in theRV data, we investigated more closely the transit-like features asso-ciated with this candidate in the
TESS light curve, at 𝑇 (cid:39) . 𝑇 (cid:39) . (Hodges1958) on the residuals of transit 1 and 4 suggests that the two residualsamples are not drawn from the same distribution (threshold level 𝛼 = .
05, statistics KS = . 𝑝 − value (cid:28) . 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 associatedwith TOI-561.03 could indeed be due to the two additional planetsinferred from the RV analysis. To check our hypothesis, we first All the 𝑇 s in this section are expressed in BJD-2457000. We used the Python version implemented in scipy.stats.ks_2samp . Figure 4.
Transit 1 ( 𝑇 (cid:39) . 𝑇 (cid:39) . TESS detrended light curve associated with the candidate TOI-561.03. The best-fitting transit model from the three-planet model photometric fit is over-plotted(black solid line). The black dots are the data points binned over 15 minutes.With respect to transit 1, the duration of transit 4 looks underestimated bythe global model, with a systematic offset in the residuals, especially in thepre-transit phase. analysed the RV dataset with a model encompassing four planets, ofwhich only .01 and .02 have period and time of transit constrainedby
TESS . In other words, we performed the same RV analysis asdescribed in Appendix B2, but without including TOI-561.03 in themodel. We repeated the analysis twice in order to disentangle theperiodicity at 78 d from its alias at 180 d, and vice versa . We usedthe posteriors of the fit to compute the expected time of transit of theouter planets. We then performed two independent fits of transit 1 and4 with
PyORBIT , following the prescriptions as specified in Section 5.We imposed a lower boundary on the period of 22 days, in order toexclude the periods that would imply a second transit of the sameplanet in the
TESS light curve, and an upper limit of 200 days. Asa counter-measure against the degeneracy between eccentricity andimpact parameter in a single-transit fit, we kept the Van Eylen et al.(2019) eccentricity prior knowing that high eccentricities for sucha compact, old system are quite unlikely (Van Eylen et al. 2019).Finally we compared the posteriors of period and time of transitfrom the photometric fit with those from radial velocities, knowingthat the former will provide extremely precise transit times, but abroad distribution in period, while RVs give us precise periods, but
MNRAS000
MNRAS000 , 1–20 (2020)
OI-561: an USP super-Earth and three mini-Neptunes Transit N o r m a li ze d po s t e r i o r N o r m a li ze d po s t e r i o r [BJD - 2457000 d]1500 1520 1540 1560 Figure 5.
Comparison between period (left panels) and 𝑇 (right panels)obtained from the RV fit and from the fit of each single transit. Top andbottom panels refer to transit 1 and 4, respectively. Each panel shows theposterior distribution of the analysed parameter, and the shaded area indicatesthe region within the 68 . little information on the transit times. The results are summarisedin Figure 5: the 25 . ± . . + . − . d signal is close to the main peak in transit 4period distribution. Moreover, Figure 5 definitely confirms that boththe conjunction times inferred from the RV fit corresponding tothe ∼
25 and ∼
78 days signals, respectively 𝑇 = + − d and 𝑇 = + − d, are consistent with the (much more precise) 𝑇 sinferred from the individual fit of transit 1 ( 𝑇 = . ± .
004 d)and 4 ( 𝑇 = . ± .
006 d) respectively. Regarding the alias at182 ± 𝑇 = ±
13 dthat is derived from our analysis is not compatible with any of thetransits in the
TESS light curve. We also note that the proportion ofthe orbital period covered by the
TESS photometry is ∼ . ∼
25 d and ∼
78 d detected in the RV data, and the 180 d signal is considered analias of the 78 d signal.Given this final configuration, hereafter we will refer to the planetswith period ∼ . ∼ . ∼
25 and ∼
78 days as planets b, c, dand e, respectively.
Given the presence of two single-transit planets in our data, a jointphotometric and RV modelling is necessary in order to characterise the orbital parameters of all members of the TOI-561 system inthe best possible way. We considered a four-planet model, with acircular orbit for the USP planet and allowing nonzero-eccentricityorbits for the others. We performed the
PyORBIT fit as specified inSection 5, running the chains for 150 000 steps, and discarding thefirst 50 000 as burn-in. We summarise the results of our best-fittingmodel in Table 5, and show the transit models, the phase folded RVs,and the global RV model in Figures 6, 7, and 8 respectively. Weobtained a robust detection of the USP planet (planet b) RV semi-amplitude ( 𝐾 b = . ± .
32 m s − ), that corresponds to a mass of 𝑀 b = . ± .
33 M ⊕ , while for the 10 . 𝐾 c = . ± .
33 m s − , corresponding to 𝑀 c = . ± .
98 M ⊕ . We point out that the here reported value of 𝐾 b and 𝑀 b is obtained from the joint photometric and RV fit. However, thefinal value of 𝐾 b and 𝑀 b that we decided to adopt (see Section 6.4 formore details) is the weighed mean between the values obtained fromthe joint fit reported in this section and from the floating chunk offsetmethod described in the next section. In addition, we inferred thepresence of two additional planets, with periods of 25 . ± .
04 days(planet d) and 77 . ± .
39 days (planet e), and robustly determinedsemi-amplitudes of 𝐾 d = . ± .
33 m s − ( 𝑀 d = . ± .
28 M ⊕ )and 𝐾 e = . ± .
41 m s − ( 𝑀 e = . ± . ⊕ ). Both planetsshow a single transit in the TESS light curve, previously attributed toa transiting planet with period ∼
16 d, whose presence has howeverbeen ruled out by our analysis. This allowed us to infer a planetaryradius of 𝑅 d = . ± .
13 R ⊕ and 𝑅 e = . ± .
11 R ⊕ for planetd and e respectively.We performed the stability analysis of our determined solution,computing the orbits for 100 Kyr with the whfast integrator (withfixed time-step of 0 . rebound package(Rein & Liu 2012; Rein & Tamayo 2015). During the integrationwe checked the dynamical stability of the solution with the MeanExponential Growth factor of Nearby Orbits (MEGNO or (cid:104) 𝑌 (cid:105) ) indi-cator developed by Cincotta & Simó (2000) and implemented within rebound by Rein & Tamayo (2016). We ran 10 simulations withinitial parameters drawn from a Gaussian distribution centred on thebest-fitting parameters and standard deviation derived in this section.All the 10 runs resulted in a MEGNO value of 2, indicating that thefamily of solutions is stable.Finally, we checked the presence of any additional signal in theRVs residuals after removing the four-planet model contribution. TheGLS periodogram showed a non-significant peak at ∼ . .
20, that is, below the 1% FAP threshold(0 . PyORBIT fit of theRVs, assuming first a four-planet model plus an additional signal, andthen a four-planet model adding a Gaussian Process (GP) regression.For the latter approach, we employed the quasi-periodic kernel asformulated by Grunblatt et al. (2015), with no priors on the GP hyper-parameters, since we could not identify any activity-related signal inthe ancillary datasets (see Section 3.3) . In both cases, the (hyper-)parameters of the additional signal did not reach convergence, whilethe results for the four transiting planets were consistent with thosereported above.Considering these results, we adopt the parameters and configurationdetermined in this section as the representative ones for the TOI-561system, with the only exception of the mass and semi-amplitude ofTOI-561 b, that we discuss in the next section. We are well aware that this is a sub-optimal use of GP regression, andthat this approach may be justified in this specific case only as an attempt toidentify additional signals. MNRAS , 1–20 (2020) G. Lacedelli et al.
Figure 6.
Top : 2-minute cadence flattened light curve of TOI-561. The transits of planet b ( 𝑃 ∼ .
45 d), c ( 𝑃 ∼ . 𝑃 ∼ . 𝑃 ∼ . Bottom : TOI-561 phase-folded 2-minute light curves over the best-fitting models (solidlines) for the four planets. The light curve residuals are shown in the bottom panel.
If the separation between the period of the planet and all the otherperiodic signals is large enough, and the RV signal has a similaror larger semi-amplitude, it is possible to determine the RV semi-amplitude for an USP planet without any assumptions about thenumber of planets in the system or the activity of the host star. Undersuch conditions, during a single night, the influence of any othersignal is much smaller than the measurement error and thus it canbe neglected. If two or more observations are gathered during thesame night and they span a large fraction of the orbital phase, theRV semi-amplitude of the USP planet can be precisely measuredby just applying nightly offsets to remove all the other signals (e.g.Hatzes et al. 2010; Howard et al. 2013; Pepe et al. 2013; Frustagliet al. 2020 for a recent example). Such an approach, also knownas floating chunk offset method (FCO; Hatzes 2014), has provenextremely reliable even in the presence of complex activity signals,as shown by Malavolta et al. (2018). In our case, the shortest, nextperiodic signal (i. e., TOI-561 c at 10 .
78 days) is (cid:39)
24 times theperiod of TOI-561 b (i. e., the USP planet at 0 .
45 days), with similarpredicted RV semi-amplitude, making this target suitable for theFCO approach. Thanks to our observational strategy (see Section2.2) we could use ten different nights for this analysis. Most notably,during two nights we managed to gather six observations spanningnearly 5 hours, i. e., more than 40% of the orbital period of TOI-561 b, at opposite orbital phases, thus providing a good coveragein phase of the RV curve. We did not include RV measurementswith an associated error greater than 2 . − (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 RVjitter as a free parameter to take into account possible short-termstellar variability and any underestimation of the errorbars. From our analysis, we obtained a RV semi-amplitude of 𝐾 p = . ± .
38 m s − , corresponding to a mass of 𝑀 p = . ± .
39 M ⊕ .The resulting RV jitter is 𝑗 < . − (84 . 𝐾 b = . ± .
35 m s − , correspondingto a mass of 𝑀 b = . ± .
36 M ⊕ . Table 5 lists the above-mentionedvalues for TOI-561 b. Our final configuration is quite different from the initial one suggestedby the
TESS automatic pipeline. However, the analyses performed onthe 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 assumingthat the two additional signals seen in the RVs were caused by twonon-transiting planets, the inner one with period of ∼
25 d and theouter one both in the case of ∼
78 d and ∼
180 d period. Sucha model, including the TOI-561.01, .02, .03 candidates plus twoadditional 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. Infact, Table 4 suggests that in this case two additional signals need
MNRAS000
MNRAS000 , 1–20 (2020)
OI-561: an USP super-Earth and three mini-Neptunes Table 5.
Final parameters of the TOI-561 system.Parameter TOI-561b TOI-561c TOI-561d TOI-561e 𝑃 (d) 0 . ± . . ± .
004 25 . ± .
04 77 . ± . 𝑇 a0 (d) 1517 . ± .
001 1527 . ± .
004 1521 . ± .
004 1538 . ± . 𝑎 / 𝑅 ★ . ± .
031 22 . ± .
26 39 . ± .
46 82 . ± . 𝑎 (AU) 0 . ± . . ± . . ± . . + . − . 𝑅 p / 𝑅 ★ . ± . . ± . . ± . . ± . 𝑅 p (R ⊕ ) 1 . ± .
066 2 . ± .
096 2 . ± .
13 2 . ± . 𝑏 . + . − . . + . − . . + . − . . + . − . 𝑖 (deg) 87 . + . − . . + . − . . + . − . . + . − . 𝑇 (hr) 1 . + . − . . + . − . . + . − . . + . − . 𝑒 . + . − . . + . − . . + . − . 𝜔 (deg) 90 (fixed) 200 + − + − ± 𝐾 b (m s − ) 1 . ± .
35 1 . ± .
33 3 . ± .
33 2 . ± . 𝑀 pb (M ⊕ ) 1 . ± .
36 5 . ± .
98 11 . ± .
28 16 . ± . 𝜌 p ( 𝜌 ⊕ ) 0 . ± .
14 0 . ± .
05 0 . ± .
14 0 . ± . 𝜌 p (g cm − ) 3 . ± . . ± . . ± . . ± . Common parameter 𝜌 ★ ( 𝜌 (cid:12) ) 1 . ± . 𝑢 . ± . 𝑢 . ± . 𝜎 cjitter , ph . + . − . 𝜎 djitter (m s − ) 1 . ± . 𝛾 e (m s − ) 79702 . ± . 𝑎 BJD
TDB -2457000. 𝑏 The here reported values of planet b correspond to the weighted mean between the values inferredfrom the floating chunk offset method ( 𝐾 b = . ± .
38 m s − , 𝑀 b = . ± .
39 M ⊕ ) and fromthe joint photometric and RV fit ( 𝐾 b = . ± .
32 m s − , 𝑀 b = . ± .
33 M ⊕ ). 𝑐 Photometric jitter term. 𝑑 Uncorrelated RV jitter term. 𝑒 RV offset. to be added to the three transiting candidates to best reproduce theRV dataset, and therefore the five-planet model should be consideredalso 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 isstrongly favoured with respect to the five-planet model in both cases,with a difference in the logarithmic Bayes factor 2 Δ ln Z (cid:29) , we used the values and standarddeviations derived from the joint photometric and RV fit, except forthe inclination of the two external planets, that we fixed to 90 ◦ .All of 10 runs yielded unstable solutions, with a close encounter oran ejection occurring within the integration time. In order to assessthe origin of the instability of the system, we tested a four-planetconfiguration following the same procedure as above, removing oneplanet each time. We found that the orbital configuration of the systemcould be stable only if we remove the candidate with period of ∼
16 d.Therefore, the stability analysis additionally confirms our determined The mass of the 16 d period planet obtained from the fit was 0 . ± .
03 M ⊕ and 1 . ± .
27 M ⊕ for the ∼
78 d and ∼
180 d external planetperiod, respectively. Obviously, when selecting the 10 samples, the mass wasconstrained to positive values.
Table 6.
Logarithmic Bayesian evidences for the models considered in Sec-tion 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 asingle transit each. Model 1 and 2 correspond the five-planet model, i. e.,including TOI-561.01, .02, .03 and the two additional RV planets (assumedin this case not to transit), in the case of an outer planet at ∼
78 d and ∼
180 dperiod respectively (see Section 6.1). All the values are expressed with re-spect to Model 0. We note that the reported errors, as obtained from the nestedsampling algorithm, are likely underestimated (Nelson et al. 2020). Model 0 Model 1 Model 2ln Z . ± . − . ± . − . ± . four-planet configuration, ruling out the presence of a ∼
16 d periodplanet.
According to our analysis, TOI-561 hosts four transiting planets,including an USP planet, a ∼ . ∼ . ∼ . TESS light curve; those transits were initially associated with acandidate planet with period of ∼
16 d, whose presence we ruledout. As a ‘lesson learned’, we would suggest that caution should be
MNRAS , 1–20 (2020) G. Lacedelli et al.
Figure 7.
Phase-folded RV fit with residuals from the joint four-planet pho-tometric and RV analysis. Planets b, c, d, and e are shown in blue, orange, redand purple, respectively. The reported errorbars include the jitter term, addedin quadrature. taken when candidate planets, detected by photometric pipelines, arebased on just two transits. In such cases, one should not hesitate toconsider alternative scenarios.TOI-561 joins the sample of 88 confirmed systems with 4 or more − − R V ( m s − ) − R e s i du a l s ( m s − ) Figure 8.
Four-planet model from the joint photometric and RV analysis. Thegrey curve is the the best-fitting model, and the blue points are the HARPS-Ndata. The residuals are shown in the bottom panel. The reported errorbarsinclude the jitter term, added in quadrature.
Figure 9.
Phase folded RVs of the ten nights used to model the RV semi-amplitude of the USP planet using the FCO approach. planets , and it is one of the few multi-planet systems with both amass and radius estimate for all the planets. Our global photometricand RV model allowed us to determine the masses and densities ofall the planets with high precision, with a significance of ∼ . 𝜎 for planet b and > 𝜎 for planets c, d and e. In Figure 10 we showthe position of TOI-561 b, c, d and e in the mass-radius diagram ofexoplanets with masses and radiii measured with a precision betterthan 30%. The comparison with the theoretical mass-radius curvesexcludes an Earth-like composition ( ∼
33% iron and 67% silicates)for all planets in the system, whose internal structure we furtheranalyse in the following sections. According to the https://exoplanetarchive.ipac.caltech.edu/.MNRAS000
33% iron and 67% silicates)for all planets in the system, whose internal structure we furtheranalyse in the following sections. According to the https://exoplanetarchive.ipac.caltech.edu/.MNRAS000 , 1–20 (2020)
OI-561: an USP super-Earth and three mini-Neptunes The density ( 𝜌 b = . ± . − ) of the USP planet is consistentwith a 50% (or even more) water composition. Such a composi-tion may be compatible with a water-world scenario, where ‘waterworlds’ are planets with massive water envelopes, in the form ofhigh pressure H O ice, comprising >
5% of the total mass. Evenassuming the higher mass value inferred with the FCO method( 𝑀 b = . ± .
39 M ⊕ , implying a density of 𝜌 b = . ± . − ),TOI-561 b would be located close to the 25% water composition the-oretical curve in the mass-radius diagram, and it would be consistentwith a rocky composition only at a confidence level greater than 2 𝜎 in both radius and mass. Given its proximity to the host star (incidentflux 𝐹 p (cid:39) 𝐹 ⊕ ), the presence of any thick H-He envelope has tobe excluded due the photo-evaporation processes that such old close-in planets are expected to suffer (e.g. Lopez 2017). Nevertheless, thepossibility of a water-world scenario is an intriguing one. An H O-dominated composition would imply that the planet formed beyondthe snow line, accreted a considerable amount of condensed water,and finally migrated inwards (Zeng et al. 2019). While the determi-nation of the precise interior composition of TOI-561 b is beyondthe scope of this work, if such an interpretation is proven trustwor-thy by future observational campaigns, TOI-561 b would supportthe hypothesis that the formation of super-Earths with a significantamount of water is indeed possible. However, an important caveatshould be considered while investigating this scenario. If TOI-561b was a water world, being more irradiated than the runaway green-house irradiation limit, the planet would present a massive and veryextended steam atmosphere. Such an atmosphere would substantiallyincrease the measured radius compared to a condensed water world(Turbet et al. 2020). Therefore, a comparison with the condensedwater-world theoretical curves should be used with caution, since inthis 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 oppositeside of the radius valley, i. e. the gap in the distribution of planetaryradii at ∼ . ⊕ (Fulton et al. 2017), with respect to all theother planets in the system. The origin of the so-called radius valleyis likely due to a transition between rocky and non-rocky planetswith extended H-He envelopes, with several physical mechanismsproposed as explanation, i.e. photoevaporation (Chen & Rogers 2016;Owen & Wu 2017; Lopez & Rice 2018; Jin & Mordasini 2018), core-powered mass loss (Ginzburg et al. 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 islocated above the radius valley and it indeed appears to require a thickH-He envelope (see next section). In the same way, the compositionsof planet d and e are consistent with the presence of a gaseousenvelope. However, the density of TOI-561 b is lower than expectedfor a planet located below the radius valley, where we mainly expectrocky compositions. Moreover, TOI-561 b is the first USP planet withsuch a low measured density (see Figure 10). We note that also theUSP 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 thanTOI-561 b, i. e., 𝜌 W47 e = . ± . − (Vanderburg et al. 2017)and 𝜌 e = . ± . − (Demory et al. 2016) respectively.Vanderburg et al. (2017) proposed the presence of water envelopes asa possible explanation for the low densities of these two planets, eventhough the inferred amount of water was smaller than the one requiredto explain TOI-561 b location in the mass-radius diagram. It shouldalso be considered that both planets are more massive than TOI-561 b,i. e., 𝑀 W47 e = . ± .
66 M ⊕ (Vanderburg et al. 2017) and 𝑀 e = . ± .
31 M ⊕ (Demory et al. 2016), thus increasing their chancesof retaining a small envelope of high-metallicity volatile materials(or water steam) that could explain their low densities (Vanderburget al. 2017). Given its smaller mass, this scenario is less probablefor TOI-561 than for WASP-47 e and 55 Cnc e, making the objecteven more peculiar. With its particular properties, this planet couldbe an intriguing case to test also other extreme planetary compositionmodels. For example, given the metal-poor alpha-enriched host star,the planet is likely to have a lighter core composition. TOI-561 c, with a density of 𝜌 c ∼ . − , is located above thethreshold of a 100% water composition, and given its position inthe mass-radius diagram we suppose the presence of a significantgaseous envelope surrounding an Earth-like iron core and a silicatemantle, and possibly a significant water layer (high-pressure ice). Ifthe inner USP planet is water-rich, there is no simple planet formationscenario in which the outer three planets are water-poor. It is simplerto assume that all four planets were formed with similar volatileabundances, and that the inner USP planet lost all of its H-He layer,plus much of its water content, while the outer planets could keepthem. Following Lopez & Fortney (2014), assuming a rocky Earth-like core and a solar composition H-He envelope, we estimate that anH-He envelope comprising ∼ .
9% of the planet mass could explainthe density of TOI-561 c, using our derived stellar and planetaryparameters.Planets TOI-561 d and e are consistent with a >
50% water com-position, a feature that may place them among the water worlds.However, such densities are also consistent with the presence of arocky core plus water mantel surrounded by a gaseous envelope. Weestimate that a H-He envelope of ∼ .
8% and ∼ .
3% of the planetmass could explain the observed planetary properties.
Our analysis shows that the orbital inclinations of planets c, d and e areall consistent within 1 𝜎 (see Table 5), and that the difference with theinclination of the USP planet is of the order of Δ 𝑖 ∼ . ◦ . Accordingto the analysis of Dai et al. (2018), when the innermost planet has 𝑎 / 𝑅 ★ <
5, the minimum mutual inclination with other planets in thesystem often reaches values up to 5 ◦ -10 ◦ , with larger period ratios( 𝑃 c / 𝑃 b > 𝑃 c / 𝑃 b ∼
24) and the value of 𝑎 b / 𝑅 ★ = .
6, the measured Δ 𝑖 ∼ . ◦ in this case is much lower thatthe expected inclination dispersion of 6 . ± . ◦ that Dai et al. (2018)inferred for systems with similar orbital configurations, indicatingthat the TOI-561 system probably evolved through a mechanism thatdid not excite the inclination of the innermost planet.We also performed a dynamical N-body simulation to check ifsignificant TTVs are expected in the TOI-561 system with our de-termined configuration. In fact, the period ratio of TOI-561 d and eindicates that the planets are close to a 3:1 commensurability, hintof a second order mean motion resonance (MMR), that may suggestthe presence of a strong dynamical interaction between these planets.Starting from the initial configuration (as reported in Table 5), wenumerically integrated the orbits using the N-body integrator ias15 within the rebound package (Rein & Liu 2012). We assumed asreference time the 𝑇 of the USP planet (see Table 5), that roughlycorresponds to the beginning of the TESS observations of TOI-561.During the integration, we computed the transit times of each planet
MNRAS , 1–20 (2020) G. Lacedelli et al. following the procedure described in Borsato et al. (2019), and wecompared the inferred transit times with the linear ephemeris in or-der to obtain the TTV signal, reported as an observed-calculateddiagram ( 𝑂 − 𝐶 , Agol & Fabrycky 2018) in Figure 11. According toour simulation, TOI-561 d and e display an anti-correlated TTV sig-nal, with a very long TTV period of ∼ ∼
13 yr), and TTVamplitudes of ∼
62 minutes (planet d) and ∼
84 minutes (planet e),calculated computing the GLS periodogram of the simulated TTVs.The anti-correlated signal demonstrates that the two planets are ex-pected to dynamically interact (Agol & Fabrycky 2018). In contrast,the predicted TTV amplitude of planet c is extremely low ( ∼ . < 𝑂 − 𝐶 peaks (or dips) in Figure 11. Accordingto our simulation, the first peak (dip) corresponds to the period be-tween March–December 2020, while the second one will be betweenJanuary–October 2026, i. e., corresponding to the time-spans be-tween ∼ ∼ 𝑇 s inferred from single transit observations, thus implying asignificant uncertainty in the TTV phase determination. Therefore,additional photometric observations are necessary to refine the linearephemeris of the planets, and consequently also the prediction of theTTV phase. Given the interesting composition of the planets in the system, wechecked if the TOI-561 planets would be accessible targets for at-mospheric characterisation through transmission spectroscopy, e.g.with the
James Webb Space Telescope ( JWST ). For all the planetsin the system, we calculated the Transmission Spectroscopy Metric(TSM, Kempton et al. 2018), which predicts the expected trans-mission spectroscopy SNR of a 10-hour observing campaign with
JWST /Near Infrared Imager and Slitless Spectrograph (NIRISS) un-der the assumptions of cloud-free atmospheres, the same atmosphericcomposition for all planets of a given type, and a fixed mass-radiusrelation. We obtained TSM values of 19, 107, 24, and 14 for planetsb, c, d, and e, respectively. According to Kempton et al. (2018) ,this classifies TOI-561 b and c as high-quality atmospheric charac-terisation targets among the TESS planetary candidates. However, itshould be noted that the TSM metric assumes rocky composition forplanets with radius < . ⊕ and according to our analysis TOI-561 b is not compatible with such a composition. The same caveatholds for planet c, for which the assumptions under which the TSMis calculated may not be totally valid (e.g. the mass obtained fromour analysis is not the same as if calculated with the Chen & Kipping(2017) mass-radius relation, that is the relation assumed in Kemp-ton et al. (2018), and that would imply a mass of 𝑀 c (cid:39) . ⊕ ).Therefore, this estimate of the atmospheric characterisation feasibil-ity should be used with caution, especially as the TSM metric has The authors suggest to select planets with TSM >
12 for 𝑅 p < . ⊕ ,TSM >
92 for 1 . ⊕ < 𝑅 p < .
75 R ⊕ , and TSM >
84 for 2 .
75 R ⊕ <𝑅 p < ⊕ . ⊕ )1 . . . . . . R a d i u s ( R ⊕ ) % F e % F e % F e % F e R o c k y % H O % H O % H O Radius valley + % H TOI-561 cTOI-561 b TOI-561 dTOI-561 e
EarthVenus F p /F ⊕ Figure 10.
Mass-radius diagram for known exoplanets with mass and ra-dius measurements more precise than 30%, colour-coded according to theirincidental flux in Earth units. The TOI-561 planets are labelled and rep-resented with coloured diamonds. The USP planets are highlighted withblack 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 con-tent predicted by collisional stripping (Marcus et al. 2010). The planetarydata are taken from the The Extrasolar Planets Encyclopaedia catalogue( http://exoplanet.eu/catalog/ ) updated to August 17, 2020. been conceived to prioritise targets for follow-up, and not to preciselydetermine the atmospheric transmission properties.
According to our analysis, TOI-561 hosts a nearly co-planar four-planet system, with an unusually low density USP super-Earth (planetb), a mini-Neptune (planet c) with a significant amount of volatilessurrounding a rocky core, and two mini-Neptunes, which are bothconsistent with a water-world scenario or with a rocky core sur-rounded by a gaseous envelope, and that are expected to show astrong, long-term TTV signal. The multi-planetary nature of TOI-561offers a unique opportunity for comparative exoplanetology. TOI-561planets may be compared with the known population of multi-planetsystems to understand their underlying distribution and occurrences,and to give insights on the formation and evolution processes ofclose-in planets, especially considering the intriguing architecture ofthe system, with the presence of a uncommonly low-density USPsuper-Earth and three mini-Neptunes on the opposite side of theradius valley.Considering the few available data (i. e., 2 transits for planet c,1 transit for planets d, e), additional observations are needed to un-equivocally confirm our solution. Further high-precision photometric(i.e. with
TESS , that will re-observe TOI-561 in sector 35 – Febru-ary/March 2021, or with the
CHEOPS satellite) and RVs observa-tions will help improving the precision on the planets parameters,both allowing for the detection of eventual TTVs and increasing thetime-span of the RV dataset, that could also unveil possible additionallong-period companions.
MNRAS000
MNRAS000 , 1–20 (2020)
OI-561: an USP super-Earth and three mini-Neptunes − − − − O - C ( m i n ) TOI-561dTOI-561e
Figure 11.
Predicted TTV signal of TOI-561 d and e assuming our best-fittingmodel (see Table 5). The planets show a strong, anti-correlated signal. Thesignals of the USP planet ( < < ACKNOWLEDGEMENTS
We thank the anonymous referee for the constructive comments andrecommendations which helped improving the paper.This paper includes data collected by the
TESS mission, which arepublicly available from the Mikulski Archive for Space Telescopes(MAST). Funding for the
TESS mission is provided by the NASAExplorer Program. Resources supporting this work were provided bythe NASA High-End Computing (HEC) Program through the NASAAdvanced 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 Galileiof the INAF at the Spanish Observatorio del Roque de los Mucha-chos of the Instituto de Astrofisica de Canarias (GTO program, andA40TAC_23 program from INAF-TAC). The HARPS-N project wasfunded 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 ItalianNational Astrophysical Institute (INAF), University of St. Andrews,Queen’s University Belfast and University of Edinburgh. Parts of thiswork have been supported by the National Aeronautics and Space Ad-ministration under grant No. NNX17AB59G issued through the Exo-planets Research Program. This research has made use of the NASAExoplanet Archive, which is operated by the California Institute ofTechnology, under contract with the National Aeronautics and SpaceAdministration under the Exoplanet Exploration Program. This workhas made use of data from the European Space Agency (ESA) mis-sion
Gaia ( ), processed bythe Gaia
Data Processing and Analysis Consortium (DPAC, ). Fund-ing for the DPAC has been provided by national institutions, in par-ticular the institutions participating in the
Gaia
Multilateral Agree-ment. This publication makes use of data products from the TwoMicron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Cen-ter/California Institute of Technology, funded by the National Aero-nautics and Space Administration and the National Science Founda-tion. This work is made possible by a grant from the John TempletonFoundation. The opinions expressed in this publication are those ofthe authors and do not necessarily reflect the views of the John Tem-pleton Foundation.GL acknowledges support by CARIPARO Foundation, according tothe agreement CARIPARO-Università degli Studi di Padova (Pratican. 2018/0098), and scholarship support by the “Soroptimist Interna-tional d’Italia” association (Cortina d’Ampezzo Club). GLa, LBo,GPi, VN, GS, and IPa acknowledge the funding support from ItalianSpace Agency (ASI) regulated by “Accordo ASI-INAF n. 2013-016-R.0 del 9 luglio 2013 e integrazione del 9 luglio 2015 CHEOPSFasi A/B/C”. DNa acknowledges the support from the French Cen-tre National d’Etudes Spatiales (CNES). AM acknowledges supportfrom the senior Kavli Institute Fellowships. ACC acknowledges sup-port from STFC consolidated grant ST/R000824/1 and UK SpaceAgency grant ST/R003203/1. ASB and MPi acknowledge financialcontribution from the ASI-INAF agreement n.2018-16-HH.0. XD isgrateful to the Branco-Weiss Fellowship for continuous support. Thisproject has received funding from the European Research Council(ERC) under the European Union’s Horizon 2020 research and in-novation program (grant agreement No. 851555). JNW thanks theHeising-Simons Foundation for support.
DATA AVAILABILITY
HARPS-N observations and data products are available throughthe Data & Analysis Center for Exoplanets (DACE) at https://dace.unige.ch/ . TESS data products can be accessed throughthe official NASA website https://heasarc.gsfc.nasa.gov/docs/tess/data-access.html .All underlying data are available either in the appendix/online sup-porting material or will be available via VizieR at CDS.
REFERENCES
Adibekyan V. Z., Sousa S. G., Santos N. C., Delgado Mena E., GonzálezHernández J. I., Israelian G., Mayor M., Khachatryan G., 2012, A&A,545, A32Agol E., Fabrycky D. C., 2018, Handbook of Exoplanets, pp 797–816Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481Bailer-Jones C. A. L., Rybizki J., Fouesneau M., Mantelet G., Andrae R.,2018, AJ, 156, 58Baranne A., et al., 1996, A&AS, 119, 373Baruteau C., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P.,Henning T., eds, Protostars and Planets VI. p. 667 ( arXiv:1312.4293 ),doi:10.2458/azu_uapress_9780816531240-ch029Baruteau C., Bai X., Mordasini C., Mollière P., 2016, Space Sci. Rev., 205,77Bensby T., Feltzing S., Oey M. S., 2014, A&A, 562, A71Borsato L., et al., 2019, MNRAS, 484, 3233Borucki W. J., et al., 2010, Science, 327, 977Buchhave L. A., et al., 2012, Nature, 486, 375Buchhave L. A., et al., 2014, Nature, 509, 593Buchner, J. et al., 2014, A&A, 564, A125Castelli F., Kurucz R. L., 2003, in Piskunov N., Weiss W. W., Gray D. F., eds,IAU Symposium Vol. 210, Modelling of Stellar Atmospheres. p. A20( arXiv:astro-ph/0405087 )Chen J., Kipping D., 2017, ApJ, 834, 17Chen H., Rogers L. A., 2016, The Astrophysical Journal, 831, 180MNRAS , 1–20 (2020) G. Lacedelli et al.
Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016,ApJ, 823, 102Cincotta P. M., Simó C., 2000, A&AS, 147, 205Claret A., 2018, A&A, 618, A20Cloutier R., et al., 2019, A&A, 621, A49Collier Cameron A., et al., 2019, MNRAS, 487, 1082Cosentino R., et al., 2012, in Society of Photo-Optical Instrumentation Engi-neers (SPIE) Conference Series. p. 1, doi:10.1117/12.925738Cosentino R., et al., 2014, in Society of Photo-Optical Instrumentation Engi-neers (SPIE) Conference Series. p. 8, doi:10.1117/12.2055813Cutri R. M., et al., 2003, VizieR Online Data Catalog, 2246Dai F., Masuda K., Winn J. N., 2018, ApJ, 864, L38Davies M. B., Adams F. C., Armitage P., Chambers J., Ford E., MorbidelliA., Raymond S. N., Veras D., 2014, in Beuther H., Klessen R. S.,Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 787( arXiv:1311.6816 ), doi:10.2458/azu_uapress_9780816531240-ch034Demory B.-O., Gillon M., Madhusudhan N., Queloz D., 2016, MNRAS, 455,2018Dotter A., 2016, ApJS, 222, 8Dotter A., Chaboyer B., Jevremović D., Kostov V., Baron E., Ferguson J. W.,2008, ApJS, 178, 89Douglas S. T., Curtis J. L., Agüeros M. A., Cargile P. A., Brewer J. M.,Meibom S., Jansen T., 2019, ApJ, 879, 100Dragomir D., et al., 2019, The Astrophysical Journal, 875, L7Dumusque X., et al., 2019, A&A, 627, A43Eastman J. D., et al., 2019, EXOFASTv2: A public, generalized, publication-quality exoplanet modeling code ( arXiv:1907.09480 )Fabrycky D. C., et al., 2014, The Astrophysical Journal, 790, 146Feroz F., Hobson M. P., 2008, MNRAS, 384, 449Feroz F., Hobson M. P., Bridges M., 2009, MNRAS, 398, 1601Feroz F., Hobson M. P., Cameron E., Pettitt A. N., 2019, The Open Journalof Astrophysics, 2, 10Figueira P., Santos N. C., Pepe F., Lovis C., Nardetto N., 2013, A&A, 557,A93Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125,306Fossati L., et al., 2017, A&A, 601, A104Fressin F., et al., 2013, The Astrophysical Journal, 766, 81Frustagli G., et al., 2020, A&A, 633, A133Fulton B. J., et al., 2017, AJ, 154, 109Gaia Collaboration et al., 2018, A&A, 616, A1Gelman A., Rubin D. B., 1992, Statistical Science, 7, 16Ginzburg S., Schlichting H. E., Sari R., 2018, Monthly Notices of the RoyalAstronomical Society, 476, 759Gomes da Silva J., Santos N. C., Bonfils X., Delfosse X., Forveille T., UdryS., 2011, A&A, 534, A30Goodman J., Weare J., 2010, Commun. Appl. Math. Comput. Sci., 5, 65Grunblatt S. K., Howard A. W., Haywood R. D., 2015, ApJ, 808, 127Günther M. N., Pozuelos F. J., Dittmann J. A., et al., 2019, Nature Astronomy,3, 1099Gupta A., Schlichting H. E., 2019, Monthly Notices of the Royal AstronomicalSociety, 487, 24Hatzes A. P., 2014, A&A, 568, A84Hatzes A. P., et al., 2010, A&A, 520, A93Helled R., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P.,Henning T., eds, Protostars and Planets VI. p. 643 ( arXiv:1311.1142 ),doi:10.2458/azu_uapress_9780816531240-ch028Hippke M., Heller R., 2019, A&A, 623, A39Hippke M., David T. J., Mulders G. D., Heller R., 2019, AJ, 158, 143Hodges J. L., 1958, Arkiv för Matematik, 3, 469Howard A. W., et al., 2013, Nature, 503, 381Isaacson H., Fischer D., 2010, ApJ, 725, 875Jenkins J. M., 2020, Kepler Data Processing Handbook, Kepler Sci-ence Document KSCI-19081-003, https://archive.stsci.edu/kepler/documents.html
Jenkins J. M., et al., 2016, in Chiozzi G., Guzman J. C., eds, Society ofPhoto-Optical Instrumentation Engineers (SPIE) Conference Series Vol.9913, Software and Cyberinfrastructure for Astronomy IV. SPIE, pp 1232 – 1251, doi:10.1117/12.2233418, https://doi.org/10.1117/12.2233418
Jin S., Mordasini C., 2018, The Astrophysical Journal, 853, 163Kass R. E., Raftery A. E., 1995, Journal of the American Statistical Associ-ation, 90, 773Kempton E. M. R., et al., 2018, PASP, 130, 114401Khan S., et al., 2019, A&A, 628, A35Kipping D. M., 2010, MNRAS, 408, 1758Kipping D. M., 2013, MNRAS, 435, 2152Kochanek C. S., et al., 2017, PASP, 129, 104502Kreidberg L., 2015, PASP, 127, 1161Lanza A. F., et al., 2018, A&A, 616, A155Latham D. W., et al., 2011, ApJ, 732, L24Lee E. J., Chiang E., 2016, The Astrophysical Journal, 817, 90Lightkurve Collaboration et al., 2018, Lightkurve: Kepler and TESStime series analysis in Python, Astrophysics Source Code Library(ascl:1812.013)Lissauer J. J., et al., 2011, The Astrophysical Journal Supplement Series, 197,8Lissauer J. J., et al., 2012, ApJ, 750, 112Lopez E. D., 2017, MNRAS, 472, 245Lopez E. D., Fortney J. J., 2014, The Astrophysical Journal, 792, 1Lopez E. D., Rice K., 2018, Monthly Notices of the Royal AstronomicalSociety, 479, 5303Lovis C., et al., 2011, preprint, ( arXiv:1107.5325 )Malavolta L., et al., 2016, A&A, 588, A118Malavolta L., et al., 2017a, AJ, 153, 224Malavolta L., Lovis C., Pepe F., Sneden C., Udry S., 2017b, MNRAS, 469,3965Malavolta L., et al., 2018, AJ, 155, 107Mamajek E. E., Hillenbrand L. A., 2008, ApJ, 687, 1264Marcus R. A., Sasselov D., Stewart S. T., Hernquist L., 2010, ApJ, 719, L45McQuillan A., Aigrain S., Mazeh T., 2013, MNRAS, 432, 1203Milbourne T. W., et al., 2019, ApJ, 874, 107Millholland S., Wang S., Laughlin G., 2017, ApJ, 849, L33Morbidelli A., Lunine J., O’Brien D., Raymond S., Walsh K., 2012, AnnualReview of Earth and Planetary Sciences, 40, 251Mortier A., Santos N. C., Sousa S. G., Fernandes J. M., Adibekyan V. Z.,Delgado Mena E., Montalto M., Israelian G., 2013, A&A, 558, A106Mortier A., Sousa S. G., Adibekyan V. Z., Brandão I. M., Santos N. C., 2014,A&A, 572, A95Morton T. D., 2015, isochrones: Stellar model grid package (ascl:1503.010)Nardetto N., Mourard D., Kervella P., Mathias P., Mérand A., Bersier D.,2006, A&A, 453, 309Nardiello D., et al., 2019, MNRAS, 490, 3806Nardiello D., et al., 2020, MNRAS, 495, 4924Nelson B. E., et al., 2020, AJ, 159, 73Noyes R. W., Hartmann L. W., Baliunas S. L., Duncan D. K., Vaughan A. H.,1984, ApJ, 279, 763Owen J. E., Wu Y., 2017, ApJ, 847, 29Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011,ApJS, 192, 3Pecaut M. J., Mamajek E. E., 2013, ApJS, 208, 9Pepe F., Mayor M., Galland F., Naef D., Queloz D., Santos N. C., Udry S.,Burnet M., 2002, A&A, 388, 632Pepe F., et al., 2013, Nature, 503, 377Planck Collaboration et al., 2018, arXiv e-prints, p. arXiv:1807.06205Raymond S. N., Kokubo E., Morbidelli A., Morishima R., Walsh K. J.,2014, in Beuther H., Klessen R. S., Dullemond C. P., HenningT., eds, Protostars and Planets VI. p. 595 ( arXiv:1312.1689 ),doi:10.2458/azu_uapress_9780816531240-ch026Rein H., Liu S. F., 2012, A&A, 537, A128Rein H., Tamayo D., 2015, MNRAS, 452, 376Rein H., Tamayo D., 2016, MNRAS, 459, 2275Ricker G. R., et al., 2014, Journal of Astronomical Telescopes, Instruments,and Systems, 1, 1Roy A., et al., 2020, AJ, 159, 161Seager S., Mallen-Ornelas G., 2003, The Astrophysical Journal, 585, 1038MNRAS000
Jin S., Mordasini C., 2018, The Astrophysical Journal, 853, 163Kass R. E., Raftery A. E., 1995, Journal of the American Statistical Associ-ation, 90, 773Kempton E. M. R., et al., 2018, PASP, 130, 114401Khan S., et al., 2019, A&A, 628, A35Kipping D. M., 2010, MNRAS, 408, 1758Kipping D. M., 2013, MNRAS, 435, 2152Kochanek C. S., et al., 2017, PASP, 129, 104502Kreidberg L., 2015, PASP, 127, 1161Lanza A. F., et al., 2018, A&A, 616, A155Latham D. W., et al., 2011, ApJ, 732, L24Lee E. J., Chiang E., 2016, The Astrophysical Journal, 817, 90Lightkurve Collaboration et al., 2018, Lightkurve: Kepler and TESStime series analysis in Python, Astrophysics Source Code Library(ascl:1812.013)Lissauer J. J., et al., 2011, The Astrophysical Journal Supplement Series, 197,8Lissauer J. J., et al., 2012, ApJ, 750, 112Lopez E. D., 2017, MNRAS, 472, 245Lopez E. D., Fortney J. J., 2014, The Astrophysical Journal, 792, 1Lopez E. D., Rice K., 2018, Monthly Notices of the Royal AstronomicalSociety, 479, 5303Lovis C., et al., 2011, preprint, ( arXiv:1107.5325 )Malavolta L., et al., 2016, A&A, 588, A118Malavolta L., et al., 2017a, AJ, 153, 224Malavolta L., Lovis C., Pepe F., Sneden C., Udry S., 2017b, MNRAS, 469,3965Malavolta L., et al., 2018, AJ, 155, 107Mamajek E. E., Hillenbrand L. A., 2008, ApJ, 687, 1264Marcus R. A., Sasselov D., Stewart S. T., Hernquist L., 2010, ApJ, 719, L45McQuillan A., Aigrain S., Mazeh T., 2013, MNRAS, 432, 1203Milbourne T. W., et al., 2019, ApJ, 874, 107Millholland S., Wang S., Laughlin G., 2017, ApJ, 849, L33Morbidelli A., Lunine J., O’Brien D., Raymond S., Walsh K., 2012, AnnualReview of Earth and Planetary Sciences, 40, 251Mortier A., Santos N. C., Sousa S. G., Fernandes J. M., Adibekyan V. Z.,Delgado Mena E., Montalto M., Israelian G., 2013, A&A, 558, A106Mortier A., Sousa S. G., Adibekyan V. Z., Brandão I. M., Santos N. C., 2014,A&A, 572, A95Morton T. D., 2015, isochrones: Stellar model grid package (ascl:1503.010)Nardetto N., Mourard D., Kervella P., Mathias P., Mérand A., Bersier D.,2006, A&A, 453, 309Nardiello D., et al., 2019, MNRAS, 490, 3806Nardiello D., et al., 2020, MNRAS, 495, 4924Nelson B. E., et al., 2020, AJ, 159, 73Noyes R. W., Hartmann L. W., Baliunas S. L., Duncan D. K., Vaughan A. H.,1984, ApJ, 279, 763Owen J. E., Wu Y., 2017, ApJ, 847, 29Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., Timmes F., 2011,ApJS, 192, 3Pecaut M. J., Mamajek E. E., 2013, ApJS, 208, 9Pepe F., Mayor M., Galland F., Naef D., Queloz D., Santos N. C., Udry S.,Burnet M., 2002, A&A, 388, 632Pepe F., et al., 2013, Nature, 503, 377Planck Collaboration et al., 2018, arXiv e-prints, p. arXiv:1807.06205Raymond S. N., Kokubo E., Morbidelli A., Morishima R., Walsh K. J.,2014, in Beuther H., Klessen R. S., Dullemond C. P., HenningT., eds, Protostars and Planets VI. p. 595 ( arXiv:1312.1689 ),doi:10.2458/azu_uapress_9780816531240-ch026Rein H., Liu S. F., 2012, A&A, 537, A128Rein H., Tamayo D., 2015, MNRAS, 452, 376Rein H., Tamayo D., 2016, MNRAS, 459, 2275Ricker G. R., et al., 2014, Journal of Astronomical Telescopes, Instruments,and Systems, 1, 1Roy A., et al., 2020, AJ, 159, 161Seager S., Mallen-Ornelas G., 2003, The Astrophysical Journal, 585, 1038MNRAS000 , 1–20 (2020)
OI-561: an USP super-Earth and three mini-Neptunes Shappee B. J., et al., 2014, ApJ, 788, 48Skilling J., 2004, in Fischer R., Preuss R., Toussaint U. V., eds, AmericanInstitute of Physics Conference Series Vol. 735, Bayesian Inference andMaximum Entropy Methods in Science and Engineering: 24th Interna-tional Workshop on Bayesian Inference and Maximum Entropy Methodsin Science and Engineering. pp 395–405, doi:10.1063/1.1835238Skilling J., 2006, Bayesian Anal., 1, 833Smith J. C., et al., 2012, PASP, 124, 1000Sneden C., 1973, ApJ, 184, 839Sousa S. G., 2014, ARES + MOOG: A Practical Overview of an EquivalentWidth (EW) Method to Derive Stellar Parameters. Springer, Cham, pp297–310, doi:10.1007/978-3-319-06956-2_26Sousa S. G., Santos N. C., Israelian G., Lovis C., Mayor M., Silva P. B., UdryS., 2011, A&A, 526, A99Sousa S. G., Santos N. C., Adibekyan V., Delgado-Mena E., Israelian G.,2015, A&A, 577, A67Speagle J. S., 2020, MNRAS, 493, 3132Stassun K. G., et al., 2018, The Astronomical Journal, 156, 102Stumpe M. C., et al., 2012, PASP, 124, 985Stumpe M. C., Smith J. C., Catanzarite J. H., Van Cleve J. E., Jenkins J. M.,Twicken J. D., Girouard F. R., 2014, PASP, 126, 100Turbet M., Bolmont E., Ehrenreich D., Gratier P., Leconte J., Selsis F., HaraN., Lovis C., 2020, A&A, 638, A41Van Eylen V., et al., 2019, AJ, 157, 61Vanderburg A., et al., 2017, AJ, 154, 237Weiss L. M., et al., 2018, AJ, 156, 254Winn J. N., 2010, arXiv e-prints, p. arXiv:1001.2010Wright E. L., et al., 2010, AJ, 140, 1868Zechmeister M., Kürster M., 2009, A&A, 496, 577Zeng L., et al., 2019, Proceedings of the National Academy of Science, 116,9723Ziegler C., Tokovinin A., Briceño C., Mang J., Law N., Mann A. W., 2020,AJ, 159, 19
APPENDIX A: PHOTOMETRIC ANALYSIS
We performed a preliminary light curve fit of the three candidateplanets found by the SPOC pipeline and our independent TLS anal-ysis, that is TOI-561.01, .02, and .03 with periods of about 10 . .
45 d, and 16 . PyORBIT ,as specified in Section 5, but assuming circular orbits for all the can-didate planets, given the uncertainty associated with the eccentricityfrom the analysis of
TESS data alone (Winn 2010). We ran the chainsfor 100 000 steps, discarding the first 20 000 as burn-in. We list theobtained parameters in Table A1 and we show the best-fitting transitmodels in Figure A1. In order to test whether our light curve flatten-ing affected the inferred parameters of the planetary candidates, wealso ran the
PyORBIT fit on the original PDCSAP light curve. For allthe candidates, the difference between the parameters of the two runswas lower than the error on the parameters themselves, indicatingthat the flattening did not significantly alter the results.We stress that, at last, our global analysis disclaimed the presenceof the planetary candidate TOI-561.03, linking the transits here asso-ciated with this candidate to single transits of two additional planetsdiscovered in the system (see Section 6).
APPENDIX B: RV ANALYSISB1 Removal of anomalous points
Before proceeding with a detailed analysis, we verified if any anoma-lous RV measurement was affecting our analysis. We followed asimilar approach to that of Cloutier et al. (2019), but slightly moresophisticated due to the presence of (possibly up to) five planetary signals. Instead of analysing the power variation of the periodogram’speaks associated with the candidate planets while removing one pointat the time, we decided to perform a full RV fit with the methodol-ogy described in Section 5, and to compare the resulting RV semi-amplitudes with those derived using the full dataset. To reduce com-putational time, we decided to remove from the dataset 5 consecutiveobservations at once (i. e., performing 17 iterations rather than 82),and then performed the leave-one-out cross-validation on those sub-sets showing deviating RV semi-amplitudes in order to identify theanomalous RV measurement. With this approach, we found out thata total of 5 RV measurements, with associated errors greater than2 . − and S/N <
35 were systematically producing a decrease inthe semi-amplitude of candidates .01 and .02 by ≈ . − . − ,and we therefore removed these points from our dataset in order toimprove the accuracy of our results, even if the total variation in RVsemi-amplitude was within the error bars. We note that these obser-vations are clearly outliers at more than 2 𝜎 in both the S/N of thespectra and the RV error distributions (see Section 2.2), which is sim-ply the consequence of having been gathered in sub-optimal weatherconditions. A much simpler sigma-clipping selection would have ledto the exclusion of the same data points. The complex approach weemployed in this work can thus be avoided in future analysis involvingHARPS-N data. B2 RV modelling and injection/retrival tests
Given the results of the frequency analysis in Section 6.1, we per-formed a
PyDE + emcee RV fit with
PyORBIT , following the method-ology as described in Section 5, and assuming the model suggestedby 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 theouter signal to be shorter or longer than 100 days, in order to disen-tangle the 78 periodicity from its alias at 180 respectively. We ranthe 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 outer-most planet, TOI-561.03 (i. e., the candidate with period of ∼ . 𝐾 (cid:46) . − , corre-sponding to a rather nonphysical mass of (cid:46) ⊕ (at 1 𝜎 ) for a planetwith 𝑅 p (cid:39) . ⊕ . We thus performed a series of injection/retrievalsimulations in order to assess the influence of the observational sam-pling and of the precision in the mass measurements of the otherplanets. In a first run, the synthetic datasets were simulated by as-suming the orbital parameters as previously determined in the RV fitsfor the candidate planets .01, .02, and the non-transiting candidates,while the RV semi-amplitude of the candidate planet at 16 d wasvaried between 0 . − and 1 . − in steps of 0 . − . Forcomputational reasons, we performed this analysis only with the 78-dsolution for the outer planet. We projected the model onto the realepochs of observation and then we added a Gaussian noise corre-sponding to the measured error plus an RV jitter of 1 . − addedin quadrature, while preserving the original value in the analysis. Webuilt 50 different noise realisations and analysed each of them withthe same methodology as before, i. e., PyDE + emcee through PyOR-BIT , but for a shorter chain length to reduce computing time. Theposteriors of each parameter were then obtained by putting togetherthe individual posterior distributions from each noise realisation. Wefinally repeated the same analysis but varying the RV semi-amplitude
10 000 steps after convergence, reached at approximately 15 000 steps.MNRAS , 1–20 (2020) G. Lacedelli et al.
Figure A1.
Top : 2-minute cadence flattened light curve of TOI-561. The transits of candidates TOI-561.02 ( 𝑃 ∼ .
45 d), .01 ( 𝑃 ∼ . 𝑃 ∼ . Bottom : TOI-561 phase-folded light curves over the best-fitting models (solid lines) for thethree planets. The grey points are the
TESS
Table A1.
Planetary parameters of the three transiting candidates from the initial light curve fitting.Parameter TOI-561.02 TOI-561.01 TOI-561.03 𝑃 (d) 0 . ± . . ± .
005 16 . + . − . 𝑇 a0 (d) 1517 . ± . . ± .
004 1521 . + . − . 𝑎 / 𝑅 ★ . ± .
030 21 . ± .
25 28 . ± . 𝑎 (AU) 0 . ± . . ± . . ± . 𝑅 p / 𝑅 ★ . ± . . ± . . ± . 𝑅 p (R ⊕ ) 1 . ± .
06 2 . ± .
10 2 . ± . 𝑏 . + . − . . ± .
12 0 . + . − . 𝑖 (deg) 86 . + . − . . + . − . . + . − . 𝑇 b14 (hr) 1 . + . − . . + . − . . ± . Common parameter 𝜌 ★ ( 𝜌 (cid:12) ) 1 . ± . 𝑢 . ± . 𝑢 . ± . 𝑎 BJD
TDB -2457000. 𝑏 Transit duration is derived from the posterior distributions using the formulasin Seager & Mallen-Ornelas (2003).MNRAS000
TDB -2457000. 𝑏 Transit duration is derived from the posterior distributions using the formulasin Seager & Mallen-Ornelas (2003).MNRAS000 , 1–20 (2020)
OI-561: an USP super-Earth and three mini-Neptunes Table B1.
Best-fitting parameters from the five-planet RV fit, assuming period boundaries of 2-100 days for the outermost planet.Parameter TOI-561.02 TOI-561.01 TOI-561.03 TOI-561.04 TOI-561.05 𝑃 (d) 0 . ± . . ± .
004 16 . ± .
008 25 . + . − . . ± . 𝑇 a0 (d) 1517 . ± . . ± .
003 1521 . ± .
004 1521 + − + − 𝑒 . + . − . . + . − . . + . − . . + . − . 𝜔 (deg) 90 (fixed) 178 ±
75 235 + − + − + − 𝐾 (m s − ) 1 . ± .
33 1 . ± . < .
37 3 . ± .
36 2 . ± . 𝑀 p (M ⊕ ) 1 . ± .
33 5 . ± . < .
27 12 . ± . . ± . Common parameter 𝜎 bjitter (m s − ) 1 . ± . 𝛾 c (m s − ) 79702 . ± . 𝑎 BJD
TDB -2457000. 𝑏 Uncorrelated jitter term. 𝑐 RV offset.
Table B2.
Best-fitting parameters from the five-planet RV fit, assuming period boundaries of 100-200 days for the outermost planet.Parameter TOI-561.02 TOI-561.01 TOI-561.03 TOI-561.04 TOI-561.05 𝑃 (d) 0 . ± . . ± .
004 16 . ± .
007 25 . ± .
19 179 . + . − . 𝑇 a0 (d) 1517 . ± . . ± .
003 1521 . ± .
004 1518 ± + − 𝑒 . + . − . . + . − . . + . − . . + . − . 𝜔 (deg) 90 (fixed) 148 + − + − + − + − 𝐾 (m s − ) 1 . ± .
32 0 . + . − . < .
54 3 . ± .
36 3 . ± . 𝑀 p (M ⊕ ) 1 . ± .
33 2 . + . − . < .
91 12 . ± . . ± . Common parameter 𝜎 bjitter (m s − ) 1 . ± . 𝛾 c (m s − ) 79703 . ± . 𝑎 BJD
TDB -2457000. 𝑏 Uncorrelated jitter term. 𝑐 RV offset. of the candidate planet .01, i. e., the closest signal in frequency spaceand the one with the most uncertain RV semi-amplitude measure-ment other than the USP candidate, by ± . − with respect tothe value of 1 . − used in the previous analysis. AFFILIATIONS Department of Physics and Astronomy, Università degli Studi diPadova, Vicolo dell’Osservatorio 3, I-35122, Padova, Italy INAF - Osservatorio Astronomico di Padova, Vicolodell’Osservatorio 5, I-35122, Padova, Italy Aix-Marseille Université, CNRS, CNES, LAM, F-13013 Marseille,France Astrophysics Group, Cavendish Laboratory, University of Cam-bridge, J.J. Thomson Avenue, Cambridge CB3 0HE, UK Kavli Institute for Cosmology, University of Cambridge, MadingleyRoad, Cambridge CB3 0HA, UK Observatoire Astronomique de l’Université de Genève, 51 Chemindes Maillettes, 1290 Versoix, Switzerland Centre for Exoplanet Science, SUPA, School of Physics andAstronomy, University of St Andrews, St Andrews KY16 9SS, UK Fundación Galileo Galilei-INAF, Rambla J. A. F. Perez, 7, E-38712,S.C. Tenerife, Spain INAF-Osservatorio Astronomico di Brera, via E. Bianchi 46, 23807 Merate (LC), Italy DTU Space, National Space Institute, Technical University ofDenmark, Elektrovej 328, DK-2800 Kgs. Lyngby, Denmark Center for Astrophysics | Harvard & Smithsonian, 60 GardenStreet, Cambridge, MA, 02138, USA INAF - Osservatorio Astrofisico di Torino, Via Osservatorio 20,I-10025, Pino Torinese, Italy NASA Ames Research Center, Moffett Field, CA, 94035, USA SUPA, Institute for Astronomy, University of Edinburgh, BlackfordHill, Edinburgh, EH9 3HJ, Scotland, UK Centre for Exoplanet Science, University of Edinburgh, Edin-burgh, EH93FD, UK Department of Astrophysical Sciences, Princeton University, 4 IvyLane, Princeton, NJ 08544, USA INAF - Osservatorio Astronomico di Roma, Via Frascati 33,00078, Monte Porzio Catone, Italy INAF - Osservatorio Astronomico di Cagliari, via della Scienza5, 09047, Selargius, Italy INAF - Osservatorio Astronomico di Palermo, Piazza delParlamento 1, I-90134 Palermo, Italy INAF - Osservatorio Astrofisico di Catania, Via S. Sofia 78,I-95123, Catania, Italy Department of Earth, Atmospheric and Planetary Sciences, andKavli Institute for Astrophysics and Space Research, MassachusettsInstitute of Technology, Cambridge, MA 02139, USA
MNRAS , 1–20 (2020) G. Lacedelli et al. Astrophysics Research Centre, School of Mathematics andPhysics, Queen’s University Belfast, Belfast, BT7 1NN, UK
This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000