# A hot mini-Neptune and a temperate, highly eccentric sub-Saturn around the bright K-dwarf TOI-2134

F Rescigno, G Hébrard, A Vanderburg, A Mann, A Mortier, S Morrell, L Buchhave, K Collins, C Mann, C Hellier, et al.

► To cite this version:

F Rescigno, G Hébrard, A Vanderburg, A Mann, A Mortier, et al. A hot mini-Neptune and a temperate, highly eccentric sub-Saturn around the bright K-dwarf TOI-2134. Monthly Notices of the Royal Astronomical Society, 2024, 527 (3), pp.5385-5407. 10.1093/mnras/stad3255. hal-04505670

HAL Id: hal-04505670

https://hal.science/hal-04505670v1

Submitted on 16 Mar 2024

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
# A hot mini-Neptune and a temperate, highly eccentric sub-Saturn around the bright K-dwarf TOI-2134

F. Rescigno<sup>1,★</sup> G. Hébrard,<sup>2,3</sup> A. Vanderburg<sup>4</sup> A. W. Mann,<sup>5</sup> A. Mortier<sup>6</sup> S. Morrell,<sup>1</sup>
L. A. Buchhave<sup>7</sup> K. A. Collins,<sup>8</sup> C. R. Mann<sup>9,10</sup> C. Hellier,<sup>11</sup> R. D. Haywood,<sup>1†</sup> R. West,<sup>12</sup>
M. Stalport,<sup>13,14</sup> N. Heidari,<sup>2</sup> D. Anderson,<sup>12</sup> C. X. Huang,<sup>15</sup> M. López-Morales,<sup>8</sup> P. Cortés-Zuleta,<sup>16</sup>
H. M. Lewis,<sup>17</sup> X. Dumusque,<sup>14</sup> I. Boisse,<sup>18</sup> P. Rowden<sup>19</sup> A. Collier Cameron<sup>19,20,21</sup> M. Deleuil,<sup>16</sup>
M. Vezie,<sup>4</sup> F. A. Pepe,<sup>14</sup> X. Delfosse,<sup>22</sup> D. Charbonneau,<sup>8</sup> K. Rice,<sup>23,24</sup> O. Demangeon,<sup>25</sup> S. N. Quinn,<sup>8</sup>
S. Udry,<sup>14</sup> T. Forveille,<sup>22</sup> J. N. Winn,<sup>26</sup> A. Sozzetti<sup>19</sup> S. Hoyer,<sup>16</sup> S. Seager,<sup>4,28,29</sup> T. G. Wilson,<sup>20,21</sup>
S. Dalal,<sup>1</sup> E. Martioli,<sup>2,30</sup> S. Striegel,<sup>31,32</sup> W. Boschin,<sup>33,34,35</sup> D. Dragomir,<sup>36</sup> A. F. Martínez Fiorenzano,<sup>33</sup>
R. Cosentino,<sup>33</sup> A. Ghedina,<sup>33</sup> L. Malavolta,<sup>37,38</sup> L. Affer,<sup>39</sup> B. S. Lakeland,<sup>1</sup> B. A. Nicholson,<sup>40,41</sup>
S. Foschino,<sup>42</sup> A. Wünsche,<sup>42</sup> K. Barkaoui,<sup>28,34,43</sup> G. Srdoc,<sup>44</sup> J. Randolph,<sup>45</sup> B. Guillet,<sup>45</sup> D. M. Conti,<sup>45</sup>
M. Ghachoui,<sup>43,46</sup> M. Gillon,<sup>43</sup> Z. Benkhaldoun,<sup>46</sup> F. J. Pozuelos,<sup>43,47</sup> M. Timmermans,<sup>43</sup> E. Girardin,<sup>48</sup>
S. Matutano,<sup>49</sup> P. Bosch-Cabot,<sup>49</sup> J. A. Muñoz<sup>50,51</sup> and R. Forés-Toribio<sup>50,51</sup>

*Affiliations are listed at the end of the paper*

Accepted 2023 October 18. Received 2023 September 20; in original form 2023 June 5

## ABSTRACT

We present the characterization of an inner mini-Neptune in a 9.2292005 ± 0.0000063 d orbit and an outer mono-transiting sub-Saturn planet in a 95.50<sup>+0.36</sup><sub>-0.25</sub> d orbit around the moderately active, bright ($m_v$ = 8.9 mag) K5V star TOI-2134. Based on our analysis of five sectors of TESS data, we determine the radii of TOI-2134b and c to be 2.69 ± 0.16 $R_⊕$ for the inner planet and 7.27 ± 0.42 $R_⊕$ for the outer one. We acquired 111 radial-velocity (RV) spectra with HARPS-N and 108 RV spectra with SOPHIE. After careful periodogram analysis, we derive masses for both planets via Gaussian Process regression: 9.13<sup>+0.78</sup><sub>-0.76</sub> $M_⊕$ for TOI-2134b and 41.89<sup>+7.69</sup><sub>-7.83</sub> $M_⊕$ for TOI-2134c. We analysed the photometric and RV data first separately, then jointly. The inner planet is a mini-Neptune with density consistent with either a water-world or a rocky core planet with a low-mass H/He envelope. The outer planet has a bulk density similar to Saturn's. The outer planet is derived to have a significant eccentricity of 0.67<sup>+0.05</sup><sub>-0.06</sub> from a combination of photometry and RVs. We compute the irradiation of TOI-2134c as 1.45 ± 0.10 times the bolometric flux received by Earth, positioning it for part of its orbit in the habitable zone of its system. We recommend further RV observations to fully constrain the orbit of TOI-2134c. With an expected Rossiter-McLaughlin (RM) effect amplitude of 7.2 ± 1.3 m s<sup>-1</sup>, we recommend TOI-2134c for follow-up RM analysis to study the spin-orbit architecture of the system. We calculate the Transmission Spectroscopy Metric, and both planets are suitable for bright-mode Near Infrared Camera (NIRCam) atmospheric characterization.

**Key words:** methods: data analysis – techniques: photometric – techniques: radial velocities – planets and satellites: detection – stars: activity – stars: individual (TOI-2134, TIC 75878355, G 204-45).

## 1 INTRODUCTION

Since the discovery of the first exoplanet circa 30 yrs ago, more than 5000 have been detected and confirmed. Radial-velocity (RV) surveys performed with instruments such as the High Accuracy Radial-velocity Planet Searcher (HARPS) coupled with the *Kepler* photometric mission started discovering a subpopulation of small

exoplanets in short (under 100 d) orbits (Mayor & Udry 2008; Fressin, Guillot & Nesta 2009; Lovis et al. 2009; Borucki et al. 2011; Batalha et al. 2013). Given their abundance in our galaxy (Chabrier et al. 2000; Winters et al. 2015), and their low mass and size, K and M dwarf stars are prime candidates for small-exoplanet searches and demographic-focused studies (Dressing & Charbonneau 2013; Crossfield et al. 2015; Astudillo-Defru et al. 2017; Pinamonti et al. 2018; Rice et al. 2019; West et al. 2019; Burt et al. 2020).

The transition point between rocky super-Earths and gaseous Neptunes is still debated (Fulton et al. 2017; Luque et al. 2021). Otegi,

<footnote>* E-mail: fr307@exeter.ac.uk</footnote>

<footnote>† STFC Ernest Rutherford Fellow</footnote>
Bouchy & Helled (2020) show that this transition range is between 5–25 M<sub>⊕</sub> and 2–3 R<sub>⊕</sub>, but several factors play into the composition of these planets. Some studies report that all planets under 1.6 R<sub>⊕</sub> must be rocky (Rogers 2015; López-Morales et al. 2016). Others give more importance to the effects of irradiation: less irradiated planets are more likely to maintain a gaseous envelope, while more irradiated ones are typically rocky (Hadden & Lithwick 2014; Jontof-Hutter et al. 2016). Owen & Adams (2019) explore how planetary magnetic fields can also decrease their mass-loss rates and therefore alter the composition of the planetary cores. A continuous effort in the detection of small planets, and in the precise characterization of their masses and sizes is therefore vital to reach a consensus on which parameters affect planetary composition.

On the other hand, our understanding of long-period planets is also lacking. The great majority of transit-detected exoplanets have periods shorter than 75 d (Jiang et al. 2019). Longer period planets are harder to detect and determining their masses can be challenging. Moreover, the baselines of most photometric surveys also limit their detection. This ‘missing’ population hampers studies of planet demographics, of planet formation, and of how planetary characteristics depend on the host star (Johnson et al. 2010; Winn 2011).

Temperate giants are located in a period valley, between 10 and 100 d, where gas planets are less frequent (Udry, Mayor & Santos 2003; Wittenmyer et al. 2010). Although more challenging to study, these cooler planets are valuable sources of information. For starters, temperate giant planets represent the middle step between the short-period Hot Jupiters and the gas giants of our own solar system. They therefore can serve as bridges between their respective formation and migration theories (Huang, Wu & Triaud 2016). The composition of giant planets depends not only on the composition of the protoplanetary disc, but also on their location at birth and migration history. Consequently, studying their metal enrichment levels can constrain the processes driving core formation and envelope enrichment (Mordasini et al. 2016; Thorngren et al. 2016). Recent studies have also shown that long-period planets are correlated to and influence the dynamical evolution of the short-period planets within their systems (Zhu & Wu 2018; Bryan et al. 2019). Moreover, theoretical models predict that the formation of inner Earth-like planets is significantly dependent on the presence of quickly accreted cold giants (Morbidelli et al. 2022). Due to their lower effective temperatures, the atmospheres of temperate giants produce entirely different molecular abundances and potentially can contain disequilibrium chemistry by-products (Fortney et al. 2020), making long-period gas planets valuable targets for atmospheric characterization. Their atmospheres are less affected by temperature-induced inflation, which in turn allows us to use cooling models of planet evolution to constrain atmospheric metallicity (Ulmer-Moll et al. 2022). Additionally, there is a clear split in the eccentricity distribution of long-period planets. They are divided into a first group of objects with significantly high eccentricities and a second group with consistently nearly circular orbits (Petrovich & Tremaine 2016). No clear cause of this bimodality has been found yet.

The numerous and highly varied scientific interests in exo-planet detection and characterization have in the years motivated many space-based missions and ground-based instruments, including the second-generation High Accuracy Radial-velocity Planet Searcher for the Northern hemisphere spectrograph (HARPS-N, Cosentino et al. 2012) and the Spectrographe pour l'Observation des Phénomènes des Intérieurs stellaires et des Exoplanàtes (SOPHIE, Perruchot et al. 2008) spectrographs. Paired with space photometric missions (e.g. Ricker et al. 2015), the combination of transit photom-

etry and RV makes the determination of precise planetary masses and radii possible. The precision of RV surveys has been steadily improving and the current uncertainty level reaches down to the tens of centimetres per second (Jurgenson et al. 2016; Thompson et al. 2016; Pepe et al. 2021), but the biggest obstacle remains stellar variability (Fischer et al. 2016; Crass et al. 2021). Great care is required when accounting for and modelling stellar activity in order to obtain accurate orbital solutions and to accurately and precisely determine planetary masses. To do so Gaussian Process (GP) regression coupled with Monte Carlo Markov Chain parameter space exploration has been implemented in this paper and its specifics will be discussed in Section 6.

In this paper, we characterize the high proper motion, bright ($m_v$ = 8.9 mag) K5-dwarf TOI-2134 and its planetary system. We detect a multitransiting mini-Neptune in a short circular orbit and an outer temperate sub-Saturn planet. We also propose these targets for Rossiter-McLaughlin effect (RM, McLaughlin 1924; Rossiter 1924; Queloz et al. 2000) follow-up and for atmospheric characterization.

This paper is structured as follows: in Section 2 we describe the photometric and spectroscopic data used in our analysis of the system. In Section 3, we characterize the host star with four independent techniques. In Section 4, we include the analysis of the stellar signals and its activity proxies to identify the stellar rotational period. In Sections 5 and 6, we fit the photometric data for transit parameters and perform a GP regression on the RV data to determine the planets' masses, radii, and orbit characteristics. Results can be found in Tables 3 and 4. In Section 7, we combine the two data sets and perform a joint photometric and RV analysis, with results in Table 5. Final results are shown in Table 6 and addressed in Section 8, together with proposed follow-ups.

## 2 DATA

### 2.1 TESS photometry

TOI-2134, also known as TIC 75 878 355 in the TESS Input Catalog (Stassun et al. 2018), was observed by NASA's *Transiting Exoplanet Survey Satellite* (TESS, Ricker et al. 2015) mission in 2-min cadence mode over five sectors (Sectors 26, 40, 52, 53, and 54) for a total of 88 431 data points between BJD 2 459 010 and 2 459 035 (2020 June 9–July 4), BJD 2 459 390 and 2 459 418 (2021 June 24–July 22), and BJD 2 459 718 and 2 459 797 (2022 May 18–August 5). The data were originally processed by the TESS Science Processing Operation Centre (SPOC) pipeline based at NASA Ames Research Center (Jenkins et al. 2016). However, Sector 40 showed strong residual systematics after the SPOC correction, so we performed our own systematics corrections of the SPOC Simple Aperture Photometry (SAP) light curves (Twicken et al. 2010; Morris et al. 2020). In particular, we modelled the systematics as a sum of moments of the spacecraft quaternion time-series (e.g. Vanderburg et al. 2019) and modelled long-term variations with a basis spline. We also included a term for variations in the background flux in our model. We performed the model fit using an analytic linear least-squares fit, excluding transits and iterating the fit several times to remove outliers. The resulting light curve was similar to the SPOC light curve (with slightly lower scatter) in most sectors, and yielded a major improvement in the problematic Sector 40.

The transit signature of a TOI-2134b candidate was initially identified in a transit search conducted by the SPOC of Sector 26 on 2020 July 24 with an adaptive, noise-compensating matched filter (Jenkins 2002; Jenkins et al. 2010). Diagnostic tests were also conducted to help make or break the planetary nature of the
<figure><img src="image_2.png" /><figcaption>**Figure 1.** TESS normalized light curve over five sectors. 14 transits of an inner planet and a mono-transit of an outer planet can be seen and are indicated by dashed lines.</figcaption></figure>

signal (Twicken et al. 2010). The transit signatures for the TOI-2134b candidate were also detected in a search of Full Frame Image (FFI) data by the Quick Look Pipeline (QLP) at MIT (Huang et al. 2020a, b) for Sector 40. A larger transit was detected by both QLP and the SPOC in searches including Sector 52. This transit was attributed to a second planetary candidate in the system, TOI-2134c. It appears to be a mono-transit and it did not re-occur in the following 75 d. The TESS Science Office (TSO) reviewed the vetting information and issued an alert on 2020 August 7 for TOI-2134b and on 2022 July 28 for TOI-2134c (Guerrero et al. 2021). The signal for the candidate TOI-2134b was repeatedly recovered as additional observations were made in Sectors 26, 40, 52, 53, and 54, and the transit signatures passed all the diagnostic tests presented in the Data Validation reports. The difference image centroiding figure and difference images for the multisector (Sectors 26–55) run for candidate TOI-2134b show that the centroid of the transit source is consistent with the target star of interest. The host star is located within $3.2 \pm 3.7$ arcsec of the source of the transit signal for candidate TOI-2134b and within $0.98 \pm 2.59$ arcsec of the source of the transit signal for candidate TOI-2134c. We flattened the light curve by simultaneously fitting transit models for the two planets along with a basis spline to model long-term variations, and then subtracting the long-term variations (a strategy similar to Vanderburg et al. 2016, except without a simultaneous systematics model; see also Pepper et al. 2020). The systematics-corrected and flattened TESS data are shown in Fig. 1. To better constrain the characteristics of the mono-transiting long-period planet candidate, we launched a ground- and space-based photometric observing campaign to catch a second transit.

## 2.2 LCOGT photometry

The Las Cumbres Observatory Global Telescope (LCOGT, Brown et al. 2013) network observed the star between BJD 2459 808 and 2459 818 (2022 August 17–27), when preliminary ephemeris prediction suggested the outer planet would re-transit.

Due to an unfortunate combination of bad weather and low visibility, only a possible egress was detected. However, the LCO 0m4 SBIG detectors are very susceptible to strong systematics and several combinations of comparison stars and aperture sizes need to be examined to assess the overall reliability of a light-curve feature, especially for ingress- or egress-only events. When using a different choice of comparison stars, a convincing egress was no longer present in the data. The apparent egress was, in fact, proven to be highly dependent on the choice of comparison star set. For this reason, we could not claim this egress as a detected transit on its own and we do not include this data in our analysis.

