# Characterization of the L 98-59 multi-planetary system with HARPS

R. Cloutier, N. Astudillo-Defru, X. Bonfils, Z. Berdiñas, G. Ricker, R.

Vanderspek, D. W. Latham, S. Seager, J. Winn, J. S. Jenkins, et al.

► To cite this version:

R. Cloutier, N. Astudillo-Defru, X. Bonfils, Z. Berdiñas, G. Ricker, et al. Characterization of the L 98-59 multi-planetary system with HARPS. Astronomy and Astrophysics - A&A, 2019, 629, pp.A111. 10.1051/0004-6361/201935957. [hal-02392912](https://hal.science/hal-02392912)

HAL Id: hal-02392912

https://hal.science/hal-02392912

Submitted on 31 Aug 2021

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
# Characterization of the L 98-59 multi-planetary system with HARPS

## Mass characterization of a hot super-Earth, a sub-Neptune, and a mass upper limit on the third planet*,**

R. Cloutier<fnref n="1" />,<fnref n="2" />,³, N. Astudillo-Defru⁴, X. Bonfils⁵, J. S. Jenkins⁶, Z. Berdiñas⁶, G. Ricker⁷, R. Vanderspek⁷,
D. W. Latham⁸, S. Seager⁹, J. Winn¹⁰, J. M. Jenkins¹¹, J. M. Almenara⁵, F. Bouchy¹², X. Delfosse⁵, M. R. Díaz⁶,
R. F. Díaz¹³,¹⁴, R. Doyon³, P. Figueira¹⁵,¹⁶, T. Forveille⁵, N. T. Kurtovic⁶, C. Lovis¹², M. Mayor¹², K. Menou<fnref n="1" />,<fnref n="2" />,
E. Morgan⁷, R. Morris¹⁷,¹¹, P. Muirhead¹⁸, F. Murgas⁵, F. Pepe¹², N. C. Santos¹⁶,¹⁹, D. Ségransan¹², J. C. Smith¹⁷,¹¹,
P. Tenenbaum¹⁷,¹¹, G. Torres⁸, S. Udry¹², M. Vezie⁷, and J. Villasenor⁷

*(Affiliations can be found after the references)*

Received 25 May 2019 / Accepted 13 August 2019

## ABSTRACT

**Aims.** L 98-59 (TIC 307210830, TOI-175) is a nearby M3 dwarf around which TESS revealed three small transiting planets (0.80, 1.35, 1.57 Earth radii) in a compact configuration with orbital periods shorter than 7.5 days. Here we aim to measure the masses of the known transiting planets in this system using precise radial velocity (RV) measurements taken with the HARPS spectrograph.

**Methods.** We considered both trained and untrained Gaussian process regression models of stellar activity, which are modeled simultaneously with the planetary signals. Our RV analysis was then supplemented with dynamical simulations to provide strong constraints on the planets' orbital eccentricities by requiring long-term stability.

**Results.** We measure the planet masses of the two outermost planets to be 2.42 ± 0.35 and 2.31 ± 0.46 Earth masses, which confirms the bulk terrestrial composition of the former and eludes to a significant radius fraction in an extended gaseous envelope for the latter. We are able to place an upper limit on the mass of the smallest, innermost planet of <1.01 Earth masses with 95% confidence. Our RV plus dynamical stability analysis places strong constraints on the orbital eccentricities and reveals that each planet's orbit likely has $e < 0.1$.

**Conclusions.** L 98-59 is likely a compact system of two rocky planets plus a third outer planet with a lower bulk density possibly indicative of the planet having retained a modest atmosphere. The system offers a unique laboratory for studies of planet formation, dynamical stability, and comparative atmospheric planetology as the two outer planets are attractive targets for atmospheric characterization through transmission spectroscopy. Continued RV monitoring will help refine the characterization of the innermost planet and potentially reveal additional planets in the system at wider separations.

**Key words.** stars: individual; L 98-59 – planetary systems – stars: low-mass – planets and satellites: terrestrial planets – techniques: radial velocities

## 1. Introduction

NASA's *Transiting Exoplanet Survey Satellite* (TESS; Ricker et al. 2015) is expected to discover thousands of new transiting planetary systems around nearby stars over ~80% of the entire sky (Sullivan et al. 2015; Ballard 2019; Barclay et al. 2018; Huang et al. 2018a). Throughout its two-year long primary mission, TESS will observe ≥200 000 targets from the TESS Input Catalog (TIC; Stassun et al. 2018) at a two minute cadence as well as many more targets within the 30-min full frame images. Indeed several confirmed planetary systems have already been uncovered by TESS within its first year of operations (Brahm et al. 2019; Jones et al. 2019; Cañas et al. 2019; Dragomir et al. 2019; Espinoza et al. 2019; Kipping et al. 2019; Kostov et al.

2019; Nielsen et al. 2019; Quinn et al. 2019; Rodriguez et al. 2019; Wang et al. 2019; Vanderspek et al. 2019) including a small number of planets that have begun to contribute to the completion of the mission's level one science requirement of delivering the masses of 50 planets smaller than $4 R_⊕$ ($\pi$ Mensae c; Gandolfi et al. 2018; Huang et al. 2018b, TOI-402.01, 02; Dumusque et al. 2019).

The nearby M3 dwarf L 98-59 (TIC 307210830¹, TOI-175², $d=10.6$ pc, Table 1) was included in the TESS Input Catalog based on its stellar parameters from the Cool Dwarf list (Muirhead et al. 2018) and so far has been observed in TESS Sector 2. Three small planetary candidates around L 98-59 (L 98-59c, TOI-175.01: $P_c = 3.69$ days, $r_{p,c} = 1.35 R_⊕$, L 98-59d, TOI-175.02: $P_d = 7.45$ days, $r_{p,d} = 1.57 R_⊕$, L 98-59b, TOI-175.03: $P_b = 2.25$ days, $r_{p,b} = 0.80 R_⊕$) were flagged by the Science Processing Operations Center Pipeline (SPOC; Jenkins et al. 2016) and subsequently passed a set of validation tests

<footnote>* Full Tables 2 and A.1 are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http://cdsarc.u-strasbg.fr/viz-bin/cat/J/A+A/629/A111</footnote>

<footnote>** Based on observations made with the HARPS instrument on the ESO 3.6 m telescope under the program IDs 198.C-0838(A), 1102.C-0339(A), 0102.C-0525, and 0102.D-0483(A) at Cerro La Silla (Chile).</footnote>

<footnote n="1">¹ TIC: TESS Input Catalog.</footnote>

<footnote n="2">² TOI: TESS Object of Interest.</footnote>
<table><caption>**Table 1.** L 98-59 stellar parameters.</caption><thead><tr><th>Parameter</th><th>Value</th><th>References</th></tr></thead><tbody><tr><td colspan="3"><em>L 98-59, TIC 307210 830, TOI-175</em></td></tr><tr><td colspan="3"><em>Astrometry</em></td></tr><tr><td>RA, α (deg)</td><td>124.532860</td><td>1,2</td></tr><tr><td>Dec, δ (deg)</td><td>-68.314466</td><td>1,2</td></tr><tr><td>RA proper motion, μ<sub>α</sub> (mas yr<sup>-1</sup>)</td><td>94.767 ± 0.054</td><td>1,2</td></tr><tr><td>Dec proper motion, μ<sub>δ</sub> (mas yr<sup>-1</sup>)</td><td>-340.470 ± 0.052</td><td>1,2</td></tr><tr><td>Parallax, ω (mas)</td><td>94.167 ± 0.028</td><td>1,2,3</td></tr><tr><td>Distance, d (pc)</td><td>10.619 ± 0.003</td><td>1,2,3</td></tr><tr><td colspan="3"><em>Photometry</em></td></tr><tr><td><em>B</em></td><td>13.289 ± 0.027</td><td>4</td></tr><tr><td><em>g'</em></td><td>12.453 ± 0.019</td><td>4</td></tr><tr><td><em>V</em></td><td>11.685 ± 0.017</td><td>4</td></tr><tr><td><em>r'</em></td><td>11.065 ± 0.044</td><td>4</td></tr><tr><td><em>G<sub>BP</sub></em></td><td>11.977 ± 0.002</td><td>1,5</td></tr><tr><td><em>G</em></td><td>10.5976 ± 0.0008</td><td>1,5</td></tr><tr><td><em>G<sub>RP</sub></em></td><td>9.472 ± 0.001</td><td>1,5</td></tr><tr><td><em>T</em></td><td>9.393</td><td>6</td></tr><tr><td><em>J</em></td><td>7.933 ± 0.027</td><td>7</td></tr><tr><td><em>H</em></td><td>7.359 ± 0.049</td><td>7</td></tr><tr><td><em>K<sub>s</sub></em></td><td>7.101 ± 0.018</td><td>7</td></tr><tr><td><em>W<sub>1</sub></em></td><td>6.935 ± 0.062</td><td>8</td></tr><tr><td><em>W<sub>2</sub></em></td><td>6.767 ± 0.021</td><td>8</td></tr><tr><td><em>W<sub>3</sub></em></td><td>6.703 ± 0.016</td><td>8</td></tr><tr><td><em>W<sub>4</sub></em></td><td>6.578 ± 0.047</td><td>8</td></tr><tr><td colspan="3"><em>Stellar parameters</em></td></tr><tr><td>Spectral type</td><td>M3V ± 1</td><td>9</td></tr><tr><td>Stellar radius, <em>R</em><sub>s</sub> (<em>R</em><sub>☉</sub>)<sup>(*)</sup></td><td>0.314 ± 0.014</td><td>3,9</td></tr><tr><td>Effective temperature, <em>T</em><sub>eff</sub> (K)</td><td>3412 ± 49</td><td>3,10</td></tr><tr><td>Stellar mass, <em>M</em><sub>s</sub> (<em>M</em><sub>☉</sub>)<sup>(†)</sup></td><td>0.312 ± 0.031</td><td>3,11</td></tr><tr><td>Surface gravity, log <em>g</em> (dex)</td><td>4.94 ± 0.06</td><td>12</td></tr><tr><td>Metallicity, [Fe/H]</td><td>-0.5 ± 0.5</td><td>9</td></tr><tr><td>log <em>R'</em><sub>HK</sub></td><td>-5.40 ± 0.11</td><td>12</td></tr><tr><td>Rotation period, <em>P</em><sub>rot</sub> (days)<sup>(‡)</sup></td><td>78 ± 13</td><td>12,13</td></tr></tbody></table>

