# A pair of sub-Neptunes transiting the bright K-dwarf TOI-1064 characterized with CHEOPS

Thomas Wilson, Elisa Goffo, Yann Alibert, Davide Gandolfi, Andrea Bonfanti, Carina Persson, Andrew Collier Cameron, Malcolm Fridlund, Luca Fossati, Judith Korth, et al.

► To cite this version:

Thomas Wilson, Elisa Goffo, Yann Alibert, Davide Gandolfi, Andrea Bonfanti, et al. A pair of sub-Neptunes transiting the bright K-dwarf TOI-1064 characterized with CHEOPS. Monthly Notices of the Royal Astronomical Society, 2022, 511, pp.1043-1071. 10.1093/mnras/stab3799. insu-03643027

HAL Id: insu-03643027

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

Submitted on 8 Mar 2023

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

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

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

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

# A pair of sub-Neptunes transiting the bright K-dwarf TOI-1064 characterized with CHEOPS

Thomas G. Wilson<sup>1,★</sup> Elisa Goffo,<sup>2,3</sup> Yann Alibert,<sup>4</sup> Davide Gandolfi<sup>1</sup>,<sup>2</sup> Andrea Bonfanti,<sup>5</sup>

Carina M. Persson,<sup>6</sup> Andrew Collier Cameron<sup>1</sup>,<sup>1</sup> Malcolm Fridlund<sup>1</sup>,<sup>6,7</sup> Luca Fossati,<sup>5</sup> Judith Korth<sup>1,8</sup>

Willy Benz,<sup>4,9</sup> Adrien Deline,<sup>10</sup> Hans-Gustav Florén,<sup>11,12</sup> Pascal Guterman,<sup>13,14</sup> Vardan Adibekyan,<sup>15</sup>

Matthew J. Hooton<sup>1</sup>,<sup>4</sup> Sergio Hoyer<sup>1</sup>,<sup>16</sup> Adrien Leleu,<sup>4,10</sup> Alexander James Mustill<sup>1</sup>,<sup>17</sup>

Sébastien Salmon,<sup>10</sup> Sérgio G. Sousa<sup>1</sup>,<sup>15</sup> Olga Suarez,<sup>18</sup> Lyu Abe,<sup>18</sup> Abdelkrim Agabi,<sup>18</sup> Roi Alonso,<sup>19,20</sup>

Guillem Anglada,<sup>21,22</sup> Joel Asquier,<sup>23</sup> Tamas Bárczy,<sup>24</sup> David Barrado Navascues,<sup>25</sup>

Susana C. C. Barros,<sup>15,26</sup> Wolfgang Baumjohann<sup>1</sup>,<sup>5</sup> Mathias Beck,<sup>10</sup> Thomas Beck,<sup>4</sup> Nicolas Billot,<sup>10</sup>

Xavier Bonfils,<sup>27</sup> Alexis Brandeker,<sup>11</sup> Christopher Broeg,<sup>4,9</sup> Edward M. Bryant<sup>1</sup>,<sup>28</sup>

Matthew R. Burleigh,<sup>29</sup> Marco Buttu,<sup>30</sup> Juan Cabrera,<sup>31</sup> Sébastien Charnoz,<sup>32</sup> David R. Ciardi,<sup>33</sup>

Ryan Cloutier<sup>1</sup>,<sup>34,†</sup> William D. Cochran,<sup>35</sup> Karen A. Collins<sup>1</sup>,<sup>34</sup> Knicole D. Colón,<sup>36</sup> Nicolas Crouzet,<sup>23</sup>

Szilard Csizmadia<sup>1</sup>,<sup>31</sup> Melvyn B. Davies,<sup>37</sup> Magali Deleuil,<sup>13</sup> Laetitia Delrez<sup>1</sup>,<sup>38,39</sup>

Olivier Demangeon,<sup>15,26</sup> Brice-Olivier Demory,<sup>9</sup> Diana Dragomir,<sup>40</sup> Georgina Dransfield,<sup>41</sup>

David Ehrenreich,<sup>10</sup> Anders Erikson,<sup>31</sup> Andrea Fortier,<sup>4,9</sup> Tianjun Gan,<sup>42</sup> Samuel Gill,<sup>28</sup> Michaël Gillon,<sup>38</sup>

Crystal L. Gnilka,<sup>43</sup> Nolan Grieves,<sup>10</sup> Sascha Grziwa,<sup>44</sup> Manuel Güdel,<sup>45</sup> Tristan Guillot,<sup>18</sup>

Jonas Haldemann<sup>1</sup>,<sup>4</sup> Kevin Heng,<sup>9,28</sup> Keith Horne,<sup>1</sup> Steve B. Howell,<sup>43</sup> Kate G. Isaak,<sup>46</sup> Jon M. Jenkins,<sup>43</sup>

Eric L. N. Jensen<sup>1</sup>,<sup>47</sup> Laszlo Kiss,<sup>48</sup> Gaia Lacedelli<sup>1</sup>,<sup>49,50</sup> Kristine Lam,<sup>31</sup> Jacques Laskar,<sup>51</sup>

David W. Latham,<sup>34</sup> Alain Lecavelier des Etangs,<sup>52</sup> Monika Lendl<sup>1</sup>,<sup>10</sup> Kathryn V. Lester,<sup>43</sup>

Alan M. Levine,<sup>53</sup> John Livingston<sup>1</sup>,<sup>54</sup> Christophe Lovis,<sup>10</sup> Rafael Luque<sup>1</sup>,<sup>55</sup> Demetrio Magrin,<sup>50</sup>

Wenceslas Marie-Sainte,<sup>30</sup> Pierre F. L. Maxted<sup>1</sup>,<sup>56</sup> Andrew W. Mayo<sup>1</sup>,<sup>57</sup> Brian McLean,<sup>58</sup>

Marko Mecina,<sup>45</sup> Djamel Mékarnia,<sup>18</sup> Valerio Nascimbeni<sup>1</sup>,<sup>50</sup> Louise D. Nielsen<sup>1</sup>,<sup>59</sup> Göran Olofsson,<sup>11</sup>

Hugh P. Osborn,<sup>9,53</sup> Hannah L. M. Osborne<sup>1</sup>,<sup>60</sup> Roland Ottensamer,<sup>45</sup> Isabella Pagano,<sup>61</sup> Enric Pallé,<sup>19,20</sup>

Gisbert Peter,<sup>62</sup> Giampaolo Piotto,<sup>49,50</sup> Don Pollacco,<sup>28</sup> Didier Queloz,<sup>10,63</sup> Roberto Ragazzoni,<sup>49,50</sup>

Nicola Rando,<sup>23</sup> Heike Rauer,<sup>31,64,65</sup> Seth Redfield<sup>1</sup>,<sup>66</sup> Ignasi Ribas,<sup>21,22</sup> George R. Ricker,<sup>53</sup>

Martin Rieder,<sup>4</sup> Nuno C. Santos,<sup>15,26</sup> Gaetano Scandariato,<sup>61</sup> François-Xavier Schmider,<sup>18</sup>

Richard P. Schwarz,<sup>67</sup> Nicholas J. Scott,<sup>43</sup> Sara Seager,<sup>53,68,69</sup> Damien Ségransan,<sup>10</sup>

Luisa Maria Serrano<sup>1</sup>,<sup>2</sup> Attila E. Simon,<sup>4</sup> Alexis M. S. Smith,<sup>31</sup> Manfred Steller,<sup>5</sup> Chris Stockdale,<sup>70</sup>

Gyula Szabó<sup>1</sup>,<sup>71,72,73</sup> Nicolas Thomas,<sup>4</sup> Eric B. Ting,<sup>43</sup> Amaury H. M. J. Triaud<sup>1</sup>,<sup>41</sup> Stéphane Udry,<sup>10</sup>

Vincent Van Eylen<sup>1</sup>,<sup>60</sup> Valérie Van Grootel,<sup>39</sup> Roland K. Vanderspek,<sup>53</sup> Valentina Viotto,<sup>50</sup>

Nicholas Walton<sup>74</sup> and Joshua N. Winn<sup>75</sup>

*Affiliations are listed at the end of the paper*

Accepted 2021 December 29. Received 2021 December 29; in original form 2021 July 12

## ABSTRACT

We report the discovery and characterization of a pair of sub-Neptunes transiting the bright K-dwarf TOI-1064 (TIC 79748331), initially detected in the *Transiting Exoplanet Survey Satellite (TESS)* photometry. To characterize the system, we performed and retrieved the *Characterising Exoplanets Satellite (CHEOPS)*, *TESS*, and ground-based photometry, the High Accuracy Radial velocity Planet Searcher (HARPS) high-resolution spectroscopy, and Gemini speckle imaging. We characterize the host

<footnote>* E-mail: tgw1@st-andrews.ac.uk</footnote>

<footnote>† Banting fellow.</footnote>
star and determine $T_{\text{eff},*} = 4734 \pm 67$ K, $R_* = 0.726 \pm 0.007 R_\odot$, and $M_* = 0.748 \pm 0.032 M_\odot$. We present a novel detrending method based on point spread function shape-change modelling and demonstrate its suitability to correct flux variations in CHEOPS data. We confirm the planetary nature of both bodies and find that TOI-1064 b has an orbital period of $P_b = 6.44387 \pm 0.00003$ d, a radius of $R_b = 2.59 \pm 0.04 R_\odot$, and a mass of $M_b = 13.5_{-1.8}^{+1.7} M_\odot$, whilst TOI-1064 c has an orbital period of $P_c = 12.22657_{-0.00004}^{+0.00005}$ d, a radius of $R_c = 2.65 \pm 0.04 R_\odot$, and a $3\sigma$ upper mass limit of 8.5 $M_\odot$. From the high-precision photometry we obtain radius uncertainties of ~1.6 per cent, allowing us to conduct internal structure and atmospheric escape modelling. TOI-1064 b is one of the densest, well-characterized sub-Neptunes, with a tenuous atmosphere that can be explained by the loss of a primordial envelope following migration through the protoplanetary disc. It is likely that TOI-1064 c has an extended atmosphere due to the tentative low density, however further radial velocities are needed to confirm this scenario and the similar radii, different masses nature of this system. The high-precision data and modelling of TOI-1064 b are important for planets in this region of mass-radius space, and it allows us to identify a trend in bulk density-stellar metallicity for massive sub-Neptunes that may hint at the formation of this population of planets.

**Key words:** techniques: photometric – techniques: radial velocities – planets and satellites: composition – planets and satellites: detection – planets and satellites: interiors – stars: individual: TOI-1064 (TIC 79748331, Gaia EDR3 6683371847364921088).

# 1 INTRODUCTION

The study of exoplanets yields improvements in our understanding of planet formation and evolution, and has done so since the discovery of the first hot Jupiter (Mayor & Queloz 1995) challenged planet formation paradigms and compelled the development of planet migration theories. Since this initial discovery and subsequent pioneering works (Charbonneau et al. 2000, 2009; Rivera et al. 2005; Léger et al. 2009; Queloz et al. 2009), thousands of planets have been detected using ground-based and space-based photometric instruments such as the Hungarian Automated Telescope Network (HATNet), Super Wide Angle Search for Planets (SuperWASP), Trans-Atlantic Exoplanet Survey (TrES), *Convection, Rotation and planetary Transits (CoRoT)*, *Kepler*, *K2*, and *Transiting Exoplanet Survey Satellite (TESS)*; Bakos et al. 2002; Alonso et al. 2004; Baglin et al. 2006; Pollacco et al. 2006; Borucki et al. 2010; Howell et al. 2014; Ricker et al. 2015), and radial velocity (RV) surveys, as in High Resolution Echelle Spectrometer (HIRES), High Accuracy Radial velocity Planet Searcher (HARPS), High Accuracy Radial velocity Planet Searcher for the Northern hemisphere (HARPS-N), Spectrographie pour l'Observation des Phénomènes des Intérieurs stellaires et des Exoplanètes (SOPHIE), and Echelle Spectrograph for Rocky Exoplanet and Stable Spectroscopic Observations (ESPRESSO; Vogt et al. 1994; Mayor et al. 2003; Bouchy & Sophie Team 2006; Cosentino et al. 2012; Pepe et al. 2021), that have subsequently paved the way for follow-up characterization of these systems and yielded further insight into planet formation and evolution.