We also attempted a TRansiting Planets and PlanetesImals Small Telescope North (Barkaoui et al. 2017) observation of the outer planet on 2022 August 22, but it was unsuccessful.

## 2.3 NEOSSat photometry

The position in the sky of TOI-2134 is such that it is not observable after late-October, which precluded the chance of a second ground-based campaign to detect a third transit of the outer planet candidate since the TESS detection. We therefore turned to space observations. TOI-2134 is outside of the CHaracterising Exoplanets Satellite (CHEOPS) field of view, but it is visible to the agile space telescope Near Earth Object Surveillance Satellite (NEOSSat, Hildebrand et al. 2004; Fox & Wiegert 2022). NEOSSat is a Canadian microsatellite orbiting the Earth in a Sun-synchronous orbit of approximately 100 min. It was originally deployed to study near-Earth satellites, but it also performs well for follow-up observations of large exoplanets transiting bright stars. It carries a 15-cm f/6 telescope, with spectral range between 350 and 1050 nm and a field of view of $0.86^\circ \times 0.86^\circ$.

NEOSSat observed TOI-2134 unevenly between BJD 2459 898 and 2459 910 (2022 November 14–26) with a 70 s cadence for a total of 3364 data points. Multiple sets of observations through the run show significant unpredictable offsets that are usually corrected with calibration on reference stars. In these orbits, however, the reference stars behave differently from each other and the correction is less precise. This is probably due to image artefacts, as the detector and readout process have quite noticeable imperfections. These high-variance orbits have been flagged in the data set and appear often enough to prevent a clear confirmation of a transit.

## 2.4 WASP photometry

TOI-2134 was also observed over 3 yr by the Wide Angle Search for Planets (WASP, Pollacco et al. 2006; Wilson et al. 2008) with coverage of about 120 nights per year. The data cover similar three-month spans between BJD 2454 580 to 2454 690 (2008 April 23–August 11), BJD 2454 941 to 2455 067 (2009 April 19–August 23), and BJD 2455 307 to 2455 432 (2010 April 20–August 23). A total of 23 097 data points were obtained and reduced with the SuperWASP pipeline (Pollacco et al. 2006). No planetary transit was detected. However, the long baseline, over three years long, allows for long-term monitoring of the stellar activity and of the rotational period of the host star, as shown in Section 4.1. All data are shown in Fig. 2.

## 2.5 HARPS-N spectroscopy

We collected a total of 111 RV observations of TOI-2134 over two seasons with the HARPS-N (Cosentino et al. 2012, 2014)
<figure><img src="image_3.png" /><figcaption>**Figure 2.** WASP normalized flux against Julian Date over the three years of coverage. All data points are plotted with error bars, and daily averages are overplotted. The predicted transits of TOI-2134c are plotted as dashed lines, while their uncertainties are plotted as shaded areas. As we address in Section 8, we did not detect any transit.</figcaption></figure>

installed on the 3.6-m Telescopio Nazionale Galileo (TNG) at the Observatorio del Roque de Los Muchachos in La Palma, Spain. HARPS-N is an updated version of HARPS at the ESO 3.6-m (Mayor et al. 2003). The spectrograph covers the wavelength range of 383–691 nm, with an average resolution $R = 115\,000$. The first 32 spectra were collected between BJD 2459417 and 2459515 (2021 July 21–October 27), and the next 79 were collected between BJD 2459638 and 2459890 (2022 February 27–November 6). All data were observed under the Guaranteed Time Observations programme with the standard observing approach of one observation per night. The average exposure time for TOI-2134 was 900s with an average signal-to-noise ratio (SNR) at 550 nm of ~100. RVs and activity indicators were extracted using the 2.3.5 version of the Data Reduction Software (DRS) adapted from the ESPRESSO pipeline (see Dumusque et al. 2021) and computed using a K6-type numerical weighted mask. The RV data show a peak-to-peak dispersion of 35 m s⁻¹, with standard rms of 7.3 m s⁻¹ and mean uncertainty of 0.7 m s⁻¹.

Several proxies are extracted by the standard DRS pipeline, including (but not limited to) the full width at half-maximum (FWHM) and the contrast of the cross-correlation function (CCF), and the S-index. The mentioned data are plotted in purple in Fig. 3. The reasoning behind the selection of plotted proxies is addressed in Section 4.2.

## 2.6 SOPHIE spectroscopy

We also obtained 113 RV observations of TOI-2134 with the SOPHIE (Perruchot et al. 2008) between BJD 2459082 and 2459894 (2020 August 20–2022 November 10). SOPHIE is a stabilized échelle spectrograph dedicated to high-precision RV measurements in optical wavelengths (387–694 nm) on the 193-cm Telescope at the Observatoire de Haute-Provence, France (Bouchy et al. 2009). We used the SOPHIE high-resolution mode (resolving power $R = 75\,000$) and the fast mode of the CCD reading. The standard stars observed at the same epochs using the same SOPHIE mode did not show significant instrumental drifts. Depending on the weather conditions, the exposure times for TOI-2134 ranged from 4.5 to 30 min (average of 11 min) and their SNR per pixel at 550 nm ranged from 21 to 77 (average of 54). Five exposures showed an SNR below 40 and were removed. The final data set therefore includes 108 epochs.

The RV data were extracted with the standard SOPHIE pipeline using CCFs (Bouchy et al. 2013) and including the CCD charge transfer inefficiency correction. The cross-correlations were made

using several numerical masks, characteristic of different types of stars. All produced similar results in terms of RV variations. We finally adopted the RVs derived using a K5-type mask, which provided the least dispersed results.

Following the method described, for example, in Pollacco et al. (2008) and Hébrard et al. (2008), we estimated and corrected for the sky background contamination (mainly due to the Moon) using the second SOPHIE fibre aperture, which is targeted 2 arcmin away from the first one pointing toward the star. We estimated that 14 of the 108 exposures were significantly polluted by sky background, each time implying a correction below 10 m s⁻¹. The final SOPHIE RVs show variations with a dispersion of 8.2 m s⁻¹ (35 m s⁻¹ peak to peak), significantly larger than their typical 2 m s⁻¹ precision. The FWHM, bisector span, and contrast of the CCF were also derived for every observation. The data are plotted in orange in Fig. 3 (for more information on proxy selection, see Section 4.2).

## 3 STELLAR CHARACTERIZATION

TOI-2134 is a bright, high-proper motion, mid-K-dwarf. As the star falls into a parameter space that is not optimal for several of the common stellar characterization pipelines, we characterized the system with multiple separate and independent methods.

### 3.1 Spectral energy distribution analysis

We estimated stellar luminosity $L_*$, effective temperature $T_{\text{eff}}$, and stellar radius $R_*$, by fitting the stellar energy distribution (SED) of TOI-2134 following the method of Mann et al. (2015), and using templates instead of the observed spectrum, as described in Mann et al. (2016). To briefly summarize, we compared available photometry [*Gaia*, 2MASS (Two-Micron All-Sky Survey), *Tycho*, and *WISE*] of the host star to a grid of flux-calibrated spectral templates from Rayner, Cushing & Vacca (2009) and Gaidos et al. (2014). We filled gaps in the spectral templates using PHOENIX BT-SETTL models from Allard et al. (2013), which also provide an estimate of $T_{\text{eff}}$. We computed the bolometric flux, $F_{\text{bol}}$, by integrating the output absolutely calibrated spectrum along wavelength. This gave us $L_*$ when combined with the *Gaia* Data Release 3 (DR3) parallax, which in turn gave us $R_*$ when combined with our estimate of $T_{\text{eff}}$ using the Stefan–Boltzmann relation. We did not correct for the offset in the *Gaia* DR3 parallax (Lindegren et al. 2021), but this effect is much smaller than the systematic uncertainties intrinsic to the rest of the analysis.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Plots of the HARPS-N and SOPHIE RV data alongside the chosen activity proxies for each data set (see Section 4.2). From the top: HARPS-N RVs, S-index, FWHM, and contrast, followed by SOPHIE RVs and their bisector span (BISS). Notice the different time axes. All error bars are plotted, but some are too small to be clearly visible in HARPS-N data. We plot only the activity proxies used in the later analysis (for more information, see Section 4.2).</figcaption></figure>

<table><caption>**Table 1.** Stellar parameters derived using the different techniques, addressed in order in Sections 3.1–3.3.</caption><thead><tr><th>Parameter</th><th>SED vs1</th><th>ARES + MOOG</th><th>SPC</th><th>SED vs2</th></tr></thead><tbody><tr><td>F<sub>bol</sub> [erg cm<sup>2</sup> s<sup>-1</sup>]</td><td>1.198 ± 0.048</td><td></td><td></td><td></td></tr><tr><td>L<sub>*</sub> [L<sub>☉</sub>]</td><td>0.192 ± 0.009</td><td></td><td></td><td>0.190<sup>+0.021</sup><sub>-0.022</sub></td></tr><tr><td>T<sub>eff</sub> [K]</td><td>4630 ± 90</td><td>4620 ± 80</td><td>4600 ± 50</td><td>4490<sup>+60</sup><sub>-70</sub></td></tr><tr><td>Radius [R<sub>☉</sub>]</td><td>0.683 ± 0.027</td><td>0.714<sup>+0.017</sup><sub>-0.028</sub></td><td></td><td>0.721<sup>+0.020</sup><sub>-0.021</sub></td></tr><tr><td>log(g) [cm s<sup>-1</sup>]</td><td></td><td>4.8 ± 0.3</td><td>4.7 ± 0.1</td><td>5.4<sup>+0.1</sup><sub>-0.5</sub></td></tr><tr><td>[Fe/H]</td><td></td><td>0.13 ± 0.04</td><td>0.09 ± 0.08</td><td>0.1<sup>a</sup></td></tr><tr><td>Mass [M<sub>☉</sub>]</td><td>0.70 ± 0.04</td><td>0.76<sup>+0.04</sup><sub>-0.02</sub></td><td></td><td>0.75<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>Microturbulence ξ<sub>t</sub> [km s<sup>-1</sup>]</td><td></td><td>0.18 ± 0.12</td><td></td><td></td></tr><tr><td>Density [ρ<sub>☉</sub>]</td><td>2.20 ± 0.63</td><td>2.11<sup>+0.06</sup><sub>-0.10</sub></td><td></td><td>1.99<sup>+0.25</sup><sub>-0.20</sub></td></tr><tr><td>Age [Gyr]</td><td></td><td>3.8<sup>+5.5</sup><sub>-2.7</sub></td><td></td><td>2<sup>a</sup></td></tr><tr><td>Distance [pc]</td><td></td><td>22.646 ± 0.015</td><td></td><td>22.657<sup>+0.006</sup><sub>-0.009</sub></td></tr></tbody></table>

*Note.*<sup>a</sup>Set as constant in the model.
More details on the uncertainties are given in Mann et al. (2015). To briefly summarize, uncertainties are incorporated as part of a Monte Carlo framework; we generate a grid of fits by sampling over the choice of template (including interpolating between templates), adjustments to the spectral shape (flux calibration uncertainties), as well as reported uncertainties in the parallax, spectra, and photometry. Two irreducible systematic effects were added separately. The first was for $T_{\text{eff}}$ and is based on comparing model-based temperatures to more empirical estimates from long-baseline optical interferometry (Mann, Gaidos & Ansdell 2013). The second was based on calibration of the zero-points and filter profiles (Mann & von Braun 2015; Maíz Apellániz & Weiler 2018). The final values are shown in Table 1 under the SED vs1 column.

As part of the analysis, we derived another estimate of $R_*$ based on the scale factor between the models and the absolutely calibrated spectrum. This scale factor is $\propto R_*^2/D_*^2$, where $D_*$ is the distance to the star. We combined it with the Gaia parallax to estimate $R_*$. This effectively is the infrared (IR)-flux method (Blackwell & Shallis 1977), and yielded $R_* = 0.700 \pm 0.028 R_\odot$, consistent with our Stefan-Boltzmann fit ($R_* = 0.683 \pm 0.027 R_\odot$).

### 3.1.1 Stellar mass from $M_{KS} - M_*$ relation

We estimated the mass of the host star using the relation between $K$ magnitude and mass, $M_{KS}$ and $M_*$, from Mann et al. (2019). This relation was built using orbits of astrometric binaries, making it empirical. Using $K_S$ photometry from the 2MASS (Skrutskie et al. 2006) and the Gaia DR3 parallax, we obtained $M_* = 0.702 \pm 0.018 M_\odot$. This $M_*$ value placed the host star at the edge of the Mann et al. (2019) relation, where errors may be underestimated due to a lack of Sun-like stars in the sample and the effects of stellar evolution. We, therefore, adopted a more realistic 5 per cent uncertainty, as shown in Table 1.

## 3.2 ARES + MOOG with isochrone fitting and SPC

We also measured stellar atmospheric parameters directly from the HARPS-N spectra. For this purpose the one-dimensional spectra were shifted to the lab frame with the DRS RVs and then co-added. The resulting spectrum had an SNR of about 600. We employed the ARES + MOOG<fnref n="1" /> method to measure the effective temperature, surface gravity, microturbulence, and iron abundance (used as a proxy for metallicity). We used the method through the FASMA<fnref n="2" /> implementation (Andreasen et al. 2017). It relies on calculating the equivalent widths of a set of isolated iron lines (taken from Tsantaki et al. 2013) and using them in the radiative transfer code MOOG (Sneden 1973) to obtain the atmospheric parameters by imposing excitation and ionization equilibrium. The stellar atmospheric models were taken from Kurucz (1993). Some iron lines were discarded as they gave equivalent-width measurements that were unreasonably large (>200 mÅ) or small (<5 mÅ). We also fixed the microturbulence following Tsantaki et al. (2013). Finally, we inflated the errors for accuracy and corrected the surface gravity following Mortier et al. (2014). The final values of $T_{\text{eff}}$, surface gravity log(g), metallicity [Fe/H], and microturbulence $\xi_t$ are shown in Table 1 under the ARES + MOOG column.

After obtaining these atmospheric parameters, we used the code ISOCHRONES (Morton 2015) to derive mass, radius, age, and distance. We ran the code four times, varying the inputs as well as the used stellar models. The common inputs for all four runs were the Gaia DR3 parallax, and the photometric magnitudes in bands B, V, J, H, and K. For two runs, we also included the effective temperature and metallicity as measured from the HARPS-N spectra. We chose not to use the spectroscopic surface gravity given its known accuracy issues (see e.g. Mortier et al. 2014). We used two stellar models (each in two runs): the Dartmouth Stellar Evolution Database (Dotter et al. 2008) and the Mesa Isochrones and Stellar Tracks (Dotter 2016). For our final results, we combined the posterior distributions of all four runs. To combine the posteriors we added them together and corrected for the sample size (as in Borsato et al. 2019). We extracted the median and 16th and 84th percentiles as the final value and its errors, as reported in Table 1.

### 3.2.1 SPC pipeline

We also derived stellar parameters using the Stellar Parameter Classification pipeline (SPC, Buchhave et al. 2012, 2014). The high SNR needed to extract precise RVs means that these spectra are more than adequate for deriving stellar parameters. We ran the SPC analysis on each individual spectrum and calculated the weighted average of the individual spectra. The weights are computed from the normalized CCF peak heights from the observed spectrum and the best-matched template (model) spectrum. Higher CCF peaks indicate a better match between the model and the observations. The normalization leads to a CCF peak height of 1 for autocorrelation. While the SNR of the observed spectra could also be used as the weighting factor, the CCF peak height better incorporates the relationship between data and model. The results are shown in Table 1 under the SPC column. We also computed $v\sin(i) < 2 \text{ km s}^{-1}$. The formal uncertainties take into account the model uncertainties, which primarily stem from systematics in the ATLAS Kurucz stellar models and degeneracies between the derived parameters when trying to compare observed spectra to model spectra (see Buchhave et al. 2012, 2014). The parameters from SPC agree well with the results from ARES + MOOG within the uncertainties.

## 3.3 Spectral energy distribution analysis with isochrone fitting