**Notes.** (*) Includes the radius uncertainty from Kostov et al. (2019). (†) We add in quadrature a fractional uncertainty of 10% based on the dispersion in stellar masses for stars with metallicities that differ from solar (Mann et al. 2019). (‡) The predicted $P_{rot}$ based on $log R'_{HK}$ and the M dwarf activity-rotation relation from Astudillo-Defru et al. (2017a). This value is consistent with the periodogram peak in the $H\alpha$ time series in Fig. 1.

**References.** (1) Gaia Collaboration (2018), (2) Lindegren et al. (2018), (3) Cloutier (2019), (4) Henden et al. (2016), (5) Evans et al. (2018), (6) Stassun et al. (2018), (7) Cutri et al. (2003), (8) Cutri et al. (2013), (9) Kostov et al. (2019), (10) Mann et al. (2015), (11) Benedict et al. (2016), (12) this work, (13) Astudillo-Defru et al. (2017a).

(Twicken et al. 2018; Li et al. 2019) prior to being published as TESS Data Alerts. Many of the properties of this multi-planet system make it of interest for radial velocity (RV) mass characterization (Cloutier et al. 2018), planetary atmospheric characterization (Kempton et al. 2018; Louie et al. 2018), and direct investigations of M dwarf planet formation, evolution, and system architectures (Lissauer et al. 2011; Fabrycky et al. 2014). As such, the system warranted an intensive follow-up

campaign presented by Kostov et al. (2019, hereafter K19) that ruled out astrophysical false positive scenarios and confirmed the planetary nature of each of the three planet candidates.

In this paper we present the results of our follow-up study to obtain precise planet masses for as many of the L 98-59 planets as possible using HARPS precision RVs. In practice we are only able to recover robust masses for the two outermost planets TOI-175.01 and 02 but we also report our derived upper limits on the mass of the smallest known planet TOI-175.03. In Sect. 2 we discuss our spectroscopic HARPS observations, in Sect. 3 we establish our model of the observed RVs before presenting our results in Sect. 4. In Sect. 5 we use the measured planet masses to perform a dynamical stability analysis of the system to provide stronger constraints on each planet's orbital eccentricity. We then conclude with a discussion in Sect. 6.

## 2. HARPS observations

### 2.1. HARPS data acquisition

Using the High Accuracy Radial velocity Planet Searcher (HARPS; Mayor et al. 2003) échelle spectrograph mounted at the 3.6 m ESO telescope at La Silla Observatory, Chile, we obtained a set of 164 spectra of L 98-59 between October 17, 2018 (BJD = 2 458 408.5) and April 28, 2019 (BJD = 2 458 601.5). The HARPS optical spectrograph at R = 115 000 is stabilized in pressure and temperature which helps enable its sub-m s⁻¹ accuracy.

Throughout the five-month observing campaign of L 98-59 we elected not to use a simultaneous wavelength calibration (i.e., on-sky calibration fiber) to prevent possible contamination of the bluer spectral orders by the calibration lamp. In the ESO programs 198.C-0838 and 1102.C-0339 (140/164 observations) the exposure time was set to 900 s, resulting in a median signal-to-noise ratio (S/N) of 41 per resolution element at 650 nm and a median measurement uncertainty of 2.08 m s⁻¹. In the ESO program 0102.C-0525 (21/164 observations) the exposure time ranged between 500 and 1800 s, with a median S/N of 49 per resolution element at 650 nm and a median measurement uncertainty of 1.61 m s⁻¹. In the ESO program 0102.D-0483 (3/164 observations) the exposure time was set to 2200 s, resulting in an average S/N of 67 per resolution element at 650 nm and an average measurement uncertainty of 1.15 m s⁻¹.

### 2.2. Radial velocity extraction

To compute the RV time series we performed a maximum likelihood analysis between a stellar template and individual spectra following Astudillo-Defru et al. (2017b). The adopted stellar template corresponded to the median of all spectra that were previously shifted to the star frame. A telluric template was derived by the median of spectra shifted to the Earth frame. For these two steps we used the stellar radial velocity derived by the HARPS Data Reduction Software (DRS; Lovis & Pepe 2007) through a cross-correlation function. We used the barycentric Earth radial velocity as computed by the DRS as well. The resulting stellar template was Doppler shifted over a window of 40 km s⁻¹ wide and centered on the average of the RVs computed by the DRS (−5.661 km s⁻¹). The telluric template was used to mask the spectral zones contaminated by telluric lines. For each RV step we computed the value of the likelihood function with the maximum of the likelihood function representing the RV of the spectrum under analysis. The process was repeated for the entire HARPS dataset and resulted in the RV time series reported in Table 2 that are used in the subsequent analysis.
<table><caption>**Table 2.** HARPS spectroscopic time series.</caption><thead><tr><th>Time<br>(BJD -<br>2 457 000)</th><th>RV<sup>(*)</sup><br>(m s<sup>-1</sup>)</th><th>σ<sub>RV</sub><br>(m s<sup>-1</sup>)</th><th>Hα<br>×10<sup>2</sup></th><th>σ<sub>Hα</sub><br>×10<sup>2</sup></th><th>Hβ<br>×10<sup>2</sup></th><th>σ<sub>Hβ</sub><br>×10<sup>2</sup></th><th>Hγ<br>×10<sup>2</sup></th><th>σ<sub>Hγ</sub><br>×10<sup>2</sup></th><th>NaD<br>×10<sup>2</sup></th><th>σ<sub>NaD</sub><br>×10<sup>2</sup></th><th>S-index<br>—</th><th>σ<sub>S</sub><br>—</th><th>FWHM<br>(km s<sup>-1</sup>)</th><th>BIS<br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1408.853661</td><td>-5678.7</td><td>2.3</td><td>6.93</td><td>0.02</td><td>5.43</td><td>0.05</td><td>12.32</td><td>0.17</td><td>0.91</td><td>0.02</td><td>0.69</td><td>0.09</td><td>3.0588</td><td>23.2798</td></tr><tr><td>1409.844622</td><td>-5678.9</td><td>2.4</td><td>6.77</td><td>0.02</td><td>5.22</td><td>0.05</td><td>11.46</td><td>0.17</td><td>0.84</td><td>0.02</td><td>0.65</td><td>0.10</td><td>3.0637</td><td>23.3301</td></tr><tr><td>1412.858886</td><td>-5679.1</td><td>2.2</td><td>7.07</td><td>0.02</td><td>5.79</td><td>0.05</td><td>12.85</td><td>0.17</td><td>0.98</td><td>0.02</td><td>0.80</td><td>0.10</td><td>3.0559</td><td>23.3108</td></tr><tr><td>1413.860798</td><td>-5676.2</td><td>1.9</td><td>6.79</td><td>0.02</td><td>5.25</td><td>0.04</td><td>11.87</td><td>0.13</td><td>0.83</td><td>0.02</td><td>0.65</td><td>0.07</td><td>3.0655</td><td>23.3935</td></tr><tr><td>1414.858984</td><td>-5674.6</td><td>2.5</td><td>6.81</td><td>0.02</td><td>5.12</td><td>0.06</td><td>11.75</td><td>0.19</td><td>0.85</td><td>0.03</td><td>0.71</td><td>0.12</td><td>3.0560</td><td>23.1147</td></tr></tbody></table>

Notes. Only the first five rows are depicted here for clarity. (*)Systemic velocity, $γ_0 = -5678.4 \pm 0.2$ m s$^{-1}$. The full table is available at the CDS.

# 3. Model setup

## 3.1. Periodogram analysis

To identify strong periodicities in our RV time series we compute the generalized Lomb-Scargle periodogram (GLSP; Zechmeister & Kürster 2009) of all spectroscopic time series derived from our HARPS spectra and of its window function (WF). The ancillary spectroscopic time series of $Hα$, $Hβ$, $Hγ$, the sodium doublet NaD, and the S-index based on the Ca H & K doublet are sensitive to chromospheric activity and may therefore be used to identify periodicities in the RV data arising from chromospheric activity sources such as plages.

The $Hα$ index was computed using the same pass-bands as Gomes da Silva et al. (2012), that is a band of 1.6 Å wide centered on 6562.8 Å and two control bands of widths 10.75 and 8.75 Å centered on 6550.87 and 6580.31 Å, respectively; the central band for $Hβ$ was limited by 4861.04–4861.60 Å, and we defined two control bands limited by 4855.04–4860.04 Å and 4862.6–4867.2 Å; for $Hγ$ we integrated over three bands bounded by 4333.60–4336.80 Å, 4340.16–4340.76 Å, and 4342.00–4344.00 Å; the central bands to calculate the NaD-index were similar as Gomes da Silva et al. (2012), namely wide of 0.5 Å and centered on 5889.95 and 5895.92 Å, but with control bands limited by 5860.0–5870.0 Å and 5904.0–5908.0 Å; the S-index was calculated following Duncan et al. (1991) and the calibration derived in Astudillo-Defru et al. (2017a, Eq. (3)) that scales the index computed from HARPS spectra to Mount Wilson.

Similarly, the full width at half maximum (FWHM) and bisector (BIS) shape parameters of the spectral CCF may be sensitive to chromospheric and/or photospheric active regions such as dark spots (Queloz et al. 2001; Desort et al. 2007). Periodicities in either the FWHM or BIS time series may therefore also allude to periodic signals arising from stellar activity. We analyzed the FWHM and BIS as derived by the HARPS DRS. The GLSP of the WF is also computed to potentially identify sources of aliasing from our time sampling.

The resulting GLSPs are shown in Fig. 1 along with their false alarm probability (FAP) curves. The FAP curves are computed via bootstrapping with replacement using $10^4$ iterations and normalizing each GLSP's power scale by its standard deviation. Although the detection of individual planetary signals at a low FAP is not required to claim a planet's detection in the RVs, the GLSP of the RVs does reveal a moderately significant peak close to the orbital period of L 98-59d (~7.45 days; FAP ~1%) plus significant peaks at the orbital period of L 98-59c (~3.69 days; FAP ≤0.1%) and centered around ~40 days (FAP ≪0.1%). The photometric rotation period of L 98-59 remains undetected in the SAP TESS light curve and the star