One such advancement has come from the precise measurements of the radii and masses of transiting exoplanets orbiting bright host stars, as knowledge of these key properties permits the determination of both empirical mass-radius curves (Weiss et al. 2013; Weiss & Marcy 2014; Chen & Kipping 2017; Otegi, Bouchy & Helled 2020) and bulk densities. Planetary densities can be converted into bulk compositional information using internal structure models (Fortney, Marley & Barnes 2007; Seager et al. 2007; Sotin, Grasset & Mocquet 2007; Valencia, Sasselov & O'Connell 2007; Zeng & Sasselov 2013; Howe, Burrows & Verne 2014; Dorn et al. 2015, 2017; Brugger et al. 2017; Zeng et al. 2019; Mousis et al. 2020; Aguichine et al. 2021). However, to break degeneracies in the modelling, both elemental stellar abundances (Santos et al. 2015; Adibekyan et al. 2021) and low uncertainties on planetary radii and masses are needed. The former can be resolved via spectroscopic follow-up to produce a combined high-resolution spectrum that can be used to extract

stellar information. High-precision photometry and RVs are needed to constrain radii and masses.

Studying multiplanet systems offers a further, unique opportunity for planetary system characterization via comparative planetology, as they have formed within the same protoplanetary disc. Therefore, observations can constrain formation and evolution models from the analysis of mutual inclinations (Steffen et al. 2010; Fang & Margot 2012; Fabrycky et al. 2014) and eccentricities (Limbach & Turner 2015; Van Eylen & Albrecht 2015; Mills et al. 2019; Van Eylen et al. 2019). Additionally, determination of orbital spacings between pairs of planets (Weiss et al. 2018; Jiang, Xie & Zhou 2020) and the assessment of correlations or discrepancies between the spacings and planetary sizes (Ciardi et al. 2013) and masses (Lissauer et al. 2011; Millholland, Wang & Laughlin 2017; Adams et al. 2020; Weiss & Petigura 2020) can lead to a better understanding on a demographic scale. Moreover, should individual planets be well characterized, the bulk internal and atmospheric compositions could be compared and, when combined the fluxes of stellar radiation at the positions of the planets, could inform formation and evolution modelling.

In this paper, we report the discovery of the two-planet system TOI-1064 using *TESS*, *CHaracterising ExoPlanets Satellite (CHEOPS)*, Las Cumbres Observatory Global Telescope (LCOGT), Next-Generation Transit Survey (NGTS), Antarctic Search for Transiting ExoPlanets (ASTEP), All Sky Automated Survey for SuperNovae (ASAS-SN), and Wide Angle Search for Planets (WASP) photometry (Sections 2.1–2.5, 2.8, and 2.9), HARPS RVs (Section 2.6), and Gemini speckle imaging (Section 2.7). We determine the stellar parameters of the host star (Section 3), and through photometric and RV analyses we precisely characterize the planetary properties of TOI-1064 b and TOI-1064 c (Section 4). Utilizing these results, we present internal structure and atmospheric modelling of the planets (Section 5), discuss important aspects of the system (Section 6), and summarize our conclusions (Section 7).

# 2 OBSERVATIONS

To fully characterize the TOI-1064 system, we collate photometric, spectroscopic, and imaging data from multiple sources detailed below.

## 2.1 TESS

The *Transiting Exoplanet Survey Satellite (TESS)* spacecraft has been conducting survey-mode observations of stars in the *TESS*
Input Catalog (TIC; Stassun et al. 2018, 2019) since launch in 2018, covering 26 sectors in the 2-yr primary mission (Ricker et al. 2015). Subsequently, *TESS* has entered its extended mission with a main goal to provide additional, shorter cadence photometry of targets observed during the primary mission. Initially listed as TIC 79748331, TOI-1064 was observed by *TESS* in Sector 13, camera 1, CCD 2, from 2019 June 19 to 2019 July 17. A gap occurred at the mid-point of the sector due to telemetry operations around spacecraft periastron passage, yielding 27.51 d of science observations. Individual frames were processed into 2 min cadence calibrated pixel files and reduced into light curves by the Science Processing Operations Center (SPOC; Jenkins et al. 2016) at NASA Ames Research Center. A transit search was conducted on the Pre-search Data Conditioning Simple Aperture Photometric (PDCSAP) light curves (Smith et al. 2012; Stumpe et al. 2012, 2014) using a wavelet-based, noise-compensating matched filter (Jenkins 2002; Jenkins et al. 2010) with two candidate planetary signals passing all diagnostic tests (Twicken et al. 2018; Li et al. 2019). The *TESS* Science Operations Center (TSO) announced these as *TESS* Objects of Interest (Guerrero et al. 2021). The two candidates were detected at signal-to-noise ratio (S/N) of 9.4 and 8.1, respectively, and no additional transiting planet signatures were identified in the residual light curves. Following the commencement of the extended mission, *TESS* observed TOI-1064 in Sector 27, camera 1, CCD 1, from 2020 July 5 to 2020 July 30. Allowing for the mid-sector observation gap, this results in 23.35 d of science observations that were subsequently processed and reduced by the SPOC into 20 s cadence photometry (in GI Cycle 3 program 3278; PI: Andrew Mayo). In total, eight transits of TOI-1064 b and four transits of TOI-1064 c were identified with the *TESS*-reported periods of the two candidates of 6.4422±0.0015 and 12.2371±0.0071 d, respectively.

We retrieved both sectors of photometry (denoted LC and FAST LC for the 2 min and 20 s cadence data, respectively) from the Mikulski Archive for Space Telescopes (MAST), by selecting the systematics-corrected PDCSAP light curves (Smith et al. 2012; Stumpe et al. 2012, 2014), and using the default quality bitmask. Finally, for the photometric analysis detailed below, we rejected data points flagged by the SPOC as being of bad quality (QUALITY > 0) and those with Not-a-Number fluxes or flux errors. This quality control yields a total of 89 642 data points from both *TESS* sectors, with the resulting light curves from Sectors 13 and 27 shown in Fig. 1.

## 2.2 CHEOPS

The *Characterising Exoplanets Satellite (CHEOPS)* spacecraft (Benz et al. 2021), launched on 2019 December 18 from Kourou, French Guiana, is an ESA small-class mission with the prime aim of observing bright ($V < 12$ mag), exoplanet-hosting stars to obtain ultrahigh-precision photometry (Broeg et al. 2013). Since launch, *CHEOPS* successfully passed In-Orbit Commissioning and was verified to achieve a photometric noise of 15 ppm per 6 h for a $V \sim 9$ mag star (Benz et al. 2021). A study of the photometric precision of *CHEOPS* has shown that the depth uncertainty of a 500 ppm transit from one *CHEOPS* observation is comparable to eight transits from *TESS* (Bonfanti et al. 2021).

The first scientific results from *CHEOPS* show the range of science that can be achieved with such a high-precision instrument (Lendl et al. 2020; Borsato et al. 2021; Delrez et al. 2021; Leleu et al. 2021; Morris et al. 2021; Szabó et al. 2021; Van Grootel et al. 2021; Maxted et al. 2022). One such study reports on the improvement in precision of exoplanet sizes in the HD 108236 system (Bonfanti et al. 2021)

that highlights a key scientific goal of the *CHEOPS* mission: the refinement of exoplanet radii to decrease bulk density uncertainties, and thereby allowing internal structure and atmospheric evolution modelling.

To better characterize and to secure the validation of both planetary candidates, we obtained six visits of TOI-1064 with the *CHEOPS* spacecraft between 2020 July 12 and 2020 August 31, as a part of the Guaranteed Time Observers programme, yielding a total of 68.38 h on target. We identify four transits of TOI-1064 b and three transits of TOI-1064 c across these runs. A breakdown of the individual visit start times and durations is detailed in Table 1. For all visits, we used an exposure time of 60 s.

As the *CHEOPS* spacecraft is in a low-Earth orbit, sections of observations are unobtainable due to the onboard rejection of images due to the level of stray light being higher than the accepted threshold, occultations of the target by the Earth, or passages through the South Atlantic Anomaly (SAA) during which no data are downlinked. These effects occur on orbital time-scales (~98.77 min) and lead to a decrease in the observational efficiency. As can be seen in Fig. 2 (and Figs B1–B5), these interruptions are apparent in the *CHEOPS* photometric data with the efficiencies for the six visits listed in Table 1.

Data for all visits were automatically processed using the latest version of the *CHEOPS* data reduction pipeline (DRP v13; Hoyer et al. 2020). The DRP undertakes image calibration, such as bias, gain, non-linearity, dark current, and flat-fielding corrections, and conducts rectifications of instrumental and environmental effects, for example cosmic ray impacts, smearing trails of nearby stars, and background variations. Subsequently, it performs aperture photometry on the corrected frames using a set of defined-radius apertures: $R = 22.5$ arcsec (RINF), 25.0 arcsec (DEFAULT), and 30.0 arcsec (RSUP), and an additional aperture that aims to optimize the radius based on contamination level and instrumental noise (ROPT). For the six observations of TOI-1064 this radius was determined to be between 15.0 and 16.0 arcsec, due to the nearby source discussed below. Furthermore, the DRP computes the contribution of nearby stars to the photometry by simulating the field of view (FoV) of the *CHEOPS* observations of the target using the *Gaia* Data Release 2 (DR2) catalogue (Gaia Collaboration 2018) as an input source list for objects' locations and brightnesses. By conducting aperture photometry on the simulated FoV with the target removed, light curves of the contamination from nearby sources are produced, as detailed in section 6.1 of Hoyer et al. (2020). As can be seen in Fig. 3, in the *CHEOPS* FoV there is a nearby object (Gaia EDR3 6683371813007224960, $\Delta G = +3.9$ mag) 24.6 arcsec away from the target that may affect the photometry of the target, and thus the contamination estimates were subtracted from the light curves of TOI-1064. The *Gaia* parallax and proper motion data indicate that this object and the target do not form a larger bound system. The right-hand panel of Fig. 3 shows the inferred FoV of the target with the contamination removed. Fig. 3 also reveals additional multiple nearby sources with stars that have $G < 17$ mag numbered. It should be noted that the sources numbered 1–5 are the same as those detected in the *TESS* FoV (Fig. 6), whereas 6–9 are too faint to be seen in the *TESS* target pixel file. The remaining objects in the FoV have $\Delta G > +7$ mag compared to TOI-1064 and thus, they will not contribute considerably to the photometry. For this study, we selected data sets that minimized the root mean square (rms) of the light curves, which for all visits were obtained with the RINF aperture. This radius minimized the contribution of nearby sources whilst ensuring that the majority of the target's point spread function (PSF) was within the aperture.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The Sector 13 (left) and Sector 27 (right) TESS light curves (with cadences of 2 min and 20 s, respectively) of TOI-1064 with the locations of the transits of planet b (green) and planet c (purple) shown. Top panel: the retrieved PDCSAP fluxes in blue with data binned every 1 h given as white circles. Middle panel: the photometry with transit and Matérn-3/2 kernel Gaussian process (GP) models used in the global analysis, as detailed in Section 4. Bottom panel: the photometry and transit model with trends modelling by the Matérn-3/2 kernel GP removed.</figcaption></figure>

<table><caption>**Table 1.** The CHEOPS observing log of TOI-1064.</caption><thead><tr><th>Visit<br>(#)</th><th>Planets</th><th>Start date<br>(UTC)</th><th>Duration<br>(h)</th><th>Data points<br>(#)</th><th>File key</th><th>Efficiency<br>(per cent)</th><th>Exp. time<br>(s)</th></tr></thead><tbody><tr><td>1</td><td>b, c</td><td>2020-07-12T10:11:36</td><td>13.64</td><td>473</td><td>CH_PR100031_TG027501_V0200</td><td>58</td><td>60</td></tr><tr><td>2</td><td>b</td><td>2020-07-31T23:57:36</td><td>7.50</td><td>290</td><td>CH_PR100031_TG029601_V0200</td><td>64</td><td>60</td></tr><tr><td>3</td><td>c</td><td>2020-08-18T12:12:36</td><td>7.52</td><td>293</td><td>CH_PR100031_TG029701_V0200</td><td>65</td><td>60</td></tr><tr><td>4</td><td>b</td><td>2020-08-20T07:04:36</td><td>7.50</td><td>270</td><td>CH_PR100031_TG029101_V0200</td><td>60</td><td>60</td></tr><tr><td>5</td><td>b</td><td>2020-08-26T18:28:17</td><td>7.99</td><td>263</td><td>CH_PR100031_TG029501_V0200</td><td>55</td><td>60</td></tr><tr><td>6</td><td>c</td><td>2020-08-30T13:04:57</td><td>24.23</td><td>830</td><td>CH_PR100031_TG027401_V0200</td><td>57</td><td>60</td></tr></tbody></table>

Because of the nature of the CHEOPS orbit and the rotating FoV, non-astrophysical, short-term photometric trends caused by a varying background, nearby contaminants, or other sources can be seen in the data. Several studies (e.g. Lendl et al. 2020; Bonfanti et al. 2021; Delrez et al. 2021; Leleu et al. 2021) have found success in removing these systematics by conducting a linear decorrelation using several basis vectors, such as background, contamination, orbital roll angle, and *x* and *y* centroid positions.

Upon inspection of the CHEOPS observations of TOI-1064, we found significant flux variations on orbital time-scales. The selection of basis vectors of concern typically involves assessing the Bayesian information criterion (BIC) upon the detrending of the data using a combination of DRP-provided vectors. However, in this study a more data-driven approach was taken to identify which basis vectors were to be used in the detrending, as is detailed in Appendix A.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The normalized light curve of the first CHEOPS visit of TOI-1064 with the locations of the transits of planet b (green) and planet c (purple) shown. Top panel: the raw data reduction pipeline (DRP)-produced fluxes. Second panel: the DRP fluxes with the fitted linear model produced from the measured point spread function (PSF) shape changes, as detailed in Appendix A. Third panel: the detrended fluxes with the nominal transit model from the global analysis (Section 4.3). Bottom panel: residuals to the fit.</figcaption></figure>

## 2.3 LCOGT

We acquired ground-based time series follow-up photometry of TOI-1064 as part of the TESS Follow-up Observing Program (TFOP)<fnref n="1" /> using the TESS Transit Finder, which is a customized version of the TAPIR software package (Jensen 2013), to schedule our transit observations.

We observed full transits of TOI-1064 b in the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS) z-short band on 2020 June 3, 2020 June 16, and 2020 August 26 from the Las Cumbres Observatory Global Telescope (LCOGT; Brown et al. 2013) 1.0 m network node at South African Astronomical Observatory (SAAO). We also observed full transits of TOI-1064 c on 2019 August 30 in Pan-STARRS z-short band from the LCOGT 1 m node at Cerro Tololo Inter-American Observatory and on 2019 October 5 in B band from the LCOGT 1 m network node at SAAO. The 0.389 arcsec pixel scale images were calibrated by the standard LCOGT BANZAI pipeline (McCully et al. 2018), and photometric data were extracted with ASTROIMAGEJ (Collins et al. 2017). The images were focused and have typical stellar PSFs with a full width at half-maximum (FWHM) of ~2 arcsec, and circular apertures with radius ~4 arcsec were used to extract the differential photometry.

## 2.4 NGTS

The Next Generation Transit Survey (NGTS; Wheatley et al. 2018) was used to observe a partial transit egress of TOI-1064 c on 2019 October 17. NGTS is a photometric facility, which is located at the ESO Paranal Observatory in Chile and consists of 12 robotic telescopes, each with a 20 cm diameter, an 8 deg² FoV, and a pixel scale of 4.97 arcsec. Each NGTS telescope is operated independently, however simultaneous observations with multiple NGTS telescopes have been shown to greatly improve the photometric precision

achieved (Bryant et al. 2020; Smith et al. 2020). TOI-1064 was observed using two NGTS telescopes with both telescopes observing with an exposure time of 10 s and using the custom NGTS filter (520–890 nm). A total of 2016 images were obtained during the observation.

The image reduction was performed using a custom photometry pipeline (Bryant et al. 2020). The source extraction and photometry in the pipeline are performed using the SEP PYTHON library (Bertin & Arnouts 1996; Barbary 2016). Comparison stars that are similar to TOI-1064 in brightness, colour, and CCD position were automatically identified using Gaia DR2 (Gaia Collaboration 2018).

## 2.5 ASTEP

We observed five transits of the TOI-1064 planets as part of the Antarctica Search for Transiting ExoPlanets (ASTEP) program (Guillot et al. 2015; Mékarnia et al. 2016). The 0.4 m ASTEP telescope is located at the French/Italian Concordia Station on the East Antarctic Plateau. It is equipped with an FLI Proline science camera with a KAF-16801E, 4096 × 4096 front-illuminated CCD. The camera has an image scale of 0.93 arcsec pixel⁻<fnref n="1" /> resulting in a 1° × 1° corrected FoV. The focal instrument dichroic plate splits the beam into a blue wavelength channel for guiding, and a non-filtered red science channel roughly matching an $R_c$ transmission curve. The telescope is automated or remotely operated when needed. Because of the extremely low data transmission rate at the Concordia Station, the data are processed on-site using an automated IDL-based pipeline, and the result is reported via email and then transferred to Europe on a server in Rome, Italy. The raw light curves of about 1000 stars are then available for deeper analysis.

Three full transits of planet TOI-1064 c were observed on 2020 June 30, 2020 August 18, and 2021 March 26. On 2021 March 14 both an egress of TOI-1064 c and a full transit of TOI-1064 b were observed during the same light curve. In all cases, the weather was fair, with temperatures ranging between -46°C and -65°C, a stable relative humidity around 50 per cent and wind speeds between 2 and 7 m s⁻<fnref n="1" />. Exposure times were chosen to be 40 s in 2020 and 50–60 s in 2021, with a read-out time of about 25 s. A 9–13 arcsec radius photometric aperture was found to give the best results.

## 2.6 HARPS

We acquired 26 high-resolution ($R = 115\,000$) spectra of TOI-1064 between 2019 September 8 and 2019 October 29 using the High Accuracy Radial velocity Planet Searcher (HARPS) spectrograph (Mayor et al. 2003) mounted at the ESO 3.6-m telescope of La Silla Observatory. The observations were carried out as part of the observing program 1102.C-0923. We set the exposure time to 1800 s, leading to an S/N per pixel at 550 nm ranging between 18 and 66, with a median of 51. We used the second fibre of the instrument to monitor the sky background and we reduced the data using the dedicated HARPS Data Reduction Software (DRS; Lovis & Pepe 2007). For each spectrum, the DRS provides also the contrast, the FWHM and the bisector inverse slope (BIS) of the cross-correlation function (CCF). We also extracted additional activity indices and spectral diagnostics, namely the Ca II H&K lines activity indicator (S-index), $H\alpha$, Na D1, and Na D2, using the code TERRA (Anglada-Escudé & Butler 2012). The 26 DRS and TERRA RV measurements and activity indicators are listed in Table 2. Time stamps are given in Barycentric Julian Date in the Barycentric Dynamical Time (BJD; TDB, from the French Temps Dynamique Barycentrique).

<footnote n="1">¹https://tess.mit.edu/followup</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The 3.3 arcmin diameter FoV of TOI-1064 with the PSF centroid indicated by the red cross and the photometric aperture given by the red circle. Left-hand panel: the FoV as observed by CHEOPS. Middle panel: DRP model of the FoV showing the contamination estimate of the background stars with the target removed. Right-hand panel: DRP model of the FoV showing TOI-1064 with contamination estimate subtracted.</figcaption></figure>

<table><caption>**Table 2.** A sample of the HARPS observing log of TOI-1064. This table is available in its entirety online.</caption><thead><tr><th>Time<br>(BJD - 245 7000)</th><th>RV<br>(km s<sup>-1</sup>)</th><th>σ<sub>RV</sub><br>(km s<sup>-1</sup>)</th><th>BIS<br>(km s<sup>-1</sup>)</th><th>FWHM<br>(km s<sup>-1</sup>)</th><th>Contrast</th><th>S-index</th><th>Hα</th><th>Na D1</th><th>Na D2</th></tr></thead><tbody><tr><td>1734.5360</td><td>21.2169</td><td>0.0024</td><td>0.0487</td><td>6.4368</td><td>46.1189</td><td>0.5954</td><td>0.6869</td><td>1.2457</td><td>0.9962</td></tr><tr><td>1738.5247</td><td>21.2254</td><td>0.0016</td><td>0.0487</td><td>6.4494</td><td>45.9029</td><td>0.5870</td><td>0.6882</td><td>1.2341</td><td>0.9905</td></tr><tr><td>1739.5374</td><td>21.2245</td><td>0.0023</td><td>0.0562</td><td>6.4469</td><td>45.7285</td><td>0.6237</td><td>0.6773</td><td>1.2366</td><td>0.9907</td></tr><tr><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr></tbody></table>

## 2.7 Gemini-South

Because of the 21 arcsec size of TESS pixels the TESS light curves may be contaminated by nearby sources. As mentioned above, for CHEOPS data sets contamination from stars in the Gaia DR2 catalogue is simulated and removed. However, should there be sources very near TOI-1064 that were not detected with Gaia, contamination may still occur. If an exoplanet host star has a spatially close companion, that companion (bound or line of sight) can create a false-positive transit signal if it is, for example, an eclipsing binary. ‘Third-light’ flux from the companion star can lead to an underestimated planetary radius if not accounted for in the transit model (Ciardi et al. 2015) and cause non-detections of small planets residing with the same exoplanetary system (Lester et al. 2021). Additionally, the discovery of close, bound companion stars, which exist in nearly one-half of FGK-type stars (Matson et al. 2018), provides crucial information toward our understanding of exoplanetary formation, dynamics, and evolution (Howell et al. 2021). Thus, to search for close-in bound companions unresolved in TESS or other ground-based follow-up observations, we obtained high-resolution speckle imaging observations of TOI-1064.

TOI-1064 was observed on 2019 September 12 using the Zorro speckle instrument on Gemini-South.<fnref n="2" /> Zorro provides simultaneous speckle imaging with a pixel scale of 0.01 arcsec in two bands (562 and 832 nm) with output data products including a reconstructed image with robust contrast limits on companion detections (e.g. Howell et al. 2016). Five sets of 1000 × 0.06 s exposures were collected and subjected to Fourier analysis in our standard reduction pipeline (see Howell et al. 2011).

## 2.8 ASAS-SN

To further assess the host star, we obtained publicly available All Sky Automated Survey for SuperNovae (ASAS-SN) V-band photometry (Shappee et al. 2014; Kochanek et al. 2017) of TOI-1064 taken over five consecutive seasons between 2014 May 4 and 2018 September 24. Upon inspection of the data, a dimming trend of roughly 0.3 mag was seen over the 4 yr of observations, with an abrupt increase (roughly 0.2 mag) in flux occurring during the final season. Therefore, as the goal of using this data set was to study shorter period variation, we rejected data after BJD 245 8300 in order to avoid photometry taken during the brightening event, and we removed the long-term trend by modelling the data set with a broad Savitzky-Golay smoothing filter and dividing the fluxes by this model. Lastly, we removed outliers by conducting a 5σ clip, which resulted in 401 data points taken on 328 epochs that can be seen in Fig. 4(a). The median flux errors for this data set are 18 ppt and thus, whilst these data are not precise enough for transit detection of the two planetary candidates around TOI-1064, they can be used to study photometric variability in the host star.

## 2.9 WASP

Additionally, TOI-1064 was observed during the SuperWASP project (Pollacco et al. 2006), with data taken between 2008 March 26 and 2014 November 10. We retrieved photometry that had been extracted following standard procedures (Pollacco et al. 2006), and detrended for systematic effects using SYSREM (see Collier Cameron et al. 2006; Mazeh, Tamuz & Zucker 2007) that preserves stellar variability and transit-like features in the data set. Flux and flux uncertainty outliers were rejected via a 5σ clipping, with the subsequent light curve shown in Fig. 4(b). This yielded a total of 66 832 data points on 387 epochs covering four seasons of observations. Similarly to the ASAS-SN data set, the median flux error for the WASP observations

<footnote n="2">²https://www.gemini.edu/sciops/instruments/alopeke-zorro/</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** (a) ASAS-SN and (b) WASP data used to study the photometric variability of TOI-1064. Left-hand panels: the full light curves (light blue) and nightly binned data (dark blue). Right-hand panels: Lomb–Scargle periodograms of the nightly binned data, with the orbital periods of planets b and c shown as vertical green and purple lines, respectively. The horizontal black dotted, dash–dot, dashed, and solid lines (from top to bottom) are the 0.1 per cent, 1 per cent, 5 per cent, and 10 per cent false alarm probabilities (FAPs), respectively. The adopted rotation period is given as red vertical marks.</figcaption></figure>

(28 ppt) means that the photometry is not precise enough to be used to detect the transits of the planetary candidates. However, due to the substantial baseline of the data, they can be used to search for photometric variability.

# 3 CHARACTERIZATION OF TOI-1064

## 3.1 Atmospheric properties and abundances

The spectral analyses of the host star of the TOI-1064 system were based on the coaddition of all the RV observations carried out with the HARPS spectrograph, detailed in Section 2.6, at a spectral resolution of 115 000. We began with applying the SPECMATCH-EMP (Yee, Petigura & von Braun 2017) software to our data. Our observed optical spectrum is compared to a spectral library of approximately 400 stars observed by Keck/HIRES with spectral classes M5 to F1. Interpolating and minimizing differences, the direct output is the effective temperature of the star $T_{\text{eff}}$, the stellar radius $R_{\star}$, and the iron abundance [Fe/H] that are found to be 4780 ± 110 K, 0.77 ± 0.10 $R_{\odot}$, and 0.05 ± 0.09 dex for TOI-1064, respectively.

We then calculate a completely synthesized spectrum, using the SME (Spectroscopy Made Easy; Valenti & Piskunov 1996; Piskunov & Valenti 2017) package version 5.22 in the fashion explained in more detail in e.g. Fridlund et al. (2017). We use as starting values the $T_{\text{eff}}$ and [Fe/H] from SPECMATCH-EMP. We model the $v \sin i_{\star}$ by fitting large numbers of narrow and unblended metal lines between 6000 and 6500 Å. By determining the depths and profiles of Na, Ca, and Fe lines, the [Na/H], [Ca/H], and [Fe/H] abundances are computed, as well as the logarithm of the stellar surface gravity, log g, from the Ca I triplet around 6200 Å. We hold the turbulent velocities, $V_{\text{mic}}$ and $V_{\text{mac}}$ fixed to 0.47 and 1.2 km s⁻¹, respectively, based on $T_{\text{eff}}$ and log g (Adibekyan et al. 2012b; Doyle et al. 2014). As a final step, our model was checked with the Na I doublet λλ5589 and 5896 Å, sensitive to both $T_{\text{eff}}$ and log g. Through these steps we arrive at the values listed in Table 3 that are in excellent agreement with the SPECMATCH-EMP model, and confirm that TOI-1064 is a K-dwarf.

We also derive stellar parameters using the ARES + MOOG tools (Sousa 2014). In particular, the equivalent widths (EWs) of iron lines are measured using the ARES code<fnref n="3" /> (Sousa et al. 2007, 2015) and the

iron abundance are then computed using Kurucz model atmospheres (Kurucz 1993) and the radiative transfer code MOOG (Sneden 1973). The values reported here are obtained via convergence of both ionization and excitation equilibria. We obtain $T_{\text{eff}} = 4636 \pm 102$ K, $\log g = 4.52 \pm 0.24$, and [Fe/H] = −0.04 ± 0.05, which are consistent with the results presented above. Because of the smaller uncertainties we adopt the SME values.

Using the stellar parameters listed in Table 3, we determine the abundances of [Mg/H] and [Si/H] to be 0.03 ± 0.06 and 0.06 ± 0.08, respectively, using the classical curve-of-growth analysis method assuming local thermodynamic equilibrium. We use the ARES v2 code (Sousa et al. 2015) to measure the EWs of the spectral lines of these elements. Then we use a grid of Kurucz model atmospheres (Kurucz 1993) and the radiative transfer code MOOG (Sneden 1973) to convert the EWs into abundances. When doing so, we closely follow the methods described in e.g. Adibekyan et al. (2012a, 2015). In Table 3, we also include the right-handed, heliocentric Galactic spatial velocities that, along with the derived [Fe/H], indicate that TOI-1064 is a member of the Galactic thin disc population.

## 3.2 Radius, mass, and age

In order to determine the stellar radius of TOI-1064 we use a modified version of the infrared flux method (IRFM; Blackwell & Shallis 1977) that allows the derivation of stellar angular diameters and effective temperatures through known relationships between these properties, and an estimate of the apparent bolometric flux, recently detailed in Schanche et al. (2020). We perform the IRFM in a Markov chain Monte Carlo (MCMC) approach in which the stellar parameters, derived via the spectral analysis detailed above, are used as priors in the construction of spectral energy distributions (SEDs) from stellar atmospheric models. The SEDs are subsequently attenuated to account for reddening and used for synthetic photometry. This is conducted by convolving the SED with the broad-band response functions for the chosen bandpasses, with the fluxes compared to the observed data to compute the apparent bolometric flux. For this study, we retrieve broad-band fluxes and uncertainties for TOI-1064 from the most recent data releases for the following bandpasses: Gaia G, $G_{\text{BP}}$, and $G_{\text{RP}}$, Two Micron All Sky Survey (2MASS) J, H, and K, and Wide-field Infrared Survey Explorer (WISE) W1 and W2 (Skrutskie et al. 2006; Wright et al. 2010; Gaia Collaboration 2021), and use the ATLAS catalogues (Castelli & Kurucz 2003) of model stellar SEDs. Internal to the MCMC, the computed posterior distributions of stellar angular diameters are converted to distributions of stellar radii using

<footnote n="3">³The last version of ARES code (ARES v2) can be downloaded at http://www.astro.up.pt/~sousasag/ares</footnote>
<table><caption>**Table 3.** Stellar properties of TOI-1064.</caption><thead><tr><th colspan="3">TOI-1064</th></tr></thead><tbody><tr><td>2MASS</td><td>J19440094-4733417</td><td></td></tr><tr><td>Gaia EDR3</td><td>6683371847364921088</td><td></td></tr><tr><td>TIC</td><td>79748331</td><td></td></tr><tr><td>UCAC2</td><td>11441398</td><td></td></tr><tr><td>Parameter</td><td>Value</td><td>Note</td></tr><tr><td>α (J2000)</td><td>19<sup>h</sup>44<sup>m</sup>00<sup>s</sup>95</td><td>1</td></tr><tr><td>δ (J2000)</td><td>-47°33'41"75</td><td>1</td></tr><tr><td>μ<sub>α</sub> (mas yr<sup>-1</sup>)</td><td>-3.543 ± 0.015</td><td>1</td></tr><tr><td>μ<sub>δ</sub> (mas yr<sup>-1</sup>)</td><td>-100.885 ± 0.012</td><td>1</td></tr><tr><td>σ (mas)</td><td>14.532 ± 0.015</td><td>1</td></tr><tr><td>d (pc)</td><td>68.81 ± 0.07</td><td>5</td></tr><tr><td>RV (km s<sup>-1</sup>)</td><td>20.7 ± 0.7</td><td>1</td></tr><tr><td>U (km s<sup>-1</sup>)</td><td>11.50 ± 0.61</td><td>5<sup>a</sup></td></tr><tr><td>V (km s<sup>-1</sup>)</td><td>-34.29 ± 0.08</td><td>5<sup>a</sup></td></tr><tr><td>W (km s<sup>-1</sup>)</td><td>-14.32 ± 0.34</td><td>5<sup>a</sup></td></tr><tr><td>V (mag)</td><td>10.95 ± 0.06</td><td>2</td></tr><tr><td>G<sub>BP</sub> (mag)</td><td>11.207 ± 0.003</td><td>1</td></tr><tr><td>G (mag)</td><td>10.645 ± 0.003</td><td>1</td></tr><tr><td>G<sub>RP</sub> (mag)</td><td>9.939 ± 0.004</td><td>1</td></tr><tr><td>J (mag)</td><td>9.10 ± 0.02</td><td>3</td></tr><tr><td>H (mag)</td><td>8.63 ± 0.04</td><td>3</td></tr><tr><td>K (mag)</td><td>8.47 ± 0.03</td><td>3</td></tr><tr><td>W1 (mag)</td><td>8.41 ± 0.02</td><td>4</td></tr><tr><td>W2 (mag)</td><td>8.48 ± 0.02</td><td>4</td></tr><tr><td>T<sub>eff</sub> (K)</td><td>4734 ± 67</td><td>5; spectroscopy</td></tr><tr><td>log g (cm s<sup>-2</sup>)</td><td>4.60 ± 0.06</td><td>5; spectroscopy</td></tr><tr><td>[Fe/H] (dex)</td><td>0.05 ± 0.08</td><td>5; spectroscopy</td></tr><tr><td>[Mg/H] (dex)</td><td>0.03 ± 0.06</td><td>5; spectroscopy</td></tr><tr><td>[Si/H] (dex)</td><td>0.06 ± 0.08</td><td>5; spectroscopy</td></tr><tr><td>[Ca/H] (dex)</td><td>0.11 ± 0.10</td><td>5; spectroscopy</td></tr><tr><td>[Na/H] (dex)</td><td>0.17 ± 0.12</td><td>5; spectroscopy</td></tr><tr><td>ν sin i<sub>α</sub> (km s<sup>-1</sup>)</td><td>2.7 ± 0.7</td><td>5; spectroscopy</td></tr><tr><td>log R'<sub>HK</sub></td><td>-4.633 ± 0.024</td><td>5; spectroscopy</td></tr><tr><td>E(B-V)</td><td>0.056 ± 0.032</td><td>5; IRFM</td></tr><tr><td>R<sub>⋆</sub> (R<sub>☉</sub>)</td><td>0.726 ± 0.007</td><td>5; IRFM</td></tr><tr><td>M<sub>⋆</sub> (M<sub>☉</sub>)</td><td>0.748 ± 0.032</td><td>5; isochrones</td></tr><tr><td>t<sub>⋆</sub> (Gyr)</td><td>9.4 ± 3.8</td><td>5; isochrones</td></tr><tr><td>L<sub>⋆</sub> (L<sub>☉</sub>)</td><td>0.238 ± 0.014</td><td>5; from R<sub>⋆</sub> and T<sub>eff</sub></td></tr><tr><td>ρ<sub>⋆</sub> (ρ<sub>☉</sub>)</td><td>1.95 ± 0.10</td><td>5; from R<sub>⋆</sub> and M<sub>⋆</sub></td></tr><tr><td>ρ<sub>⋆</sub> (g cm<sup>-3</sup>)</td><td>2.76 ± 0.14</td><td>5; from R<sub>⋆</sub> and M<sub>⋆</sub></td></tr></tbody></table>

Note. 1 – Gaia Collaboration (2021); 2 – Zacharias et al. (2013); 3 – Skrutskie et al. (2006); 4 – Wright et al. (2010); 5 – this work.

<sup>a</sup>Calculated via the right-handed, heliocentric Galactic spatial velocity formulation of Johnson & Soderblom (1987) using the proper motions, parallax, and RV from Gaia Collaboration (2021).

the Gaia Early Data Release 3 (EDR3) parallax (Gaia Collaboration 2021), from which we obtain the stellar radius of TOI-1064 and E(B-V) to be R<sub>⋆</sub> = 0.726 ± 0.007 R<sub>☉</sub> and E(B-V) = 0.056 ± 0.032, respectively. The distance to TOI-1064 is calculated in the IRFM using the Gaia EDR3 parallax with the parallax offset of Lindegren et al. (2021) applied. These stellar parameters are reported in Table 3.

By adopting T<sub>eff</sub>, [Fe/H], and R<sub>⋆</sub> as input parameters, we infer the stellar age t<sub>⋆</sub> and mass M<sub>⋆</sub> through stellar evolutionary models. To make our analysis more robust, we employ two different techniques each applied to a different set of stellar isochrones and tracks. The first technique derives t<sub>⋆</sub> and M<sub>⋆</sub> using the isochrone placement method described in Bonfanti et al. (2015) and Bonfanti, Ortolani & Nascimbeni (2016), which interpolates the input parameters within pre-computed grids of isochrones and tracks generated by the

PARSEC<sup>4</sup> v1.2S code (Marigo et al. 2017). The second technique, instead, directly fits the input parameters in the CLES (Code Liègeois d'Évolution Stellar; Scuflaire et al. 2008) code to then retrieve t<sub>⋆</sub> and M<sub>⋆</sub> following a Levenberg-Marquardt minimization scheme as described in Salmon et al. (2021). Once the two pairs of age and mass values are derived, first their consistency is checked through a χ<sup>2</sup> test, and then their probability distributions are combined together to provide the final t<sub>⋆</sub> and M<sub>⋆</sub> values with errors at the 1σ level (see Table 3). Specific details about the statistical derivation of both t<sub>⋆</sub> and M<sub>⋆</sub> may be found in Bonfanti et al. (2021).

From the stellar radius and mass determined via the IRFM and isochrone placement techniques, we obtain log g = 4.59 ± 0.02 that is in agreement with the value derived from the spectral analysis.

The derived stellar radius was checked with the PYTHON code ARIADNE (described in e.g. Acton et al. 2020) that fitted broad-band photometry to the PHOENIX v2 (Husser et al. 2013), BTSETTL (Allard, Homeier & Freytag 2012), Castelli & Kurucz (2003), and Kurucz (1979) atmospheric model grids, utilizing data in the following bandpasses Gaia EDR3 G, G<sub>BP</sub>, and G<sub>RP</sub>, 2MASS J, H, and K, WISE W1 and W2, and the Johnson B and V magnitudes from the AAVSO Photometric All Sky Survey (APASS). We used SME values for T<sub>eff</sub>, log g, and [Fe/H] as priors to the model. The final radius is computed with Bayesian model averaging. We obtain a stellar radius of 0.723 ± 0.007 R<sub>☉</sub>. Combining this radius with the surface gravity, we obtain a mass of 0.762 ± 0.050 M<sub>☉</sub>. Both values are in excellent agreement with the above derived mass and radius.

## 3.3 Stellar variability

Stellar activity can contribute strong signals that are apparent in RV observations and can hinder the detection and characterization of small exoplanets via precise measurements of planetary masses (Haywood et al. 2014; Rajpaul et al. 2015; Mortier et al. 2016; Dumusque et al. 2017; Faria et al. 2020). Therefore, there have been recent efforts made to mitigate the effect of stellar activity on RV observations of exoplanets (de Beurs et al. 2020; Collier Cameron et al. 2021), which will be discussed in greater detail in Section 4.2.

In order to properly account for the stellar activity we need to measure the stellar rotation period, which can be done via inspection of the RV and photometric data of TOI-1064. As can be seen in the Lomb–Scargle periodograms of standard stellar activity indicators of our HARPS observations (Fig. C1), no significant peaks potentially related to the stellar rotation period are apparent. Therefore, we assess the long-baseline ASAS-SN and WASP light curves to search for photometric variability. First, for both data sets we performed a nightly binning that results in rms scatters of 28 and 16 ppt for ASAS-SN and WASP, respectively. Subsequently, we produced Lomb–Scargle periodograms of both binned light curves and find significant peaks at 25.9 and 26.6 d as shown in Figs 4(a) and (b). To independently confirm this variability period we applied a Gaussian process (GP) regression with a quasi-periodic kernel to both ASAS-SN and WASP light curves separately. This kernel is chosen as multiple previous studies have found that it accurately represents flux modulation from stellar activity (Haywood et al. 2014; Dubber et al. 2019; Mortier et al. 2020). Using the JULIET PYTHON package (Espinoza, Kossakowski & Brahm 2019) with unconstrained priors, we obtain median and 1σ uncertainties for the rotation period of the kernel to be 27.0 ± 4.3 and 26.6 ± 1.0 d, for ASAS-SN and WASP, respectively. Because of the higher S/N of the WASP observations,

<footnote n="4"><sup>4</sup>PAdova and TRieste Stellar Evolution Code.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** The 5σ contrast curves of Gemini/Zorro speckle high-resolution images at 562 nm (blue) and 832 nm (red), with insets showing the central region of the images centred on TOI-1064.</figcaption></figure>

we adopt the value obtained from that data set and interpret this signal to be the stellar rotation period. These values agree well with the rotation period of 26.9 ± 1.6 d derived from the mean HARPS value of log R<sub>HK</sub> following the empirical relations in Mamajek & Hillenbrand (2008). It is worth noting that the derived rotation periods are distinct from the average lunar synodic period of approximately 29.53 d.

# 4 VALIDATING AND FITTING THE SYSTEM

## 4.1 High-resolution imaging analysis and planet validations

To establish if there is a potentially contaminating source nearby to TOI-1064 we analysed high-resolution images obtained with the Gemini/Zorro instrument by determining 5σ contrast curves in both bandpasses. Fig. 5 shows our final contrast curves and the reconstructed speckle images. We find that TOI-1064 is a single star from 20 mas out to 1.2 arcsec with no companion brighter than 5–8 mag below that of the target star beyond 100 mas. At the distance calculated in Section 3, this corresponds to the absence of a main-sequence star at spatial limits of 1.36–82 au. This isolation is supported by the Southern Astrophysical Research (SOAR)/High-Resolution Camera (HRCam) observations that find no companion within 1 arcsec at a 4.5 mag contrast limit (Ziegler et al. 2021).

Therefore, we conclude that the CHEOPS photometry is not contaminated by nearby sources undetected by Gaia, whereas all other sources are corrected for by subtracting the contamination estimate described in Section 2.2. However, as the 21 arcsec TESS pixels are substantially larger than those of CHEOPS, the nearby source with $ΔG = +3.9$ mag seen in Fig. 3 may affect the TESS photometry. Thus, we analysed the centroid position of the TESS observations in order to assess if this source affects the photometry.

First, we inspected the TESS target pixel files of both sectors to ascertain which nearby sources may affect the observations. As can be seen in the top panels of Fig. 6 there is only one $G < 15$ mag source with the core of its PSF within the pipeline aperture photometry mask. The wings of the PSFs of additional nearby sources may also fall within the aperture photometry mask; however, as the core of the PSFs and the majority of the flux from these stars are outside of the photometric mask, we conclude that these objects do not

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Top panels: example TESS target pixel files covering 3.5 × 3.5 arcmin² for Sectors 13 (left) and 27 (right) showing TOI-1064 (red star) and nearby stars with $G < 15$ mag (numbered, black stars). The TESS pipeline aperture photometry mask is shown as green squares with the average centroid position given as a red plus. Bottom panels: 35 × 35 arcsec² regions for Sectors 13 (left) and 27 (right) showing all nearby stars with $G < 20$ mag, including source 1 seen in the TESS target pixel files. In-transit offsets to out-of-transit positions for planets b and c are shown in green and purple, respectively.</figcaption></figure>

contaminate the TESS photometry of TOI-1064. To assess if the nearby object is affecting the data, we computed the average in-transit centroid positions for all transits of both planets in the system and subsequently determined the offsets between these values and the average out-of-transit centroid positions. The bottom panels of Fig. 6 show that the in-transit data are obtained on target with an average offset in Sector 13 of 0.1 ± 1.3 and 0.0 ± 1.5 arcsec, and in Sector 27 of 0.2 ± 3.6 and 0.2 ± 3.6 arcsec for planet candidates 01 and 02, respectively. Thus, as the singular $G < 15$ mag object in the aperture mask is 22 arcsec (roughly 15σ away) and 29 arcsec (roughly 8σ away) in Sectors 13 and 27, respectively, we conclude that it does not affect the TESS photometry of TOI-1064 and that the observed transits can be attributed to the target star. This finding is consistent with the TESS SPOC difference image centroiding results that constrained the source of the transits to within 4 arcsec of the target using data from both sectors. Lastly, it should be noted that whilst there are additional nearby objects, their comparatively greater $G$-band magnitudes ($ΔG > +9$ mag) mean that they do not contaminate the TESS photometry substantially.

Previous work has noted that multiplanet systems have a very low probability of being false positives (Lissauer et al. 2012) and thus we may consider these candidates to be verified based on the multiplicity of the system. However, to further confirm both planet candidates we utilized the statistical validation tool, TRICERATOPS (Giacalone et al. 2021), that uses stellar and transit parameters, transit photometry, and high-resolution speckle imaging in order to determine the false positive probability (FPP) of planetary candidates. This is done by first querying the TIC to calculate the flux contribution of nearby stars in a given aperture and determining which sources are bright enough
to feasibly host a body with a given transit depth. The TESS and any additional light curves are subsequently fitted with transit and eclipsing binary models to determine the probability in a Bayesian manner of a transiting planet, an eclipsing binary, and an eclipsing binary on twice the orbital period around the target or nearby star in the case of no unresolved companion, around the primary or secondary star in the case of an unresolved bound companion, and around the target or background star in the case of an unresolved background star, resulting in 18 scenarios. From these probabilities, the overall FPP of a transiting planet is computed. Using our stellar mass, radius, effective temperature, and parallax values from Table 3, and all photometric data we find that both TOI-1064 b and TOI-1064 c have FPP values <1 per cent, and thus confirm the presence of both planets.

## 4.2 PSF and CCF shape monitoring using SCALPELS

### 4.2.1 PSF

CHEOPS photometry is known to suffer from at least two types of systematic error arising from thermal effects that are found to be correlated with the output of temperature sensors in the telescope structure (Morris et al. 2021). It has been established that when the sunlight illumination pattern on the spacecraft changes, thermal flexure of the telescope structure causes subtle changes in the shape of the point spread function (PSF). This in turn causes a change in the fraction of the flux entering the pupil that falls within the fixed-radius photometric aperture in the focal plane used for signal extraction. The first type of systematic is a recognized secular 'ramp' effect as the telescope structure settles into a new equilibrium following a change of pointing direction (Maxted et al. 2022; Morris et al. 2021), and the second type of systematic is shorter term periodic fluctuation modulated on the 98.77 min orbital period as the spacecraft passes in and out of the Earth's shadow. Such 'ramps' have been seen in space-based photometry of other telescopes (Deming et al. 2006; Berta et al. 2012; Demory et al. 2015). The details of the interplay between spacecraft illumination, thermal changes in the telescope structure, and the response of the PSF shape cannot be modelled directly. We can, however, use a simple linear unsupervised machine-learning approach to establish correlations between changes in the shape of the PSF and changes in the encircled fraction of the total stellar flux.

Conceptually the problem is similar to the systematic errors produced in RV measurements by stellar-activity-driven changes in the shapes of spectral lines, and the solution we adopt here is modelled on the SCALPELS algorithm developed by Collier Cameron et al. (2021) for separating shape-driven RV offsets from genuine stellar Doppler shifts, and is detailed in Appendix A.

In this study, we applied our novel method to the CHEOPS photometry and the CCF-based SCALPELS (Collier Cameron et al. 2021) to the HARPS data in order to model flux modulation due to PSF shape changes and RV variation due to stellar activity. For the CHEOPS light curves, we used the aforementioned method on the DRP produced fluxes (Section 2.2) for the six visits separately, with the number of principal components, $\theta$, chosen by the leave-one-out cross-validation (LOOCV) method to be 24, 17, 23, 10, 26, and 42, respectively, corresponding to 5 per cent, 6 per cent, 8 per cent, 4 per cent, 10 per cent, and 5 per cent, of the vectors produced by the principal component analysis. Subsequently, we decorrelated the CHEOPS data sets against the selected vectors using the linear regression method, masking the in-transit fluxes. The linear models produced by the regression, and used to decorrelate the light curves,

<figure><img src="image_9.png" /><figcaption>**Figure 7.** The HARPS time series of TOI-1064, offset for clarity, showing the observed, median-subtracted RVs in blue, the CCF shape-based RVs in orange, the CCF shift-based RVs in green, the fit of two Keplerian orbits to the cleaned RVs in red, and the residuals to the fit in purple.</figcaption></figure>

are shown in the bottom panels of Fig. 2 and Figs B1–B5. In order to be conservative in the subsequent global analysis detailed below and preserve the uncertainty of the decorrelation fits, the errors on the linear models were estimated from one thousand samples drawn from the posterior distributions of the regression coefficients and added in quadrature with the flux errors of the detrended photometry for each data set. These decorrelated light curves yield 3 h noise estimates of 58.4, 51.9, 49.5, 53.9, 45.8, and 53.0 ppm, respectively, and were used in our joint analysis below.

### 4.2.2 CCF

Following the stellar variability analysis (see Section 3.3), it was noted that the TESS-derived 12.2 d orbital period for TOI-1064 c is close to the 13.3 d first harmonic of the stellar rotation period. As stellar activity signals can impede the detection of exoplanets in RV data (Haywood et al. 2014; Rajpaul et al. 2015; Mortier et al. 2016; Dumusque et al. 2017), we employed the cross-correlation function (CCF)-based SCALPELS (Collier Cameron et al. 2021) to separate any stellar activity signals from planet-induced RV variations via the re-extraction of RVs from the 26 HARPS observations (see Section 2.6). For this data set, the LOOCV approach selected three principal components (12 per cent) that well model the CCF shape variations and hence stellar activity.

As reported in Collier Cameron et al. (2021), by conducting a joint fit of the stellar activity and planetary RVs, orbital signals injected into solar RV data can be more reliably retrieved using SCALPELS. Therefore, as a preliminary study of the HARPS data, we performed simultaneous detrending using the CCF-based SCALPELS method and a fit to two Keplerian orbits utilizing the approach presented in the appendix of the aforementioned paper. Using the TESS-derived values as priors and assuming circular orbits, we detected TOI-1064 b with a semi-amplitude, $K_b = 4.7 \pm 0.5$ m s$^{-1}$, and tentatively found TOI-1064 c with a semi-amplitude, $K_c = 1.8 \pm 0.6$ m s$^{-1}$. The median-subtracted raw RV time series, along with the SCALPELS-determined stellar activity signal, the resulting planetary RV component, and the orbital fit and residuals, is shown in Fig. 7.

The SCALPELS produced vectors represent variations in observed RVs due to stellar activity, and so can be used to potentially remove such variation via a linear regression. Thus, the three identified components were subsequently used to detrend the HARPS RVs in our global analysis simultaneously with the RV fitting.
## 4.3 Joint photometric and radial velocity analysis

In order to determine the properties of TOI-1064 b and TOI-1064 c, we conducted a global analysis of the TESS, CHEOPS, LCOGT, NGTS, and ASTEP transit photometry, and the HARPS RV data. Prior to the joint fitting, we carried out additional preliminary checks to further assess the data and ascertain a complete model for use in the analysis.

### 4.3.1 Preliminary analyses

As can be seen from Fig. 2 and Figs B1–B5, the linear models produced by the PSF-based SCALPELS provide an excellent reproduction of the flux variation seen in the CHEOPS light curves, and indicate that the flux modulation typically seen on CHEOPS orbit time-scales can be well accounted for via modelling of the PSF shape changes alone. To check whether additional detrended basis vectors were needed to model CHEOPS roll angle flux variation, we used the PYCHEOPS<fnref n="5" /> PYTHON package (Maxted et al. 2022) to evaluate the detrended data produced by the method outlined above. For each visit, we performed simultaneous transit fitting and detrending for all combinations of standard basis vectors used in the decorrelation of CHEOPS data (i.e. background, contamination, smear, x and y centroid positions, and first-, second-, and third-order harmonics of the roll angle), and assessed the reported Bayes factors for the considered basis vectors. We found that for all visits, the fitting favoured not including additional detrending parameters, and thus for the joint fit no further basis vectors are used.

We also performed a preliminary analysis of the HARPS data utilizing the RV fitting capabilities of RADVEL within JULIET (Fulton et al. 2018; Espinoza et al. 2019) in order to determine whether there are any long-term trends seen in the data that should be accounted for. Using the stellar and planetary priors outlined below, we conducted two fits with the inclusion or not of a RV slope and intercept. By comparing the nested sampling-produced Bayes evidences, we found that excluding a long-term RV trend is favoured, and therefore, it is omitted from the global analysis below.

Lastly, we conducted a stability analysis of the system to determine dynamically plausible regions of the eccentricity and argument of periastron planes for both planets that can be used to constrain priors on these parameters in the global fit. We consider an architecture unstable if the trajectories of planets would result in intersecting orbits calculated as satisfying the following:

$$ \left| \frac{p_1 - p_2}{\sqrt{(p_1 e_2)^2 + (p_2 e_1)^2 - 2p_1 p_2 e_1 e_2 \cos(\omega_1 - \omega_2)}} \right| < 1, \quad (1) $$

where $p_j = a_j(1 - e_j)$ for $j = 1, 2$. Using this formulation, we perform 100 000 calculations randomly sampling the $\sqrt{e\cos\omega}$ and $\sqrt{e\sin\omega}$ parameter space between -1 and 1 for both planets, ensuring that the eccentricities remain below 1, whilst fixing the nominal values for all remaining parameters. In total, we find that 23.2 per cent of scenarios remain stable with the remaining architectures resulting in orbit crossing events that occur more frequently with larger values of the eccentricity and argument of periastron components. Therefore, in order to provide bounds for the eccentricity and argument of periastron priors and to avoid the subsequent global fit from returning unstable results, cut-offs for these priors are needed. In our simulations we find that bounds of $\pm 0.5$ for both components on both planets result in 96.8 per cent of the returned orbits being stable. Thus, this constraint provides a good compromise and assurance

that the fitted eccentricities and arguments of periastron yield a dynamically stable scenario and are used in the joint fit below.

### 4.3.2 Global fit

For the global analysis of the TOI-1064 system, we used the JULIET package (Espinoza et al. 2019) that employs the BATMAN code (Kreidberg 2015) for transit photometry fitting and RADVEL (Fulton et al. 2018) for the modelling of RVs. To explore the parameter space of the fit and estimate the Bayesian posteriors, we used nested sampling algorithms provided in the DYNESTY package (Speagle 2020) that computes the Bayesian evidences allowing for robust model comparison. The parametrization of the fit, and the priors used, is listed in Table D1 and outlined below.

(i) The orbital period, $P$, and the mid-transit time, $T_0$, for both planets were taken as uniform priors with bounds equal to the preliminary TESS values minus and plus three times the corresponding uncertainty.

(ii) The ($r_1, r_2$) parametrization as introduced by Espinoza (2018) was used, as it permits the physically plausible area of the ($b, p$) parameter space to be efficiently investigated through the introduced random variates $r_1$ and $r_2$, where $b$ is the transit impact parameter and $p$ is the planet-to-star radius ratio. We set uniform priors on both $r_1$ and $r_2$ of $U(0, 1)$ for both planets.

(iii) The $\sqrt{e\cos\omega}$ and $\sqrt{e\sin\omega}$ parametrization (Eastman, Gaudi & Agol 2013), where $e$ is the eccentricity and $\omega$ is the argument of periastron, were used, with uniform priors on $\sqrt{e\cos\omega}$ and $\sqrt{e\sin\omega}$ with bounds of $U(-0.5, 0.5)$ taken from the preliminary analysis detailed above on both components for both planets.

(iv) The RV semi-amplitude, $K$, for both planets was taken as uniform priors, $U(0, 20)$ m s⁻¹.

(v) The stellar density, $\rho_*$, parametrization was used as, when combined with the orbital periods, it provides scaled semimajor axes for both planets that are anchored to a single common value rather than setting priors on separate scaled semimajor axes. We set a normal prior on $\rho_*$ using the value from Table 3.

(vi) The quadratic limb-darkening coefficients parametrized in the ($q_1, q_2$) plane (Kipping 2013) were used, with the coefficients for each bandpass (TESS, CHEOPS, LCOGT PanSTARRS $z_s$ and Bessel $B$, NGTS, and ASTEP Cousins $R$) calculated using the 1d module of PYCHEOPS (Maxted et al. 2022). Normal priors were set for all limb-darkening coefficients centred on the determined values with a $1\sigma$ uncertainty of 0.1 taken.

In addition to the modelling of planetary signals in transit photometry and RV data, and SCALPELS basis vectors to simulate the stellar activity in the RVs, we included further parameters to model instrumental and astrophysical noise in the form of GPs. For all transit photometry data sets, we used GPs with a Matérn-3/2 kernel against time in order to model long-term correlated noise, utilizing the CELERITE package (Foreman-Mackey et al. 2017) within JULIET. An example of a fitted GP to the TESS fluxes can be seen in Fig. 1. Moreover, for both photometry and RV data we fitted jitter terms in order to reflect any extra noise not previously accounted for, which is subsequently summed in quadrature with the uncertainties on the data. Lastly, we fitted the HARPS RVs with a zero-point RV offset, $\gamma_*$.

### 4.3.3 Results

We report the fitted results of our global analysis for both planets in the TOI-1064 system in Table 4. The detrended and phase-folded TESS and CHEOPS photometry along with the fitted transit models

<footnote n="5">⁵https://github.com/pmaxted/pycheops</footnote>
<table><caption>**Table 4.** Fitted and derived parameter values for TOI-1064 b and TOI-1064 c based on the joint fit to the photometric and RV data detailed in Section 4.3.</caption><thead><tr><th>Parameter (unit)</th><th>TOI-1064 b</th><th>TOI-1064 c</th></tr></thead><tbody><tr><td colspan="3" style="text-align: center;">Fitted parameters</td></tr><tr><td>P (d)</td><td>6.443868 ± 0.000025</td><td>12.226574<sup>+0.000046</sup><sub>-0.000043</sub></td></tr><tr><td>T<sub>0</sub> (BJD - 245 7000)</td><td>2036.85340<sup>+0.000060</sup><sub>-0.000053</sub></td><td>2043.51289<sup>+0.000069</sup><sub>-0.000067</sub></td></tr><tr><td>r<sub>1</sub><sup>a</sup></td><td>0.818<sup>+0.020</sup><sub>-0.015</sub></td><td>0.835<sup>+0.016</sup><sub>-0.017</sub></td></tr><tr><td>r<sub>2</sub><sup>a</sup></td><td>0.03267<sup>+0.00044</sup><sub>-0.00041</sub></td><td>0.03347<sup>+0.00042</sup><sub>-0.00041</sub></td></tr><tr><td>√e cosω</td><td>-0.01<sup>+0.08</sup><sub>-0.09</sub></td><td>0.25<sup>+0.13</sup><sub>-0.14</sub></td></tr><tr><td>√e sinω</td><td>-0.20<sup>+0.11</sup><sub>-0.08</sub></td><td>0.06 ± 0.14</td></tr><tr><td>K (m s<sup>-1</sup>)</td><td>5.62<sup>+0.67</sup><sub>-0.75</sub></td><td>0.85<sup>+0.59</sup><sub>-0.67</sub></td></tr><tr><td>γ<sub>*</sub> (m s<sup>-1</sup>)</td><td></td><td>21 214.5 ± 0.7</td></tr><tr><td>ρ<sub>*</sub> (kg m<sup>-3</sup>)</td><td></td><td>2711<sup>+75</sup><sub>-60</sub></td></tr><tr><td>ρ<sub>*</sub> (ρ<sub>☉</sub>)</td><td></td><td>1.92<sup>+0.05</sup><sub>-0.04</sub></td></tr><tr><td colspan="3" style="text-align: center;">Derived parameters</td></tr><tr><td>δ<sub>tr</sub> (ppm)</td><td>1067<sup>+29</sup><sub>-27</sub></td><td>1120<sup>+28</sup><sub>-27</sub></td></tr><tr><td>R<sub>p</sub>/R<sub>*</sub></td><td>0.03267<sup>+0.00044</sup><sub>-0.00041</sub></td><td>0.03347<sup>+0.00042</sup><sub>-0.00041</sub></td></tr><tr><td>R<sub>*</sub>/a</td><td>0.05488<sup>+0.00097</sup><sub>-0.00092</sub></td><td>0.03581<sup>+0.00063</sup><sub>-0.00060</sub></td></tr><tr><td>R<sub>p</sub>/a</td><td>0.001793<sup>+0.000040</sup><sub>-0.000038</sub></td><td>0.001199<sup>+0.000026</sup><sub>-0.000025</sub></td></tr><tr><td>R<sub>p</sub> (R<sub>☉</sub>)</td><td>2.587<sup>+0.043</sup><sub>-0.042</sub></td><td>2.651<sup>+0.043</sup><sub>-0.042</sub></td></tr><tr><td>a (au)</td><td>0.06152<sup>+0.00086</sup><sub>-0.00089</sub></td><td>0.09429<sup>+0.00132</sup><sub>-0.00136</sub></td></tr><tr><td>t<sub>14</sub> (h)</td><td>1.98<sup>+0.07</sup><sub>-0.08</sub></td><td>2.37 ± 0.10</td></tr><tr><td>b</td><td>0.728<sup>+0.029</sup><sub>-0.023</sub></td><td>0.753<sup>+0.024</sup><sub>-0.026</sub></td></tr><tr><td>i (°)</td><td>87.709<sup>+0.083</sup><sub>-0.097</sub></td><td>88.455<sup>+0.058</sup><sub>-0.056</sub></td></tr><tr><td>e</td><td>0.047<sup>+0.038</sup><sub>-0.030</sub></td><td>0.088<sup>+0.081</sup><sub>-0.064</sub></td></tr><tr><td>ω (°)</td><td>120<sup>+49</sup><sub>-89</sub></td><td>25<sup>+22</sup><sub>-16</sub></td></tr><tr><td>S<sub>p</sub> (S<sub>⊕</sub>)</td><td>62.9<sup>+4.2</sup><sub>-4.0</sub></td><td>26.8<sup>+1.8</sup><sub>-1.7</sub></td></tr><tr><td>T<sub>eq</sub> (K)</td><td>784 ± 13</td><td>634 ± 10</td></tr><tr><td>M<sub>p</sub> (M<sub>☉</sub>)</td><td>13.5<sup>+1.7</sup><sub>-1.8</sub></td><td>2.5<sup>+1.8</sup><sub>-2.0</sub> (<8.5)<sup>b</sup></td></tr><tr><td>ρ<sub>p</sub> (ρ<sub>⊕</sub>)</td><td>0.78<sup>+0.10</sup><sub>-0.11</sub></td><td>0.14<sup>+0.09</sup><sub>-0.11</sub> (<4.6)<sup>b</sup></td></tr><tr><td>g<sub>p</sub> (m s<sup>-2</sup>)</td><td>19.7<sup>+2.5</sup><sub>-2.7</sub></td><td>3.5<sup>+2.4</sup><sub>-2.8</sub> (<11.8)<sup>b</sup></td></tr></tbody></table>

for planets b and c are shown Fig. 8, whilst the detrended LCOGT, NGTS, and ASTEP data, and fitted transit models of both planets are reported in Figs 9 and 10. The SCALPELS-corrected HARPS RV time series with the fit of two Keplerian orbits is presented in Fig. 11, and the orbital period phase-folded RVs and the planetary models for TOI-1064 b and TOI-1064 c are shown in Fig. 11. The posterior distributions for the main fitted parameters for both planets are given in Fig. D1 with the posterior values for the fitted limb-darkening coefficients and noise terms reported in Table D2.

From our analysis, we detect TOI-1064 b and TOI-1064 c in the combined TESS, CHEOPS, LCOGT, NGTS, and ASTEP transit photometry at 38.2σ and 40.3σ, respectively. The fitted depths and derived stellar radius yield planetary radii of $R_{p, b} = 2.587_{-0.042}^{+0.043}$ R<sub>⊕</sub> and $R_{p, c} = 2.651_{-0.042}^{+0.043}$ R<sub>⊕</sub>. We report a 8.0σ detection of TOI-1064 b in the HARPS RVs, and a 1.4σ signal of TOI-1064 c, that results in $M_{p, b} = 13.5_{-1.8}^{+1.7}$ M<sub>⊕</sub> and a $M_{p, c}$ 3σ upper limit of 8.5 M<sub>⊕</sub> (nominal value = $2.5_{-2.0}^{+1.8}$ M<sub>⊕</sub>, to be discussed in Section 4.3.4).

Combining the radius and mass of TOI-1064 b gives a bulk density of $\rho_{p, b} = 0.78_{-0.11}^{+0.10}$ $\rho_{\oplus}$ (4.28<sup>+0.57</sup><sub>-0.61</sub> g cm<sup>-3</sup>). We determine the orbital

period of planet TOI-1064 b to be $P_b = 6.443868 \pm 0.000025$ d, and the corresponding semimajor axis of $a_b = 0.06152_{-0.00089}^{+0.00086}$ au. At this distance TOI-1064 b receives $62.9_{-4.0}^{+4.2}$ times Earth's insolation and has a zero Bond albedo equilibrium temperature of 784 ± 13 K. For TOI-1064 c, the subtle RV signal together with the precise radius yields a $\rho_{p, c}$ 3σ upper limit of 2.56 g cm<sup>-3</sup> (0.46 $\rho_{\oplus}$; nominal value = $0.14_{-0.11}^{+0.09}$ $\rho_{\oplus}$). The fitted orbital period of $P_c = 12.226573_{-0.00043}^{+0.00046}$ d gives a semimajor axis of $a_c = 0.09429_{-0.00134}^{+0.00132}$ au, which results in the stellar irradiance of TOI-1064 c being $26.8_{-1.7}^{+1.8}$ that of Earth with a zero Bond albedo equilibrium temperature of 634 ± 10 K.

Placing TOI-1064 b and TOI-1064 c on a mass-radius diagram, as shown in Fig. 12 alongside well-characterized planets with masses and radii known to better than 20 per cent, we see that planet TOI-1064 b is one of the smallest planets known with a mass above 10 M<sub>⊕</sub>. Conversely, the tentative signal of planet TOI-1064 c places it at the lower end of the mass distribution for planets of this size. However, further RV follow-up observations are needed to confirm the mass precisely and to determine if TOI-1064 c remains one of the lowest density sub-Neptunes known.

Given the refined nature of the radii of the two planets, with uncertainties of 1.63 per cent and 1.59 per cent, we are able to characterize the planets further and conduct internal structure and atmospheric escape modelling of the bodies, as detailed below. Additionally, we find the timing uncertainties for individual transits of TOI-1064 b to be between 2 and 10 min, and between 4 and 13 min for TOI-1064 c.

In addition to the information about the planets in the TOI-1064 system, from our analysis we determined the density and RV of the host star. We find a stellar density of $\rho_{\star} = 1.92_{-0.04}^{+0.05}$ $\rho_{\odot}$, which is in agreement with the value derived from IRFM and isochrone placement methods reported in Table 3. From the HARPS data, we obtain a $\gamma_{\star}$ value of 21.2 km s<sup>-1</sup> that agrees with the Gaia value of 20.7 km s<sup>-1</sup> (Gaia Collaboration 2021). To assess the robustness of the fitted eccentricities we conducted two further analyses of the data. First, we carry out a fit of the complete transit photometry and RV data set setting the eccentricities of both planets to 0 whilst using the same priors as detailed above for the remaining parameters. We find a difference in log evidences between this fit and our global analysis of $\Delta \ln Z = 17.8$ in favour of the non-zero eccentricity solution, indicating a decisive preference for this result (Kass & Raftery 1995; Gordon & Trotta 2007). Second, we fit the transit photometry using the same priors as previously set on the appropriate parameters to evaluate the contribution of the RVs. We find that the derived eccentricities from the transit photometry alone ($e_b = 0.38_{-0.20}^{+0.05}$ and $e_c = 0.33 \pm 0.07$) agree within 2σ with the lower eccentricities listed in Table 4 that were derived from the RV and photometric data. Because of the compact nature of the system, we conclude that the lower, non-zero eccentricity solution from our global analysis is likely correct.

### 4.3.4 RV model comparison

To give confidence in the subtle RV signal of TOI-1064 c found by the global analysis, we conduct two additional fits of the complete transit photometry and RV data set. In these supplementary analyses we set the priors of the majority of the parameters as detailed above, and either remove planet TOI-1064 c from the RV model or fix the semi-amplitude of TOI-1064 c to 0. By comparing the log Bayesian evidences between these analyses and our previous fit, we find $\Delta \ln Zs$ of 55.6 and 26.1, respectively, in favour of our

<footnote><sup>a</sup>The r<sub>1</sub> and r<sub>2</sub> fitting parameters used in the JULIET package (Espinoza et al. 2019) are the result of a physically plausible parametrization of the (b, p) planet as detailed in Espinoza (2018).</footnote>

<footnote><sup>b</sup>3σ upper limit.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Top panels: detrended TESS (left) and CHEOPS (right) photometry phase-folded to the orbital periods of planets b and c (offset by -0.015 for clarity). Unbinned data are shown in grey, with fluxes binned every 15 min given as white circles. The best-fitting transit models for TOI-1064 b and TOI-1064 c are presented in green and purple, respectively. Bottom panels: residuals to the fits shown in the top panels with the same colouration.</figcaption></figure>

global analysis presented in Section 4.3.2. This indicates decisive preferences of using a non-massless prior to model the RV signal of TOI-1064 c (Kass & Raftery 1995; Gordon & Trotta 2007), and provides evidence for the validity of the tenuous signal seen in the HARPS RVs.

It should be noted that the semi-amplitudes and masses for both planets derived from the global analysis are subtly different to the values retrieved by the CCF-based SCALPELS algorithm. As presented in the Sections 4.2 and 4.3.3, TOI-1064 b and TOI-1064 c are more and less massive by 0.4σ and 0.8σ, respectively, in the nested sampling joint fit. Importantly, the SCALPELS-fitted semi-amplitude uncertainties are lower that allows for a more confident detection of planet TOI-1064 c in the HARPS RVs. As these differences can have a significant effect on the bulk density, and internal structure and

atmospheric escape modelling, we conducted two fits of the RV data using RADVEL within JULIET (Fulton et al. 2018; Espinoza et al. 2019) and the nested sampling algorithms in the DYNESTY package (Speagle 2020) to reconcile this difference. For these analyses, we set the semi-amplitude priors to the values and uncertainties produced by either the SCALPELS or our global fit, with wide uniform priors on P, T₀, √eicosω, and √esinω identical to the values set out in Section 4.3.2.

To assess if one model (priors taken from the fit produced by the SCALPELS) is preferred over the other (priors taken from our global analysis), we consult the log Bayesian evidences reported by the nested sampling and compute an odds ratio between the two models. We find a difference in log evidences of Δ lnZ = 0.49 and thus an odds ratio of 1.63 in favour of the model produced by our global fit. Following the standard reference levels (Kass & Raftery
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Top panel: detrended LCOGT (top) and ASTEP (bottom) photometry of TOI-1064 b (offset by -0.01 for clarity). Unbinned data are shown in grey, with fluxes binned every 15 min given as white circles. The best-fitting transit model for TOI-1064 b is presented in green. Bottom panel: residuals to the fits shown in the top panel with the same colouration.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Top panel: detrended LCOGT (top and upper middle), NGTS (lower middle), and ASTEP (bottom) photometry of planet TOI-1064 c (offset by -0.01 for clarity). Unbinned data are shown in grey, with fluxes binned every 15 min given as white circles. The best-fitting transit model for TOI-1064 c is presented in purple. Bottom panel: residuals to the fits shown in the top panel with the same colouration.</figcaption></figure>

1995; Gordon & Trotta 2007), this translates to a weak and non-conclusive preference. However, given the more comprehensive data sets analysed in our global fit, we adopt the values presented in Table 4 as the nominal values for this system.

# 5 CHARACTERIZATION OF THE SYSTEM

## 5.1 Internal structure

We analysed the internal structure of the two planets in the TOI-1064 system using the method employed by Leleu et al. (2021) for TOI-178. The method is based on a global Bayesian model that fits the observed properties of the star (mass, radius, age, effective temperature, and the photospheric abundances [Si/Fe] and [Mg/Fe]) and planets (planet-star radius ratio, the RV semi-amplitude, and the orbital period).

In terms of the forward model, we assume a fully differentiated planet, consisting of a core composed of Fe and S, a mantle composed of Si, Mg, Fe, and O, a pure water layer, and an H and He layer. The temperature profile is adiabatic, and the equations of state (EoS) used for these calculations are taken from Hakim et al. (2018) and Fei et al. (2016) for the core materials, from Sotin et al. (2007) for the mantle materials, and from Haldemann et al. (2020) for water. The thickness of the gas envelope is computed using the semi-analytical model of Lopez & Fortney (2014).

In our model, we assume the following priors: the logarithms of the gas-to-solid ratios in planets have uniform distributions; the mass fractions of the planetary cores, mantles, and overlying water layers have uniform positive priors except that the mass fractions of water are limited to a maximum value of 0.5. The bulk Si/Fe and Mg/Fe mole ratios in the planet are assumed to be equal to the values deter-
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Top panels: the SCALPELS-corrected, offset-subtracted, HARPS RV time series plotted in blue with the jitter term added in quadrature to the RV uncertainty in the error bars. The best-fitting model is plotted in black with the model uncertainty shown in grey, and the residuals to the solution are given in the lower subpanel. Middle panels: the RVs phase folded to the orbital period of TOI-1064 b with the best-fitting model and corresponding uncertainty plotted in green and light green, respectively, and the residuals to the fit given in the lower subpanel. Lower panels: same as the middle panels, but for planet TOI-1064 c. All data shown include the jitter noise term added in quadrature with the RV error.</figcaption></figure>

mined for the atmosphere of the star, given in Table 3. Recent work by Adibekyan et al. (2021) has, however, found that whilst the abundances of planets and host stars are correlated, the relation may not be one-to-one. Finally, we note that the solid and gas parts of the planets are computed independently, which means that we do not include the compression effect of the planetary envelope on its core. This is justified given the small masses of the gas envelopes (see below). The results of our modelling for both planets are shown in Fig. 13.

## 5.2 Atmospheric evolution

We modelled the atmospheric evolution of the two planets using a modified version of the algorithm presented by Kubyshkina et al. (2019a,b). This Bayesian tool requires input parameters that are both stellar ($M_*$, $t_*$, and the present-day rotation period $P_{rot,*}$) and planetary (semimajor axes $a$ and masses $M_p$). The stellar rotational period $P_{rot}$ is assumed as a proxy for the stellar high-energy emission, which plays a significant role in controlling the atmospheric escape rate, and is modelled over time as a broken power law with a variable exponent within the first 2 Gyr (Tu et al. 2015), and afterwards $P_{rot} \propto t^{0.566}$ following Mamajek & Hillenbrand (2008). We translated $P_{rot}$ into the stellar X-ray and extreme ultraviolet (XUV) luminosities using the scaling relations in Wright et al. (2011), Sanz-Forcada et al. (2011), and McDonald, Kreidberg & Lopez (2019). Lastly, the remaining input parameters include the planetary equilibrium temperature $T_{eq}$, radius $R_p$, mass $M_p$, atmospheric mass $M_{atm}$, and orbital semimajor axis $a$, with the evolution of $T_{eq}$ over time due to changes in stellar bolometric luminosity $L_{bol,*}$ calculated by interpolating within grids of stellar evolutionary models (Choi et al. 2016; Dotter 2016).

The main model hypotheses are that orbital semimajor axes are fixed to the present-day value and that planetary atmospheres are hydrogen dominated. The free parameters of the algorithm are the exponent of the stellar rotation period evolution power law and the initial atmospheric mass fractions of the planets $f_{at}^{start}$. Millions of planetary evolutionary tracks are generated within a Bayesian context employing an MCMC scheme (Cubillos et al. 2017). After rejecting those that do not fulfil the constraint imposed by the present-day atmospheric content derived in Section 5.1, we obtained the posterior distributions of the free parameters. Further details about the tool may be found in Delrez et al. (2021). The results of the joint evolution of both planet TOI-1064 b and TOI-1064 c are shown in Figs E1 and E2.

## 6 DISCUSSION

From our analysis, we clearly detect TOI-1064 b in the *TESS*, *CHEOPS*, LCOGT, and ASTEP photometry and the HARPS RVs. TOI-1064 c is confidently detected in the *TESS*, *CHEOPS*, LCOGT, NGTS, and ASTEP photometry, but does not register a significant signal in the HARPS RV data. The photometric periods and HARPS RVs secure a mass of $13.5 \pm 2.0 M_⊕$ for planet TOI-1064 b, and a $3\sigma$ upper limit of $8.5 M_⊕$ for planet TOI-1064 c. The proximity of the orbital period to the first harmonic of the stellar rotation period hinders a more confident detection in the RVs. However, with more data over an extended baseline SCALPELS should be able to separate the two signals due to a more apparent shift in the CCFs, whilst periodograms may also be able to find two peaks due to a higher frequency resolution.

### 6.1 Mass-radius comparisons and bulk density-metallicity correlation

Figs 12 and 14 highlight that, although the radii of TOI-1064 b and TOI-1064 c are similar (in agreement with the ‘peas in a pod’ scenario; Weiss et al. 2018), the current mass values are likely significantly different and thus, the bulk planetary densities are considerably distinct. TOI-1064 b is one of the smallest, and therefore densest, well-characterized sub-Neptunes with a mass greater $10 M_⊕$ known. This places planet TOI-1064 b amongst a small family of dense ($4.0–4.3 \text{ g cm}^{-3}$), warm ($690–810$ K), and mildly irradiated ($32–71 S_⊕$) sub-Neptunes orbiting around K-dwarf stars that also
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Mass–radius diagram for exoplanets with mass and radius measurements more precise than 20 per cent coloured by the zero Bond albedo equilibrium temperature with data taken from TEPCat (Southworth 2011). TOI-1064 b and TOI-1064 c are highlighted via black bordered circles, with other planets discussed in Section 6 shown as stars, diamonds, and squares. The black solid line represents the upper bound of the forbidden region according to collisional mantle stripping via impact (Marcus et al. 2010). The remaining solid lines depict theoretical mass–radius curves for solid planets with various bulk compositions (Zeng & Sasselov 2013), and the dashed lines represent mass–radius curves for an Earth-like core enveloped with a 1 per cent mass fraction H₂ atmosphere with various equilibrium temperatures (Zeng et al. 2019). Venus, Earth, Uranus, and Neptune are indicated in blue.</figcaption></figure>

includes HIP 116454 b (Vanderburg et al. 2015) and Kepler-48 c (Steffen et al. 2013; Marcy et al. 2014). The other planets in this parameter space orbit stars that are more metal poor ([M/H] = −0.16 ± 0.08) and metal rich ([Fe/H] = 0.17 ± 0.07) than TOI-1064 ([Fe/H] = 0.05 ± 0.08), respectively. There are two further planets in this region of the mass–radius diagram that both are similarly warm and irradiated, but are more dense (5.0–5.6 g cm⁻³); K2-110 b (Osborn et al. 2017) and K2-180 b (Korth et al. 2019). Interestingly they both orbit metal-poor K-dwarfs with [Fe/H] = −0.34 ± 0.03 and −0.65 ± 0.10 for K2-110 b and K2-180 b, respectively. These four planets are highlighted in Figs 12 and 14 as black bordered stars. This may hint that metallicity could affect the bulk density of sub-Neptunes in this parameter space via differing formation conditions.

To test this apparent correlation for massive sub-Neptunes, we select all planets with masses greater than 10 M⊕ in our well-characterized sample (uncertainties on the radii and masses less than 20 per cent) and compared their bulk densities against host star metallicities as seen in Fig. 15. Whilst large uncertainties on the bulk density of some exoplanets, especially at larger values, results in a scatter to the data there appear to be a negative correlation. We quantify this trend for our sample using a Bayesian correlation tool, that aims characterize the strength of the correlation between

two parameters in a Bayesian framework (Figueira et al. 2016), and reports a correlation distribution that represent values akin to a Spearman’s rank value. For this sample, we find the peak of the correlation posterior distribution to be −0.26, with 95 per cent lower and upper bounds of −0.60 and 0.11. Interestingly, if we exclude highly irradiated planets (S > 200 S⊕) we retrieve a peak density–metallicity correlation of −0.45, with 95 per cent lower and upper bounds of −0.78 and −0.09. If this trend of less dense sub-Neptunes orbiting metal-rich stars is indicative of formation conditions it could suggest that metal-rich stars form planets that accrete bigger envelopes or have more massive cores and are able to retain their atmospheres (Owen & Murray-Clay 2018). This correlation could also imply that planets around metal-rich stars have metal-rich atmospheres that have reduced photoevaporation-driven atmospheric mass-loss rates (Owen & Jackson 2012). This correlation is similar to that seen in giant planets in which more massive planets were found around metal-rich stars (Guillot et al. 2006). This formation picture could be supported by the strengthening of the correlation with the removal of highly irradiated planets that may have an observed high density due to evolution processes such as atmospheric stripping.

Thus, as there are overlapping mass–radius curves in this region, accurate characterization of additional dense sub-Neptunes
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Corner plot showing the main internal structure parameters of TOI-1064 b (left) and TOI-1064 c (right). Shown are the mass fraction of the inner core, the mass fraction of water, the Si and Mg mole fractions in the mantle, the Fe mole fraction in the inner core, and the mass of gas (log scale). The values on top of each column are the mean and 5 per cent and 95 per cent quantiles. Next to each corner plot, we also show illustrations of both planets representing the radius fractions of the high-Z components (inner core and mantle; brown), water layer (blue), and gas envelope (grey), corresponding to the medians of the posterior distributions, with the dashed line outer rings representing the uncertainty on the total radius.</figcaption></figure>

over a range of metallicities should help clarify this trend and inform whether these planets are water dominated or terrestrial with a H₂ atmosphere, and indicate the processes that sculpt these planets.

We note that the similarly sized and irradiated planet HD 119130 b is denser and orbits a more massive solar-like star than the dense sub-Neptunes discussed above (Luque et al. 2019). Lastly, two planets (Kepler-411 b and K2-66 b; Wang et al. 2014; Crossfield et al. 2016; Sinukoff et al. 2017; Sun et al. 2019) have similar radii to TOI-1064 b, but larger masses and substantially higher irradiation (>200 S<sub>⊕</sub>) meaning that these ultradense bodies are likely separate from the warm, dense sub-Neptunes mentioned previously that include TOI-1064 b. However, these objects may represent the lower mass end of a stripped planet core family that includes the more massive, but equally dense TOI-849 b (Armstrong et al. 2020). These three planets are highlighted in Figs 12 and 14 as black bordered diamonds.

From the substantial constraint on the radius of TOI-1064 c there exists a band of possible bulk densities from the large uncertainty on the mass. Taken at the nominal value, planet TOI-1064 c would be one of the least dense sub-Neptunes known and could represent a planet akin to Kepler-307 c (Xie 2014; Jontof-Hutter et al. 2016), albeit with a more extended atmosphere. If this is the case, the TOI-1064 planets would resemble TOI-178 c and TOI-178 d both in density difference and orbital period ratio (Leleu et al. 2021). However, if the mass is confirmed to be nearer the upper limit, the low density and mild irradiation would place TOI-1064 c closer to the outer planets of Kepler-80 (Kepler-80 b and Kepler-80 c; Xie 2014; MacDonald et al. 2016) or GJ 1214 b (Charbonneau et al. 2009), and in between planets of similar density and higher (HIP 41378 b and TOI-125 c; Vanderburg et al. 2016; Quinn et al. 2019; Santerne et al. 2019; Nielsen et al. 2020) and lower irradiation (Kepler-26 b and Kepler-

26 c, TOI-270 c, and LTT 3780 c; Steffen et al. 2012; Jontof-Hutter et al. 2016; Günther et al. 2019; Cloutier et al. 2020; Nowak et al. 2020; Van Eylen et al. 2021). These eight planets are highlighted in Figs 12 and 14 as black bordered squares.

## 6.2 Internal structure and atmospheric escape

The internal structure models show that TOI-1064 b has a small gas envelope (log M<sub>gas</sub> between −9.2 and −1.0 – all given values are the 5 per cent or 95 per cent quantiles). It has a relatively small inner core and a large mass fraction of water (comprised between 0.01 and 0.26, and 0.06 and 0.49 of the mass of the solid body, respectively), with the remainder of the mass in the Si- and Mg-dominant mantle. It should be noted that in our atmospheric escape modelling we assume the presence of an accreted primordial atmosphere. However, the internal structure modelling results highlight that the tenuous atmosphere of TOI-1064 b means that the mass of the high-Z elements (i.e. the inner core, mantle, and volatile layers) is essentially equal to the total mass determined from the RV fitting. This value, and other modelled properties, could thus be similar to the properties of the high-Z core prior to the accretion of an atmosphere and may inform formation models (Pollack et al. 1996; Ikoma & Hori 2012). Further high-precision photometry and internal structure modelling of additional ultradense sub-Neptunes (such as Kepler-411 b and K2-66 b; Wang et al. 2014; Crossfield et al. 2016; Sinukoff et al. 2017; Sun et al. 2019) may help shed light on the formation of these bodies.

In contrast, planet TOI-1064 c likely has a substantial gas envelope (log M<sub>gas</sub> between −2.5 and −1.7, corresponding to mass fractions of up to 3 per cent). The presence of a non-negligible gas envelope results in the mass fraction of water being unconstrained as seen in Fig. 13.
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Mass–density diagram for exoplanets with mass and radius measurements more precise than 20 per cent coloured by the host star metallicity with data taken from TEPCat (Southworth 2011). TOI-1064 b and TOI-1064 c are highlighted via black bordered circles, with other planets discussed in Section 6 shown as stars, diamonds, and squares. The black solid line represents the lower bound of the forbidden region according to collisional mantle stripping via impact (Marcus et al. 2010). The remaining solid lines depict theoretical mass–density curves for solid planets with various bulk compositions (Zeng & Sasselov 2013) and the dashed lines represent mass–density curves for an Earth-like core enveloped with a 1 per cent mass fraction H₂ atmosphere with various equilibrium temperatures (Zeng et al. 2019). Venus, Earth, Uranus, and Neptune are indicated in blue.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 15.** The bulk density versus stellar metallicity for planets with masses above 10 M<sub>⊕</sub> and radius and mass uncertainties less than 20 per cent, shown in Figs 12 and 14. Black circles indicate highly irradiated planets (S > 200 S<sub>⊕</sub>), with green circles showing mildly irradiated bodies. TOI-1064 b is highlighted by a black bordered star, with Uranus and Neptune indicated in blue.</figcaption></figure>

The results of the atmospheric escape modelling based on fitting simultaneously both planets, respect the priors well (see Figs E1 and E2). These results suggest that the rotation of TOI-1064 evolved as observed in the majority of late-type stars of similar mass. We also find that planet TOI-1064 b has lost (almost) the entirety of its H/He envelope, assuming it accreted one, at some time in the past, as indicated by the flat posterior on the initial atmospheric mass fraction (Fig. E2, left-hand column). As a consequence of the main assumptions of the atmospheric evolution algorithm (namely that planets are assumed to have formed and spent their whole evolution at the same distance from the star), planet TOI-1064 b is considered to have migrated to its current position while still embedded in the protoplanetary disc, as is thought to have happened to the Kepler-223 (Mills et al. 2016), TRAPPIST-1 (Luger et al. 2017), and TOI-178 planetary systems (Leleu et al. 2021). This migration scenario is favoured with respect to *in situ* formation as, to form TOI-1064 b at 0.06 au, the protoplanetary disc mass must have been substantially enhanced compared to solar values (Schlichting 2014) and the presence of a substantial volatile layer, as determined by our internal structure modelling, that must have been accreted beyond the ice line. To test this migration scenario, atmospheric studies of TOI-1064 b could be undertaken to determine if the tenuous envelope is a
secondary, high molecular weight atmosphere (Fortney et al. 2013), similar to that inferred at $\pi$ Men c (García Muñoz et al. 2021), that could be composed of steam (Marounina & Rogers 2020; Mousis et al. 2020), given the large volatile fraction. Another possibility is that the thin layer is primordial in nature and is the remnant of an atmosphere that has undergone escape. This picture would support recent modelling that has suggested that highly irradiated planets with a high water mass fraction can lose their H/He envelopes and result in observed sub-Neptunes (Aguichine et al. 2021).

Indeed, due to the overlap of multiple mass–radius curves in the parameter space around TOI-1064 b, future work assessing the atmospheric properties of high-density sub-Neptunes will help inform internal structure and atmospheric modelling in this region.

As seen in the right-hand column of Fig. E2, the posterior of the initial atmospheric mass fraction of TOI-1064 c is also flat. This distribution could suggest that planet TOI-1064 c also migrated through the protoplanetary disc having accreted a larger primordial envelope at an exterior orbital distance prior to migrating inwards and losing part of its envelope. However, it is also possible that the large uncertainty on the planet mass means that the atmospheric mass-loss rate of the atmosphere of planet TOI-1064 c is unconstrained and thus the initial atmospheric mass fraction cannot be modelled.

We derived the current atmospheric mass-loss rate by inserting the system parameters and the current stellar high-energy emission at the distance of the planet derived from the atmospheric evolution framework ($4.3 \times 10^4 \text{ erg cm}^{-2} \text{ s}^{-1}$) into the interpolation routine of Kubyshkina et al. (2018) obtaining a value of $6.11 \times 10^{11} \text{ g s}^{-1}$. With such a high mass-loss rate, the planet is expected to lose its H/He atmosphere as estimated by the planetary structure models within just 8 Myr, which implies that the planetary mass is likely to be higher than the best fit. With a mass of $4.3 M_\oplus$, corresponding to the best-fitting planetary mass plus $1\sigma$, we obtain that the atmosphere would require 100 Myr to escape. At the upper mass limit of $8.5 M_\oplus$, instead, the atmosphere would escape within about 400 Myr. We remark that these calculations consider the atmospheric mass fraction computed for the best-fitting planetary mass. If the planet were more massive, then the planetary surface gravity would be higher, leading to longer atmospheric-escape time-scales. This strengthens the possibility that planet TOI-1064 c has likely gone through, and possibly is still going through, significant migration. Utilizing a formulation of the Jeans escape criteria (Jeans 1925) that include Roche lobe effects (Erkaev et al. 2007) as detailed in equations (7)–(9) of Fossati et al. (2017) and setting $\tilde{\lambda}^*$ to the recommended value of 25, we compute the minimum mass of planet TOI-1064 c required to retain a hydrogen atmosphere to be ~6.1 $M_\oplus$. Whilst these analyses imply that the true mass of TOI-1064 c is near the upper mass limit, in order to accurately constrain the mass of this planet and to advance our understanding of the formation and evolution of this planetary system, further RV observations are needed.

Lastly, our atmospheric escape modelling permits us to calculate the integrated X-ray flux of the planets over the lifetime of the system. Even under the assumption that both planets in the TOI-1064 system formed at their current orbital distances, the integrated X-ray fluxes we determine ($F_X = 1.9 \times 10^{20}$ and $8.1 \times 10^{19}$ erg cm$^{-2}$, respectively) are lower than amount of flux required to strip a sub-Neptune-sized planet around a K-dwarf via photoevaporation alone (McDonald et al. 2019). Therefore, as we think that migration may have occurred, hence lowering the received X-ray flux, and that our internal structure modelling finds that planet TOI-1064 b has a tenuous atmosphere, we conclude that core-powered mass-loss also plays an important role in sculpting exoplanet atmospheres.

## 6.3 Planetary architecture and formation mechanisms

The planets of the TOI-1064 system have a transit-chord-length ratio, $\xi = 1.03 \pm 0.06$, indicative of coplanar orbits (Steffen et al. 2010; Fabrycky et al. 2014) that is common for multiplanet systems (Fang & Margot 2012). The eccentricities for both planets are also typical of compact, multiplanet systems (Limbach & Turner 2015; Van Eylen & Albrecht 2015; Mills et al. 2019; Van Eylen et al. 2019). The coplanar and roughly circular orbits support formation at larger semimajor axes and subsequent migration through a disc, either protoplanetary or planetesimal (Terquem & Papaloizou 2007). Currently, the planets are in orbits that are short or ‘narrow’ of the 2:1 mean motion resonance (MMR; 1.90:1) that is somewhat rare (Fabrycky et al. 2014) as planets have been observed to occur more frequently beyond MMRs in so-called ‘pile-ups’ (Batygin & Morbidelli 2013; Silburt & Rein 2015; Winn & Fabrycky 2015). Whilst a precise mass value for planet TOI-1064 c still needs to be obtained, our $3\sigma$ upper limit of $8.5 M_\oplus$ yields one of the largest mass ratios among pairs of small planets with similar radii. Thus, TOI-1064 likely joins other ‘different twin’ systems that have a significant mass, and therefore density, disparity such as Kepler-107 (Bonomo et al. 2019), KOI-1599 (Panichi, Migaszewski & Goździewski 2019), and TOI-125 (Nielsen et al. 2020), although TOI-1064 would be only the second system with such planets near an MMR. If this mass imbalance is verified, the migration scenario for TOI-1064 may support the energy optimization modelling of ‘peas in a pod’ systems presented in Adams et al. (2020), which proposes that compact systems of similar-mass planets can form in situ, whereas non-comparable mass planets should form in the outer planetary systems where the critical mass (equation 22 in Adams et al. 2020) is small and there can be large mass disparities.

Subsequent to the discovery of the radius valley (Fulton et al. 2017; Owen & Wu 2017), many studies have aimed at understanding the underlying physical processes with observationally well-characterized planets playing a key role (Cloutier et al. 2020; Van Eylen et al. 2021). As noted above, planet TOI-1064 b has a thin gaseous envelope likely having undergone atmospheric loss in the past. Using the quantification of the radius valley in period–radius space from Van Eylen et al. (2018), we find that TOI-1064 b lies on the upper edge of the valley that is supported by our observations and internal structure modelling.

Recent work using archival data has shown that there appears to be a gap in the period–mass distribution for small, short-period planets (Armstrong et al. 2019). This valley is reported to decrease from $18.9 M_\oplus$ at a rate of $-0.9 M_\oplus \text{ d}^{-1}$ and extend over a width of $4.0 M_\oplus$. With a mass of $13.5_{-1.8}^{+1.7} M_\oplus$, TOI-1064 b falls within this valley. Although the underlying process is currently unknown, the presence of planet TOI-1064 b in the gap can help refine the gradient and width of the valley, and under the assumption of the migration of TOI-1064 b, it could point towards a zero-torque location mechanism for shaping the gap, although other processes have been proposed (Lyra, Paardekooper & Mac Low 2010; Bitsch et al. 2015; Morbidelli & Raymond 2016).

## 7 SUMMARY

In this study, we present the multiplanet TOI-1064 system identified by TESS. We report the discovery of the pair of sub-Neptunes, TOI-1064 b and TOI-1064 c, and characterize both planets using high-precision space-based photometry from TESS and CHEOPS, and ground-based photometry from LCOGT, NGTS, and ASTEP, and HARPS RVs. We also present a novel method for detrending
photometric data via the monitoring of PSF shape changes and successfully apply it to the CHEOPS photometry.

From our analyses find that TOI-1064 is a bright K-dwarf with a stellar rotation period of ~26 d that hosts two warm sub-Neptunes, TOI-1064 b and TOI-1064 c, on orbital periods of ~6.44 and ~12.23 d. Interestingly, although this pair of sub-Neptunes has similar radii, the masses are substantially different resulting in a density disparity. Placing planets TOI-1064 b and TOI-1064 c on a mass-radius diagram, we find that TOI-1064 b is one of the smallest planets with a mass greater than $10 M_⊕$. It thus joins a small group of dense, warm and mildly irradiated sub-Neptunes around old K-dwarfs. TOI-1064 c, on the other hand, is potentially one of the least dense small planets known, however further RV observations are needed to confirm this.

Upon analysis of well-characterized, massive sub-Neptunes, including TOI-1064 b, we find a negative correlation between bulk density and host star metallicity. This trend is strengthened when highly irradiated planets are removed from the sample, potentially indicating that this correlation is a result of the formation of planets rather than subsequent evolutionary processes.

Our internal structure analysis shows that TOI-1064 b has a tenuous atmosphere and a large volatile component by mass fraction. The mass fractions of the component layers and the total refractory mass of TOI-1064 b may help shed light on the formation mechanisms for such dense bodies. Whilst there are similarities in the mass fractions for the high-Z components of both planets, TOI-1064 c likely has an extended gaseous envelope. Our atmospheric escape evolution modelling shows that TOI-1064 b likely lost its primordial atmosphere after migration through the protoplanetary disc. TOI-1064 c may have also lost some of its envelope subsequent to migration, however understanding the origin and evolution of this planetary system requires first the acquisition of additional RV data to better measure the mass of planet TOI-1064 c.

Future atmospheric studies may shed light on whether the atmosphere of TOI-1064 b is primordial or secondary in nature, and may potentially constrain formation pathways of TOI-1064 c after the acquisition of further RV data.

## ACKNOWLEDGEMENTS

This study is based on observations made with ESO Telescopes at the La Silla Observatory under program ID 1102.C-0923. The authors would like to thank the anonymous referee for helpful comments that improved the paper. CHEOPS is an ESA mission in partnership with Switzerland with important contributions to the payload and the ground segment from Austria, Belgium, France, Germany, Hungary, Italy, Portugal, Spain, Sweden, and the UK. The CHEOPS Consortium would like to gratefully acknowledge the support received by all the agencies, offices, universities, and industries involved. Their flexibility and willingness to explore new approaches were essential to the success of this mission. Funding for the TESS mission is provided by NASA Science Mission Directorate. We acknowledge the use of public TESS data from pipelines at the TESS Science Office and at the TESS Science Processing Operations Center. This research has made use of the Exoplanet Follow-up Observation Program website, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center for the production of the SPOC data products. Part of this work was done using data taken by

KESPRINT, an international consortium devoted to the characterization and research of exoplanets discovered with space-based missions (http://www.kesprint.science/). This work makes use of observations from the LCOGT network. Part of the LCOGT telescope time was granted by NOIRLab through the Mid-Scale Innovations Program (MSIP). MSIP is funded by NSF. This study is based on data collected under the NGTS project at the ESO La Silla Paranal Observatory. The NGTS facility is operated by the consortium institutes with support from the UK Science and Technology Facilities Council (STFC) projects ST/M001962/1 and ST/S002642/1. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This work makes use of observations from the ASTEP telescope. ASTEP benefited from the support of the French and Italian polar agencies IPEV and PNRA in the framework of the Concordia station program and from Idex UCAJEDI (ANR-15-IDEX-01). Some of the observations in the paper made use of the High-Resolution Imaging instrument Zorro. Zorro was funded by the NASA Exoplanet Exploration Program and built at the NASA Ames Research Center by Steve B. Howell, Nic Scott, Elliott P. Horch, and Emmett Quigley. Zorro was mounted on the Gemini-South telescope of the international Gemini Observatory, a program of NSF's OIR Lab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation on behalf of the Gemini partnership: the National Science Foundation (USA), National Research Council Canada (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea). TGW, ACC, and KH acknowledge support from STFC consolidated grant numbers ST/R000824/1 and ST/V000861/1, and UKSA grant ST/R003203/1. YA and MJH acknowledge the support of the Swiss National Fund under grant 200020-172746. DG and LMS gratefully acknowledge financial support from the CRT Foundation under grant no. 2018.2323 'Gaseous or rocky? Unveiling the nature of small worlds'. DG, MF, XB, SC, and JL acknowledge their roles as ESA-appointed CHEOPS science team members. CMP, MF, JK, and AJM gratefully acknowledge the support of the Swedish National Space Agency (SNSA; DNR 65/19, 174/18, 2020-00104, and Career grant 120/19C). ADe and DE acknowledge support from the European Research Council (ERC) under the European Union's Horizon 2020 Framework Programme (project FOUR ACES; grant agreement no. 724427). ADe, ALe, and HO acknowledge support from the Swiss National Centre for Competence in Research 'PlanetS' and the Swiss National Science Foundation (SNSF). SH gratefully acknowledges CNES funding through the grant 837319. SES have received funding from the European Research Council (ERC) under the European Union's Horizon 2020 Framework Programme (grant agreement no. 833925, project STAREX). SGS acknowledges support from FCT through FCT contract no. CEECIND/00826/2018 and POPH/FSE (EC). We acknowledge support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grants ESP2016-80435-C2-1-R, ESP2016-80435-C2-2-R, PGC2018-098153-B-C33, PGC2018-098153-B-C31, ESP2017-87676-C5-1-R, MDM-2017-0737 Unidad de Excelencia Maria de Maeztu-Centro de Astrobiología (INTA-CSIC), and the support of the Generalitat de Catalunya/CERCA
programme. The MOC activities have been supported by the ESA contract no. 4000124370. SCCB and VA acknowledge support from FCT through FCT contracts no. IF/01312/2014/CP1215/CT0004 and IF/00650/2015/CP1273/CT0001, respectively. ABr was supported by the SNSA. This project was supported by the CNES. LD is an F.R.S.-FNRS Postdoctoral Researcher. The Belgian participation to CHEOPS has been supported by the Belgian Federal Science Policy Office (BELSPO) in the framework of the PRODEX Program, and by the University of Liège through an ARC grant for Concerted Research Actions financed by the Wallonia-Brussels Federation. This work was supported by FCT – Fundação para a Ciência e a Tecnologia – through national funds and by FEDER through COMPETE2020 – Programa Operacional Competitividade e Internacionalização by these grants: UID/FIS/04434/2019, UIDB/04434/2020, UIDP/04434/2020, PTDC/FIS-AST/32113/2017 & POCI-01-0145-FEDER-032113, PTDC/FIS-AST/28953/2017 & POCI-01-0145-FEDER-028953, PTDC/FIS-AST/28987/2017, and POCI-01-0145-FEDER-028987, ODSD is supported in the form of work contract (DL 57/2016/CP1364/CT0004) funded by national funds through FCT. B-OD acknowledges support from the Swiss National Science Foundation (PP00P2-190080). DD acknowledges support from the TESS Guest Investigator Program grant 80NSSC19K1727 and NASA Exoplanet Research Program grant 18-2XRP18.2-0136. MG is an F.R.S.-FNRS Senior Research Associate. KGI is the ESA CHEOPS Project Scientist and is responsible for the ESA CHEOPS Guest Observers Programme. She does not participate in, or contribute to, the definition of the Guaranteed Time Programme of the CHEOPS mission through which observations described in this paper have been taken, nor to any aspect of target selection for the programme. GL acknowledges support by CARIPARO Foundation, according to the agreement CARIPARO-Università degli Studi di Padova (Pratica no. 2018/0098). This work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche. ML acknowledges support from the Swiss National Science Foundation under grant no. PCEFP2_194576. PFLM acknowledges support from STFC research grant number ST/M001040/1. LDN thanks the Swiss National Science Foundation for support under Early Postdoc. Mobility grant P2GEP2_200044. This work was also partially supported by a grant from the Simons Foundation (PI: Queloz, grant number 327127). IR acknowledges support from the Spanish Ministry of Science and Innovation and the European Regional Development Fund through grant PGC2018-098153-B-C33, as well as the support of the Generalitat de Catalunya/CERCA programme. This project has been supported by the Hungarian National Research, Development and Innovation Office (NKFIH) grant K-125015, the MTA-ELTE Lendület Milky Way Research Group, and the City of Szombathely under agreement no. 67.177-21/2016. This research received funding from the European Research Council (ERC) under the European Union's Horizon 2020 Framework Programme (grant agreement no. 803193/BEBOP), and from the Science and Technology Facilities Council (STFC; grant no. ST/S00193X/1). VVG is an F.R.S-FNRS Research Associate.

## DATA AVAILABILITY

The data underlying this paper will be made available in the CHEOPS mission archive (https://cheops.unige.ch/archive browser/). This paper includes data collected by the TESS mission, which is publicly available from the Mikulski Archive for Space Telescopes (MAST)

at the Space Telescope Science Institute (STScI) (https://mast.stsci.edu).

## REFERENCES

Acton J. S. et al., 2020, MNRAS, 498, 3115

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

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

Adibekyan V. Z., Delgado Mena E., Sousa S. G., Santos N. C., Israelian G., González Hernández J. I., Mayor M., Hakobyan A. A., 2012b, A&A, 547, A36

Adibekyan V. et al., 2015, A&A, 583, A94

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

Aguichine A., Mousis O., Deleuil M., Marcq E., 2021, ApJ, 914, 84

Allard F., Homeier D., Freytag B., 2012, Philos. Trans. R. Soc. Lond. Ser. A, 370, 2765

Alonso R. et al., 2004, ApJ, 613, L153

Anglada-Escudé G., Butler R. P., 2012, ApJS, 200, 15

Armstrong D. J., Meru F., Bayliss D., Kennedy G. M., Veras D., 2019, ApJ, 880, L1

Armstrong D. J. et al., 2020, Nature, 583, 39

Baglin A., Auvergne M., Barge P., Deleuil M., Catala C., Michel E., Weiss W., COROT Team, 2006, in Fridlund M., Baglin A., Lochard J., Conroy L., eds, The CoRoT Mission Pre-Launch Status - Stellar Seismology and Planet Finding (ESA SP-1306). ESA, Noordwijk, p. 33

Bakos G. Á., Lázár J., Papp I., Sári P., Green E. M., 2002, PASP, 114, 974

Barbary K., 2016, J. Open Source Softw., 1, 58

Batygin K., Morbidelli A., 2013, AJ, 145, 1

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

Berta Z. K. et al., 2012, ApJ, 747, 35

Bertin E., Arnouts S., 1996, A&AS, 117, 393

Bitsch B., Johansen A., Lambrechts M., Morbidelli A., 2015, A&A, 575, A28

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

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

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

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

Bonomo A. S. et al., 2019, Nat. Astron., 3, 416

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

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

Bouchy F., Sophie Team, 2006, in Arnold L., Bouchy F., Moutou C., eds, Tenth Anniversary of 51 Peg-b: Status of and prospects for Hot Jupiter Studies. Frontier Group, Paris, p. 319

Broeg C. et al., 2013, EPJ Web Conf., 47, 03005

Brown T. M. et al., 2013, PASP, 125, 1031

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

Bryant E. M. et al., 2020, MNRAS, 494, 5872

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

Celisse A., 2014, Ann. Stat., 42, 1879

Charbonneau D., Brown T. M., Latham D. W., Mayor M., 2000, ApJ, 529, L45

Charbonneau D. et al., 2009, Nature, 462, 891

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

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

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

Ciardi D. R., Beichman C. A., Horch E. P., Howell S. B., 2015, ApJ, 805, 16

Cloutier R. et al., 2020, AJ, 160, 3

Collier Cameron A. et al., 2006, MNRAS, 373, 799

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

Collins K. A., Kielkopf J. F., Stassun K. G., Hessman F. V., 2017, AJ, 153, 77
Cosentino R. et al., 2012, in McLean I. S., Ramsay S. K., Takami H., eds, Proc. SPIE Vol. 8446, Ground-Based and Airborne Instrumentation for Astronomy IV. SPIE, Bellingham, p. 84461V

Crossfield I. J. M. et al., 2016, ApJS, 226, 7

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

de Beurs Z. L. et al., 2020, preprint (arXiv:2011.00003)

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

Deming D., Harrington J., Seager S., Richardson L. J., 2006, ApJ, 644, 560

Demory B.-O. et al., 2015, MNRAS, 450, 2043

Dorn C., Khan A., Heng K., Connolly J. A. D., Alibert Y., Benz W., Tackley P., 2015, A&A, 577, A83

Dorn C., Venturini J., Khan A., Heng K., Alibert Y., Helled R., Rivoldini A., Benz W., 2017, A&A, 597, A37

Dotter A., 2016, ApJS, 222, 8

Doyle A. P., Davies G. R., Smalley B., Chaplin W. J., Elsworth Y., 2014, MNRAS, 444, 3592

Dubber S. C. et al., 2019, MNRAS, 490, 5103

Dumusque X. et al., 2017, A&A, 598, A133

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

Erkaev N. V., Kulikov Y. N., Lammer H., Selsis F., Langmayr D., Jaritz G. F., Biernat H. K., 2007, A&A, 472, 329

Espinoza N., 2018, Res. Notes Am. Astron. Soc., 2, 209

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

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

Fang J., Margot J.-L., 2012, ApJ, 761, 92

Faria J. P. et al., 2020, A&A, 635, A13

Fei Y., Murphy C., Shibazaki Y., Shahar A., Huang H., 2016, Geophys. Res. Lett., 43, 6837

Figueira P., Faria J. P., Adibekyan V. Z., Oshagh M., Santos N. C., 2016, Origins Life Evolution Biosphere, 46, 385

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

Fortney J. J., Marley M. S., Barnes J. W., 2007, ApJ, 659, 1661

Fortney J. J., Mordasini C., Nettelmann N., Kempton E. M. R., Greene T. P., Zahnle K., 2013, ApJ, 775, 80

Fossati L. et al., 2017, A&A, 598, A90

Fridlund M. et al., 2017, A&A, 604, A16

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

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

Gaia Collaboration, 2018, A&A, 616, A1

Gaia Collaboration, 2021, A&A, 649, A1

García Muñoz A., Fossati L., Youngblood A., Nettelmann N., Gandolfi D., Cabrera J., Rauer H., 2021, ApJ, 907, L36

Giacalone S. et al., 2021, AJ, 161, 24

Gordon C., Trotta R., 2007, MNRAS, 382, 1859

Guerrero N. M. et al., 2021, ApJS, 254, 39

Guillot T., Santos N. C., Pont F., Iro N., Melo C., Ribas I., 2006, A&A, 453, L21

Guillot T. et al., 2015, Astron. Nachr., 336, 638

Günther M. N. et al., 2019, Nat. Astron., 3, 1099

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

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

Haywood R. D. et al., 2014, MNRAS, 443, 2517

Howe A. R., Burrows A., Verne W., 2014, ApJ, 787, 173

Howell S. B., Everett M. E., Sherry W., Horch E., Ciardi D. R., 2011, AJ, 142, 19

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

Howell S. B., Everett M. E., Horch E. P., Winters J. G., Hirsch L., Nusdeo D., Scott N. J., 2016, ApJ, 829, L2

Howell S. B., Scott N. J., Matson R. A., Everett M. E., Furlan E., Gnilka C. L., Ciardi D. R., Lester K. V., 2021, Frontiers Astron. Space Sci., 8, 10

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

Husser T. O., Wende-von Berg S., Dreizler S., Homeier D., Reiners A., Barman T., Hauschildt P. H., 2013, A&A, 553, A6

Ikoma M., Hori Y., 2012, ApJ, 753, 66

Jeans J., 1925, The Dynamical Theory of Gases. Cambridge Univ. Press, Cambridge

Jenkins J. M., 2002, ApJ, 575, 493

Jenkins J. M. et al., 2010, in Radziwill N. M., Bridger A., eds, Proc. SPIE Vol. 7740, Software and Cyberinfrastructure for Astronomy. SPIE, Bellingham, p. 77400D

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

Jensen E., 2013, Astrophysics Source Code Library, record ascl:1306.007

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

Johnson D. R. H., Soderblom D. R., 1987, AJ, 93, 864

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

Jontof-Hutter D. et al., 2016, ApJ, 820, 39

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

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

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

Korth J. et al., 2019, MNRAS, 482, 1807

Kreidberg L., 2015, PASP, 127, 1161

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

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

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

Kurucz R. L., 1979, ApJS, 40, 1

Kurucz R. L., 1993, SYNTHE Spectrum Synthesis Programs and Line Data. Smithsonian Astrophysical Observatory, Cambridge, MA (Kurucz CD-ROM No. 18)

Léger A. et al., 2009, A&A, 506, 287

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

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

Lester K., Howell S., Ciardi D., Furlan E., Gnilka C., Matson R., Scott N., 2021, Am. Astron. Soc. Meeting Abstr., #237, 344.01

Li J., Tenenbaum P., Twicken J. D., Burke C. J., Jenkins J. M., Quintana E. V., Rowe J. F., Seader S. E., 2019, PASP, 131, 024506

Limbach M. A., Turner E. L., 2015, Proc. Natl. Acad. Sci., 112, 20

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

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

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

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

Lovis C., Pepe F., 2007, A&A, 468, 1115

Luger R. et al., 2017, Nat. Astron., 1, 0129

Luque R. et al., 2019, A&A, 623, A114

Lyra W., Paardekooper S.-J., Mac Low M.-M., 2010, ApJ, 715, L68

MacDonald M. G. et al., 2016, AJ, 152, 105

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

Marcus R. A., Sasselov D., Hernquist L., Stewart S. T., 2010, ApJ, 712, L73

Marcy G. W. et al., 2014, ApJS, 210, 20

Marigo P. et al., 2017, ApJ, 835, 77

Marounina N., Rogers L. A., 2020, ApJ, 890, 107

Matson R. A., Howell S. B., Horch E. P., Everett M. E., 2018, AJ, 156, 31

Maxted P. F. L. et al., 2022, MNRAS, in press

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

Mayor M. et al., 2003, Messenger, 114, 20

Mazeh T., Tamuz O., Zucker S., 2007, in Afonso C., Weldrake D., Henning T., eds, ASP Conf. Ser. Vol. 366, Transiting Extrasolar Planets Workshop. Astron. Soc. Pac., San Francisco, p. 119

McCully C., Volgenau N. H., Harbeck D.-R., Lister T. A., Saunders E. S., Turner M. L., Siiverd R. J., Bowman M., 2018, in Guzman J. C., Ibsen J., eds, Proc. SPIE Vol. 10707, Software and Cyberinfrastructure for Astronomy V. SPIE, Bellingham, p. 107070K

McDonald G. D., Kreidberg L., Lopez E., 2019, ApJ, 876, 22

Mékarnia D. et al., 2016, MNRAS, 463, 45

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

Mills S. M., Fabrycky D. C., Migaszewski C., Ford E. B., Petigura E., Isaacson H., 2016, Nature, 533, 509

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

Morbidelli A., Raymond S. N., 2016, J. Geophys. Res.: Planets, 121, 1962
Morris B. M. et al., 2021, A&A, 653, A173

Mortier A. et al., 2016, A&A, 585, A135

Mortier A. et al., 2020, MNRAS, 499, 5004

Mousis O., Deleuil M., Aguichine A., Marcq E., Naar J., Aguirre L. A., Brugger B., Gonçalves T., 2020, ApJ, 896, L22

Murdoch K. A., Hearnshaw J. B., Clark M., 1993, ApJ, 413, 349

Nielsen L. D. et al., 2020, MNRAS, 492, 5399

Nowak G. et al., 2020, A&A, 642, A173

Osborn H. P. et al., 2017, A&A, 604, A19

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

Owen J. E., Jackson A. P., 2012, MNRAS, 425, 2931

Owen J. E., Murray-Clay R., 2018, MNRAS, 480, 2206

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

Panichi F., Migaszewski C., Goździewski K., 2019, MNRAS, 485, 4601

Pepe F. et al., 2021, A&A, 645, A96

Piskunov N., Valenti J. A., 2017, A&A, 597, A16

Pollacco D. et al., 2006, Ap&SS, 304, 253

Pollack J. B., Hubickyj O., Bodenheimer P., Lissauer J. J., Podolak M., Greenzweig Y., 1996, Icarus, 124, 62

Queloz D. et al., 2009, A&A, 506, 303

Quinn S. N. et al., 2019, AJ, 158, 177

Rajpaul V., Aigrain S., Osborne M. A., Reece S., Roberts S., 2015, MNRAS, 452, 2269

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

Rivera E. J. et al., 2005, ApJ, 634, 625

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

Santerne A. et al., 2019, preprint (arXiv:1911.07355)

Santos N. C. et al., 2015, A&A, 580, L13

Sanz-Forcada J., Micela G., Ribas I., Pollock A. M. T., Eiroa C., Velasco A., Solano E., García-Álvarez D., 2011, A&A, 532, A6

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

Schlichting H. E., 2014, ApJ, 795, L15

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

Seager S., Kuchner M., Hier-Majumder C. A., Militzer B., 2007, ApJ, 669, 1279

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

Silburt A., Rein H., 2015, MNRAS, 453, 4089

Sinukoff E. et al., 2017, AJ, 153, 271

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

Smith A. M. S. et al., 2020, Astron. Nachr., 341, 273

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

Sneden C. A., 1973, PhD thesis, The University of Texas, Austin, TX

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

Sousa S. G., 2014, in Niemczura E., Smalley B., Pych W., eds, Determination of Atmospheric Parameters of B-, A-, F- and G-Type Stars. Springer, Cham, Switzerland, p. 297

Sousa S. G., Santos N. C., Israelian G., Mayor M., Monteiro M. J. P. F. G., 2007, A&A, 469, 783

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

Southworth J., 2011, MNRAS, 417, 2166

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

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

Stassun K. G. et al., 2019, AJ, 158, 138

Steffen J. H. et al., 2010, ApJ, 725, 1226

Steffen J. H. et al., 2012, MNRAS, 421, 2342

Steffen J. H. et al., 2013, MNRAS, 428, 1077

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

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

Sun L., Ioannidis P., Gu S., Schmitt J. H. M. M., Wang X., Kouwenhoven M. B. N., 2019, A&A, 624, A15

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

Terquem C., Papaloizou J. C. B., 2007, ApJ, 654, 1110

Tu L., Johnstone C. P., Güdel M., Lammer H., 2015, A&A, 577, L3

Twicken J. D. et al., 2018, PASP, 130, 064502

Valencia D., Sasselov D. D., O'Connell R. J., 2007, ApJ, 656, 545

Valenti J. A., Piskunov N., 1996, A&AS, 118, 595

Van Eylen V., Albrecht S., 2015, ApJ, 808, 126

Van Eylen V., Agentoft C., Lundkvist M. S., Kjeldsen H., Owen J. E., Fulton B. J., Petigura E., Snellen I., 2018, MNRAS, 479, 4786

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

Van Eylen V. et al., 2021, MNRAS, 507, 2154

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

Vanderburg A. et al., 2015, ApJ, 800, 59

Vanderburg A. et al., 2016, ApJ, 827, L10

Vogt S. S. et al., 1994, in Crawford D. L., Craine E. R., eds, Proc. SPIE Vol. 2198, Instrumentation in Astronomy VIII. SPIE, Bellingham, p. 362

Wang J., Xie J.-W., Barclay T., Fischer D. A., 2014, ApJ, 783, 4

Weiss L. M., Marcy G. W., 2014, ApJ, 783, L6

Weiss L. M., Petigura E. A., 2020, ApJ, 893, L1

Weiss L. M. et al., 2013, ApJ, 768, 14

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

Wheatley P. J. et al., 2018, MNRAS, 475, 4476

Winn J. N., Fabrycky D. C., 2015, ARA&A, 53, 409

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

Wright N. J., Drake J. J., Mamajek E. E., Henry G. W., 2011, ApJ, 743, 48

Xie J.-W., 2014, ApJS, 210, 25

Yee S. W., Petigura E. A., von Braun K., 2017, ApJ, 836, 77

Zacharias N., Finch C. T., Girard T. M., Henden A., Bartlett J. L., Monet D. G., Zacharias M. I., 2013, AJ, 145, 44

Zechmeister M., Kürster M., Endl M., 2009, A&A, 505, 859

Zeng L., Sasselov D., 2013, PASP, 125, 227

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

Ziegler C., Tokovinin A., Latiolais M., Briceño C., Law N., Mann A. W., 2021, AJ, 162, 192

## SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

### Wilson2021_Table2ext.txt

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

## APPENDIX A: A DESCRIPTION OF THE PSF-BASED SCALPELS METHOD

Here we describe the PSF-based SCALPELS method developed to model PSF shape changes in CHEOPS observations. We make use of the fact that the autocorrelation function (ACF) of the PSF is, to a high degree of approximation, invariant to the small shifts imposed by spacecraft pointing jitter. The PSF is a positive definite function of position on the CCD. We apply a simple autocorrelation procedure, imposing a sequence of shifts by integer-pixel offsets in x and y covering the image, and comultiplying the shifted and unshifted images. The resulting ACF is unwrapped into a one-dimensional array and normalized by its own mean that is given as $A(\delta r_i)$, where $(\delta r_i)$ is the aforementioned shift.

From this point on, the procedure is as described by Collier Cameron et al. (2021), i.e. we compute the singular value decomposition (SVD) of the time sequence of ACFs computed from the images following the standard procedure:

$$A(\delta r_i, t_j) = \langle A(\delta r_i, t_j) \rangle + U(t_j) \cdot \text{diag}(S) \cdot V(\delta r_i), \quad (A1)$$

where $U$, $S$, and $V$ are the components of the SVD. The diagonal matrix $S$ lists the singular values (eigenvalues) of the principal components in decreasing order. The columns of $U(t_j)$ define an
orthonormal basis in the time domain with each column comprising the time-domain coefficients of the corresponding eigenvector $V(\delta r_i)$.

We mask out images affected by cosmic rays, saturation, and other unwanted chance effects by identifying strong outliers in the time-domain coefficients $U$ of the eigenfunctions $V$. We use LOOCV (Celisse 2014) to separate the subset of the columns of $U$ whose variations are driven by genuine PSF changes, from those that simply describe noise. The LOOCV method systematically leaves single data points (or in this study, rows of the ACF) out of the analysis and aims to reconstruct the missing data. If this prediction accurately matches the removed data we consider, this row to represent PSF changes, whereas if the reconstruction poorly reproduces the data, it is interpreted as noise.

The response of the photometric data time series $d_{\text{obs}}$ to the set of normal modes of PSF shape variation described by the columns of $U$ is found by applying the projection operator $U \cdot U^T$ to the light curve. Subtracting this response from the light curve is equivalent to projecting the light curve into the orthogonal complement of $U$, via the linear operator $P_{\perp} = (I - U \cdot U^T)$.

This procedure is only valid if the target flux is already corrected for intrinsic stellar variability, because the target variations may have a non-zero projection into the null space that would distort the aperture correction. In reality we are interested in modelling the intrinsic variations of the target with a time-dependent model light curve $m(\theta)$, where $\theta$ is the vector of parameters defining the light-curve model.

The likelihood of the data, $d_{\text{obs}}$, given the light-curve model parameters, $\theta$, is found by subtracting the model from the data and projecting into the orthogonal complement of the null space:

$$ 
\begin{aligned}
\ln \mathcal{L} &= -\frac{1}{2} [P_{\perp} \cdot (d_{\text{obs}} - m(\theta))]^T \cdot \Sigma^{-1} \cdot [P_{\perp} \cdot (d_{\text{obs}} - m(\theta))] \\
&= -\frac{1}{2} \ln |\Sigma| - \frac{N}{2} \ln(2\pi) \\
&= -\frac{1}{2} [P_{\perp} \cdot (d_{\text{obs}} - m(\theta))]^T \cdot \Sigma^{-1} \cdot [P_{\perp} \cdot (d_{\text{obs}} - m(\theta))] \\
&\quad -\frac{1}{2} \left( \sum_{j=1}^{N} \ln \Sigma_{jj} \right) - \frac{N}{2} \ln(2\pi),
\end{aligned}
\quad (A2) $$

where the covariance matrices are $\Sigma_{jj}$ and thus $\sigma_j^2 = \text{var}(d_j)$ are the variances of the individual photometric measurements $d_j$.

Used with MCMC sampling, this allows us to determine the posterior probability distribution of the model parameters $\theta$, while simultaneously correcting for aperture losses arising from PSF shape changes. For the optimal set of model parameters $\theta_*$, the magnitude correction to the light curve is

$$ d_{\parallel} = P_{\parallel} \cdot (d_{\text{obs}} - m(\theta_*)), \quad (A3) $$

where $P_{\parallel} = U \cdot U^T$ is the projection operator into the null space. From this approach, the corrected flux light curve is then $f_{\text{corr}} = d_{\text{obs}} - d_{\parallel}$.

## APPENDIX B: RAW CHEOPS LIGHT CURVES

In addition to the raw data and fitted linear model of the first CHEOPS visit of TOI-1064 seen in Fig. 2, here for completeness we present the

<figure><img src="image_18.png" /><figcaption>**Figure B1.** The normalized light curve of the second CHEOPS visit of TOI-1064 with the location of the transit of planet TOI-1064 b (green) shown. Top panel: the raw DRP-produced fluxes. Second panel: the DRP fluxes with the fitted linear model produced from the measured PSF shape changes, as detailed in Appendix A. Third panel: the detrended fluxes with the nominal transit model from the global analysis (Section 4.3). Bottom panel: residuals to the fit.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure B2.** Same as Fig. B1, but for the third CHEOPS visit with the location of the transit of planet TOI-1064 c (purple) shown.</figcaption></figure>

raw light curves and linear models for the remaining five CHEOPS observations (Figs B1–B5). The raw data plotted in the top panels are the CHEOPS DRP-produced photometry outlined in Section 2.2, with the linear models in the bottom panels constructed following the principal component analysis of the ACF of the CHEOPS PSF as described in Appendix A.
<figure><img src="image_20.png" /><figcaption>**Figure B3.** Same as Fig. B1, but for the fourth CHEOPS visit with the location of the transit of planet TOI-1064 b (green) shown.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure B4.** Same as Fig. B1, but for the fifth CHEOPS visit with the location of the transit of planet TOI-1064 b (green) shown.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure B5.** Same as Fig. B1, but for the sixth CHEOPS visit with the location of the transit of planet TOI-1064 c (purple) shown.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure C1.** Generalized Lomb–Scargle periodograms of the HARPS RVs and spectral activity indicators from the HARPS DRS and TERRA. The horizontal dashed line marks the 0.1 per cent FAP level. The vertical dashed red lines mark the orbital frequencies of the two transiting planets ($f_b = 0.155$ d$^{-1}$ and $f_c = 0.082$ d$^{-1}$). The vertical dashed blue lines mark the frequencies of the stellar signal at ~0.038 d$^{-1}$ and its first and second harmonics. Upper panel: HARPS DRS RVs. Second panel: RV residuals following the subtraction of the signal of TOI-1064 b. Remaining panels: the activity indicators of TOI-1064: the BIS, the FWHM, and the contrast of the CCF, the Ca II H&K lines activity indicator (S-index), $H\alpha$, Na D1, and Na D2.</figcaption></figure>

## APPENDIX C: PERIODOGRAM ANALYSIS OF THE HARPS MEASUREMENTS

We performed a frequency analysis of the HARPS RV measurements and spectral diagnostics in order to search for the Doppler reflex motion induced by the two transiting planets TOI-1064 b and TOI-1064 c, and look for possible additional signals that might be associated with stellar activity or other orbiting bodies in the system.
Fig. C1 displays the generalized Lomb–Scargle periodograms (GLS; Zechmeister, Kürster & Endl 2009) of the HARPS DRS RVs and activity indicators extracted with DRS and with TERRA. The horizontal black dashed lines mark the power corresponding to the 0.1 per cent false alarm probability (FAP). We estimated the FAP following the bootstrap method described in Murdoch, Hearnshaw & Clark (1993), i.e. by computing the GLS periodogram of 10⁶ mock time series obtained by randomly shuffling the RV measurements and their uncertainties whilst keeping the time stamps fixed. We defined the FAP as the fraction of those periodograms whose highest power exceeds the observed power of $f_b$ in the periodogram of the original HARPS data over the frequency range 0–0.3 d⁻¹.

The GLS periodogram of the HARPS DRS RVs (Fig. C1, upper panel) shows its highest peak at $f_b = 0.155$ d⁻¹, i.e. the orbital frequency of the inner transiting planet initially detected in the TESS light curve (right dashed red line). Although this peak is not significant (FAP > 0.1 per cent) within the frequency range 0–0.3 d⁻¹, the occurrence of this peak at a known frequency provides strong evidence that this signal is related to planet TOI-1064 b. We estimated the FAP for a peak at the orbital frequency of TOI-1064 b and found the probability that the periodogram of white/red noise can display a peak at the orbital frequency of the inner planet, with a power higher than the observed power, is only 0.06 per cent, making the peak at $f_b$ significant. The peak at $f_b$ has no counterpart in the periodograms of the activity indicators, suggesting that the signal detected in Doppler data is induced by an orbiting companion and confirms the presence of the inner transiting planet seen in both TESS and CHEOPS data with a period of about 6.44 d.

We subtracted the Doppler signal of TOI-1064 b from the HARPS RVs by fitting a circular-orbit model, fixing period and time of

first transit as derived from the analysis of the light curves. The periodogram of the RV residuals (Fig. C1, second panel) displays a non-significant peak at ~8.7 d (FAP > 0.1 per cent). This non-significant peak is also observed in the periodogram of the CCF BIS, FWHM, contrast, and $H\alpha$. However, observing the same peak in different unrelated activity periodograms provides evidence that it could be real. The rotation period of the star – as inferred from the WASP light curve – is about 26.6 d. The peak at 8.7 d would then be near the second harmonic of this period. Active regions separated by about 120° in longitude could account for this signal. In addition, the FWHM, the contrast, and the S-index show powers at frequencies lower than 0.04 d⁻¹ (25 d), which might be due to activity coupled to stellar rotation. Unfortunately, although seen in TESS, CHEOPS, LCOGT, NGTS, and ASTEP photometry, the outer planet at 12.23 d is undetected in these HARPS RVs periodograms. Therefore, the use of the SCALPELS algorithm to potentially separate the stellar and planetary signals is needed; see Section 4.2.

## APPENDIX D: JOINT PHOTOMETRY AND RV FIT PRIORS AND POSTERIOR DISTRIBUTIONS

Here we present a table of the priors used in the global analysis of the TOI-1064 system (Tables D1 and D2), and corner plots of the posterior distributions of the main fitted transit and RV parameters for planets TOI-1064 b and TOI-1064 c, and the stellar density (Fig. D1). The corner plots show a slight correlation between $r_1$ and $\sqrt{e\sin\omega}$ that might indicate a tension between the fitting the ingresses and egresses of the transit photometry to obtain the impact parameters and two Keplerian orbits to obtain the eccentricities of the planets. This could manifest as the bimodality around 0 for the $\sqrt{e\sin\omega}$ component of TOI-1064 c.

<figure><img src="image_24.png" /><figcaption>**Figure D1.** The posterior distributions of the fitted transit and RV parameters for (a) planet TOI-1064 b and (b) planet TOI-1064 c, including the orbital period ($P$), the mid-transit time ($T_0$), the ($r_1, r_2$) parametrization of the impact parameter ($b$) and planet-to-star radius ratio ($p$) plane, the $\sqrt{e\cos\omega}$ and $\sqrt{e\sin\omega}$ parametrization, where $e$ is the eccentricity and $\omega$ is the argument of periastron, the semi-amplitude ($K$), and the stellar density ($\rho_*$).</figcaption></figure>
<table><caption>**Table D1.** Priors used in the joint photometry and RV fit. Uniform priors are represented by $U(a, b)$, with lower and upper bounds of $a$ and $b$, whereas $N(\mu, \sigma)$ indicates a normal (Gaussian) prior with mean, $\mu$, and standard deviation, $\sigma$.</caption><thead><tr><th>Parameter (unit)</th><th>TOI-1064 b</th><th>TOI-1064 c</th></tr></thead><tbody><tr><td>P (d)</td><td>U(6.437785, 6.446653)</td><td>U(12.215945, 12.258251)</td></tr><tr><td>T<sub>0</sub> (BJD - 245 7000)</td><td>U(2036.849036, 2036.865344)</td><td>U(2043.489928, 2043.519832)</td></tr><tr><td>r<sub>1</sub></td><td>U(0, 1)</td><td>U(0, 1)</td></tr><tr><td>r<sub>2</sub></td><td>U(0, 1)</td><td>U(0, 1)</td></tr><tr><td>√e cosω</td><td>U(-0.5, 0.5)</td><td>U(-0.5, 0.5)</td></tr><tr><td>√e sinω</td><td>U(-0.5, 0.5)</td><td>U(-0.5, 0.5)</td></tr><tr><td>K (m s<sup>-1</sup>)</td><td>U(0, 20)</td><td>U(0, 20)</td></tr><tr><td>γ<sub>*</sub> (m s<sup>-1</sup>)</td><td></td><td>U(21 186.8, 21 243.0)</td></tr><tr><td>ρ<sub>*</sub> (kg m<sup>-3</sup>)</td><td></td><td>N(2760, 140)</td></tr><tr><td>q<sub>1</sub>,TESS</td><td></td><td>N(0.32, 0.1)</td></tr><tr><td>q<sub>2</sub>,TESS</td><td></td><td>N(0.53, 0.1)</td></tr><tr><td>q<sub>1</sub>,CHEOPS</td><td></td><td>N(0.33, 0.1)</td></tr><tr><td>q<sub>2</sub>,CHEOPS</td><td></td><td>N(0.44, 0.1)</td></tr><tr><td>q<sub>1</sub>,LCOGT<sub>zs</sub></td><td></td><td>N(0.34, 0.1)</td></tr><tr><td>q<sub>2</sub>,LCOGT<sub>zs</sub></td><td></td><td>N(0.60, 0.1)</td></tr><tr><td>q<sub>1</sub>,LCOGT<sub>B</sub></td><td></td><td>N(0.36, 0.1)</td></tr><tr><td>q<sub>2</sub>,LCOGT<sub>B</sub></td><td></td><td>N(0.26, 0.1)</td></tr><tr><td>q<sub>1</sub>,NGTS</td><td></td><td>N(0.32, 0.1)</td></tr><tr><td>q<sub>2</sub>,NGTS</td><td></td><td>N(0.47, 0.1)</td></tr><tr><td>q<sub>1</sub>,ASTEP<sub>Rc</sub></td><td></td><td>N(0.30, 0.1)</td></tr><tr><td>q<sub>2</sub>,ASTEP<sub>Rc</sub></td><td></td><td>N(0.43, 0.1)</td></tr></tbody></table>

<table><caption>**Table D2.** Fitted light curve and RV parameter posterior values for the joint fit to the data detailed in Section 4.3.</caption><caption>Fitted light curve and RV parameters</caption><tbody><tr><td>q<sub>1</sub>,TESS</td><td>0.336<sup>+0.070</sup><sub>-0.068</sub></td></tr><tr><td>q<sub>2</sub>,TESS</td><td>0.464<sup>+0.076</sup><sub>-0.068</sub></td></tr><tr><td>σ<sub>jitter,TESS</sub></td><td>0.0007652<sup>+0.0000094</sup><sub>-0.0000095</sub></td></tr><tr><td>q<sub>1</sub>,CHEOPS</td><td>0.265<sup>+0.047</sup><sub>-0.050</sub></td></tr><tr><td>q<sub>2</sub>,CHEOPS</td><td>0.521<sup>+0.066</sup><sub>-0.058</sub></td></tr><tr><td>σ<sub>jitter,CHEOPS</sub></td><td>0.0000012<sup>+0.0001557</sup><sub>-0.0000012</sub></td></tr><tr><td>q<sub>1</sub>,LCOGT<sub>zs</sub></td><td>0.343<sup>+0.081</sup><sub>-0.075</sub></td></tr><tr><td>q<sub>2</sub>,LCOGT<sub>zs</sub></td><td>0.592<sup>+0.075</sup><sub>-0.079</sub></td></tr><tr><td>σ<sub>jitter,LCOGT<sub>zs</sub></sub></td><td>0.00099960<sup>+0.00000027</sup><sub>-0.00000039</sub></td></tr><tr><td>q<sub>1</sub>,LCOGT<sub>B</sub></td><td>0.934<sup>+0.275</sup><sub>-0.222</sub></td></tr><tr><td>q<sub>2</sub>,LCOGT<sub>B</sub></td><td>0.721<sup>+0.289</sup><sub>-0.284</sub></td></tr><tr><td>σ<sub>jitter,LCOGT<sub>B</sub></sub></td><td>0.000711<sup>+0.000073</sup><sub>-0.000069</sub></td></tr><tr><td>q<sub>1</sub>,NGTS</td><td>0.300<sup>+0.072</sup><sub>-0.074</sub></td></tr><tr><td>q<sub>2</sub>,NGTS</td><td>0.505<sup>+0.074</sup><sub>-0.070</sub></td></tr><tr><td>σ<sub>jitter,NGTS</sub></td><td>0.000012<sup>+0.000094</sup><sub>-0.000011</sub></td></tr><tr><td>q<sub>1</sub>,ASTEP<sub>Rc</sub></td><td>0.352<sup>+0.062</sup><sub>-0.053</sub></td></tr><tr><td>q<sub>2</sub>,ASTEP<sub>Rc</sub></td><td>0.428<sup>+0.074</sup><sub>-0.075</sub></td></tr><tr><td>σ<sub>jitter,ASTEP<sub>Rc</sub></sub></td><td>0.000851<sup>+0.000042</sup><sub>-0.000043</sub></td></tr><tr><td>σ<sub>jitter,HARPS</sub> (m s<sup>-1</sup>)</td><td>4.01<sup>+0.58</sup><sub>-0.49</sub></td></tr></tbody></table>

## APPENDIX E: ATMOSPHERIC EVOLUTION POSTERIOR DISTRIBUTIONS

Here we present the posterior distribution plots of results of our atmospheric evolution modelling of the TOI-1064 system (Figs E1 and E2).
<figure><img src="image_25.png" /><figcaption>**Figure E1.** Posterior distributions as inferred from the atmospheric evolution code described in Section 5.2. Purple areas highlight the highest probability density (HPD) intervals at the 68 per cent level. Red lines are the normal priors imposed on the parameters, while the background histogram plotted with the thin black line in the leftmost panel is the distribution of the stellar rotation period observed in open clusters of 150 Myr (Johnstone et al. 2015) considering only stars whose masses deviate from $M_*$ by less than 0.1 $M_\odot$.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Figure E2.** Posterior distributions of the planetary masses (upper panels) and initial atmospheric mass fractions (bottom panels). Purple areas highlight the highest probability density (HPD) intervals at the 68 per cent level. The red lines represent the priors imposed on the parameters, which are uniform in the case of $f_{at}^{start}$ (free parameters), while the light blue lines, spiking near $f=0$, in the bottom panels represent the present-day atmospheric mass fractions. The atmospheric content is basically negligible for planet TOI-1064 b, while it is barely visible, though non-negligible, for planet TOI-1064 c.</figcaption></figure>

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

<footnote n="2">²Dipartimento di Fisica, Universita degli Studi di Torino, via Pietro Giuria 1, I-10125 Torino, Italy</footnote>

<footnote n="3">³Thüringer Landessternwarte Tautenburg, Sternwarte 5, D-07778 Tautenburg, Germany</footnote>

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

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

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

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

<footnote n="8">⁸Department of Space, Earth and Environment, Astronomy and Plasma Physics, Chalmers University of Technology, SE-41296 Gothenburg, Sweden</footnote>

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

<footnote n="10">¹⁰Observatoire Astronomique de l'Université de Genève, Chemin Pegasi 51, CH-1290 Versoix, Switzerland</footnote>

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

<footnote n="12">¹²Department of Astronomy, Stockholm University, SE-10691 Stockholm, Sweden</footnote>

<footnote n="13">¹³Aix Marseille Univ, CNRS, CNES, LAM, F-13013 Marseille, France</footnote>

<footnote n="14">¹⁴Division Technique INSU, BP 330, F-83507 La Seyne cedex, France</footnote>

<footnote n="15">¹⁵Instituto de Astrofísica e Ciencias do Espaco, Universidade do Porto, CAUP, Rua das Estrelas, P-4150-762 Porto, Portugal</footnote>
<footnote n="16">¹⁶Laboratoire d'Astrophysique de Marseille, 38 Rue Frédéric Joliot Curie, F-13013 Marseille, France</footnote>

<footnote n="17">¹⁷Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University, Box 43, SE-22100 Lund, Sweden</footnote>

<footnote n="18">¹⁸Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, CS 34229, F-06304 Nice Cedex 4, France</footnote>

<footnote n="19">¹⁹Instituto de Astrofísica de Canarias, E-38200 La Laguna, Tenerife, Spain</footnote>

<footnote n="20">²⁰Departamento de Astrofísica, Universidad de La Laguna, E-38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="21">²¹Institut de Ciencies de l'Espai (ICE, CSIC), Campus UAB, Can Magrans s/n, E-08193 Bellaterra, Spain</footnote>

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

<footnote n="23">²³European Space Agency (ESA), European Space Research and Technology Centre (ESTEC), Keplerlaan 1, NL-2201 AZ Noordwijk, the Netherlands</footnote>

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

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

<footnote n="26">²⁶Departamento de Física e Astronomia, Faculdade de Ciências, Universi- dade do Porto, Rua do Campo Alegre, P-4169-007 Porto, Portugal</footnote>

<footnote n="27">²⁷Université Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

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

<footnote n="29">²⁹Department of Physics and Astronomy, University of Leicester, Leicester LE1 7RH, UK</footnote>

<footnote n="30">³⁰Concordia Station, Dome C, Antarctica</footnote>

<footnote n="31">³¹Institute of Planetary Research, German Aerospace Center (DLR), Ruther- fordstrasse 2, D-12489 Berlin, Germany</footnote>

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

<footnote n="33">³³NASA Exoplanet Science Institute, Caltech/IPAC, 1200 E. California Ave, Pasadena, CA 91125, USA</footnote>

<footnote n="34">³⁴Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="35">³⁵Center for Planetary Systems Habitability and McDonald Observatory, The University of Texas at Austin, Austin, TX 78712, USA</footnote>

<footnote n="36">³⁶NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA</footnote>

<footnote n="37">³⁷Centre for Mathematical Sciences, Lund University, Box 118, SE-22100 Lund, Sweden</footnote>

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

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

<footnote n="40">⁴⁰Department of Physics and Astronomy, University of New Mexico, 1919 Lomas Blvd NE, Albuquerque, NM 87131, USA</footnote>

<footnote n="41">⁴¹School of Physics and Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK</footnote>

<footnote n="42">⁴²Department of Astronomy and Tsinghua Centre for Astrophysics, Tsinghua University, Beijing 100084, China</footnote>

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

<footnote n="44">⁴⁴Rhenish Institute for Environmental Research, Department of Planetary Research, University of Cologne, Aachener Strasse 209, D-50931 Cologne, Germany</footnote>

<footnote n="45">⁴⁵Department of Astrophysics, University of Vienna, Türkenschanzstrasse 17, A-1180 Vienna, Austria</footnote>

<footnote n="46">⁴⁶Science and Operations Department – Science Division (SCI-SC), Direc- torate of Science, European Space Agency (ESA), European Space Research</footnote>

<footnote>and Technology Centre (ESTEC), Keplerlaan 1, NL-2201 AZ Noordwijk, the Netherlands</footnote>

<footnote n="47">⁴⁷Department of Physics and Astronomy, Swarthmore College, Swarthmore, PA 19081, USA</footnote>

<footnote n="48">⁴⁸Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, Konkoly Thege Miklós ut 15-17, 1121 Budapest, Hungary</footnote>

<footnote n="49">⁴⁹Dipartimento di Fisica e Astronomia 'Galileo Galilei', Universita degli Studi di Padova, Vicolo dell'Osservatorio 3, I-35122 Padova, Italy</footnote>

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

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

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

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

<footnote n="54">⁵⁴Department of Astronomy, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan</footnote>

<footnote n="55">⁵⁵Instituto de Astrofísica de Andalucía (IAA-CSIC), Glorieta de la Astronomía s/n, E-18008 Granada, Spain</footnote>

<footnote n="56">⁵⁶Astrophysics Group, Keele University, Staffordshire ST5 5BG, UK</footnote>

<footnote n="57">⁵⁷Astronomy Department, University of California, Berkeley, CA 94720, USA</footnote>

<footnote n="58">⁵⁸Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="59">⁵⁹Department of Physics, University of Oxford, Oxford OX1 3RH, UK</footnote>

<footnote n="60">⁶⁰Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK</footnote>

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

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

<footnote n="63">⁶³Cavendish Laboratory, JJ Thomson Avenue, Cambridge CB3 0HE, UK</footnote>

<footnote n="64">⁶⁴Center for Astronomy and Astrophysics, Technical University Berlin, Hardenberstrasse 36, D-10623 Berlin, Germany</footnote>

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

<footnote n="66">⁶⁶Astronomy Department and Van Vleck Observatory, Wesleyan University, Middletown, CT 06459, USA</footnote>

<footnote n="67">⁶⁷Patashnick Voorheesville Observatory, Voorheesville, NY 12186, USA</footnote>

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

<footnote n="69">⁶⁹Department of Aeronautics and Astronautics, MIT, 77 Massachusetts Avenue, Cambridge, MA 02139, USA</footnote>

<footnote n="70">⁷⁰Hazelwood Observatory, Australia</footnote>

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

<footnote n="72">⁷²MTA-ELTE Exoplanet Research Group, Szent Imre h. u. 112, 9700 Szom- bathely, Hungary</footnote>

<footnote n="73">⁷³MTA-ELTE Lendület Milky Way Research Group, Szent Imre h. u. 112, 9700 Szombathely, Hungary</footnote>

<footnote n="74">⁷⁴Institute of Astronomy, University of Cambridge, Madingley Road, Cam- bridge CB3 0HA, UK</footnote>

<footnote n="75">⁷⁵Department of Astrophysical Sciences, Princeton University, 4 Ivy Ln, Princeton, NJ 08544, USA</footnote>

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