We have also computed an estimate of $R_*$ and $T_{\text{eff}}$ using the SED fitting method presented in Morrell & Naylor (2019, 2020). This method compares multiband photometry placed across the stellar SED with synthetic photometry, generated from the BT-SETTL Cosmological Impact of the First STars (CIFIST) (Allard, Homeier & Freytag 2012) atmosphere grid, and diluted using the distances of Bailer-Jones et al. (2021). By best matching the area beneath the SED and the overall shape of the SED, we determined the luminosity $L_{\text{SED}}$ and temperature $T_{\text{SED}}$ respectively – which together unambiguously define $R_*$. Unlike the method presented in Section 3.1, which makes use of spectroscopic templates for the measurement of $T_{\text{eff}}$, this method self-consistently measures both $T_{\text{eff}}$ and $R_*$ using only photometry and distances, effectively providing an alternate measure of temperature to the other methods.

For this fitting, we used the $G_{\text{BP}}$ and $G_{\text{RP}}$ bands from Gaia DR3 (Gaia Collaboration 2016, 2023), the J, H, and K bands from 2MASS (Skrutskie et al. 2006), and the W1, W2, and W3 bands from All Wide-field Infrared Survey Explorer (AllWISE) (Wright et al. 2010). As with Morrell & Naylor (2019), we adopted a floor value of 0.01 mag,

<footnote n="1">¹ARESV2: http://www.astro.up.pt/~sousasag/ares/;
MOOG 2017: http://www.as.utexas.edu/~chris/moog.html</footnote>

<footnote n="2">²FASMA: http://www.iastro.pt/fasma/index.html</footnote>
corresponding to about 1 per cent, for the photometric uncertainty for all bands. The parameters resulting from our fitting are shown in the SED vs2 column in Table 1. At first glance, $T_{\text{eff}}$ and $R_*$ from this method are inconsistent with the other determinations. However, the resulting $L_*$ from these parameters is consistent with that described in Section 3.1, supporting the validity of both sets of parameters. We considered the possibility of extinction contributing to the aforementioned difference, however the star is close enough that extinction should be negligible. Moreover the extinction required to match the results of Section 3.1 is 0.1, which is too large to be probable. Furthermore, the measurement of $R_*$ using this method is consistent with the secondary, IR flux-based method determination from Section 3.1. From our study, the $G_{\text{BP}}$ and $G_{\text{RP}}$ bands appear to be sampling a redder SED than the bands at longer wavelengths, resulting in a cooler measured $T_{\text{SED}}$. Given that the photometric data were not contemporaneous, with the visible and IR photometry being 5–10 yr separated, it is possible for the observed SED to have changed over this intervening period. Though, as we can find no quality issues or physical reason for this discrepancy, the fitting for our parameter determinations for this section did employ the $G_{\text{BP}}$ and $G_{\text{RP}}$ bands.

We then determined the stellar mass $M_*$ using the PAdova and TRieste Stellar Evolution Code (PARSEC) 1.2S isochrones (Bressan et al. 2012; Chen et al. 2014; Tang et al. 2014; Chen et al. 2015; Marigo et al. 2017; Pastorelli et al. 2019). We used CMD $3.7^3$ to generate evolution tracks at a metallicity of $[\text{M/H}] = 0.1$, which is in line with the value determined in Section 3.2. Given that the ARES + MOOG age estimation places the star on the main sequence, we interpolated the 2 Gyr isochrone to estimate the $M_*$ at our measured $L_*$ and its uncertainty bounds, also shown in Table 1. We note that, due to not having access to the posterior for distance and instead just assuming it to be Gaussian, the uncertainty bounds for luminosity, mass, and stellar density from this method are likely to be overestimated.

Overall, all analysis agree with each other within their uncertainties. For the scope of this work, we characterized TOI-2134 via the mean of all the computed values weighted by the inverse of their errors, as compiled in Table 2. Their uncertainties are computed as the standard deviation between measurements in each method, to avoid improper averaging down of systematic effects.

## 4 STELLAR ACTIVITY SIGNAL

We conducted a thorough preliminary analysis of the available data in order to search for and to best characterize the stellar activity-induced signals in both the photometric and the spectroscopic observations.

To begin with, the projected rotational velocity $v\sin(i)$ of TOI-2134 was determined to be $<2 \text{ km s}^{-1}$ from the HARPS-N spectra (as mentioned in the previous section), and $1.5 \pm 1.0 \text{ km s}^{-1}$ from the SOPHIE CCFs (following the method in Boisse et al. 2010). No more precise measurement could be derived from the spectra. We therefore calculated a minimum stellar rotation period $P_{\text{rot, min}}$ associated to the lower maximum limit of $v\sin(i)$ as:

$$P_{\text{rot,min}} = \frac{2\pi R_*}{v\sin(i)} \approx 23 \text{ d.} \quad (1)$$

Using the method described in Noyes, Weiss & Vaughan (1984), we computed the average $\log R'_{\text{HK}}$ to be $-4.83 \pm 0.45$ from the S-index measurements taken by HARPS-N. There was significant scatter in the S-index measurements which degraded the quality of

<table><caption>**Table 2.** Stellar parameters of TOI-2134.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Source</th></tr></thead><tbody><tr><td>Name</td><td>TOI-2134<br>TIC 75 878 355<br>G204-45</td><td>TESS project<sup>a</sup><br>Stassun et al. (2019)<br>Giclas, Burnham &<br>Thomas (1979)</td></tr><tr><td>RA [h:m:s]</td><td>18:07:44.52</td><td>Gaia Collaboration<br>(2020)</td></tr><tr><td>Dec. [d:m:s]</td><td>+ 39:04:22.54</td><td>Gaia Collaboration<br>(2020)</td></tr><tr><td>Spectral type</td><td>K5V</td><td>Stephenson (1986)</td></tr><tr><td><i>m</i><sub>V</sub> [mag]</td><td>8.933 ± 0.003</td><td>TESS project<sup>a</sup></td></tr><tr><td><i>m</i><sub>J</sub> [mag]</td><td>6.776 ± 0.023</td><td>TESS project<sup>a</sup></td></tr><tr><td><i>m</i><sub>K</sub> [mag]</td><td>6.091 ± 0.017</td><td>TESS project<sup>a</sup></td></tr><tr><td>(<i>B</i> − <i>V</i>) [mag]</td><td>1.192 ± 0.033</td><td>TESS project<sup>a</sup></td></tr><tr><td>Parallax [mas]</td><td>44.1087 ± 0.0144</td><td>Gaia Collaboration<br>(2020)</td></tr><tr><td>Distance [pc]</td><td>22.655 ± 0.007</td><td>This work</td></tr><tr><td>Proper motion<br>[mas yr<sup>-1</sup>]</td><td>288.257 ± 0.016</td><td>Gaia Collaboration<br>(2020)</td></tr><tr><td><i>L</i><sub>*</sub> [L<sub>☉</sub>]</td><td>0.192 ± 0.008</td><td>This work</td></tr><tr><td><i>F</i><sub>bol</sub> [erg cm<sup>2</sup> s<sup>-1</sup>]</td><td>1.198 ± 0.048</td><td>This work</td></tr><tr><td><i>T</i><sub>eff</sub> [K]</td><td>4580 ± 50</td><td>This work</td></tr><tr><td>log(g) [cm s<sup>-1</sup>]</td><td>4.8 ± 0.3</td><td>This work</td></tr><tr><td>[Fe/H]</td><td>0.12 ± 0.02</td><td>This work</td></tr><tr><td>Mass [M<sub>☉</sub>]</td><td>0.744 ± 0.027</td><td>This work</td></tr><tr><td>Radius [R<sub>☉</sub>]</td><td>0.709 ± 0.017</td><td>This work</td></tr><tr><td>Density [ρ<sub>☉</sub>]</td><td>2.09 ± 0.10</td><td>This work</td></tr><tr><td>Age [Gyr]</td><td>3.8<sup>+5.5</sup><sub>-2.7</sub></td><td>This work</td></tr><tr><td>vsin(i) [km s<sup>-1</sup>]</td><td>0.78 ± 0.09</td><td>This work</td></tr><tr><td>&lt; log R'<sub>HK</sub> &gt;</td><td>-4.83 ± 0.45</td><td>This work</td></tr><tr><td><i>P</i><sub>rot</sub> [d]</td><td>45.78<sup>+5.56</sup><sub>-5.31</sub></td><td>This work</td></tr></tbody></table>

Note.<sup>a</sup> See ExoFOP: https://exofop.ipac.caltech.edu/tess/target.php?id = 75878355

the results, but the empirical relations of Noyes et al. (1984) yielded a stellar rotation period of $\sim 42$ d.

To better identify the stellar rotational period, we performed a periodogram analysis.

### 4.1 Photometry

We computed the Bayesian Generalised Lomb–Scargle (BGLS) periodograms (Mortier et al. 2015) for both the WASP and the TESS photometric data, shown, respectively, in green and blue in the first and second rows of Fig. 4. The same periodograms in frequency space, alongside their window functions are shown in Fig. 5. The TESS data showed a forest of peaks at $\sim 9.2$ d (highlighted by a black dashed line), which is generated by the repeated transits of the inner planet. As expected given the detection of no transits due to lower precision, the WASP periodogram had no power around this period. It instead showed two significant forests of peaks centred around $\sim 29$ and $\sim 58$ d (shown as blue bands in Fig. 4), which were originally attributed to the stellar rotational period, but could also be generated by the moon cycle. To further investigate this, we also plotted the BGLS periodograms of each yearly season of WASP, as shown in the first row of Fig. 4 as blue, red, and purple dashed lines. The BGLS periodograms of the two later years also presented a significant peak at 58 d, but the 2008 data did not. Instead, its most significant peak was at 29 d. A peak at $\sim 29$ d was also present in 2010, but not in 2009. While some of the discrepancies could be attributed to differing coverage, these result hinted at either a different lunar contribution over the different seasons, or at evolving surface inhomogeneities structure trends over the years, possibly

<footnote n="3">³CMD 3.7: http://stev.oapd.inaf.it/cgi-bin/cmd</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Set of BGLS periodograms of the acquired data plotted as period versus logarithmic probability normalized to 1. From the top, WASP photometry in the solid green with yearly seasons in blue for 2008, red for 2009, and purple for 2010 as dashed lines, *TESS* photometry, HARPS-N RVs, HARPS-N activity proxies (FWHM, *S*-index, and contrast in, respectively, blue, red, and green), SOPHIE RVs, SOPHIE activity proxy (BISS), and the combined SOPHIE and HARPS-N RVs. The dashed vertical dashed lines represent the periods of the two planets at 9.2 and 95 d. The shaded vertical bands indicate the possible stellar rotational signals at 29, 48, and 58 d.</figcaption></figure>

related to a stellar magnetic cycle. After alias analysis, we found that the 29 d forest of peaks in the full periodogram can be explained as the extended aliases generated by the 365 d period. The WASP data span over ~850 d. SOPHIE RVs (taken 10 yr later) also cover a similar stretch of time. Therefore, assuming these signals are stellar, we can expect the structure of surface inhomogeneities that allow us to detect stellar rotational period in periodogoram analyses to also evolve during the three years of RV data. This evolution could be the reason behind the difficulties constraining the stellar rotational period in the further RV analyses.

## 4.2 Radial-velocity data and proxies

We conducted a full periodogram analysis of the spectroscopic data. The last five rows of Fig. 4 show the BGLS periodograms of, in order, the HARPS-N RVs, the HARPS-N derived proxies (FWHM, *S*-index, and contrast), the SOPHIE RVs, the SOPHIE derived activity indicator (bisector span, or BISS), and the combined RV data. We were able to combine the RVs with a simple offset, as they are derived from similar wavelength windows and therefore are probing the same section of the stellar spectra. The same periodograms in frequency space, alongside their window functions, are once again shown in Fig. 5.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Same set of BGLS periodograms as Fig. 4 in frequency space. The window functions for each data set are also included. The dashed vertical lines indicate the periods of the two planet candidates. The shaded vertical band shows the stellar rotational period.</figcaption></figure>

Although the star showed significant variation in the activity indicators, and the average log R'HK also classified the star as moderately active, both sets of RVs had little to no correlation to their

activity indicators. The specific reason for this lack of correlation is ultimately beyond the scope of this paper, as the activity indicators were only used as a starting point to the analysis, but we propose some
possible origins. As a first most likely option, the Keplerian signals introduced by the planets in the system are large enough to 'muddle' the correlation to activity indicators. In this case, the RV amplitude of the stellar activity computed in the next sections is shown to be comparable to the amplitude of the RV signals generated by the planets. It is likely that these signals are significant enough to prevent a clean correlation between RVs and activity indicators (which only map the variations induced by stellar activity). To test this, we also computed the correlation between the activity indicators and the RVs after subtracting the best-fitting Keplerian models computed in Section 8. While the correlation did improve by a factor of 2, they still remained low. So other reasons may be considered. As an example, the stellar rotation axis inclination angle with respect to the observer can influence the strength of this correlation, weakening it for unfavourable line of sights: as the stellar rotational axis becomes parallel to the observer line of sight, the signal from active regions coming in and out of view becomes less rotationally modulated. At the same time, in late K-dwarfs convective redshift may in some cases prevail against blueshift. This can happen either due to an opacity effect (like in M-type stars), or if most of the photospheric absorption lines used for RV measurements form in regions of convective overshoot (Norris et al. 2017). Costes et al. (2021) note that a possible explanation for low correlation between RVs and activity proxies, as is the case for our target, is that the convective blue- and redshifts are 'cancelling' one another. The possibility of a temporal lag (Collier Cameron et al. 2019) between the RVs and the proxies was also considered, but a visual inspection of their time-series did not strongly support this possibility.

For our analysis, we nevertheless selected and plotted the indicators with the strongest correlation to their RVs. For HARPS-N, we selected the S-index, the FWHM, and the contrast. Their Spearman's rank correlation coefficients with the RVs were computed to be 0.15, 0.11, and -0.12-, respectively. For SOPHIE, we selected only the bisector span, with correlation coefficient of -0.16, as the FWHM and contrast seem to be affected by instrumental systematics.

While the BGLS periodograms of the RV data sets did not show clear peaks for the inner planet, there was a strong periodic signal at ~95 d (shown as a black dashed line) shared between HARPS-N and SOPHIE RVs that was not present in any of the HARPS-N stellar activity proxies. The SOPHIE bisector does have a peak at ~100 d, but its normalized logarithmic probability is comparable to most other peaks in the periodogram and therefore does not have a strong relevance. This preliminary analysis suggested a period of ~95 d for the mono-transiting planet detected by TESS. This signal presented minorly relevant yearly aliases at 129 and 75 d in both the HARPS-N only and the combined data, which could be easily discarded in the analysis of the periodogram. No statistically significant yearly aliases arise for the 95 d signal in the SOPHIE data. The only major peak of both HARPS-N RVs and of all its activity indicators was centred around 48 d (shown as a blue band). In the HARPS-N data, we could also see some of the yearly aliases of this signal, at 42 and 38 d. This peaks were only moderately relevant and could be easily identified. No such signal can be found in either SOPHIE RVs or its indicator. On the other hand, SOPHIE data presented a minor peak at ~24 d, half of the HARPS-N value. This disagreement could be due to the different sampling and observing strategies between the two observatories. Further alias analysis showed that 24 d was also a yearly alias of 48 d. The 48 d period, although not in perfect agreement, is compatible with the longer modulation in the WASP data, especially given the fact that the data in each season only span just more than twice this period.

<figure><img src="image_7.png" /><figcaption>**Figure 6.** $\ell$1 periodograms of from top to bottom HARPS-N, SOPHIE, and combined RVs. The periods of the major identified signals are highlighted.</figcaption></figure>

To further analyse the signals within the spectroscopic data sets, we have also included an $\ell$1 periodogram$^4$ analysis with correlated noise (Hara et al. 2017; Hara & Mari 2021), as shown in Fig. 6. This periodogram formulation was first devised to overcome the distortions in the residuals that arise when fitting planets one by one, and can help isolate the most relevant signals in a data set. Once again, HARPS-N and SOPHIE RVs on their own, as well as their combination, all showed a clear peak at ~95 d. Similarly, the $\ell$1 periodograms of HARPS-N and SOPHIE both also peaked at ~9.2 d. The $\ell$1 periodogram is also able to isolate the signal of the inner planet in the combined RV data set. Regarding the possible stellar rotation period, HARPS-N data again showed a clear modulation at ~48 d, while the strongest peak in SOPHIE not attributed to planetary signals was at half that value. The $\ell$1 periodograms have therefore re-confirmed the previous results from the BGLS analysis and have allowed for a clearer understanding of the SOPHIE data.