exhibits a negligible rotational broadening ($v \sin i < 1.9$ km s$^{-1}$; K19) indicative of L 98-59 being a largely inactive, old M dwarf with a likely rotation period ≥10 days and a correspondingly low amplitude of photometric variability (Newton et al. 2016). However, the strongest periodic signal as seen in any activity sensitive time series in Fig. 1 is a feature centered around ~80 days in the $Hα$ GLSP (FAP ≪0.1%). This signal is consistent with the expected $P_{rot} = 78 \pm 13$ days based on the star's value of $\log R'_{HK} = -5.4 \pm 0.11$ and using the M dwarf magnetic activity-rotation relation from Astudillo-Defru et al. (2017a). If this periodic signal is indeed due to stellar rotation at $P_{rot} \sim 80$ days then it could explain the ~40 day signal in the RV GLSP as being the first harmonic of $P_{rot}$. Nearly all of the remaining activity sensitive time series have GLSPs that are consistent with noise (i.e., FAP ≥10%) implying that no strong periodic signals are resolved in those time series although two notable exceptions exist. The first is the peak at ≥100 days in the NaD time series that likely arising from the excess power seen in the WF at short frequencies around 1/100 days$^{-1}$. The second is a less significant peak that is intermediate between 40 and 80 days and persists in each of the $Hβ$, $Hγ$, S-index, and FWHM GLSPs. The origin of this weak, intermediate peak is unknown but may also be related to stellar rotation as the ~80 day $Hα$ peak is posited to be.

A periodic signal from the remaining planet L 98-59b (~2.25 days) is not seen at a low FAP in the GLSP of the raw RVs (Fig. 1). However our iterative periodogram analysis will ultimately reveal the presence of this signal, as well as the increased strength of the L 98-59d planetary signal (~7.45 days), upon the joint modeling of the planets with RV stellar activity in Sect. 4.

## 3.2. Stellar activity

Stellar activity on M dwarfs predominantly arises from active regions in the stellar photosphere and chromosphere (Lindegren & Dravins 2003). The resulting RV signal is modulated by stellar rotation as active regions traverse the visible stellar disk and disrupt its symmetry thus creating a temporarily correlated RV variation that can mask or even mimic planetary signals under certain circumstances (Vanderburg et al. 2016). In this study, we will use the $Hα$ activity indicator, which is unaffected by planet-induced Doppler shifts, to inform our stellar activity model.

Following its successful application to activity modeling in M dwarf planetary systems (e.g., Astudillo-Defru et al. 2017c; Cloutier et al. 2017, 2019; Bonfils et al. 2018; Ment et al. 2019), we adopt a semi-parametric Gaussian process (GP) regression model of RV stellar activity to simultaneously model activity with the RV planetary signals. Given a parameterization of
<figure><img src="image_2.png" /><figcaption>Fig. 1. Generalized Lomb-Scargle periodograms of HARPS spectroscopic time series. Left column: GLSPs of our HARPS RV time series, its window function, and $H\alpha$, $H\beta$, $H\gamma$, sodium doublet, S-index, full width half maximum, and bisector activity indicators. The vertical dotted lines highlight the orbital frequencies (i.e., inverse periods) of the three known transiting planets plus the posited rotation period of L 98-59 at $P_{\text{rot}} \sim 80$ days and its first harmonic $P_{\text{rot}}/2$. Right column: false alarm probabilities computed from bootstrapping with replacement.</figcaption></figure>

the temporal covariance structure in our time series, the semi-parametric nature of the GP activity treatment is well-suited to modeling a stochastic physical process like stellar activity without requiring a deterministic functional form. Here we assume that the apparent stellar activity signal seen in the RV data at ~40 days has a manifestation in the $H\alpha$ data at ~80 days and whose temporal covariance structure is quasi-periodic as it is rotationally modulated and yet it is not purely periodic due to evolution in the active region sizes, contrasts, and spatial distribution over multiple rotation cycles (Giles et al. 2017). The corresponding covariance kernel function of our GP model,

trained on the $H\alpha$ time series, is

$$k(t_i, t_j) = a^2 \exp \left[ -\frac{(t_i - t_j)^2}{2\lambda^2} - \Gamma^2 \sin^2 \left( \frac{\pi|t_i - t_j|}{P_{\text{rot}}} \right) \right] \quad (1)$$

and is parameterized by the following GP hyperparameters: the covariance amplitude $a$, an exponential decay time scale $\lambda$, a coherence parameter $\Gamma$, and the strong periodic signal seen in the $H\alpha$ GLSP (Fig. 1) which we attribute to the L 98-59 rotation period $P_{\text{rot}}$.

We sample the posterior probability density function (PDF) of the logarithmic GP hyperparameters by running the $H\alpha$ time
<table><caption>**Table 3.** L 98-59 model parameter priors.</caption><thead><tr><th>Parameter</th><th>Prior</th></tr></thead><tbody><tr><td colspan="2"><em>Hα training model</em></td></tr><tr><td>ln covariance amplitude, ln a</td><td>U(-10, -2)</td></tr><tr><td>ln exponential time scale, ln λ</td><td>U(1, 20)</td></tr><tr><td>ln coherence, ln Γ</td><td>U(-10, 10)</td></tr><tr><td>ln rotation period, ln P<sub>rot</sub></td><td>U(2, 5)</td></tr><tr><td>Additive jitter, s</td><td>J(10<sup>-6</sup>, 10<sup>-3</sup>)</td></tr><tr><td colspan="2"><em>RV model</em></td></tr><tr><td>ln covariance amplitude, ln a</td><td>U(-5, 5)</td></tr><tr><td>ln exponential time scale, ln λ</td><td>p(ln λ|Hα)</td></tr><tr><td>ln coherence, ln Γ</td><td>p(ln Γ|Hα)</td></tr><tr><td>ln periodic time scale, ln P<sub>GP</sub></td><td>p(ln [P<sub>rot</sub>/2]|Hα)</td></tr><tr><td>Additive jitter, s (m s<sup>-1</sup>)</td><td>J(10<sup>-2</sup>, 10)</td></tr><tr><td>Systemic velocity, γ<sub>0</sub> (m s<sup>-1</sup>)</td><td>U(-10, 10)</td></tr><tr><td colspan="2"><em>L 98-59b (TOI-175.03)</em></td></tr><tr><td>Orbital period, P<sub>b</sub> (days)</td><td>N(2.2532, 3 × 10<sup>-4</sup>)</td></tr><tr><td>Time of mid-transit, T<sub>0,b</sub></td><td></td></tr><tr><td>[BJD - 2 457 000]</td><td>N(1366.1708, 1 × 10<sup>-4</sup>)</td></tr><tr><td>Semi-amplitude, K<sub>b</sub> (m s<sup>-1</sup>)</td><td>J(0.1, 10)</td></tr><tr><td>h<sub>b</sub> = √e<sub>b</sub> cos ω<sub>b</sub><sup>(*)</sup></td><td>U(-1, 1)</td></tr><tr><td>k<sub>b</sub> = √e<sub>b</sub> sin ω<sub>b</sub><sup>(*)</sup></td><td>U(-1, 1)</td></tr><tr><td colspan="2"><em>L 98-59c (TOI-175.01)</em></td></tr><tr><td>Orbital period, P<sub>c</sub> (days)</td><td>N(3.6904, 2.5 × 10<sup>-4</sup>)</td></tr><tr><td>Time of mid-transit, T<sub>0,c</sub></td><td></td></tr><tr><td>[BJD - 2 457 000]</td><td>N(1367.2751, 6 × 10<sup>-4</sup>)</td></tr><tr><td>Semi-amplitude, K<sub>c</sub> (m s<sup>-1</sup>)</td><td>J(0.1, 10)</td></tr><tr><td>h<sub>c</sub> = √e<sub>c</sub> cos ω<sub>c</sub><sup>(*)</sup></td><td>U(-1, 1)</td></tr><tr><td>k<sub>c</sub> = √e<sub>c</sub> sin ω<sub>c</sub><sup>(*)</sup></td><td>U(-1, 1)</td></tr><tr><td colspan="2"><em>L 98-59d (TOI-175.02)</em></td></tr><tr><td>Orbital period, P<sub>d</sub> (days)</td><td>N(7.4513, 7 × 10<sup>-4</sup>)</td></tr><tr><td>Time of mid-transit, T<sub>0,d</sub></td><td></td></tr><tr><td>[BJD - 2 457 000]</td><td>N(1362.7375, 8 × 10<sup>-4</sup>)</td></tr><tr><td>Semi-amplitude, K<sub>d</sub> (m s<sup>-1</sup>)</td><td>J(0.1, 10)</td></tr><tr><td>h<sub>d</sub> = √e<sub>d</sub> cos ω<sub>d</sub><sup>(*)</sup></td><td>U(-1, 1)</td></tr><tr><td>k<sub>d</sub> = √e<sub>d</sub> sin ω<sub>d</sub><sup>(*)</sup></td><td>U(-1, 1)</td></tr></tbody></table>

**Notes.** (*)We also require that the corresponding eccentricity value $e_i < 1$.

series through a Markov chain Monte Carlo (MCMC) simulation using the emcee ensemble sampler package (Foreman-Mackey et al. 2013). The prior PDFs on each of the GP hyperparameters are reported in Table 3. The ln likelihood function used to sample their joint posterior PDF is given by

$$ \ln \mathcal{L} = - \frac{1}{2} (\mathbf{y}^{\mathrm{T}} \cdot \mathbf{K} \cdot \mathbf{y} + \ln \det \mathbf{K} + N \ln 2\pi), \quad (2) $$

where $\mathbf{y}$ is the vector of $N H\alpha$ measurements taken at times $t = \{t_1, t_2, \dots, t_N\}$ and the $N \times N$ covariance matrix $\mathbf{K}$ is given by