Finally, to probe the coherence of these signals, we plotted the Stacked BGLS periodograms (Mortier et al. 2015; Mortier & Collier Cameron 2017) of the three sets of RV data in Fig. 7. The Stacked BGLS periodogram was developed to better identify the signals that are generated by stellar activity. Planetary signals are coherent in nature, meaning their probability should consistently increase with increasing number of observations. Signals produced by stellar activity are incoherent, meaning that their probability will change and oscillate. Fig. 7 clearly showed that the signals indicated by the blue vertical lines (respectively, 24 and 48 d, as identified by the $\ell$1 periodograms) were incoherent. They therefore could not

<footnote n="4">⁴Available at https://github.com/nathanchara/11periodogram</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** From the top, stacked BGLS periodograms of HARPS-N, SOPHIE, and combined RVs. The blue dashed blue lines identify 48 d (and 25 d for SOPHIE data). The dashed grey lines show 9.23 and 95 d, the proposed periods of the two planets.</figcaption></figure>

be attributed to planets and were more likely generated by stellar activity. The 9.2 and 95 d signals, highlighted by the grey dashed lines, showed more coherent trends. With the exception of a 1 d alias, no other major signals could be identified.

## 5 TRANSIT PHOTOMETRY

We then performed an analysis to determine the best-fitting transit parameters and uncertainties for the two planet candidates orbiting TOI-2134. We modelled the TESS photometry (after systematics correction and flattening as described in Section 2.1) with Mandel & Agol (2002) transit models. Our model included three parameters describing the host star (its mean density, and both linear and quadratic $q_1$ and $q_2$ limb-darkening coefficient parametrizations sampled following Kipping 2013). The inner planet TOI-2134b was described by six parameters (its orbital period, time of transit, orbital inclination, the logarithm of the planet/star radius ratio $\log R_p/R_\star$, and combinations of the eccentricity and argument of periastron of the planet $\sqrt{e} \cos \omega_p$ and $\sqrt{e} \sin \omega_p$, which will be further explained in Section 6). The transit of the outer planet TOI-2134c was described by four parameters (time of transit, transit duration, impact parameter, and the logarithm of the planet/star radius ratio). Finally, we included two parameters characterizing the data set itself (a constant flux offset and the white noise level).

### 5.1 Selection of priors

We imposed an informative Gaussian prior on the stellar density based on our analysis of the stellar parameters. All other parameters

were bound by uniform priors. We restricted the inclination of planet b to be less than 90° and the impact parameter of planet c to be greater than 0 (to avoid the degeneracy for transit configurations with inclinations greater 90°). We restricted $\sqrt{e} \cos \omega_p$ and $\sqrt{e} \sin \omega_p$ to be in the interval $[-1, 1]$ (as necessary as per their definition), and the impact parameters (in the case of TOI-2134b after conversion from inclination) to be in the range $[0, 1 + R_p/R_\star]$ (requiring the planets transit the star). $\log R_p/R_\star$ was allowed to vary in the range $[-\infty, 0]$ (planets must be smaller than the host star), and $q_1$ and $q_2$ in the range $[0, 1]$ following Kipping (2013). All other parameters with uniform priors were allowed to explore the range $[-\infty, \infty]$.$^5$

### 5.2 Transit results

We explored the parameter space using a Markov Chain Monte Carlo (MCMC) algorithm with a Differential Evolution sampler (Ter Braak 2006). We simultaneously evolved 100 chains for 100 000 steps each, discarding the first 30 000 as burn-in. We assessed convergence by calculating the Gelman–Rubin statistic and found values less than 1.006 for all parameters. Our best-fitting models are phase-folded and plotted in Fig. 8 and the results of our planetary fit are given in Table 3. We chose to initially not derive eccentricity, angle of periastron, and period for the outer planet candidate, given the mono-transit. Those parameters will be extracted in a second step we discuss in Section 6.2. The multiple transits of the inner planet allow us to precisely measure its period and planet-to-star radius ratio. The radius ratio of TOI-2134c is also constrained to over 100$\sigma$.

## 6 RADIAL-VELOCITY ANALYSIS

To analyse the RVs, we used the new code MAGPy_RV.$^6$ MAGPy_RV is a pipeline for GP regression with an affine invariant MCMC parameter space searching algorithm (as defined in Foreman-Mackey et al. 2013).

GPs have been extensively employed in astrophysical literature to successfully model stellar activity-induced variations and instrumental noise in both RV and photometric measurements (e.g. Haywood et al. 2014; Rajpaul et al. 2015; Faria et al. 2016; Serrano et al. 2018; Barros et al. 2020).

We modelled the RV data as a combination of two planetary signals in the form of Keplerians (for the two transiting objects), and the stellar activity in the form of a quasi-periodic kernel. We selected the quasi-periodic kernel defined in Haywood et al. (2014) with the inclusion of a white noise ‘jitter’ term, in the form

$$k(t_n, t_m) = \theta_1^2 \exp \left[ -\frac{|t_n - t_m|^2}{\theta_2^2} - \frac{\sin^2 \left( \frac{\pi |t_n - t_m|}{\theta_3} \right)}{\theta_4^2} \right] + \delta_{n,m} \beta^2, (2)$$

in which $t_n$ and $t_m$ are two data points, the four hyperparameters $\theta_s$ are in order the maximum amplitude, the time-scale over which the quasi-periodicity evolves, the period of the periodic variation (mapping the stellar rotation), and the ‘smoothness’ of the fit (its amount of high-frequency structure) also often referred to as the harmonic complexity. The ‘jitter’ term is represented by the delta function and $\beta$ can be thought of as the contribution to the RVs from the precision on the spectrograph.

<footnote>$^5$To be precise PYTHON defines its minimum and maximum float values to specific numbers, so these are actually uniform priors between $[-1.7976931348623157 \times 10^{308}, 1.7976931348623157 \times 10^{308}]$.</footnote>

<footnote>$^6$Available at https://github.com/frescigno/magpy_rv</footnote>
While eccentricity *e*, and planetary angle of periastron $\omega_p$ were used within the Keplerian model, when iterating in the MCMC algorithm we instead took steps in a different set of variables $S_k$ and $C_k$, defined as

$$ 
\begin{aligned}
S_k &= \sqrt{e} \sin \omega_p, \\
C_k &= \sqrt{e} \cos \omega_p.
\end{aligned} 
\qquad (3) 
$$

As explained in Eastman, Gaudi & Agol (2013), this reparametrization avoids a boundary condition at zero eccentricity, allowing for a better sampling around zero while maintaining the overall prior flat over eccentricity.

The Keplerian models also depended on time of periastron passage $t_p$, rather than the time of transit $t_0$, derived by transit photometry. However, the two variables are linked via the following equation

$$ t_p = t_0 - \frac{P}{2\pi} \cdot [E_{tr} - e \cdot \sin(E_{tr})], \qquad (4) $$

in which *P* is the orbital period of the considered planet, *e* its eccentricity, and the eccentric anomaly $E_{tr}$ is computed from the argument of periastron and the eccentricity as

$$ E_{tr} = 2 \arctan \left[ \sqrt{\frac{1-e}{1+e}} \cdot \tan \left( \frac{\pi - 2\omega_p}{4} \right) \right]. \qquad (5) $$

We conducted our investigation on the combined HARPS-N and SOPHIE data set, as well as on the two data sets separately. Once again we were able to combine the two RV data sets with a simple offset parameter and could use a single GP to describe both because they have comparable jitters and they are derived by similar spectral windows in the optical range. Therefore, they are expected to map the same physical processes and to be sensitive to Doppler shift in the same way.

## 6.1 Selections of priors

In this section, we describe the choices of priors for the analysis of the RV data. The same priors are used for all three analyses. They are also summarized in Table 4.

Starting with the Keplerians, we imposed a strict $1\sigma$ Gaussian prior on the orbital period of the inner transiting planet, $P_b$, derived from the posterior distribution of the same variable in the transit photometry analysis. Similarly, we imposed a strict Gaussian prior to the time of periastron passage, $t_p$, inflating the $\sigma$ to account for the uncertainties in the eccentricity of the planet. The period of the outer planet was bound by a uniform prior between [75,150], derived from the minimum period allowed by consecutive TESS photometry and the information derived from the periodogram analysis. Given the inability to derive a period from transit photometry, the time of periastron passage of the outer planet $t_{p,c}$ was bound by a uniform prior in the range [2459678.5, 2459773.5], determined by the preliminary $P_c$ from the periodograms. $S_k$ and $C_k$ for both planets are also bound by uniform priors in the range [-1, 1]. The SOPHIE-HARPS-N offset was allowed to vary only in the [-5,5] m s⁻¹ interval. The rest of the parameters are left with wide positive (larger than zero) uniform priors.

Regarding the kernel hyperparameters, we applied a strict Gaussian prior to $\theta_4$ (the 'smoothness' of the fit) centred on $0.5 \pm 0.05$, as recommended by Jeffers & Keller (2009). This choice is grounded in the fact that even highly complex active-region distributions average out to just two or three large active regions per rotation. We set a wide Gaussian prior on the stellar rotation period $\theta_2$ derived from

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Phase-folded TESS light curves of TOI-2134b and c. Faint points are individual TESS two-minute cadence measurements, bold darker points are data binned in orbital phase, and the line curves are the best-fitting transit models. The error bars on the binned points are smaller than the symbols. For the transit of TOI-2134c, we artificially offset the out-of-transit flux measurements for improved visibility.</figcaption></figure>

<table><caption>**Table 3.** Results and uncertainties of the planetary parameters for the photometry analysis described in Section 5.</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Radius ratio (R<sub>b</sub>/R<sub>*</sub>)</td><td>0.03475 ± 0.00038</td></tr><tr><td>Orbital period P<sub>b</sub> [d]</td><td>9.2292005 ± 0.0000063</td></tr><tr><td>Time of transit t<sub>0, b</sub> [BJD]</td><td>2459407.54493 ± 0.00027</td></tr><tr><td>Orbital inclination i<sub>b</sub> [deg]</td><td>89.49 ± 0.37</td></tr><tr><td>Transit impact parameter b<sub>b</sub></td><td>0.21 ± 0.14</td></tr><tr><td>Radius ratio (R<sub>c</sub>/R<sub>*</sub>)</td><td>0.09404 ± 0.00078</td></tr><tr><td>Time of transit t<sub>0, c</sub> [BJD]</td><td>2459718.96939 ± 0.00020</td></tr><tr><td>Transit impact parameter b<sub>c</sub></td><td>0.464 ± 0.042</td></tr></tbody></table>
<table><caption>**Table 4.** Results from the three GP regression analysis. We include the priors applied to each parameter. In order the HARPS-N RVs only, and the SOPHIE RVs only results, followed by the combined HARPS-N and SOPHIE data results (used for all further analysis). We abbreviated uniform priors as *U*, Gaussian priors as *G*, and Jeffreys' priors as *J*. We only show the results obtained for the high-eccentricity case, as addressed in Section 6.2.</caption><thead><tr><th>Parameter</th><th>Prior</th><th>HARPS-N RVs</th><th>SOPHIE RVs</th><th>Combined RVs</th></tr></thead><tbody><tr><td>GP amplitude θ₁ [m s<sup>-1</sup>]</td><td>U[0, 20]</td><td>4.24<sup>+0.81</sup><sub>-0.59</sub></td><td>5.52<sup>+0.98</sup><sub>-0.68</sub></td><td>5.52<sup>+0.67</sup><sub>-0.71</sub></td></tr><tr><td>GP time-scale θ₂ [d]</td><td>J[0, 100]</td><td>31.84<sup>+9.93</sup><sub>-10.36</sub></td><td>10.15<sup>+22.51</sup><sub>-7.99</sub></td><td>25.05<sup>+8.48</sup><sub>-8.53</sub></td></tr><tr><td>GP period θ₃ [d]</td><td>G[48, 10]</td><td>45.85<sup>+4.89</sup><sub>-4.84</sub></td><td>38.89<sup>+13.99</sup><sub>14.17</sub></td><td>45.78<sup>+5.56</sup><sub>-5.31</sub></td></tr><tr><td>GP smoothness θ₄</td><td>G[0.5, 0.05]</td><td>0.48<sup>+0.05</sup><sub>-0.05</sub></td><td>0.48<sup>+0.06</sup><sub>-0.05</sub></td><td>0.48<sup>+0.05</sup><sub>-0.05</sub></td></tr><tr><td>Jitter [m s<sup>-1</sup>]</td><td>U[0, 2]</td><td>0.69<sup>+0.13</sup><sub>-0.12</sub></td><td>0.82<sup>+0.24</sup><sub>22</sub></td><td>0.91<sup>+0.14</sup><sub>-0.13</sub></td></tr><tr><td>SOPHIE—HARPS-N offset [m s<sup>-1</sup>]</td><td>U[-5, 5]</td><td></td><td></td><td>2.30<sup>+0.44</sup><sub>-0.45</sub></td></tr><tr><td>Orbital period P<sub>b</sub> [d]</td><td>G[9.2292004, 0.0000063]</td><td>9.22923<sup>+0.00004</sup><sub>-0.00003</sub></td><td>9.2292<sup>+0.0002</sup><sub>-0.0001</sub></td><td>9.22923<sup>+0.00004</sup><sub>-0.00004</sub></td></tr><tr><td>RV amplitude K<sub>b</sub> [m s<sup>-1</sup>]</td><td>U[0, 20]</td><td>3.01<sup>+0.32</sup><sub>-0.32</sub></td><td>4.13<sup>+0.84</sup><sub>-0.87</sub></td><td>3.40<sup>+0.28</sup><sub>-0.29</sub></td></tr><tr><td>S<sub>k, b</sub></td><td>U[-1, 1]</td><td>-0.04<sup>+0.06</sup><sub>-0.08</sub></td><td>0.21<sup>+0.08</sup><sub>-0.09</sub></td><td>-0.07<sup>+0.06</sup><sub>-0.09</sub></td></tr><tr><td>C<sub>k, b</sub></td><td>U[-1, 1]</td><td>0.22<sup>+0.06</sup><sub>-0.09</sub></td><td>0.21<sup>+0.07</sup><sub>-0.09</sub></td><td>0.21<sup>+0.06</sup><sub>-0.08</sub></td></tr><tr><td>Time of periastron t<sub>peri, b</sub> [BJD]</td><td>G[2459408.22, 0.50]</td><td>2459407.71<sup>+0.46</sup><sub>-0.38</sub></td><td>2459407.55<sup>+0.43</sup><sub>-0.44</sub></td><td>2459407.89<sup>+1.52</sup><sub>-1.23</sub></td></tr><tr><td>Orbital period P<sub>c</sub> [d]</td><td>U[75, 150]</td><td>94.71<sup>+1.17</sup><sub>-1.11</sub></td><td>94.86<sup>+1.13</sup><sub>-0.83</sub></td><td>95.50<sup>+0.36</sup><sub>-0.25</sub></td></tr><tr><td>RV amplitude K<sub>c</sub> [m s<sup>-1</sup>]</td><td>U[0, 20]</td><td>11.92<sup>+1.82</sup><sub>-1.82</sub></td><td>10.28<sup>+2.99</sup><sub>-2.94</sub></td><td>9.74<sup>+1.60</sup><sub>-1.63</sub></td></tr><tr><td>S<sub>k, c</sub></td><td>U[-1, 1]</td><td>-0.65<sup>+0.11</sup><sub>-0.07</sub></td><td>0.69<sup>+0.29</sup><sub>-0.10</sub></td><td>-0.57<sup>+0.19</sup><sub>-0.13</sub></td></tr><tr><td>C<sub>k, c</sub></td><td>U[-1, 1]</td><td>0.42<sup>+0.13</sup><sub>-0.21</sub></td><td>0.41<sup>+0.66</sup><sub>-0.32</sub></td><td>0.59<sup>+0.13</sup><sub>-0.16</sub></td></tr><tr><td>Time of periastron t<sub>peri, c</sub> [BJD]</td><td>U[2459678.5, 2459773.5]</td><td>2459724.33<sup>+3.27</sup><sub>-2.53</sub></td><td>2459731.05<sup>+4.07</sup><sub>-8.86</sub></td><td>2459721.20<sup>+1.52</sup><sub>-1.23</sub></td></tr></tbody></table>

the periodogram analysis centred in 48 d with $\sigma = 10$ d, as wide as the forest of peaks in the WASP BGLS periodogram. The evolution time-scale $\theta_3$ is bound by a wide Jeffreys' prior. A Jeffreys' prior is a uniform, uninformed prior that is invariant under reparametrization of the given parameter vector. It is less informative than a uniform prior when the scale and range of the considered parameter is not known, as it corresponds to a uniform probability density in logarithmic frequency. A wide positive (larger than zero) uniform prior is applied to the GP amplitude $\theta_1$, and the jitter is only allowed to vary in the interval [0,2] m s$^{-1}$.

## 6.2 The eccentricity of TOI-2134c

Initial analysis of the RV data showed a significant trimodality in the distribution of the eccentricity of the outer 95 d-orbit planet, $e_c$. After further investigation we found that multiple fully converged models with different outer planet eccentricities existed. The RVs allowed for eccentricities of TOI-2134c equal to $0.002^{+0.029}_{-0.002}$, $0.45 \pm 0.05$, and $0.67^{+0.05}_{-0.06}$. All the models agreed within their uncertainties for most other parameters. Significantly large eccentricities have been detected before for temperate gas planets (as mentioned in Introduction) and stability can be reached within this system, so we could not a priori exclude any of the models. The stellar rotational period derived from the analysis is close to half the period of TOI-2134c. We therefore postulated that an interaction between the fit of the Keplerian model and the stellar activity-induced signal by the GP could be the reason behind the multiple models. While the flexibility of GPs are what makes them valuable tools to model stellar activity, we believe that in this case this flexibility allowed the Keplerian to take different accepted forms, while absorbing any 'left-over' signal into the activity model. To further compare the final likelihoods of the three solutions, we computed the corrected Akaike Information Criterion, AICc, (Sugiura 1978) for all converged models:

$$ \text{AICc} = \text{AIC} + 2 \left( \frac{N_{\text{free}}(1 + N_{\text{free}})}{N_{\text{data}} - N_{\text{free}} + 1} \right), \quad (6) $$

where $N_{\text{free}}$ is the number of free parameters and $N_{\text{data}}$ is the number of data points. The original AIC (Akaike 1983) is calculated as

$$ \text{AIC} = -2 \ln \mathcal{L} + 2N_{\text{free}}, \quad (7) $$

where $\ln \mathcal{L}$ is the logarithmic-likelihood maximized after the MCMC analysis. The larger the AICc the less likely the model. The AICcs of the combined (HARPS-N + SOPHIE) RV data for the low-, medium-, and high-eccentricity models were, respectively, 1224.0 and 1195.7 and 1196.7. As a further check, and to test whether this system would significantly benefit from a simpler analysis, we also computed the Keplerian-only best-fitting model to the data. For this analysis, we only included the planetary model with a jitter term and no stellar activity or GP component. This last model struggled to converge and its AICc was 1253.2. This analysis led us to strongly disfavour the Keplerian-only and the circular-orbit models (with AICc difference from the best model larger than 7). However, the AICc values for the medium- and the high-eccentricity cases were similar enough that no single model was significantly favoured and no significant statistical preference could be reached.

We then turned to the obtained photometric data. We estimated the orbital period of the singly transiting planet candidate using only the *TESS* light curve, following the procedure of Vanderburg et al. (2018). This method does not take into consideration the results from RV, and derives the planetary period directly from the photometric mono-transit. We extracted the impact parameter $b_c$, planet-star radius ratio $R_c/R_*$, and total transit duration of the single transit candidate from the MCMC posteriors from our two-planet transit fit, and solved for the orbital period assuming the stellar parameters reported in this paper and an eccentricity probability distribution from Kipping (2014). We also imposed the constraint that a second transit was not observed by *TESS*, which requires the orbital period be longer than about 75 d. We found that the short duration of the transit and minimum period allowed by *TESS* rule out circular orbits for this planet with periastron passage happening near the time of transit (as expected from geometric arguments), as the RV model comparison also had found. We then estimated the eccentricity $e_c$
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Posterior distribution corner of the eccentricity and the argument of periastron $\omega_p$ of the outer planet c derived after MCMC model optimization on the deep mono-transit present in the *TESS* data, as explained in Section 6.2. Most notable, the eccentricity of TOI-2134c converges to a high ~0.7 value.</figcaption></figure>

and argument of periastron $\omega_{p,c}$ required to reproduce the transit data, assuming the orbital period larger than 75 d. The posterior probability distributions of $e_c$ and $\omega_{p,c}$ are shown in Fig. 9. The eccentricity is required to be high (~0.7), and the argument of periastron is broadly to happen near the conjunction of the orbit of the planet.

This eccentricity value derived from transit photometry was then used to constrain the RVs. Given the high-eccentricity preference, we added a Gaussian prior centred in 0.7 with a $\sigma$ of 0.1 to $e_c$. In this paper, we chose to only report the high-eccentricity RV models for the HARPS-N, SOPHIE, and the combined RV data set consistent with the results from photometry.

## 6.3 RV results

A summary of the final results of our RV analyses can be found in Table 4. For this MCMC analysis, we simultaneously evolved 100 chains for 100 000 iterations each, discarding a burn-in phase of 20 000 steps. We assessed the health and convergence of the chains by computing the Gelman–Rubin statistic and all parameters reached values under the 1.1 convergence cut. As mentioned in the previous section we tested a series of models. For each set of HARPS-N only, SOPHIE only and combined RVs we evolved Keplerian-only models with no stellar activity (which overall struggled to converge or did not converge), forced circular-orbit models, medium-eccentricity models, and finally high-eccentricity models bound with an eccentricity prior derived by the photometry analysis. In this paper, we only present the last set.

The HARPS-N only data can constrain the amplitude and period of the inner TOI-2134b better than the SOPHIE data can, but conversely the SOPHIE RVs are able to better identify the signal of the outer planet, especially its period. A combined analysis allows us to more robustly constrain both planets with a single model. Since all three of the GP regression models fully converged and reached final values consistently within $1\sigma$ of

each other, we only discuss the results of the combined RV analysis.

The periods of the two planets are well defined. Their RV amplitudes are constrained to $12\sigma$ for planet b and to $6\sigma$ for planet c. The MCMC struggles to constrain the stellar activity evolution time-scale $\theta_2$, as expected from the low correlation with activity indicators and the weak overall rotational modulation (see Section 4.2). The stellar rotation period is derived to be $45.78^{+5.56}_{-5.31}$ d.

# 7 JOINT PHOTOMETRY AND RV ANALYSIS

Finally, we also modelled the *TESS* photometry and the RV data jointly, to more robustly test whether the high-eccentricity model was still favoured. This more complex analysis allowed for simultaneous modelling of the orbital solutions for both planets. We once again used the code MAGPy_RV,<fnref n="7" /> which for joint photometry analysis includes transit modelling with the python package BATMAN (Kreidberg 2015).

We modelled the RVs similarly to Section 6, as two Keplerian signals for the planet candidates with a quasi-periodic kernel describing the stellar activity and an offset parameter to match the zero-line of the HARPS-N and the SOPHIE data sets. For the *TESS* data, we described the transits of both planets with six parameters each (period, time of transit, $S_k$, $C_k$, planet to stellar radius ratio, and orbital inclination). Our photometric model also included five parameters to describe the host star (its mean density, $q_1$, $q_2$, photometric jitter and offset). In this analysis, we are jointly modelling the periods, time of transits, eccentricity, and angle of periastron of each planet.

## 7.1 Selection of priors

We imposed the similar priors on the GP hyperparameters as described in Section 6.1: Gaussian priors on the stellar rotational period and the harmonic complexity, uniform priors on amplitude and RV jitter, and a Jeffreys’ prior on the evolution time-scale. The RV offset between SOPHIE and HARPS-N data was also similarly bound by a uniform prior between [-5,5]. The period of the inner planet, $P_b$, was bound by a Gaussian prior centred on 9.2 d with $\sigma$ of 0.2 d derived from preliminary transit analysis. The time of transit $t_{tr,b}$ was also similarly bound by a Gaussian prior. The period of the outer planet, $P_c$, was bound by a uniform prior between [75, 150] d, as it was in the original RV analysis. The RV amplitude of both planets were as before bound between [0,20] m s⁻¹. $S_k$ and $C_k$ of both planets were only allowed to vary in the interval [-1,1] by definition. For the photometry, the stellar density was bound by a Gaussian prior centred on the derived density in Section 3 with $\sigma$ equal to its uncertainty. We allowed both planet-to-star radius ratios, $R_b/R_*$ and $R_c/R_*$, to only vary between [0,1] (we expect the planets to be smaller then the star), $q_1$ and $q_2$ between [0,1] as per their definition, and we required both inclinations $i$ to be less than 90°. All other priors were flat uninformative priors.

## 7.2 Joint analysis results

We simultaneously evolved 100 chains for 100 000 iterations each, discarding once again a burn-in phase of 20 000 steps and we tested for convergence with the Gelman–Rubin statistic. The results of our combined analysis are listed in Table 5. All parameters agree

<footnote n="7">⁷This version of MAGPy_RV is not yet public.</footnote>
<table><caption>**Table 5.** Results and uncertainties of the planetary parameters for the joint photometry and RV analysis described in Section 7.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>GP amplitude θ<sub>1</sub> [m s<sup>-1</sup>]</td><td>4.59<sup>+1.38</sup><sub>-1.29</sub></td></tr><tr><td>GP time-scale θ<sub>2</sub> [d]</td><td>28.01<sup>+21.31</sup><sub>-22.15</sub></td></tr><tr><td>GP period θ<sub>3</sub> [d]</td><td>53.87<sup>+3.14</sup><sub>-3.02</sub></td></tr><tr><td>GP smoothness θ<sub>4</sub></td><td>0.44<sup>+0.08</sup><sub>-0.06</sub></td></tr><tr><td>Jitter [m s<sup>-1</sup>]</td><td>0.85<sup>+0.95</sup><sub>-0.59</sub></td></tr><tr><td>SOPHIE HARPS-N Offset [m s<sup>-1</sup>]</td><td>2.64<sup>+0.09</sup><sub>-0.06</sub></td></tr><tr><td>Orbital period P<sub>b</sub> [d]</td><td>9.229209<sup>+0.000006</sup><sub>-0.000004</sub></td></tr><tr><td>Radius ratio (R<sub>b</sub>/R<sub>*</sub>)</td><td>0.02 ± 0.01</td></tr><tr><td>Orbital inclination i<sub>b</sub> [deg]</td><td>89.91<sup>+0.05</sup><sub>-0.06</sub></td></tr><tr><td>RV Amplitude K<sub>b</sub> [m s<sup>-1</sup>]</td><td>3.51<sup>+0.33</sup><sub>-0.41</sub></td></tr><tr><td>Eccentricity e<sub>b</sub></td><td>0.05<sup>+0.03</sup><sub>-0.03</sub></td></tr><tr><td>Argument of periastron ω<sub>p, b</sub> [rad]</td><td>-0.75<sup>+0.47</sup><sub>-0.88</sub></td></tr><tr><td>Time of periastron t<sub>p, b</sub> [BJD]</td><td>2459407.82<sup>+0.09</sup><sub>-0.06</sub></td></tr><tr><td>Orbital period P<sub>c</sub> [d]</td><td>94.98<sup>+0.95</sup><sub>-1.02</sub></td></tr><tr><td>Radius ratio (R<sub>c</sub>/R<sub>*</sub>)</td><td>0.09 ± 0.01</td></tr><tr><td>Orbital inclination i<sub>b</sub> [deg]</td><td>89.91<sup>+0.02</sup><sub>-0.03</sub></td></tr><tr><td>RV Amplitude K<sub>c</sub> [m s<sup>-1</sup>]</td><td>9.83<sup>+0.85</sup><sub>-0.89</sub></td></tr><tr><td>Eccentricity e<sub>c</sub></td><td>0.62<sup>+0.09</sup><sub>-0.02</sub></td></tr><tr><td>Argument of periastron ω<sub>p, c</sub> [rad]</td><td>1.41<sup>+0.49</sup><sub>-0.48</sub></td></tr><tr><td>Time of periastron t<sub>p, c</sub> [BJD]</td><td>2459432.39<sup>+3.11</sup><sub>-3.01</sub></td></tr></tbody></table>

within $1\sigma$ uncertainty with the results from the previous less complex transit and RV analyses, shown in Tables 3 and 4. These results once again confirmed the high-eccentricity model for the outer planet TOI-2134c.

Overall, we were able fully recover both planet candidates and their periods. Their RV amplitudes were constrained to $10\sigma$ for the inner planet and $11\sigma$ for the outer one. The joint photometry and RV analysis is minorly less effective in the retrieval of the RV signal of inner planet than the RV data on their own, but it performed better for TOI-2134c. Once again, the stellar activity evolution time-scale is not very well-constrained. The stellar rotational period was here derived to be slightly longer ($54.27^{+3.27}_{-3.23}$), but it is still consistent with the previous analysis. Both planet radius ratios were fully retrieved to $2\sigma$ and $9\sigma$ for TOI-2134 b and c, respectively.

## 8 RESULTS AND DISCUSSION

The results of the joint photometry and RV analysis fully agree within their $1\sigma$ uncertainties with the results from the separate transit and RV analyses. While the joint method successfully retrieved and characterized both planet candidates, from here on, we chose to use the results from the less complex, separated analyses undertaken in Sections 5 and 6. All the final results are compiled in Table 6. In Fig. 10, we plot the combined SOPHIE and HARPS-N data set alongside the complete best-fitting model in grey, as well as the GP-predicted activity as a black dashed line. Fig. 11 shows the phase-folded, best-fitting Keplerian orbital models, after subtracting the stellar activity-induced signal modelled by the GP, and their residuals.

As a result of our investigation, we establish the presence of an inner planet TOI-2134b, and an outer planet TOI-2134c. All derived planetary characteristics are listed in Table 6. Fig. 12 shows the two planets in a mass-radius diagram.