$$ \mathbf{K}_{ij} = k(t_i, t_j) + \delta_{ij}(\sigma_{H\alpha}(t_i)^2 + s^2). \quad (3) $$

The inclusion of the Kronecker delta $\delta_{ij}$ term adds the $H\alpha$ measurement uncertainties $\sigma_{H\alpha}$ to the diagonal elements of $\mathbf{K}$ and includes an additive jitter factor $s$. Hence the full set parameters sampled during the $H\alpha$ training phase is $\{\ln a, \ln \lambda, \ln \Gamma, \ln P_{\mathrm{rot}}, s\}$.

The underlying physical process of active regions in the stellar chromosphere and photosphere rotating in and out of view at $P_{\mathrm{rot}}$ is responsible for the observed temporal variations in the $H\alpha$ time series. In the following analysis we assume that this process also has a manifestation in the observed RVs such that we can use the constraints on the GP hyperparameters from training on $H\alpha$ to inform our RV model of stellar activity. In particular, the posterior PDFs of the covariance parameters $\{\ln \lambda, \ln \Gamma, \ln P_{\mathrm{rot}}\}$ will be used throughout our modeling of the stellar RVs to derive self-consistent activity and planetary solutions with minimal contamination of the latter by the former as a result of training.

### 3.3. Radial velocity model

Following the training of the GP activity model on the $H\alpha$ time series we can proceed with modeling the RVs. Our RV model contains four physical components from stellar activity plus the three known transiting planets around L 98-59. The RV GP activity model features the same covariance function as was adopted during training (Eq. (1)) and therefore contains the five hyperparameters $\{a, \lambda, \Gamma, P_{GP}, s\}$ where $a$ and $s$ are unique to the RVs whereas the priors on the remaining hyperparameters $\lambda, \Gamma$ and $P_{GP}$ are constrained by their joint posterior PDF from training. Recall that the apparent rotation signal in $H\alpha$ at $P_{\mathrm{rot}} \sim 80$ days appears to be manifested at its first harmonic ($P_{\mathrm{rot}}/2$) in the RVs at $\sim 40$ days (see Fig. 1). As such, we modify the marginalized posterior on $P_{\mathrm{rot}}$ from training by rescaling $P_{GP} \rightarrow P_{\mathrm{rot}}/2$ and use the modified PDF as a prior on $P_{GP}$ in our RV model.

The three planetary signals in our RV model are treated as independent Keplerian orbital solutions. This simplification neglects any gravitational interactions between the planets and makes the sampling of the planetary parameter posterior PDFs much more computationally tractable by negating the need to run dynamical simulations at every step in the MCMC chains. We can justify this simplification by first noting that the planets do not appear to exhibit significant transit timing variations (TTVs) at the level of precision for which such TTVs would be resolvable with the TESS photometric precision ($\sim 5.1$ min for L 98-59b; K19). Furthermore, K19 performed long-term dynamical simulations of the system by considering the maximum a-posteriori (MAP) planet mass predictions and their $+1\sigma$ values for each of the known L 98-59 planets. The planet mass predictions were based on the planets' measured radii and the use of the forecaster tool (Chen & Kipping 2017)<fnref n="3" />. Combining each planet's predicted mass with their osculating orbital elements, and assuming initially circular orbits, the orbital eccentricities of the three planets remained nearly circular (i.e., $\le 0.006$) after one million orbits of the outermost planet (i.e., $\sim 21$ thousand years). In contrast, K19 reported that half of the simulations with initial eccentricities of 0.1 became unstable. These results support the notion that the orbits of the L 98-59 planets are nearly circular, a result that itself is consistent with other compact multi-planet systems exhibiting low eccentricities ($\le 0.05$; Hadden & Lithwick 2014; Van Eylen & Albrecht 2015).

As a back of the envelope calculation, we compare the RV semi-amplitudes $K$ under circular orbits to orbits with $e = 0.1$ for each of the three L 98-59 planets. Assuming the MAP predicted planet masses, the difference in $K$ between $e = 0$ and $e = 0.1$ is $\le 1$ cm $s^{-1}$ for all planets. Even for eccentricities of $e = 0.23$, for which an orbit crossing would occur for either planet pair, the

<footnote n="3">³ The MAP predicted planet masses for L 98-59b, c and d are 0.3, 2.0, and 2.3 $M_{\oplus}$ respectively. Their MAP $+1\sigma$ predicted masses are 0.5, 3.6, and 4.2 $M_{\oplus}$.</footnote>
difference in RV semi-amplitudes between that and the circular orbit scenario is $\le 5$ cm s$^{-1}$. These discrepancies between the circular and maximally elliptical system architectures are well below the typical HARPS RV measurement uncertainty of 2.06 m s$^{-1}$ such that differences in the amplitudes of planet-induced stellar RV signals are negligible. Furthermore, our RV observations only span ~6 months which is not enough time to allow for significant dynamical evolution of the planetary orbits away from their osculating that are measured in Sect. 4. For these reasons we expect the difference between the superposition of Keplerian planet solutions and N-body integrations, for which small non-zero eccentricities would develop, to be negligible. By adopting the simplification of Keplerian orbits, in our MCMC we are effectively sampling the orbital parameters of each planet's osculating orbit rather than tracking the time evolution of those orbital parameters due to mutual planetary interactions. This simplification holds given that the dynamical variations in those orbits are small compared to the level of precision of our data.

Our complete RV model therefore includes the five GP hyperparameters of stellar activity, the L 98-59 systemic velocity $\gamma_0$ plus five Keplerian parameters for each known planet. Namely, each planet's orbital period $P_i$ ($i$ in the planet index; $i = b, c, d$), time of mid-transit $T_{0,i}$, RV semi-amplitude $K_i$, $h_i = \sqrt{e_i} \cos \omega_i$, and $k_i = \sqrt{e_i} \sin \omega_i$ (Ford 2006) where $e_i$ and $\omega_i$ are the planet's orbital eccentricity and argument of periastron respectively. Our complete RV model therefore contains 21 parameters $\{\ln a, \ln \lambda, \ln \Gamma, \ln P_{GP}, s, \gamma_0, P_b, T_{0,b}, K_b, h_b, k_b, P_c, T_{0,c}, K_c, h_c, k_c, P_d, T_{0,d}, K_d, h_d, k_d\}$. We adopt Gaussian priors on each planet's orbital period and time of mid-transit based on the results of their transit light curve analysis (K19) as those data have much more constraining power on the planet ephemerides than do the RVs alone. We adopt broad uninformative priors on the remaining Keplerian parameters which are reported in Table 3.

## 4. Results

The resulting joint and marginalized posterior PDFs from our MCMC analysis are depicted in Fig. A.1. Unless stated otherwise, point estimates of each parameter correspond to their MAP values and are reported in Table 4 along with their uncertainties from the 16th and 84th percentiles of their marginalized posterior PDF.

The marginalized posterior PDF of the $\ln P_{GP}$ hyperparameter does not have a well-defined solution that was expected from the training phase. Furthermore, the $\ln \Gamma$ posterior PDF is highly asymmetric and clearly favours larger values than were favoured by the $H\alpha$ time series. Although we are not principally interested in the values of the RV GP hyperparameters, their values might exhibit a direct effect on the planetary parameters. To investigate this we also consider an alternative model consisting of the three planets plus an untrained GP activity model. This analysis is carried out identically to when using the trained GP except that the priors on $\{\ln \lambda, \ln \Gamma, \ln P_{GP}\}$ are modified to the following uninformative priors: $U(1, 20)$, $U(-10, 10)$, and $U(0, 5)$ respectively. The resulting point estimates of the model parameters are also reported in Table 4 and we find that all Keplerian planet parameters are consistent at the $1\sigma$ level between the two models considered.

We also estimated the Bayesian evidence $Z$ for each RV model featuring a trained and untrained GP activity component respectively. The evidences were computed using the Perrakis et al. (2014) estimator and the marginalized posterior PDFs from our MCMC analyses as importance samplers. This evidence

estimator has been shown to result in quantitatively similar results to other more robust but computationally expensive methods (e.g., nested samplers; Nelson et al. 2018). The resulting Bayes factor, or evidence ratio, between competing models containing a trained and untrained GP activity component is 0.3 thus indicating that inferences resulting from either model are nearly equivalent. Following the consistency of the two models considered we opt to focus on the results from the trained RV model in the subsequent analysis and discussion.

Figure 2 depicts each of the MAP components of our RV model along with its corresponding GLSP. The stellar activity component (i.e., the RVs less the three Keplerian solutions) has a maximum amplitude of ~6 m s$^{-1}$ and is dominated by the clear periodicity at ~40 days that was seen in the GLSP of the raw RVs. Removal of the mean GP activity model mitigates that signal at ~40 days. The planetary RV components from L 98-59c and d are each dominated by their known orbital periods with some aliasing at shorter orbital periods that are consequently mitigated once the planet's MAP Keplerian solution is subtracted off. The Keplerian model of the remaining planet L 98-59b has a small median RV semi-amplitude of $K_b = 0.48$ m s$^{-1}$ making its orbital period only slightly resolved in its GLSP. Lastly, the GLSP of the RV residuals is consistent with noise indicating that we have modeled all major sources of RV variation in our HARPS time series.

Only the two outer planets have RV semi-amplitude “detections” in that their measured values are robustly >0 m s$^{-1}$. These values are $K_c = 2.21 \pm 0.28$ m s$^{-1}$ and $K_d = 1.67 \pm 0.31$ m s$^{-1}$ and represent 7.9 and $5.4\sigma$ detections respectively. The marginalized posterior PDF of $K_b$, corresponding to the smallest and innermost planet in our model, has a median value of $0.48$ m s$^{-1}$ but is consistent with $0$ m s$^{-1}$ therefore resulting in a non-detection of $K_b$ in our dataset. Instead we are only able to place an upper limit on $K_b$ of <1.06 at 95% confidence. The phase-folded RVs are shown in Fig. 3 along with the MAP L 98-59c and d Keplerian orbital solutions and the median L 98-59b Keplerian orbit. The periodic modulation from L 98-59c and d are clearly discernible. Meanwhile the median value of $K_b$ corresponds to a $\le 1.5\sigma$ detection and is not discernible in the phase-folded RVs.