<table><caption>**Table 6.** System parameters for the TOI-2134 system. The transit and RV parameters are derived in Sections 5 and 6. Derived parameters are addressed in Section 8 and its subsections alongside the necessary assumptions.</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td colspan="2">GP regression – modelled activity parameters</td></tr><tr><td>GP amplitude θ<sub>1</sub> [m s<sup>-1</sup>]</td><td>5.52<sup>+0.67</sup><sub>-0.71</sub></td></tr><tr><td>GP time-scale θ<sub>2</sub> [d]</td><td>25.05<sup>+8.48</sup><sub>-8.53</sub></td></tr><tr><td>GP period θ<sub>3</sub> [d]</td><td>45.78<sup>+5.56</sup><sub>-5.31</sub></td></tr><tr><td>GP smoothness θ<sub>4</sub></td><td>0.48<sup>+0.05</sup><sub>-0.05</sub></td></tr><tr><td>Jitter [m s<sup>-1</sup>]</td><td>0.91<sup>+0.14</sup><sub>-0.13</sub></td></tr><tr><td>SOPHIE HARPS-N Offset [m s<sup>-1</sup>]</td><td>2.30<sup>+0.44</sup><sub>-0.45</sub></td></tr><tr><td colspan="2">TOI-2134 b</td></tr><tr><td colspan="2">Transit and RV parameters</td></tr><tr><td>Orbital period P<sub>b</sub> [d]</td><td>9.2292005 ± 0.0000063</td></tr><tr><td>Time of transit t<sub>0, b</sub> [BJD]</td><td>2459407.54493 ± 0.00027</td></tr><tr><td>Radius ratio (R<sub>b</sub>/R<sub>*</sub>)</td><td>0.03475 ± 0.00038</td></tr><tr><td>Orbital inclination i<sub>b</sub> [deg]</td><td>89.49 ± 0.37</td></tr><tr><td>Transit impact parameter b<sub>b</sub></td><td>0.21 ± 0.14</td></tr><tr><td>Transit duration τ<sub>b</sub> [h]</td><td>2.995 ± 0.047</td></tr><tr><td>RV Amplitude K<sub>b</sub> [m s<sup>-1</sup>]</td><td>3.40<sup>+0.28</sup><sub>-0.29</sub></td></tr><tr><td>Eccentricity e<sub>b</sub></td><td>0.06<sup>+0.03</sup><sub>-0.04</sub></td></tr><tr><td>Argument of periastron ω<sub>p, b</sub> [rad]</td><td>1.91<sup>+0.32</sup><sub>-0.34</sub></td></tr><tr><td>Time of periastron t<sub>p, b</sub> [BJD]</td><td>2459407.89<sup>+0.45</sup><sub>-0.49</sub></td></tr><tr><td colspan="2">Derived parameters</td></tr><tr><td>Radius R<sub>c</sub> [R⊕]</td><td>2.69 ± 0.16</td></tr><tr><td>Mass M<sub>c</sub> [M⊕]</td><td>9.13<sup>+0.78</sup><sub>-0.76</sub></td></tr><tr><td>Density ρ<sub>b</sub> [kg m<sup>-3</sup>]</td><td>2607 ± 516</td></tr><tr><td>Density ρ<sub>c</sub> [kg m<sup>-3</sup>]</td><td>0.47 ± 0.09</td></tr><tr><td>Scaled semi-major axis (a<sub>b</sub>/R<sub>*</sub>)</td><td>23.66 ± 0.52</td></tr><tr><td>Semi-major axis a<sub>b</sub> [au]</td><td>0.0780 ± 0.0009</td></tr><tr><td>Incident flux F<sub>inc, b</sub> [F<sub>inc, ⊕</sub>]</td><td>32 ± 2</td></tr><tr><td>Equilibrium temperature T<sub>eq, b</sub> [K]</td><td>666 ± 8</td></tr><tr><td colspan="2">TOI-2134 c</td></tr><tr><td colspan="2">Transit and RV parameters</td></tr><tr><td>Orbital period P<sub>c</sub> [d]</td><td>95.50<sup>+0.36</sup><sub>-0.25</sub></td></tr><tr><td>Time of transit t<sub>0, c</sub> [BJD]</td><td>2459718.96939 ± 0.00020</td></tr><tr><td>Radius ratio (R<sub>c</sub>/R<sub>*</sub>)</td><td>0.09404 ± 0.00078</td></tr><tr><td>Transit impact parameter b<sub>c</sub></td><td>0.464 ± 0.042</td></tr><tr><td>Transit duration τ<sub>c, c</sub> [h]</td><td>5.267 ± 0.028</td></tr><tr><td>RV Amplitude K<sub>c</sub> [m s<sup>-1</sup>]</td><td>9.74<sup>+1.60</sup><sub>-1.63</sub></td></tr><tr><td>Eccentricity e<sub>c</sub></td><td>0.67<sup>+0.05</sup><sub>-0.06</sub></td></tr><tr><td>Argument of periastron ω<sub>p, c</sub> [rad]</td><td>0.232<sup>+0.22</sup><sub>-0.32</sub></td></tr><tr><td>Time of periastron t<sub>p, c</sub> [BJD]</td><td>2459721.20<sup>+1.52</sup><sub>-1.23</sub></td></tr><tr><td colspan="2">Derived parameters</td></tr><tr><td>Radius R<sub>c</sub> [R⊕]</td><td>7.27 ± 0.42</td></tr><tr><td>Mass M<sub>c</sub> [M⊕]</td><td>41.89<sup>+7.69</sup><sub>-7.83</sub></td></tr><tr><td>Density ρ<sub>c</sub> [kg m<sup>-3</sup>]</td><td>599 ± 152</td></tr><tr><td>Density ρ<sub>c</sub> [kg m<sup>-3</sup>]</td><td>0.11 ± 0.03</td></tr><tr><td>Scaled semi-major axis (a<sub>c</sub>/R<sub>*</sub>)</td><td>112 ± 2</td></tr><tr><td>Semi-major axis a<sub>c</sub> [au]</td><td>0.371 ± 0.004</td></tr><tr><td>Incident Flux F<sub>inc, c</sub> [F<sub>inc, ⊕</sub>]</td><td>1.4 ± 0.1</td></tr><tr><td>Equilibrium temperature T<sub>eq, c</sub> [K]</td><td>306 ± 4</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Combined SOPHIE (orange) and HARPS-N (purple) RV data plotted with error bars (HARPS-N error bars are too small to be clearly visible). The complete model, which includes two Keplerians and the predicted activity, is plotted in as a solid, with its uncertainties as the shaded area. The dashed darker line represents the GP activity prediction only. On the bottom, the residuals between the data (in the corresponding colour) and the complete model are plotted.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Phase-folded activity model-subtracted plots for the inner (top) and outer (bottom) planets. In orange are the SOPHIE RVs and in purple the HARPS-N ones with respective error bars (some HARPS-N error bars may be too small to be visible). The Keplerian model is plotted as a solid line, with the residuals shown on the bottom. The phase has also been extended on both sides.</figcaption></figure>

We computed for the inner planet TOI-2134b a mass $M_b$ of 9.13$^{+0.78}_{-0.76}$ M$_⊕$ and a radius of 2.69 ± 0.16 R$_⊕$, for an orbital period of 9.2292005 ± 0.0000063 d. Combining mass and radius yielded a bulk density of 0.47 ± 0.09 $\rho_⊕$. In the mass–radius diagram TOI-2134b falls in a parameter space significantly degenerate in composition. Planet b could be a 100 per cent water-planet (Zeng, Sasselov & Jacobsen 2016). At the same time, it could also have a rocky core, a water (or other heavy volatile elements) layer and a low-mass H/He

envelope. Overall, it is not possible to distinguish the composition of planet b without additional information. For more information about the atmospheric characteristics of TOI-2134b, see Zhang et al. (2023).

The outer planet TOI-2134c has mass $M_c$ of 41.89$^{+7.69}_{-7.83}$ M$_⊕$ and a radius of 7.27 ± 0.42 R$_⊕$, for a period of 95.50$^{+0.36}_{-0.25}$ d. The derived mass from the RVs and radius from photometry are well-matched and further justify the association of the mono-transit and the detected RV signal. The presence of a third planet with similar mass to TOI-2134c that could instead explain the transit would have been detected in the RV. The bulk density of TOI-2134c is calculated to be 0.11 ± 0.03 $\rho_⊕$ (similar to the density of Saturn). It can therefore be considered a long-orbit mini-Saturn. Given its derived period, we also went back to the other photometric data and computed when transits would have occurred. The derived transit times are plotted in Fig. 2 as black dashed lines, and their uncertainty windows as grey shaded regions. TOI-2134c transited 5 times over the 3 yr of WASP coverage, but none of those transits was originally detected. The possible explanation for this is twofold. On one hand, WASP is a ground instrument and therefore only observes during dark hours; given the transit duration of ~5 h, the event could have easily been missed. At the same time, the precision of the WASP data fluctuates significantly and a 0.01 flux deficit (as it is for TOI-2134c) is often too shallow for WASP to reliably detect.

## 8.1 System orbital stability

As a preliminary test of the stability of the system given the high eccentricity of TOI-2134c, we calculated the radius of the Hill Sphere (Hamilton & Burns 1992) of the outer planet and compared it to the closest approach distance between the two planets. If the orbit of the inner TOI-2134b at any point falls within the Hill Sphere of TOI-2134c, we expect the two bodies to gravitationally interact enough to de-stabilize their orbits. If a body of mass m is orbiting a larger body of mass M at semi-major axis a with an eccentricity e, the Hill Radius $R_{Hill}$ of the smaller body can be approximated to be

$$R_{Hill} \approx a(1 - e)^3 \sqrt{\frac{m}{3M}} \qquad (8)$$
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Mass–radius diagram with zoom-in for sub-Neptunian planets. The data are taken from the EU Exoplanet catalogue: http://exoplanet.eu/catalog/ on 2023 Feb 17. The solid line shows the mass–radius relation developed by Chen & Kipping (2017), with its categorization of Terran ($M < 2M_⊕$), Neptunian ($2M_⊕ < M < 0.4M_⊕$) and Jovian worlds ($M > 0.4M_⊕$). The zoomed-in plot includes composition lines taken from Zeng, Sasselov & Jacobsen (2016), and the Radius Valley band. Solar system planets are included for scale.</figcaption></figure>

For planet c, we computed an $R_{Hill,c}$ of $0.0047 \pm 0.0008$ au. The closest approach between the outer and inner planets is $0.048 \pm 0.026$ au. Therefore, the orbit of planet b at no point intersects with the Hill Sphere of TOI-2134c. To further assess the stability of the system under the high-eccentricity $e_c$ model, we also computed the chaos map in the neighbourhood of the best-fitting solution to the high-eccentricity model, shown in Fig. 13. We created a grid of $81 \times 81$ system configurations that vary between each other based on period $P_c$ and eccentricity $e_c$. All other parameters were fixed to their values derived from the MCMC best-fitting estimation. Each system defined a unique set of initial conditions that was then used for 50 kyr numerical integrations with REBOUND<fnref n="8" /> (Rein & Liu 2012) with the 15th order adaptive time-step integrator IAS15 (Rein & Spiegel 2015). We also included in our analysis the correction from general relativity implemented in the REBOUND extension REBOUNDX<fnref n="9" /> by Tamayo et al. (2020). After the simulations, we computed the Numerical Analysis of Fundamental Frequencies (NAFF, Laskar 1990, 1993). The NAFF indicator informs about the amount of chaos in a planetary orbit by precisely estimating its main frequency via a technique called frequency analysis (Laskar 1988). The main frequency of a planetary orbit corresponds to the mean-motion, which does not drift over time in non-chaotic dynamics, but does drift if the system is chaotic. Therefore, we apply frequency analysis on the two halves of each simulation, and for each planetary orbit, to estimate the amount of drift in the mean-motions. Weakly chaotic (hence stable) orbits should only show small differences in mean motions between

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Chaos map for the outer planet TOI-2134c. The period $P_c$ and eccentricity $e_c$ are explored on a $81 \times 81$ grid of different system configurations. After numerical integrations the NAFF indicator is computed and plotted as a colour scale. Blue regions (with more negative NAFF values) correspond to weakly chaotic, therefore more stable, planetary systems, while red areas (larger NAFF) refer to strongly chaotic systems, and hence more unstable. The best-fitting system position in this space together with its $1\sigma$ uncertainties indicate that both stable and unstable solutions are compatible with our high-eccentricity fit.</figcaption></figure>

the two integration halves. In this work we consider as the NAFF of the system the logarithmic maximum value of this drift, defined as

$$ \text{NAFF} = \max_i \left[ \log_{10} \frac{\Delta n_i}{n_0} \right], \qquad (9) $$

<footnote n="8">⁸REBOUND is an open-source software package dedicated to N-body integrations: http://rebound.readthedocs.org</footnote>

<footnote n="9">⁹Available at https://reboundx.readthedocs.io</footnote>
in which *i* refers to the chosen planet, $\Delta n_i$ is the difference in the mean-motion of planet *i* between its estimation on the first and second halves of the integrations, and $n_0$ is the initial mean motion of that planet *i*. In Fig. 13, blue regions have lower NAFF, and are weakly chaotic. Red regions correspond to systems that undergo strong chaos, and likely lead to rapid instability.<fnref n="10" /> White regions refer to those systems which had an escape or a close encounter between two bodies, and for which the simulation was stopped. We also overplot the area of $1\sigma$ limit uncertainties on the estimates of $P_c$ and $e_c$. Inside the subsequent square, we find that both chaotic and regular systems can exist. In other words, the high-eccentricity model is not incompatible with the system stability.

## 8.2 Planetary incident flux and equilibrium temperature

The incident flux of a planet $F_{inc}$ is computed from stellar luminosity $L_\star$ and planetary semi-major axis *a* with the following formula:

$$F_{inc} = \frac{L_{\star}}{4\pi a^2} = \frac{4\pi R_{\star}^2 \sigma_{SB} T_{eff}^4}{4\pi a^2}, \quad (10)$$

where $T_{eff}$ and $R_\star$ are the stellar effective temperature and radius and $\sigma_{SB}$ is the Stefan–Boltzmann constant. We can express this same incident flux in Earth units as:

$$\frac{F_{inc}}{F_{inc, \oplus}} = \left(\frac{T_{eff}}{T_{\odot}}\right)^4 \left(\frac{R_{\star}}{R_{\odot}}\right)^2 \left(\frac{1}{a}\right)^2, \quad (11)$$

in which $T_\odot$ and $R_\odot$ are the solar effective temperature and radius and *a* is expressed in au. Given semi-major axes $a_b$ and $a_c$ of 0.0780 ± 0.0009 and 0.371 ± 0.004 au, respectively, we computed incident fluxes of 33 ± 2 and 1.4 ± 0.1 $F_{inc, \oplus}$ for planet *b* and *c*.

The planets' equilibrium temperatures $T_{eq}$ can be derived as

$$T_{eq} = T_{eff} \sqrt{\frac{R_{\star}}{2a}} [f(1 - A_B)]^{1/4}, \quad (12)$$

where $A_B$ is the Bond albedo of the considered planet and *f* represents the effectiveness of atmospheric circulation. Assuming isotropic re-emission and a uniform equilibrium temperature over the entire planet (therefore *f* = 1), an upper limit on $T_{eq}$ can be derived from equation (12) by setting $A_B = 0$. We, therefore, calculated the upper limit of the equilibrium temperature of planet *b* to be 666 ± 8 K, and of planet *c* to be 305 ± 4 K.

From this analysis, the upper limit of the equilibrium temperature of the sub-Saturn object would be compatible with liquid water. Planet *c* is a gas giant, but could host potentially temperate rocky moons. However, the orbit of TOI-2134c is highly eccentric and the distance of the planet from the star changes significantly during its orbit, as shown in purple in Fig. 14. The boundaries of the habitable zone (HZ) of the system, $r_{HZ, \star}$, can be derived from the solar luminosity $L_\odot$ and the stellar luminosity as:

$$\frac{L_{\odot}}{r_{HZ, \odot}} = \frac{L_{\star}}{r_{HZ, \star}}, \quad (13)$$

where $r_{HZ, \odot}$ is the radius of the boundaries of the solar HZ. The boundaries in this paper were determined following the two models for narrow and empirical HZs described in Kopparapu et al. (2014). The narrow HZ is bound by an inner Runaway Greenhouse limit and an outer Maximum Greenhouse limit. The boundaries of the empirical HZ are defined by the Recent Venus and Early Mars limits.

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Depiction of the configuration of the TOI-2134 system. We include the inner planet with a circular orbit of 9.2292004 ± 0.0000063 d in blue, and the outer planet with an eccentric ($e_c = 0.67^{+0.05}_{-0.06}$) orbit of 95.50$^{+0.36}_{-0.25}$ d in purple. Their uncertainties are depicted as lighter orbits. The HZ boundaries are indicated as green shaded regions: the empirical HZ is plotted in lighter green, while the narrow HZ is overplotted in darker green. The boundaries are computed as described in Section 8.2 based on results from Kopparapu et al. (2014).</figcaption></figure>

The narrow and empirical HZs for the TOI-2134 system are shown in Fig. 14 respectively in dark and light green. As Fig. 14 clearly shows, TOI-2134c only spends less than half of its orbit within the HZ boundaries. In fact, we also computed the incident flux and upper limit of the equilibrium temperature planet *c* at periastron to be 13 ± 4 $F_{inc, \oplus}$ and 533 ± 8 K.

## 8.3 Suggested follow-up observations

### 8.3.1 Long-term RV observations and transit detection for TOI-2134c

This system would foremost benefit from long-term RV observations to better constrain the period and eccentricity of the outer planet. Both HARPS-N and SOPHIE plan on continuing observing the star sporadically. A second photometric observing campaign aimed at detecting another transit of the outer planet candidate would also be valuable. In the current mission plan, TESS will re-observe TOI-2134 in Sectors 74, 79, and 80 in 2024. A transit of planet *c* should occur in Sector 80 (2024 late-June to early-July). Given the brightness of TOI-2134 and the larger radius ratio between planet *c* and its host star, transits of the outer planet can also be observed with ground-based telescopes. Another firm detection of a transit would re-confirm its period and further inform the eccentricity model choice. We include a list of the times of transit between the original detection and the end of 2025 in Table 7. The uncertainties on the times of transit $\sigma_{tr}$ increase with increasing number of 'missed' transits as:

$$\sigma_{tr} = \sqrt{(n\sigma_p)^2 + \sigma_{t0}^2} \approx n\sigma_p, \quad (14)$$