The planetary masses corresponding to the RV semi-amplitudes measured with our data are $m_{p,b} < 1.01 M_\oplus$ (at 95% confidence), $m_{p,c} = 2.42_{-0.34}^{+0.35} M_\oplus$, and $m_{p,d} = 2.31_{-0.45}^{+0.46} M_\oplus$. K19 measure planet-star radius ratios ($r_{p,b}/R_s = 0.0234 \pm 0.0009$, $r_{p,c}/R_s = 0.0396 \pm 0.0010$, and $r_{p,d}/R_s = 0.0462 \pm 0.0029$) from which they derive planetary radii of $r_{p,b} = 0.80 \pm 0.05 R_\oplus$, $r_{p,c} = 1.35 \pm 0.07 R_\oplus$, and $r_{p,d} = 1.57 \pm 0.14 R_\oplus$. The measured planetary masses and radii result in constraints on the planets' bulk densities of $\rho_{p,b} < 12.7$ g cm$^{-3}$, $\rho_{p,c} = 5.4_{-1.0}^{+1.3}$ g cm$^{-3}$, $\rho_{p,d} = 3.3_{-0.9}^{+1.3}$ g cm$^{-3}$. Our mass measurements allow the L 98-59 planets to be added to the planetary mass-radius plane in Fig. 4. In doing so we confirm that the bulk composition of the middle planet L 98-59c is consistent with terrestrial and is therefore a bona-fide super-Earth whose interior appears to be dominated by silicate plus an iron core. The bulk composition of the outermost planet L 98-59d is inconsistent with a terrestrial composition. This sub-Neptune planet must instead contain either a significant fraction of its size in water or have retained a significant gaseous atmosphere. Although the non-detection of $K_b$ prevents a precise value of its bulk density from being derived, the close proximity of L 98-59b to its host star and its intermediate size between that of Earth and Mars are evidence for its terrestrial nature (Owen & Wu 2013, 2017; Jin et al. 2014; Lopez & Fortney 2014; Chen & Rogers 2016; Lopez & Rice 2018). Further constraints on $\rho_{p,b}$ may be realized as we note that the upper limit on
<table><caption>**Table 4.** Model parameters for the L 98-59 three-planet system.</caption><thead><tr><th></th><th>Three-planet model plus GP activity model trained on Hα</th><th>Three-planet model plus an untrained GP activity model</th></tr></thead><tbody><tr><td>In Bayesian evidence, ln Z</td><td>-411.5</td><td>-410.3</td></tr><tr><td>Systemic velocity, γ₀ (m s⁻¹)</td><td>-5678.4 ± 0.2</td><td>-5678.3 ± 0.3</td></tr><tr><td colspan="3"><strong>GP hyperparameters</strong></td></tr><tr><td>Covariance amplitude, a (m s⁻¹)</td><td>7.6<sup>+4.8</sup><sub>-2.6</sub></td><td>6.9<sup>+4.0</sup><sub>-2.3</sub></td></tr><tr><td>Exponential timescale, λ (days)</td><td>770<sup>+3308</sup><sub>-632</sub></td><td>1452<sup>+2819</sup><sub>-1359</sub></td></tr><tr><td>Coherence, Γ</td><td>4.5<sup>+1.9</sup><sub>-2.3</sub></td><td>7.5<sup>+10.1</sup><sub>-4.9</sub></td></tr><tr><td>Periodic timescale, P<sub>GP</sub> (days)</td><td>51.4<sup>+1.3</sup><sub>-24.6</sub></td><td>51.2<sup>+6.5</sup><sub>-24.3</sub></td></tr><tr><td>Additive jitter, s (m s⁻¹)</td><td>0.06<sup>+0.18</sup><sub>-0.05</sub></td><td>0.08<sup>+0.19</sup><sub>-0.07</sub></td></tr><tr><td colspan="3"><strong>Measured parameters</strong></td></tr><tr><td colspan="3"><strong>L 98-59b (TOI-175.03)</strong></td></tr><tr><td>Orbital period, P<sub>b</sub> (days)</td><td>2.2531 ± 0.0004</td><td></td></tr><tr><td>Time of mid-transit, T<sub>0,b</sub> [BJD-2 457 000]</td><td>1366.1708 ± 0.0001</td><td></td></tr><tr><td>Semi-amplitude, K<sub>b</sub> (m s⁻¹)<sup>*</sup></td><td>&lt;1.06</td><td>&lt;1.07</td></tr><tr><td>h<sub>b</sub> = √e<sub>b</sub> cos ω<sub>b</sub></td><td>0.03 ± 0.22</td><td>-0.07 ± 0.24</td></tr><tr><td>k<sub>b</sub> = √e<sub>b</sub> sin ω<sub>b</sub></td><td>-0.01 ± 0.22</td><td>-0.06 ± 0.26</td></tr><tr><td colspan="3"><strong>L 98-59c (TOI-175.01)</strong></td></tr><tr><td>Orbital period, P<sub>c</sub> (days)</td><td>3.6904 ± 0.0003</td><td></td></tr><tr><td>Time of mid-transit, T<sub>0,c</sub> [BJD-2 457 000]</td><td>1367.2752 ± 0.0006</td><td></td></tr><tr><td>Semi-amplitude, K<sub>c</sub> (m s⁻¹)</td><td>2.21 ± 0.28</td><td>2.13 ± 0.32</td></tr><tr><td>h<sub>c</sub> = √e<sub>c</sub> cos ω<sub>c</sub></td><td>0.08 ± 0.15</td><td>-0.05 ± 0.16</td></tr><tr><td>k<sub>c</sub> = √e<sub>c</sub> sin ω<sub>c</sub></td><td>0.01 ± 0.19</td><td>0.04 ± 0.22</td></tr><tr><td colspan="3"><strong>L 98-59d (TOI-175.02)</strong></td></tr><tr><td>Orbital period, P<sub>d</sub> (days)</td><td>7.4512 ± 0.0007</td><td></td></tr><tr><td>Time of mid-transit, T<sub>0,d</sub> [BJD-2 457 000]</td><td>1362.7376 ± 0.0009</td><td></td></tr><tr><td>Semi-amplitude, K<sub>d</sub> (m s⁻¹)</td><td>1.67 ± 0.31</td><td>1.59 ± 0.38</td></tr><tr><td>h<sub>d</sub> = √e<sub>d</sub> cos ω<sub>d</sub></td><td>-0.04 ± 0.20</td><td>0.02 ± 0.18</td></tr><tr><td>k<sub>d</sub> = √e<sub>d</sub> sin ω<sub>d</sub></td><td>0.08 ± 0.25</td><td>-0.01 ± 0.23</td></tr><tr><td colspan="3"><strong>Derived parameters</strong></td></tr><tr><td colspan="3"><strong>L 98-59b (TOI-175.03)</strong></td></tr><tr><td>Semi-major axis, a<sub>b</sub> (AU)</td><td>0.02282 ± 0.00008</td><td></td></tr><tr><td>Equilibrium temperature, T<sub>eq,b</sub> (K)</td><td></td><td></td></tr><tr><td>Bond albedo = 0</td><td>610 ± 13</td><td></td></tr><tr><td>Bond albedo = 0.3</td><td>558 ± 12</td><td></td></tr><tr><td>Planet radius, r<sub>p,b</sub> (R<sub>⊕</sub>)<sup>(*)</sup></td><td>0.80 ± 0.05</td><td></td></tr><tr><td>Planet mass, m<sub>p,b</sub> (M<sub>⊕</sub>)<sup>(*)</sup></td><td>&lt;1.01</td><td>&lt;1.01</td></tr><tr><td>Bulk density, ρ<sub>c</sub> (g cm⁻³)<sup>(*)</sup></td><td>&lt;12.7</td><td>&lt;13.1</td></tr><tr><td>Surface gravity, g<sub>b</sub> (m s⁻²)<sup>(*)</sup></td><td>&lt;16.1</td><td>&lt;17.1</td></tr><tr><td>Escape velocity, v<sub>esc,b</sub> (km s⁻¹)<sup>(*)</sup></td><td>&lt;14.3</td><td>&lt;14.7</td></tr><tr><td>Eccentricity, e<sub>b</sub><sup>(‡)</sup></td><td>&lt;0.12</td><td>&lt;0.15</td></tr><tr><td colspan="3"><strong>L 98-59c (TOI-175.01)</strong></td></tr><tr><td>Semi-major axis, a<sub>c</sub> (AU)</td><td>0.0317 ± 0.0001</td><td></td></tr><tr><td>Equilibrium temperature, T<sub>eq,c</sub> (K)</td><td></td><td></td></tr><tr><td>Bond albedo = 0</td><td>517 ± 11</td><td></td></tr><tr><td>Bond albedo = 0.3</td><td>473 ± 10</td><td></td></tr><tr><td>Planet radius, r<sub>p,c</sub> (R<sub>⊕</sub>)<sup>(‡)</sup></td><td>1.35 ± 0.07</td><td></td></tr><tr><td>Planet mass, m<sub>p,c</sub> (M<sub>⊕</sub>)</td><td>2.42<sup>+0.35</sup><sub>-0.34</sub></td><td>2.36<sup>+0.36</sup><sub>-0.35</sub></td></tr><tr><td>Bulk density, ρ<sub>c</sub> (g cm⁻³)</td><td>5.4<sup>+1.3</sup><sub>-1.0</sub></td><td>5.3<sup>+1.3</sup><sub>-1.1</sub></td></tr><tr><td>Surface gravity, g<sub>d</sub> (m s⁻²)</td><td>13.0<sup>+2.5</sup><sub>-2.1</sub></td><td>12.8<sup>+2.5</sup><sub>-2.2</sub></td></tr><tr><td>Escape velocity, v<sub>esc,c</sub> (km s⁻¹)</td><td>15.0 ± 1.1</td><td>14.8 ± 1.2</td></tr><tr><td>Eccentricity, e<sub>c</sub><sup>(‡)</sup></td><td>&lt; 0.07</td><td>&lt; 0.07</td></tr><tr><td colspan="3"><strong>L 98-59d (TOI-175.02)</strong></td></tr><tr><td>Semi-major axis, a<sub>d</sub> (AU)</td><td>0.0506 ± 0.0002</td><td></td></tr><tr><td>Equilibrium temperature, T<sub>eq,d</sub> (K)</td><td></td><td></td></tr><tr><td>Bond albedo = 0</td><td>409 ± 8</td><td></td></tr><tr><td>Bond albedo = 0.3</td><td>374 ± 8</td><td></td></tr><tr><td>Planet radius, r<sub>p,d</sub> (R<sub>⊕</sub>)<sup>(‡)</sup></td><td>1.57 ± 0.14</td><td></td></tr><tr><td>Planet mass, m<sub>p,d</sub> (M<sub>⊕</sub>)</td><td>2.31<sup>+0.46</sup><sub>-0.45</sub></td><td>2.20<sup>+0.48</sup><sub>-0.47</sub></td></tr><tr><td>Bulk density, ρ<sub>d</sub> (g cm⁻³)</td><td>3.3<sup>+1.3</sup><sub>-0.9</sub></td><td>3.1<sup>+1.3</sup><sub>-0.9</sub></td></tr><tr><td>Surface gravity, g<sub>d</sub> (m s⁻²)</td><td>9.1<sup>+2.7</sup><sub>-2.2</sub></td><td>8.7<sup>+2.7</sup><sub>-2.2</sub></td></tr><tr><td>Escape velocity, v<sub>esc,d</sub> (km s⁻¹)</td><td>13.6 ± 1.5</td><td>13.2<sup>+1.5</sup><sub>-1.6</sub></td></tr><tr><td>Eccentricity, e<sub>d</sub><sup>(‡)</sup></td><td>&lt;0.09</td><td>&lt;0.08</td></tr></tbody></table>