<footnote n="10">¹⁰We refer to Stalport et al. (2022) for details on the link between NAFF and orbital stability.</footnote>
<table><caption>**Table 7.** List of times of transits of TOI-2134c between the detected mono-transit and the end of 2025. The uncertainty on the dates computed as shown in Section 8.3.1. The transit that should be observed by TESS in Sector 80 is highlighted in bold.</caption><thead><tr><th>BJD</th><th>UT date (yyyy-mm-dd)</th><th>UT time (hh:mm:ss)</th></tr></thead><tbody><tr><td>2459814.5 ± 0.3</td><td>2022-08-22</td><td>23:20:35</td></tr><tr><td>2459910.0 ± 0.6</td><td>2022-11-26</td><td>11:25:12</td></tr><tr><td>2460005.5 ± 0.9</td><td>2023-03-01</td><td>23:29:53</td></tr><tr><td>2460101.0 ± 1.2</td><td>2023-06-05</td><td>11:34:34</td></tr><tr><td>2460196.5 ± 1.5</td><td>2023-09-08</td><td>23:39:11</td></tr><tr><td>2460292.0 ± 1.8</td><td>2023-12-13</td><td>11:43:52</td></tr><tr><td>2460387.4 ± 2.1</td><td>2024-03-17</td><td>23:48:29</td></tr><tr><td><b>2460483.0 ± 2.4</b></td><td><b>2024-06-21</b></td><td><b>11:53:10</b></td></tr><tr><td>2460578.5 ± 2.7</td><td>2024-09-24</td><td>23:57:50</td></tr><tr><td>2460674.0 ± 3.0</td><td>2024-12-29</td><td>12:02:28</td></tr><tr><td>2460769.5 ± 3.3</td><td>2025-04-04</td><td>00:07:08</td></tr><tr><td>2460865.0 ± 3.6</td><td>2025-07-08</td><td>12:11:46</td></tr><tr><td>2460960.5 ± 3.9</td><td>2025-10-12</td><td>00:16:26</td></tr></tbody></table>

in which *n* is the epoch since the observed transit, and $\sigma_p$ and $\sigma_{i0}$ are the uncertainties on, respectively, the period of the planet and its observed transit time.

### 8.3.2 Rossiter–McLaughlin analysis

Given the presence of both the inner mini-Neptune and the outer temperate sub-Saturn, (once the period of the outer planet is better defined with follow up RV observations or a second transit detection), TOI-2134 and its planets are scientifically valuable targets for follow-up RM (McLaughlin 1924; Rossiter 1924) analysis to determine the spin-orbit alignment of the system. The RM amplitude $K_{RM}$ can be computed as

$$K_{RM} = 52.8\,\text{ms}^{-1} \frac{v\sin(i)}{5\,\text{kms}^{-1}} \left(\frac{R_{\text{pl}}}{R_{\text{J}}}\right)^2 \left(\frac{R_{\star}}{R_{\odot}}\right)^{-2}, \quad (15)$$

in which $R_{pl}$ and $R_*$ are the radius of the considered transiting planet and the radius of the star. Instead of using a maximum limit for $v\sin(i)$, we recomputed it starting from the derived stellar rotational period to be $0.78 \pm 0.09 \, \text{km s}^{-1}$. Since both TOI-2134b and c transit, we computed the minimum expected RM amplitude for both: $K_{RM, b} = 0.98 \pm 0.17 \, \text{m s}^{-1}$ and $K_{RM, c} = 7.2 \pm 1.2 \, \text{m s}^{-1}$. Although the longer transit duration can be an obstacle, RM observations of temperate gas giants as TOI-2134c are valuable to further our understanding of planet migration. A significant fraction of hot giants are shown to have orbits that are misaligned with the rotational axis of their star (Winn et al. 2010; Albrecht et al. 2012). The origin of such misalignment is still unclear, but a leading hypothesis is that high-eccentricity migration tilts the orbit of the planet away from its initial plane via dynamical interactions (e.g. Fabrycky & Tremaine 2007; Ford & Rasio 2008; Petrovich 2015). Unlike hot giants, it is significantly more challenging to form temperate gas planets via high-eccentricity migration (Dong, Katz & Socrates 2013), and it is even less likely in the case of this system due to the presence of an inner small planet. Therefore, if high-eccentricity migration is in fact the driving factor behind the misalignment, the majority of temperate giants should have orbits aligned to spin of their star. However, given their lower transit probabilities, there are only few RM observations of temperate giants. Whether the aim is the whole transit or just observing the ingress or egress in a shorter summer

night, the temperate sub-Saturn planet c has a large peak-to-peak amplitude ($7.2 \pm 1.3 \, \text{m s}^{-1}$) that makes it easily observable. With a more firmly constrained eccentricity model, TOI-2134c will be a great candidate for RM follow-up.

### 8.3.3 Transmission spectroscopy

We also discussed the suitability of TOI-2134b and c for follow-up atmospheric characterization via transmission spectroscopy. Kemp- ton et al. (2018) developed an analytic metric to estimate the expected SNR of transmission-spectroscopy observations based on the strength of the spectral features and the brightness of the star, the Transmission Spectroscopy Metric, or TSM. It can be computed as:

$$\text{TSM} = \epsilon \cdot \frac{R_{\text{pl}}^3 T_{\text{eq}}}{M_{\text{pl}} R_*^2} \cdot 10^{-m_J/5}, \quad (16)$$

in which $R_{pl}$ and $M_{pl}$ are the radius and mass of the considered planet in Earth radii and masses, $R_*$ is the stellar radius in solar radii, $T_{eq}$ is the equilibrium temperature of the planet computed at zero albedo and full day-night heat redistribution (as in Section 8.2), and $m_J$ is the apparent magnitude of the host star in the J band. The term $\epsilon$ is a normalization factor to give one-to-one scaling to the JWST/NIRISS 10-h simulated observations described in Louie et al. (2018). This scaling factor also absorbs the unit conversion factors so that the parameters can be in natural units. The term $\epsilon$ changes depending on the radius of the planet, and is equal to 1.26 for TOI-2134b, and 1.15 for TOI-2134c. We computed a $TSM_b = 172 \pm 42$ and a $TSM_c = 243 \pm 54$. The TSMs of both planets are therefore considered well above the suggested cut-offs for their size bin. It is however important to note that the TSM was developed for targeted JWST effort and therefore it is not optimized for stars with $m_J < 9$ mag, as brighter stars require the bright readout mode and have substantially lower duty cycles. Given its brightness, TOI-2134 currently is observable without saturation by the JWST with NIRCam in its bright mode, with similar observational strategies as the ones successfully proposed by Dr Hu for 55 Cancri (program ID: 1952) and by Dr Deming for HD 189733b (program ID: 1633). Moreover, higher efficiency read modes for JWST observations are being investigated (Batalha et al. 2018) and future dedicated missions such as Ariel, and the ground-based Extremely Large Telescopes (ELTs) are suitable for brighter targets such as TOI-2134 (Houllé et al. 2021; Danielski et al. 2022).

## 9 SUMMARY AND CONCLUSIONS

In this work, we presented the photometric light curves of five TESS sectors and of three years of WASP monitoring, alongside 219 high-precision RV measurements obtained with HARPS-N and SOPHIE of the star TOI-2134. We characterized the star with multiple independent techniques and we studied its periodograms to better understand its stellar activity signals. We then performed a transit photometry analysis on the photometric data and a GP regression analysis on the RV data to constrain the radii and masses of the planets in the system. To test the statistical strength of the derived model, we also completed a joint analysis of the photometric and the RV data. The resulting planetary parameters fully agree within $1\sigma$ uncertainties with the results of the previous investigations. We therefore selected to focus on the results of the less complex, separate analyses for our discussion. As a result, we reached the following conclusions:
(i) We characterize the new multitransiting inner planet TOI-2134b in a 9.2292004 ± 0.0000063 d orbit with $M_b = 9.13_{-0.76}^{+0.78}$ M$_⊕$ (12σ detection) and $R_b = 2.69_{-0.16}^{+0.11}$ R$_⊕$. Its bulk density ($\rho_b = 0.47 \pm 0.09 \rho_⊕$) identifies the planet as either a water-world or a mini-Neptune with a rocky core and a low-mass H/He envelope. We computed the upper limit of the equilibrium temperature of the planet to be 666 ± 8 K.

(ii) We also constrain a second mono-transiting planet TOI-2134c with $M_c = 41.89_{-7.83}^{+7.69}$ M$_⊕$ (5σ detection) and $R_c = 7.27 \pm 0.42$ R$_⊕$ in a 95.50$^{+0.36}_{-0.25}$ d orbit, with an upper limit of the equilibrium temperature of 306 ± 4 K. Its bulk density ($\rho_c = 0.11 \pm 0.03 \rho_⊕$) is similar to Saturn's.

(iii) After GP regression, we find three possible orbital architectures for the outer TOI-2134c that model the RV data, one with low eccentricity (0.0002$^{+0.0025}_{-0.0002}$), one with medium eccentricity (0.45 ± 0.05), and one with high eccentricity (0.67$^{+0.05}_{-0.06}$). While we were able to disfavour the circular orbit case, the AICc values of the latter two solutions are comparable, therefore statistically there is no preference. We noted that in all models the rotational period of the star is half the orbital period of the outer TOI-2134c. We postulated that fitting interactions between the Keplerian model for the planet, and the activity-induced signal that the GP is extrapolating are the reason behind the multiple fully converged solutions. The flexibility of the GP allows the Keplerian to take different accepted forms while the GP model absorbs the residual signal and attributes it to stellar activity. As described in Section 6.2, further analysis of the photometry data showed that, given the derived orbital period for planet c, its transit duration time was too short to allow circular orbits. In fact, the mono-transit in the TESS data strongly prefers the high-eccentricity case. To further strengthen our results, we also undertook joint modelling of the photometric and the RV data. This investigation yielded a single converged state with an $e_c = 0.61_{-0.03}^{+0.08}$. In this paper, we therefore chose to present the high-eccentricity model of the separate, less complex RV only analysis and to use it for all further analysis. We also tested the stability of the system given these results and reached the conclusion that the high-eccentricity model is not incompatible with a stable system.

(iv) Since the mass-radius parameter space planet TOI-2134c resides in is not well populated and in order to better constrain its period and eccentricity, we recommend further RV observations and a second photometric observing campaign to detect another transit. To further characterize the architecture of the system we also recommend RM follow-up observations. We compute the expected RM amplitude of the temperate sub-Saturn TOI-2134c as $7.2 \pm 1.2$ m s$^{-1}$, making it accessible to ground-based instruments.

(v) We also compute the TSM of both planets of the system for possible follow-up atmospheric characterization via transmission spectroscopy. Although the projected SNRs place the planets well above the recommended cuts, TOI-2134 is close to the bright limits of most instruments on JWST, and is currently only observable with NIRCam in its bright mode. Future missions such as Ariel or ground-based transition spectroscopy will be suited for brighter target such as TOI-2134.

## ACKNOWLEDGEMENTS

This analysis is based on observations made with the Italian TNG operated on the island of La Palma by the Fundación Galileo Galilei of INAF (Istituto Nazionale di Astrofisica) at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias.

FR is funded by the University of Exeter's College of Engineering, Maths and Physical Sciences, UK.

The HARPS-N project was funded by the Prodex Program of the Swiss Space Office (SSO), the Harvard University Origin of Life Initiative (HUOLI), the Scottish Universities Physics Alliance (SUPA), the University of Geneva, the Smithsonian Astrophysical Observatory (SAO), the Italian National Astrophysical Institute (INAF), University of St. Andrews, Queen's University Belfast, and University of Edinburgh.

This work has been supported by the National Aeronautics and Space Administration under grant no. NNX17AB59G, issued through the Exoplanets Research Program.

This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40.182901 and 51NF40.205606.

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement SCORE no. 851555).

ACC and TGW acknowledge support from STFC consolidated grant numbers ST/R000824/1 and ST/V000861/1, and UKSA grant number ST/R003203/1.

FPE and CLO would like to acknowledge the Swiss National Science Foundation (SNSF) for supporting research with HARPS-N through the SNSF grants nos 140649, 152721, 166227, and 184618. The HARPS-N Instrument Project was partially funded through the Swiss ESA-PRODEX Programme.

Funding for the TESS mission is provided by NASA's Science Mission Directorate. KAC acknowledges support from the TESS mission via subaward s3449 from MIT.

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.

This paper includes data collected by the TESS mission that are publicly available from the Mikulski Archive for Space Telescopes (MAST).

We thank the Observatoire de Haute-Provence (CNRS) staff for its support in collecting SOPHIE data. This work was supported by the Programme National de Planétologie (PNP) of CNRS/INSU and CNES.

The postdoctoral fellowship of KB is funded by F.R.S.-FNRS grant T.0109.20 and by the Francqui Foundation.

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.

We acknowledge the use of public TESS data from pipelines at the TESS Science Office and at the TESS Science Processing Operations Center.

DD acknowledges support from the TESS Guest Investigator Program grants 80NSSC21K0108 and 80NSSC22K0185, and from the NASA Exoplanet Research Program grant 18-2XRP18.2-0136.

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.

RDH is funded by the UK Science and Technology Facilities Council (STFC)'s Ernest Rutherford Fellowship (grant no. ST/V004735/1).

SD is funded by the UK Science and Technology Facilities Council (grant no. ST/V004735/1).

BSL is funded by a UK Science and Technology Facilities Council (STFC) studentship (ST/V506679/1).
The postdoctoral fellowship of KB is funded by F.R.S.-FNRS grant T.0109.20 and by The Francqui Foundation.

XD and TF acknowledge funding by the French National Research Agency in the framework of the Investissements d'Avenir program (ANR-15-IDEX-02), through the funding of the 'Origin of Life' project of the Grenoble-Alpes University.

s financial support from the grant CEX2021-001131-S funded by MCIN/AEI/ 10.13039/501100011033.

DATA AVAILABILITY

The observational data presented in this publication are openly available. The TESS data are available at: https://mast.stsci.edu/portal/Mashup/Clients/Mast/Portal.html. The mentioned NEOSSat data can be found at: https://open.canada.ca/data/en/dataset/9ae3e718-8b6d-40b7-8aa4-858f00e84b30. The LOCGT data can be found at: https://exofop.ipac.caltech.edu/tess/target.php?id = 75878355. The WASP data cannot be found online, but can be made available on request. The HARPS-N and SOPHIE RV data, alongside their mentioned activity proxies are included as Supporting Information.

REFERENCES

Akaike H., 1983, Inter. Stat. Inst., 44, 277

Albrecht S. et al., 2012, ApJ, 757, 18

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

Allard F., Homeier D., Freytag B., Schaffenberger W., Rajpurohit A. S., 2013, Mem. Soc. Astron. Ital. Suppl., 24, 128

Andreasen D. T. et al., 2017, A&A, 600, A69

Astudillo-Defru N. et al., 2017, A&A, 605, L11

Bailer-Jones C. A. L., Rybizki J., Fouesneau M., Demleitner M., Andrae R., 2021, AJ, 161, 147

Barkaoui K., Gillon M., Benkhaldoun Z., Emmanuel J., Elhalkouj T., Daassou A., Burdanov A., Delrez L., 2017, J. Phys. Conf. Ser., 869, 012073

Barros S. C. C., Demangeon O., Díaz R. F., Cabrera J., Santos N. C., Faria J. P., Pereira F., 2020, A&A, 634, A75

Batalha N. M. et al., 2013, ApJS, 204, 24

Batalha N. E., Lewis N. K., Line M. R., Valenti J., Stevenson K., 2018, ApJ, 856, L34

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

Boisse I. et al., 2010, A&A, 523, A88

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

Borucki W. J. et al., 2011, ApJ, 736, 19

Bouchy F. et al., 2009, A&A, 505, 853

Bouchy F., Díaz R. F., Hébrard G., Arnold L., Boisse I., Delfosse X., Perruchot S., Santerne A., 2013, A&A, 549, A49

Bressan A., Marigo P., Girardi L., Salasnich B., Dal Cero C., Rubele S., Nanni A., 2012, MNRAS, 427, 127

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

Bryan M. L., Knutson H. A., Lee E. J., Fulton B. J., Batygin K., Ngo H., Meshkat T., 2019, AJ, 157, 52

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

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

Burt J. A. et al., 2020, AJ, 160, 153

Chabrier G., Baraffe I., Allard F., Hauschildt P., 2000, ApJ, 542, 464

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

Chen Y., Girardi L., Bressan A., Marigo P., Barbieri M., Kong X., 2014, MNRAS, 444, 2525

Chen Y., Bressan A., Girardi L., Marigo P., Kong X., Lanza A., 2015, MNRAS, 452, 1068

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

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

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

Costes J. C. et al., 2021, MNRAS, 505, 830

Crass J. et al., 2021, preprint (arXiv:2107.14291)

Crossfield I. J. M. et al., 2015, ApJ, 804, 10

Danielski C. et al., 2022, Exp. Astron., 53, 473

Dong S., Katz B., Socrates A., 2013, ApJ, 762, L26

Dotter A., 2016, ApJS, 222, 8

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

Dressing C. D., Charbonneau D., 2013, ApJ, 767, 95

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

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

Fabrycky D., Tremaine S., 2007, ApJ, 669, 1298

Faria J. P., Haywood R. D., Brewer B. J., Figueira P., Oshagh M., Santerne A., Santos N. C., 2016, A&A, 588, A31

Fischer D. A. et al., 2016, PASP, 128, 066001

Ford E. B., Rasio F. A., 2008, ApJ, 686, 621

Foreman-Mackey D. et al., 2013, Astrophysics Source Code Library, record(ascl:1303.002)

Fortney J. J., Visscher C., Marley M. S., Hood C. E., Line M. R., Thorngren D. P., Freedman R. S., Lupu R., 2020, AJ, 160, 288

Fox C., Wiegert P., 2022, BAAS, 40, 102.25

Fressin F., Guillot T., Nesta L., 2009, A&A, 504, 605

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

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

Gaia Collaboration, 2016, A&A, 595, A1

Gaia Collaboration, 2023, A&A, 674, A1

Gaidos E. et al., 2014, MNRAS, 443, 2561

Giclas H. L., Burnham R. J., Thomas N. G., 1979, Lowell Obs. Bull., 8, 145

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

Hadden S., Lithwick Y., 2014, ApJ, 787, 80

Hamilton D. P., Burns J. A., 1992, Icarus, 96, 43

Hara N. C., Mari A. R., 2021, Astrophysics Source Code Library, record(ascl:2112.024)

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

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

Hébrard G. et al., 2008, A&A, 488, 763

Hildebrand A. R. et al., 2004, Earth Moon Planets, 95, 33

Houllé M. et al., 2021, A&A, 652, A67

Huang C., Wu Y., Triaud A. H. M. J., 2016, ApJ, 825, 98

Huang C. X. et al., 2020a, Res. Notes Am. Astron. Soc., 4, 204

Huang C. X. et al., 2020b, Res. Notes Am. Astron. Soc., 4, 206

Jeffers S. V., Keller C. U., 2009, in Stempels E. ed., AIP Conf. Proc. Vol. 1094, 15th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun. Am. Inst. Phys., New York, p. 664

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

Jenkins J. M. et al., 2010, ApJ, 713, L87

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

Jiang J. H., Ji X., Cowan N., Hu R., Zhu Z., 2019, AJ, 158, 96

Johnson J. A., Aller K. M., Howard A. W., Crepp J. R., 2010, PASP, 122, 905

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

Jurgenson C., Fischer D., McCracken T., Sawyer D., Szymkowiak A., Davis A., Muller G., Santoro F., 2016, in Evans C. J., Simard L., Takami H. eds, Proc. SPIE Conf. Ser. Vol. 9908, Ground-based and Airborne Instrumentation for Astronomy VI. SPIE, Bellingham, p. 99086T

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

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

Kipping D. M., 2014, MNRAS, 444, 2263

Kopparapu R. K., Ramirez R. M., SchottelKotte J., Kasting J. F., Domagal-Goldman S., Eymet V., 2014, ApJ, 787, L29

Kreidberg L., 2015, PASP, 127, 1161

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

Laskar J., 1988, A&A, 198, 341

Laskar J., 1990, Icarus, 88, 266
Laskar J., 1993, *Phys. D Nonlinear Phenom.*, 67, 257

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

López-Morales M. et al., 2016, AJ, 152, 204

Louie D. R., Deming D., Albert L., Bouma L. G., Bean J., Lopez-Morales M., 2018, PASP, 130, 044401

Lovis C., Mayor M., Bouchy F., Pepe F., Queloz D., Udry S., Benz W., Mordasini C., 2009, in Pont F., Sasselov D., Holman M. J. eds, *Transiting Planets, Vol. 253*. Cambridge University Press, Cambridge, p. 502

Luque R. et al., 2021, A&A, 645, A41

Mańcz Apellániz J., Weiler M., 2018, A&A, 619, A180

Mandel K., Agol E., 2002, ApJ, 580, L171

Mann A. W., von Braun K., 2015, PASP, 127, 102

Mann A. W., Gaidos E., Ansdell M., 2013, ApJ, 779, 188

Mann A. W., Feiden G. A., Gaidos E., Boyajian T., von Braun K., 2015, ApJ, 804, 64

Mann A. W. et al., 2016, AJ, 152, 61

Mann A. W. et al., 2019, ApJ, 871, 63

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

Mayor M., Udry S., 2008, *Phys. Scr. T*, 1310, 014010

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

McLaughlin D. B., 1924, ApJ, 60, 22

Morbidelli A., Baillié K., Batygin K., Charnoz S., Guillot T., Rubie D. C., Kleine T., 2022, Nat. Astron., 6, 72

Mordasini C., van Boekel R., Mollière P., Henning T., Benneke B., 2016, ApJ, 832, 41

Morrell S., Naylor T., 2019, MNRAS, 489, 2615

Morrell S., Naylor T., 2020, MNRAS, 498, 5048

Morris R. L., Twicken J. D., Smith J. C., Clarke B. D., Jenkins J. M., Bryson S. T., Girouard F., Klaus T. C., 2020, Kepler Data Processing Handbook: Photometric Analysis, Kepler Science Document KSCI-19081-003

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

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

Mortier A., Faria J. P., Correia C. M., Santerne A., Santos N. C., 2015, A&A, 573, A101

Morton T. D., 2015, Astrophysics Source Code Library, record (asci:1503.010)

Norris C. M., Beeck B., Unruh Y. C., Solanki S. K., Krivova N. A., Yeo K. L., 2017, A&A, 605, A45

Noyes R. W., Weiss N. O., Vaughan A. H., 1984, ApJ, 287, 769

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

Owen J. E., Adams F. C., 2019, MNRAS, 490, 15

Pastorelli G. et al., 2019, MNRAS, 485, 5666

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

Pepper J. et al., 2020, AJ, 159, 243

Perruchot S. et al., 2008, in McLean I. S., Casali M. M. eds, *Proc. SPIE Conf. Ser. Vol. 7014*, Ground-based and Airborne Instrumentation for Astronomy II. SPIE, Bellingham, p. 70140J

Petrovich C., 2015, ApJ, 805, 75

Petrovich C., Tremaine S., 2016, ApJ, 829, 132

Pinamonti M. et al., 2018, A&A, 617, A104

Pollacco D. L. et al., 2006, PASP, 118, 1407

Pollacco D. et al., 2008, MNRAS, 385, 1576

Queloz D., Eggenberger A., Mayor M., Perrier C., Beuzit J. L., Naef D., Sivan J. P., Udry S., 2000, A&A, 359, L13

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

Rayner J. T., Cushing M. C., Vacca W. D., 2009, ApJS, 185, 289

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

Rein H., Spiegel D. S., 2015, MNRAS, 446, 1424

Rice K. et al., 2019, MNRAS, 484, 3731

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

Rogers L. A., 2015, ApJ, 801, 41

Rossiter R. A., 1924, ApJ, 60, 15

Serrano L. M., Barros S. C. C., Oshagh M., Santos N. C., Faria J. P., Demangeon O., Sousa S. G., Lendl M., 2018, A&A, 611, A8

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

Sneden C., 1973, ApJ, 184, 839

Stalport M., Delisle J. B., Udry S., Matthews E. C., Bourrier V., Leleu A., 2022, A&A, 664, A53

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

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

Stephenson C. B., 1986, AJ, 92, 139

Sugiura N., 1978, Commun. Stat. – Theory Methods, 7, 13

Tamayo D., Rein H., Shi P., Hernandez D. M., 2020, MNRAS, 491, 2885

Tang J., Bressan A., Rosenfield P., Slemer A., Marigo P., Girardi L., Bianchi L., 2014, MNRAS, 445, 4287

Ter Braak C. J., 2006, Stat. Comput., 16, 239

Thompson S. J. et al., 2016, in Evans C. J., Simard L., Takami H. eds, *Proc. SPIE Conf. Ser. Vol. 9908*, Ground-based and Airborne Instrumentation for Astronomy VI. SPIE, Bellingham, p. 99086F

Thorngren D. P., Fortney J. J., Murray-Clay R. A., Lopez E. D., 2016, ApJ, 831, 64

Tsantaki M., Sousa S. G., Adibekyan V. Z., Santos N. C., Mortier A., Israeli G., 2013, A&A, 555, A150

Twicken J. D., Clarke B. D., Bryson S. T., Tenenbaum P., Wu H., Jenkins J. M., Girouard F., Klaus T. C., 2010, Proc. SPIE, 7740, 774023

Udry S., Mayor M., Santos N. C., 2003, A&A, 407, 369

Ulmer-Moll S. et al., 2022, A&A, 666, A46

Vanderburg A. et al., 2016, ApJS, 222, 14

Vanderburg A. et al., 2018, AJ, 156, 46

Vanderburg A. et al., 2019, ApJ, 881, L19

West R. G. et al., 2019, MNRAS, 486, 5094

Wilson D. M. et al., 2008, ApJ, 675, L113

Winn J. N., 2011, in Sozzetti A., Lattanzi M. G., Boss A. P. eds, *Proc. IAU Symp. 276*, The Astrophysics of Planetary Systems: Formation, Structure, and Dynamical Evolution. Kluwer, Dordrecht, p. 230

Winn J. N., Fabrycky D., Albrecht S., Johnson J. A., 2010, ApJ, 718, L145

Winters J. G. et al., 2015, AJ, 149, 5

Wittenmyer R. A., O'Toole S. J., Jones H. R. A., Tinney C. G., Butler R. P., Carter B. D., Bailey J., 2010, ApJ, 722, 1854

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

Zeng L., Sasselov D. D., Jacobsen S. B., 2016, ApJ, 819, 127

Zhang M., Dai F., Bean J. L., Knutson H. A., Rescigno F., 2023, ApJ, 953, L25

Zhu W., Wu Y., 2018, AJ, 156, 92

## SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

### suppl.data

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.

<footnote n="1">¹Department of Astrophysics, University of Exeter, Stocker Rd, Exeter EX4 4QL, UK</footnote>

<footnote n="2">²Institut d'Astrophysique de Paris, UMR7095 CNRS, Université Pierre and Marie Curie, 98bis Boulevard Arago, F-75014 Paris, France</footnote>

<footnote n="3">³Observatoire de Haute-Provence, CNRS, Université d'Aix-Marseille, 04870 Saint-Michel-l'Observatoire, France</footnote>

<footnote n="4">⁴Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA</footnote>

<footnote n="5">⁵Department of Physics and Astronomy, The University of North Carolina at Chapel Hill, Chapel Hill, NC 27599, USA</footnote>

<footnote n="6">⁶Department of Astrophysics, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK</footnote>

<footnote n="7">⁷DTU Space, National Space Institute, Technical University of Denmark, Elektrovej 328, DK-2800 Kgs. Lyngby, Denmark</footnote>
<footnote n="8">⁸Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="9">⁹Trottier Institute for Research on Exoplanets (iREx)), Université de Montréal, Complexe des Sciences, Montréal, QC H3C 3J7, Canada</footnote>

<footnote n="10">¹⁰Département de Physique, Université de Montréal, 1375 Avenue Thérèse-Lavoie-Roux, Montreal, QC, H2V 0B3, Canada</footnote>

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

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

<footnote n="13">¹³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="14">¹⁴Observatoire de Geneève, Université de Geneève, Chemin de Pegasi, 51, CH-1290 Versoix, Switzerland</footnote>

<footnote n="15">¹⁵Centre for Astrophysics, University of Southern Queensland, West Street, Toowoomba, QLD 4350 Australia</footnote>

<footnote n="16">¹⁶Laboratoire d'Astrophysique de Marseille, Université de Provence, UMR6110 CNRS, 38 rue F. Joliot Curie, F-13388 Marseille cedex 13, France</footnote>

<footnote n="17">¹⁷Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="18">¹⁸Laboratoire d'Astrophysique de Marseille, Aix Marseille Université, CNRS, CNES, 13388 Marseille CEDEX 13, France</footnote>

<footnote n="19">¹⁹Royal Astronomical Society, Burlington House, Piccadilly, London W1J 0BQ, UK</footnote>

<footnote n="20">²⁰SUPA, School of Physics and Astronomy, University of St Andrews, North Haugh, St Andrews, KY169SS, UK</footnote>

<footnote n="21">²¹Centre for Exoplanet Science, University of St Andrews, North Haugh, St Andrews, KY169SS, UK</footnote>

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

<footnote n="23">²³SUPA, Institute for Astronomy, Royal Observatory, University of Edinburgh, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="24">²⁴Centre for Exoplanet Science, University of Edinburgh, Edinburgh EH9 3FD, UK</footnote>

<footnote n="25">²⁵Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150–762 Porto, Portugal</footnote>

<footnote n="26">²⁶Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA</footnote>

<footnote n="27">²⁷INAF – Osservatorio Astrofisico di Torino, Strada Osservatorio, 201-10025 Pino Torinese (TO), Italy</footnote>

<footnote n="28">²⁸Department of Earth, Atmospheric and Planetary Science, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA</footnote>

<footnote n="29">²⁹Department of Aeronautics and Astronautics, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA</footnote>

<footnote n="30">³⁰Laboratório Nacional de Astrofísica, Rua Estados Unidos 154, 37504–364, Itajubá—MG, Brazil</footnote>

<footnote n="31">³¹SETI Institute, Mountain View, CA 94043, USA</footnote>

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

<footnote n="33">³³Fundación Galileo Galilei – INAF (Telescopio Nazionale Galileo), Rambla J. A. F. Perez 7, E-38712 Breña Baja (La Palma), Canary Islands, Spain</footnote>

<footnote n="34">³⁴Instituto de Astrofísica de Canarias, C/Vía Láctea s/n, E-38205 La Laguna (Tenerife), Canary Islands, Spain</footnote>

<footnote n="35">³⁵Departamento de Astrofísica, Universidad de La Laguna, Av. del Astrofísico Francisco Sánchez s/n, E-38205 La Laguna (Tenerife), Canary Islands, Spain</footnote>

<footnote n="36">³⁶Department of Physics and Astronomy, University of New Mexico, 210 Yale Blvd NE, Albuquerque, NM 87106, USA</footnote>

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

<footnote n="38">³⁸INAF – Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, I-35122 Padova, Italy</footnote>

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

<footnote n="40">⁴⁰Centre for Astrophysics, University of Southern Queensland, West Street, Toowoomba, QLD 4350, Australia</footnote>

<footnote n="41">⁴¹Sub-department of Astrophysics, University of Oxford, Keble Rd, OX13RH, Oxford, UK</footnote>

<footnote n="42">⁴²Observatoire des Baronnies Provençales, F-05150 Moydans, France</footnote>

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

<footnote n="44">⁴⁴Kotizarović Observatory, Sarsoni 90, 51216 Viskovo, Croatia</footnote>

<footnote n="45">⁴⁵American Association of Variable Star Observers, 185 Alewife Brook Parkway, Cambridge, MA 02138, USA</footnote>

<footnote n="46">⁴⁶Oukaimeden Observatory, High Energy Physics and Astrophysics Laboratory, Cadi Ayyad University, Abdelkrim Khattabi Avenue, Marrakech, 40000, Morocco</footnote>

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

<footnote n="48">⁴⁸Grand Pra Observatory, 1984 Les Haudàres, Switzerland</footnote>

<footnote n="49">⁴⁹Observatori Astronómico Albanyá, Girona, Spain</footnote>

<footnote n="50">⁵⁰Departamento de Astronomía y Astrofísica, Universidad de Valencia, E-46100 Burjassot, Valencia, Spain</footnote>

<footnote n="51">⁵¹Observatorio Astronómico, Universidad de Valencia, E-46980 Paterna, Valencia, Spain</footnote>

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