Notes. <sup>(*)</sup>Upper limit given by the 95% confidence interval. <sup>(‡)</sup>Planetary radii from K19. <sup>(‡)</sup>Upper limit given by the 95% confidence interval derived from the joint RV and dynamical stability analyses.
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Time series of each physical RV component in our L 98-59 models. **Left column:** raw RVs (**top panel**), RV activity (**second panel**), three L 98-59 planets (**third, fourth, and fifth panels**), and RV residuals (**bottom panel**). The RV activity model depicted is the mean GP function along with its ±1σ uncertainty in the surrounding shaded region. The L 98-59c and d planet curves are their MAP Keplerian orbital solutions while the L 98-59b curve is its median Keplerian orbit. **Right column:** the GLSP corresponding to each RV component.</figcaption></figure>

$\rho_{p,b} < 12.7$ g cm$^{-3}$ from the 95% confidence interval of the $\rho_{p,b}$ marginalized posterior exceeds the bulk density of a pure iron ball the size of L 98-59b (12.2 g cm$^{-3}$; Zeng & Sasselov 2013). This implies that the true RV semi-amplitude of L 98-59b is likely $\le 1$ m s$^{-1}$ and that the detection of $K_b$ will require much more stringent RV follow-up with an instrument whose performance on L 98-59 is similar to or better than HARPS, such as ESPRESSO (Pepe et al. 2010).

**5. Dynamical stability and eccentricity constraints**

The presence of three planets in a compact configuration around L 98-59 provides a unique opportunity to provide additional

constraints on the planets' orbital eccentricities using stability criteria to limit the range of permissible eccentricities. K19 showed through dynamical simulations (assuming MAP planet mass predictions from Chen & Kipping 2017) that for initially circular orbits the system can be long-lived but as the initial eccentricities were increased to just 0.1, many of their simulated planetary systems became unstable in $\le 20$ 000 yr. Using the planetary mass measurements and upper limits derived in this paper we can use dynamical simulations to constrain each planet's eccentricity given that the system must remain stable for at least the duration of the simulation.

We proceed with deriving the fraction of stable systems as a function of each planet's orbital eccentricity by simulating $10^4$
<figure><img src="image_4.png" /><figcaption>Fig. 3. Phase-folded RVs for each known L 98-59 planet. Each set of RVs has been corrected for stellar activity and the two planets not depicted in its panel. Only the two outermost planets are detected with semi-amplitudes that are inconsistent with 0 m s⁻¹.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Planetary mass and radius plane for small planets. The three L 98-59 planets are highlighted in the planetary mass and radius space along with a set of selected exoplanets, plus the Earth and Venus for comparison. The solid lines depict theoretical mass-radius curves from two component models of fully differentiated planetary interiors with fractional compositions by mass in water (H₂O), silicate (MgSiO₃), and/or iron (Fe) (Zeng & Sasselov 2013). Each model's interior composition is annotated above its respective curve.</figcaption></figure>

realizations of the L 98-59 planetary system and integrating each system forward in time using the WHFast symplectic integrator (Rein & Tamayo 2015) within the open-source REBOUND N-body package (Rein & Liu 2012). In each realization, the stellar mass is drawn from N(0.312, 0.031) M☉ which in turn prescribes each planet's initial semi-major axis when combined with its orbital period that are drawn from their marginalized posterior PDF from Sect. 4. K19 noted that despite having a period ratio of 2.02, the two outer planets are likely just wide of a resonant configuration such that we do not attempt to force the outer planet pair to converge towards a mean motion resonance in our dynamical simulations. Similarly to the orbital periods, each planet's mass and orbital phase at t = 0 are drawn from their marginalized posterior PDFs from Sect. 4. Orbital inclinations are drawn from the approximately Gaussian posterior PDFs reported in Table 2 of K19. The argument of periastron and longitude of the ascending node for each planet are both drawn from U(0, 2π). Lastly, the orbital eccentricities of the planets in each realization are treated as free parameters and are drawn from U(0, 0.3) where the upper eccentricity limit was chosen as any orbit that is initialized with e ≥ 0.3 will undergo an immediate orbit crossing in less than one orbital timescale.

Each simulated planetary system is integrated forward in time until one of the following stopping conditions is reached:

1. Any pair of planets come within one mutual Hill radius:

$$R_{\text{Hill}} = \left( \frac{m_{p,i} + m_{p,i+1}}{3M_s} \right)^{1/3} \frac{a_i + a_{i+1}}{2} \quad (4)$$

2. Any planet that travels beyond the imposed maximum barycentric distance of 0.2 AU (~4aₚ,ₐ).

3. The integration reaches its stopping time of 10⁶ orbits of the outermost planet; ~2 × 10⁴ yr.

Simulations that are halted because of either of the former two stopping criteria are flagged as unstable systems and the corresponding initial eccentricities are ruled out due to instability. The remaining systems that survive until the end of the simulation are deemed stable. We note that due to the short duration of the simulations performed here compared to the expected age of L 98-59 (>1 Gyr; K19), these simulations are not intended to provide a detailed overview of the system's long-term stability but instead are used solely for the purpose of constraining the planetary eccentricities beyond that which can be measured by the RV data alone.

The fraction of stable systems as a function of each planet's initial orbital eccentricity is shown in Fig. 5. The strong stability constraints on each planet's eccentricity are evident as the majority of the three-parameter space exhibits a stability fraction that is consistent with zero. In particular, for the two more massive planets in the system (i.e., L 98-59c and d), planetary systems for which eₐ or eₐ ≥ 0.1 have a stability fraction of <1%. Conversely, a small fraction of planetary systems with 0.1 ≤ eₐ ≤ 0.2 can remain stable as the median mass of L 98-59b is only ~16% of either of the other two planets and therefore its eccentricity has a reduced effect on the overall stability of the system (Barnes & Greenberg 2006).

The large fraction of initial eccentricity values that result in an unstable orbital configuration as seen in Fig. 5, provides constraints on the planets' orbital eccentricity values. If initially we ignore the fractional stability criteria derived from our dynamical simulations, we can derive posterior PDFs of the orbital eccentricities eᵢ from the PDFs of hᵢ and kᵢ obtained from our MCMC analysis (Fig. A.1). The resulting eᵢ posteriors are
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** L 98-59 stability maps. The stability fraction of the L 98-59 three planet system as functions of the planet's initial orbital eccentricities. The 2D maps depict the fraction of stable systems computed from a set of N-body integrations using the planetary masses measured in this study. The 1D histograms depict the number of stable and unstable systems as a function of each planet's eccentricity separately and marginalized over all other dynamical parameters. The histograms are depicted on a logarithmic scale. The annotated numbers indicate each histogram bin's stability fraction in percentages.</figcaption></figure>

depicted in Fig. 6 and represent our measurements of $e_i$ from the RV data alone. Next, we treat the stability fraction as a function of each $e_i$ as an additional prior on $e_i$ and resample the $e_i$ posterior PDFs according to the stability fraction. That is, for each sample from the joint {$e_b, e_c, e_d$} posterior derived from MCMC, the probability that that sample is retained is given by the stability fraction of simulated planetary systems with those eccentricity values ±0.02. In this way high eccentricity values that cannot be ruled out by the MCMC analysis alone are frequently rejected because they often result in an unstable orbital configuration.

The resampled $e_i$ posteriors that account for system stability are compared to the MCMC only results in Fig. 6. The distinct narrowing of each $e_i$ posterior after including the stability criteria indicates that the joint RV + stability data provide the strongest constraints on the orbital eccentricities of the planets in the compact L 98-59 system. From each set of $e_i$ posteriors we derive eccentricity upper limits at 95% confidence. The resulting upper limits from the RV data alone are $e_b < 0.53$, $e_c < 0.19$, and $e_d < 0.31$. By comparison, the inclusion of the stability criteria results in drastically improved upper limits of $e_b < 0.12$, $e_c < 0.07$, and $e_d < 0.09$. These measurements confirm that the planets in the L 98-59 compact planetary system all likely have eccentricities ≲0.1, a result that is consistent with similarly compact systems exhibiting low eccentricities (Hadden & Lithwick 2014; Van Eylen & Albrecht 2015) and the low dispersion in mutual inclinations in the system (Δi ~ 0.4°; K19).

# 6. Discussion and conclusions

In this study we conducted an intensive HARPS RV follow-up campaign of the L 98-59 multi-planet system to characterize the

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Orbital eccentricity marginalized posterior PDFs of each L 98-59 planet. The broader red PDF for each planet corresponds to the MCMC only results representing the eccentricity constraint from the RV data alone. The shallower blue PDFs combine the MCMC results with an additional prior from stability and therefore provides a stronger constraint on each planet's eccentricity. The shaded regions highlight the 95% confidence intervals. The PDFs are plotted on a logarithmic scale for improved visibility.</figcaption></figure>

masses of its three known transiting planets (K19). We measure planet masses of the two outermost planets of $m_{p,c} = 2.42_{-0.34}^{+0.35}$ $M_\oplus$ and $m_{p,d} = 2.31_{-0.45}^{+0.46} M_\oplus$ and derive an upper limit on the mass of the innermost planet of $m_{p,b} < 1.01$ at 95% confidence. The resulting bulk density of the super-Earth L 98-59c is consistent with a bulk terrestrial composition (see Fig. 4) while the bulk density of the outer sub-Neptune L 98-59d is inconsistent with being solely terrestrial and requires either a significant size fraction in water or in an extended gaseous envelope. Although the mass of the inner planet L 98-59b is not robustly measured by our data, its small size and small orbital separation place the innermost planet interior to the photoevaporation valley (Owen & Wu 2013, 2017; Jin et al. 2014; Lopez & Fortney 2014; Chen & Rogers 2016; Lopez & Rice 2018) thus providing supporting evidence for its terrestrial nature as well. Confirmation of the terrestrial nature of L 98-59b will likely require *O*(500) additional RVs with a similar level of precision as our HARPS RVs to measure the semi-amplitude of L 98-59b at $3\sigma$ (Cloutier et al. 2018) given its expected semi-amplitude of $0.32 \text{ m s}^{-1}$ (Chen & Kipping 2017).

With the precise RV planet masses presented in this study, L 98-59c and d add to the growing list of planets to directly contribute to the completion of the TESS level one science
requirement of delivering the masses of fifty planets smaller than $4 R_⊕$. Furthermore, at 1.35 and 1.57 $R_⊕$ respectively, L 98-59c and d are among the smallest TESS planets to have precisely measured masses via RV follow-up observations.

The nearby L 98-59 system of three small planets in a compact configuration within 7.5 days presents an ideal opportunity for comparative atmospheric planetology. To quantify the feasibility of detecting atmospheric signatures from the L 98-59 planets in transmission using an observatory like JWST, we compute the transmission spectroscopy metric from Kempton et al. (2018) using the MAP planet parameters measured in this study and assuming cloud-free atmospheres. We find that $TSM_b > 14.6$, $TSM_c = 23.6 \pm 5.2$, and $TSM_d = 212 \pm 76$. For L 98-59c and d this amounts to ~0.8 – 1.3 and ~6 – 13 times that of GJ 1132b (Dittmann et al. 2017; Bonfils et al. 2018), the previously “best” prospect for the atmospheric characterization of a terrestrial-sized exoplanet from the pre-TESS era (Morley et al. 2017). Hence both L 98-59c and L 98-59d belong to an important set of targets from the TESS mission<fnref n="4" /> that are extremely promising for the atmospheric characterization of hot small exoplanets through transmission spectroscopy.

Similarly, we compute the emission spectroscopy metric from Kempton et al. (2018) assuming that the planet daytime temperatures are equivalent to their equilibrium temperatures assuming an Earth-like albedo (i.e., $A = 0.3$). We find $ESM_b = 2.2 \pm 0.4$, $ESM_c = 3.4 \pm 0.6$, $ESM_d = 1.5 \pm 0.4$ such that the L 98-59 planets are somewhat less favorable for emission spectroscopy characterization compared to GJ 1132b ($ESM_{GJ 1132b} = 3.6 \pm 0.5$) although L 98-59c still represents a viable target for such observations. The disfavourability of the L 98-59 planets compared to GJ 1132b is largely due to the larger stellar radius of L 98-59. Nevertheless, the close proximity of L 98-59 ($d = 10.6$ pc) continues to make the super-Earth L 98-59c a viable candidates for the characterization of a hot terrestrial exoplanet's atmosphere in emission if it can first be demonstrated on a more favourable target such as LHS 3844b (Vanderspek et al. 2019).

Regarding the prospect of the direct detection of the L 98-59 planets in reflected light using near-IR imagers on-board the next generation of extremely large telescopes, the planets' small angular separations ($\theta_b = 0.0021''$, $\theta_c = 0.0030''$, $\theta_d = 0.0048''$) make them difficult to resolve despite having modest planet-star contrasts (0.6, 1.0, 0.5 × 10⁻⁶ respectively). Thus despite its close proximity, the orbital architecture of the known planets around L 98-59 is likely too compact for any of the planets to be directly imaginable within the next decade.

Lastly, we emphasize that L 98-59 is slated to be re-observed by TESS within sectors 5, 8, 9, 10, 11, and 12. The extended baseline beyond the single 27 day field from sector 2 will provide opportunities to improve the orbital ephemerides and radii of the known planets and to continue to search for TTVs. If detected, TTV measurements could enable independent measurements of the planet masses for direct comparison to the RV results presented herein. The extended observational baseline may also enable the detection of additional planets at long orbital periods although the RVs presented herein do not show any significant evidence for such planets.

*Acknowledgements*. R.C. is supported in part by the Natural Sciences and Engineering Research Council of Canada and acknowledges that this work was performed on land traditionally inhabited by the Wendat, the Anishnaabeg, Haudenosaunee, Metis, and the Mississaugas of the New Credit First Nation.

N.A.-D. acknowledges the support of FONDECYT project 3180063. X.B. acknowledges funding from the European Research Council under the ERC Grant Agreement n. 337591-ExTrA. J.S.J acknowledges support by FONDECYT grant 1161218 and partial support from CONICYT project Basal AFB-170002. This work is supported 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 Univ. Grenoble-Alpes. Z.B. acknowledges CONICYT-FONDECYT/Chile Postdoctorado 3180405. We thank the Swiss National Science Foundation (SNSF) and the Geneva University for their continuous support to our exoplanet researches. This work has been in particular carried out in the frame of the National Centre for Competence in Research "PlanetS" supported by SNSF. N.C.S was supported by FCT/MCTES through national funds and by FEDER – Fundo Europeu de Desenvolvimento Regional through COMPETE2020 – Programa Operacional Competitiva e Internacionalização by these grants: UID/FIS/04434/2019; PTDC/FIS-AST/32113/2017 & POCI-01-0145-FEDER-032113; PTDC/FIS-AST/28953/2017 & POCI-01-0145-FEDER-028953. Funding for the TESS mission is provided by NASA's Science Mission directorate. We acknowledge the use of public TESS Alert data from pipelines at the TESS Science Office and at the TESS Science Processing Operations Center. This research has made use of the Exoplanet Follow-up Observation Program website, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center for the production of the SPOC data products. This paper includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST).

## References

Astudillo-Defru, N., Delfosse, X., Bonfils, X., et al. 2017a, A&A, 600, A13

Astudillo-Defru, N., Forveille, T., Bonfils, X., et al. 2017b, A&A, 602, A88

Astudillo-Defru, N., Díaz, R. F., Bonfils, X., et al. 2017c, A&A, 605, L11

Ballard, S. 2019, AJ, 157, 113

Barclay, T., Pepper, J., & Quintana, E. V. 2018, ApJS, 239, 2

Barnes, R., & Greenberg, R. 2006, ApJ, 647, L163

Benedict, G. F., Henry, T. J., Franz, O. G., et al. 2016, AJ, 152, 141

Bonfils, X., Almenara, J. M., Cloutier, R., et al. 2018, A&A, 618, A142

Brahm, R., Espinoza, N., Jordán, A., et al. 2019, AJ, 158, 45

Cañas, C. I., Stefansson, G., Monson, A. J., et al. 2019, ApJ, 877, L29

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

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

Cloutier, R. 2019, AJ, 158, 81

Cloutier, R., Astudillo-Defru, N., Doyon, R., et al. 2017, A&A, 608, A35

Cloutier, R., Doyon, R., Bouchy, F., & Hébrard, G. 2018, AJ, 156, 82

Cloutier, R., Astudillo-Defru, N., Doyon, R., et al. 2019, A&A, 621, A49

Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, 2MASS All Sky Catalog of point sources

Cutri, R. M., Wright, E. L., Conrow, T., et al. 2013, Explanatory Supplement to the AllWISE Data Release Products, Tech. rep

Desort, M., Lagrange, A.-M., Galland, F., Udry, S., & Mayor, M. 2007, A&A, 473, 983

Dittmann, J. A., Irwin, J. M., Charbonneau, D., et al. 2017, Nature, 544, 333

Dragomir, D., Teske, J., Günther, M. N., et al. 2019, ApJ, 875, L7

Dumusque, X., Turner, O., Dorn, C., et al. 2019, A&A, 627, A43

Duncan, D. K., Vaughan, A. H., Wilson, O. C., et al. 1991, ApJS, 76, 383

Espinoza, N., Brahm, R., Henning, T., et al. 2019, A&A, 627, A43

Evans, D. W., Riello, M., De Angeli, F., et al. 2018, A&A, 616, A4

Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ, 790, 146

Ford, E. B. 2006, ApJ, 642, 505

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

Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1

Gandolfi, D., Barragán, O., Livingston, J. H., et al. 2018, A&A, 619, L10

Giles, H. A. C., Collier Cameron, A., & Haywood, R. D. 2017, MNRAS, 472, 1618

Gomes da Silva, J., Santos, N. C., Bonfils, X., et al. 2012, A&A, 541, A9

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

Henden, A. A., Templeton, M., Terrell, D., et al. 2016, VizieR Online Data Catalog: II/336

Huang, C. X., Shporer, A., Dragomir, D., et al. 2018a, AJ, submitted [arXiv:1807.11129]

Huang, C. X., Burt, J., Vanderburg, A., et al. 2018b, ApJ, 868, L39

Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, in Software and Cyberinfrastructure for Astronomy IV, Proc. SPIE, 9913, 99133E

<footnote n="4">⁴ Also of note are LTT 1445Ab (Winters et al. 2019) and GJ 357b (Luque et al. 2019).</footnote>
Jin, S., Mordasini, C., Parmentier, V., et al. 2014, ApJ, 795, 65

Jones, M. I., Brahm, R., Espinoza, N., et al. 2019, A&A, 625, A16

Kempton, E. M. R., Bean, J. L., Louie, D. R., et al. 2018, PASP, 130, 114401

Kipping, D., Nesvorný, D., Hartman, J., et al. 2019, MNRAS, 486, 4980

Kostov, V. B., Schlieder, J. E., Barclay, T., et al. 2019, AJ, 158, 32

Li, J., Tenenbaum, P., Twicken, J. D., et al. 2019, PASP, 131, 024506

Lindegren, L., & Dravins, D. 2003, A&A, 401, 1185

Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2

Lissauer, J. J., Ragozzine, D., Fabrycky, D. C., et al. 2011, ApJS, 197, 8

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

Lopez, E. D., & Rice, K. 2018, MNRAS, 479, 4353

Louie, D. R., Deming, D., Albert, L., et al. 2018, PASP, 130, 044401

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

Luque, R., Pallé, E., Kossakowski, D., et al. 2019, A&A, 628, A39

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

Mann, A. W., Dupuy, T., Kraus, A. L., et al. 2019, ApJ, 871, 63

Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20

Ment, K., Dittmann, J. A., Astudillo-Defru, N., et al. 2019, AJ, 157, 32

Morley, C. V., Kreidberg, L., Rustamkulov, Z., Robinson, T., & Fortney, J. J. 2017, ApJ, 850, 121

Muirhead, P. S., Dressing, C. D., Mann, A. W., et al. 2018, AJ, 155, 180

Nelson, B. E., Ford, E. B., Buchner, J., et al. 2018, AAS J., submitted [arXiv:1806.04683]

Newton, E. R., Irwin, J., Charbonneau, D., et al. 2016, ApJ, 821, 93

Nielsen, L. D., Bouchy, F., Turner, O., et al. 2019, A&A, 623, A100

Owen, J. E., & Wu, Y. 2013, ApJ, 775, 105

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

Pepe, F. A., Cristiani, S., Rebolo Lopez, R., et al. 2010, in Ground-based and Airborne Instrumentation for Astronomy III, Proc. SPIE, 7735, 77350F

Perrakis, K., Ntzoufras, I., & Tsionas, E. G. 2014, Comput. Stat. Data Anal., 77, 54

Queloz, D., Henry, G. W., Sivan, J. P., et al. 2001, A&A, 379, 279

Quinn, S. N., Becker, J. C., Rodriguez, J. E., et al. 2019, AAS J., submitted [arXiv:1901.09092]

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

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

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

Rodriguez, J. E., Quinn, S. N., Huang, C. X., et al. 2019, AJ, 157, 191

Stassun, K. G., Oelkers, R. J., Pepper, J., et al. 2018, AJ, 156, 102

Sullivan, P. W., Winn, J. N., Berta-Thompson, Z. K., et al. 2015, ApJ, 809, 77

Twicken, J. D., Catanzarite, J. H., Clarke, B. D., et al. 2018, PASP, 130, 064502

Vanderburg, A., Plavchan, P., Johnson, J. A., et al. 2016, MNRAS, 459, 3565

Vanderspek, R., Huang, C. X., Vanderburg, A., et al. 2019, ApJ, 871, L24

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

Wang, S., Jones, M., Shporer, A., et al. 2019, AJ, 157, 51

Winters, J. G., Medina, A. A., Irwin, J. M., et al. 2019, AJ, accepted [arXiv:1906.10147]

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

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

¹ Department of Astronomy & Astrophysics, University of Toronto, 50 St. George Street, M5S 3H4, Toronto, ON, Canada
e-mail: cloutier@astro.utoronto.ca

² Centre for Planetary Sciences, Department of Physical & Environmental Sciences, University of Toronto Scarborough, 1265 Military Trail, M1C 1A4, Toronto, ON, Canada

³ Institut de Recherche sur les Exoplanètes, Département de Physique, Université de Montréal, Montréal QC, H3C 3J7, Canada

⁴ Departamento de Matemática y Física Aplicadas, Universidad Católica de la Santísima Concepción, Alonso de Rivera 2850, Concepción, Chile

⁵ Université Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France

⁶ Departamento de Astronomía, Universidad de Chile, Camino El Observatorio 1515, Las Condes, Santiago, Chile

⁷ Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

⁸ Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA

⁹ Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA 02139, USA

¹⁰ Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

¹¹ NASA Ames Research Center, Moffett Field, CA 94035, USA

¹² Observatoire de Genève, Université de Genève, 51 ch. des Maillettes, 1290 Sauverny, Switzerland

¹³ Universidad de Buenos Aires, Facultad de Ciencias Exactas y Naturales. Buenos Aires, Argentina

¹⁴ CONICET – Universidad de Buenos Aires. Instituto de Astronomía y Física del Espacio (IAFE). Buenos Aires, Argentina

¹⁵ European Southern Observatory, Alonso de Córdova 3107, Vitacura, Región Metropolitana, Chile

¹⁶ Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150-762 Porto, Portugal

¹⁷ SETI Institute, Mountain View, CA 94043, USA

¹⁸ Department of Astronomy & Institute for Astrophysical Research, Boston University, 725 Commonwealth Avenue, Boston, MA 02215, USA

¹⁹ Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, Portugal
**Appendix A: Joint and marginalized posterior**
**PDFs from the RV analysis**

The joint and marginalized posteriors PDFs from the RV anal-
ysis presented in Sect. 4 is depicted in Fig. A.1. These results

correspond to the three planet model that includes a GP activity
model trained on the *Hα* time series. Point estimates from each
parameter's marginalized posterior are reported in Table 4.

For the purpose of propagating our posteriors in future
studies, we also include 10⁴ posterior samples in Table A.1.

<figure><img src="image_8.png" /><figcaption>**Fig. A.1.** Marginalized and joint posterior probability density functions of the model parameters from the RV analysis. The adopted RV model includes a trained GP activity model (ln a, ln λ, ln Γ, ln P<sub>GP</sub>, s) plus the star’s systemic velocity (γ<sub>0</sub>) and three Keplerian planet solutions (P<sub>i</sub>, T<sub>0,i</sub>, K<sub>i</sub>, h<sub>i</sub>, k<sub>i</sub>) for i = b, c, d).</figcaption></figure>
<table><caption>Table A.1. 10⁻⁴ posterior samples from the RV analysis presented in Sect. 4.</caption><thead><tr><th>ln a<br>(m s<sup>-1</sup>)</th><th>ln l<br>(days)</th><th>ln Γ<br>(days)</th><th>ln P<sub>GP</sub><br>(days)</th><th>s<br>(m s<sup>-1</sup>)</th><th>γ<sub>0</sub><br>(m s<sup>-1</sup>)</th><th>P<sub>b</sub><br>(days)</th><th>T<sub>0,b</sub><br>(BJD -<br>2 457 000)</th><th>K<sub>b</sub><br>(m s<sup>-1</sup>)</th><th>h<sub>b</sub></th><th>k<sub>b</sub></th><th>P<sub>c</sub><br>(days)</th><th>T<sub>0,c</sub><br>(BJD -<br>2 457 000)</th><th>K<sub>c</sub><br>(m s<sup>-1</sup>)</th><th>h<sub>c</sub></th><th>k<sub>c</sub></th><th>P<sub>d</sub><br>(days)</th><th>T<sub>0,d</sub><br>(BJD -<br>2 457 000)</th><th>K<sub>d</sub><br>(m s<sup>-1</sup>)</th><th>h<sub>d</sub></th><th>k<sub>d</sub></th></tr></thead><tbody><tr><td>1.587511</td><td>3.260021</td><td>-1.111042</td><td>2.741766</td><td>0.000210</td><td>0.080813</td><td>2.252934</td><td>1366.170773</td><td>0.053808</td><td>0.314354</td><td>0.460166</td><td>3.690616</td><td>1367.274264</td><td>2.635292</td><td>-0.145689</td><td>0.178376</td><td>7.452218</td><td>1362.737059</td><td>1.455353</td><td>-0.031272</td><td>0.034064</td></tr><tr><td>1.315478</td><td>6.064731</td><td>1.929885</td><td>3.242569</td><td>0.137362</td><td>0.110952</td><td>2.252928</td><td>1366.170718</td><td>0.139298</td><td>0.154179</td><td>0.226071</td><td>3.690254</td><td>1367.275423</td><td>2.533262</td><td>-0.250329</td><td>0.063469</td><td>7.452395</td><td>1362.737471</td><td>1.159128</td><td>-0.112440</td><td>0.265805</td></tr><tr><td>2.092124</td><td>6.628766</td><td>0.457813</td><td>3.352636</td><td>0.013836</td><td>0.123912</td><td>2.252987</td><td>1366.170525</td><td>0.430868</td><td>0.266768</td><td>0.294228</td><td>3.690091</td><td>1367.277172</td><td>2.016934</td><td>0.105155</td><td>0.032139</td><td>7.451449</td><td>1362.738421</td><td>1.819288</td><td>-0.159901</td><td>0.235015</td></tr><tr><td>1.941346</td><td>5.929036</td><td>1.524704</td><td>3.303236</td><td>0.107918</td><td>0.391627</td><td>2.253225</td><td>1366.170705</td><td>0.564593</td><td>-0.001245</td><td>-0.003462</td><td>3.690240</td><td>1367.274308</td><td>2.555389</td><td>-0.163492</td><td>0.350411</td><td>7.451894</td><td>1362.737876</td><td>1.947260</td><td>0.059668</td><td>0.172321</td></tr><tr><td>1.736321</td><td>4.501731</td><td>1.822428</td><td>3.603358</td><td>0.004258</td><td>0.173040</td><td>2.253674</td><td>1366.170820</td><td>0.459478</td><td>-0.352986</td><td>0.037728</td><td>3.690085</td><td>1367.275824</td><td>2.265025</td><td>-0.037000</td><td>0.441254</td><td>7.452289</td><td>1362.739038</td><td>1.263433</td><td>-0.344930</td><td>0.200037</td></tr></tbody></table>

Notes. The full table is available at the CDS.