Issue 
A&A
Volume 629, September 2019



Article Number  A111  
Number of page(s)  14  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/201935957  
Published online  13 September 2019 
Characterization of the L 9859 multiplanetary system with HARPS
Mass characterization of a hot superEarth, a subNeptune, and a mass upper limit on the third planet^{★,}^{★★}
^{1}
Department of Astronomy & Astrophysics, University of Toronto,
50 St. George Street,
M5S 3H4, Toronto,
ON,
Canada
email: cloutier@astro.utoronto.ca
^{2}
Centre for Planetary Sciences, Department of Physical & Environmental Sciences, University of Toronto Scarborough,
1265 Military Trail,
M1C 1A4,
Toronto,
ON,
Canada
^{3}
Institut de Recherche sur les Exoplanètes, Département de Physique, Université de Montréal,
Montréal QC,
H3C 3J7,
Canada
^{4}
Departamento de Matemática y Física Aplicadas, Universidad Católica de la Santísima Concepción,
Alonso de Rivera 2850, Concepción,
Chile
^{5}
Université Grenoble Alpes, CNRS, IPAG,
38000 Grenoble,
France
^{6}
Departamento de Astronomía, Universidad de Chile,
Camino El Observatorio 1515,
Las Condes, Santiago,
Chile
^{7}
Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology,
Cambridge,
MA 02139,
USA
^{8}
HarvardSmithsonian Center for Astrophysics,
60 Garden Street,
Cambridge,
MA 02138,
USA
^{9}
Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology,
Cambridge,
MA 02139,
USA
^{10}
Department of Astrophysical Sciences, Princeton University,
Princeton,
NJ 08544,
USA
^{11}
NASA Ames Research Center,
Moffett Field,
CA 94035,
USA
^{12}
Observatoire de Genève, Université de Genève,
51 ch. des Maillettes, 1290 Sauverny,
Switzerland
^{13}
Universidad de Buenos Aires,
Facultad de Ciencias Exactas y Naturales. Buenos Aires,
Argentina
^{14}
CONICET – Universidad de Buenos Aires. Instituto de Astronomía y Física del Espacio (IAFE). Buenos Aires,
Argentina
^{15}
European Southern Observatory,
Alonso de Córdova 3107,
Vitacura, Región Metropolitana,
Chile
^{16}
Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas,
4150762 Porto,
Portugal
^{17}
SETI Institute,
Mountain View,
CA 94043,
USA
^{18}
Department of Astronomy & Institute for Astrophysical Research, Boston University,
725 Commonwealth Avenue,
Boston,
MA 02215,
USA
^{19}
Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto,
Portugal
Received:
25
May
2019
Accepted:
13
August
2019
Aims. L 9859 (TIC 307210830, TOI175) 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 longterm 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 9859 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 9859 / planetary systems / stars: lowmass / planets and satellites: terrestrial planets / techniques: radial velocities
Full Tables 2 and A.1 are only available at the CDS via anonymous ftp to cdsarc.ustrasbg.fr (130.79.128.5) or via http://cdsarc.ustrasbg.fr/vizbin/cat/J/A+A/629/A111
© ESO 2019
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 twoyear long primary mission, TESS will observe ≳200 000 targets fromthe TESS Input Catalog (TIC; Stassun et al. 2018) at a two minute cadence as well as many more targets within the 30min 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_{⊕} (π Mensae c; Gandolfi et al. 2018; Huang et al. 2018b, TOI402.01, 02; Dumusque et al. 2019).
The nearby M3 dwarf L 9859 (TIC 307210830^{1}, TOI175^{2}, 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 9859 (L 9859c, TOI175.01: P_{c} = 3.69 days, r_{p,c} = 1.35 R_{⊕}. L 9859d, TOI175.02: P_{d} = 7.45 days, r_{p,d} = 1.57 R_{⊕}. L 9859b, TOI175.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 (Twicken et al. 2018; Li et al. 2019) prior to being published as TESS Data Alerts. Many of the properties of this multiplanetsystem make it of interest for radial velocity (RV) mass characterization (Cloutier et al. 2018), planetary atmosphericcharacterization (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 followup 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 followup study to obtain precise planet masses for as many of the L 9859 planets as possible using HARPS precision RVs. In practice we are only able to recover robust masses for the two outermost planets TOI175.01and 02 but we also report our derived upper limits on the mass of the smallest known planet TOI175.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.
L 9859 stellar parameters.
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 9859 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 subm s^{−1} accuracy.
Throughout the fivemonth observing campaign of L 9859 we elected not to use a simultaneous wavelength calibration (i.e., onsky calibration fiber) to prevent possible contamination of the bluer spectral orders by the calibration lamp. In the ESO programs 198.C0838 and 1102.C0339 (140/164 observations) the exposure time was set to 900 s, resulting in a median signaltonoise ratio (S/N) of 41 per resolution element at 650 nm and a median measurement uncertainty of 2.08 m s^{−1}. In the ESO program 0102.C0525 (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^{−1}. In the ESO program 0102.D0483 (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^{−1}.
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 AstudilloDefru 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 crosscorrelation 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^{−1} wide and centered on the average of the RVs computed by the DRS (−5.661 km s^{−1}). 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.
HARPS spectroscopic time series.
3 Model setup
3.1 Periodogram analysis
To identify strong periodicities in our RV time series we compute the generalized LombScargle 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 Sindex based on the Ca H & K doublet are sensitive to chromospheric activity and may therefore beused to identify periodicities in the RV data arising from chromospheric activity sources such as plages.
The Hα index was computed using the same passbands 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 NaDindex 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 Sindex was calculated following Duncan et al. (1991) and the calibration derived in AstudilloDefru 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 9859d (~ 7.45 days; FAP ~1%) plus significant peaks at the orbital period of L 9859c (~3.69 days; FAP ≲ 0.1%) and centered around ~40 days (FAP ≪ 0.1%). The photometric rotation period of L 9859 remains undetected in the SAP TESS light curve and the star exhibits a negligible rotational broadening (vsini < 1.9 km s^{−1}; K19) indicative of L 9859 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 ± 13 days based on the star’s value of and using the M dwarf magnetic activityrotation relation from AstudilloDefru et al. (2017a). If this periodic signal is indeed due to stellar rotation at P_{rot} ~ 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γ, Sindex, 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 9859b (~ 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 9859d planetary signal (~ 7.45 days), upon the joint modeling of the planets with RV stellar activity in Sect. 4.
Fig. 1
Generalized LombScargle periodograms of HARPS spectroscopic time series. Left column: GLSPs of our HARPS RV time series, its window function, and Hα, Hβ, Hγ, sodium doublet, Sindex, 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 9859 at P_{rot} ~ 80 days and its first harmonic P_{rot}∕2. Right column: false alarm probabilities computed from bootstrapping with replacement. 
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 temporally 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 planetinduced Doppler shifts, to inform our stellar activity model.
Following its successful application to activity modeling in M dwarf planetary systems (e.g., AstudilloDefru et al. 2017c; Cloutier et al. 2017, 2019; Bonfils et al. 2018; Ment et al. 2019), we adopt a semiparametric Gaussian process (GP) regression model of RV stellar activity to simultaneously model activity with the RV planetary signals. Given a parameterization of the temporal covariance structure in our time series, the semiparametric nature of the GP activity treatment is wellsuited 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α data at ~80 days and whose temporal covariance structure is quasiperiodic 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α time series, is (1)
and is parameterized by the following GP hyperparameters: the covariance amplitude a, an exponential decay time scale λ, a coherence parameter Γ, and the strong periodic signal seen in the Hα GLSP (Fig. 1) which we attribute to the L 9859 rotation period P_{rot}.
We sample the posterior probability density function (PDF) of the logarithmic GP hyperparameters by running the Hα time series through a Markov chain Monte Carlo (MCMC) simulation using the emcee ensemble sampler package (ForemanMackey 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 (2)
where y is the vector of N Hα measurements taken at times t = {t_{1}, t_{2}, …, t_{N}} and the N × N covariance matrix K is given by (3)
The inclusion of the Kronecker delta δ_{ij} term adds the Hα measurement uncertainties σ_{Hα} to the diagonal elements of K and includes an additive jitter factor s. Hence the full set parameters sampled during the Hα training phase is {lna, lnλ, lnΓ, lnP_{rot}, s}.
The underlying physical process of active regions in the stellar chromosphere and photosphere rotating in and out of view at P_{rot} is responsible for the observed temporal variations in the Hα 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α to inform our RV model of stellar activity. In particular, the posterior PDFs of the covariance parameters {ln λ, lnΓ, lnP_{rot}} will be used throughout our modeling of the stellar RVs to derive selfconsistent activity and planetary solutions with minimal contamination of the latter by the former as a result of training.
L 9859 model parameter priors.
3.3 Radial velocity model
Following the training of the GP activity model on the Hα 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 9859. The RV GP activity model features the same covariance function as was adopted during training (Eq. (1)) and therefore contains the five hyperparameters {a, λ, Γ, P_{GP}, s} where a and s are unique to the RVs whereas the priors on the remaining hyperparameters λ, Γ and P_{GP} are constrained by their joint posterior PDF from training. Recall that the apparent rotation signal in Hα at P_{rot} ~ 80 days appears to be manifested at its first harmonic (P_{rot}∕2) in the RVs at ~40 days (see Fig. 1). As such, we modify the marginalized posterior on P_{rot} from training by rescaling P_{GP} → P_{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 anygravitational 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 (~ 5.1 min for L 9859b; K19). Furthermore, K19 performed longterm dynamical simulations of the system by considering the maximum aposteriori (MAP) planet mass predictions and their + 1σ values for each of the known L9859 planets. The planet mass predictions were based on the planets’ measured radii and the use of the forecaster tool (Chen & Kipping 2017)^{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., ≲ 0.006) after one million orbits of the outermost planet (i.e., ~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 9859 planets are nearly circular, a result that itself is consistent with other compact multiplanet systems exhibiting low eccentricities (≲ 0.05; Hadden & Lithwick 2014; Van Eylen & Albrecht 2015).
As a back of the envelope calculation, we compare the RV semiamplitudes K under circular orbits to orbits with e = 0.1 for of each the three L 9859 planets. Assuming the MAP predicted planet masses, the difference in K between e = 0 and e = 0.1 is ≲ 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 difference in RV semiamplitudes between that and the circular orbit scenario is ≲ 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 planetinduced 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 Nbody integrations, for which small nonzero 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 9859 systemic velocity γ_{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 midtransit T_{0,i}, RV semiamplitude K_{i}, h_{i} = , and (Ford 2006) where e_{i} and ω_{i} are the planet’s orbital eccentricity and argument of periastron respectively. Our complete RV model therefore contains 21 parameters {ln a, lnλ, lnΓ, lnP_{GP}, s, γ_{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 midtransit 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 lnP_{GP} hyperparameter does not have a welldefined solution that was expected from the training phase. Furthermore, the ln Γ posterior PDF is highly asymmetric and clearly favours larger values than were favoured by the Hα 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 λ, lnΓ, lnP_{GP}} are modified to the following uninformative priors: , , and 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σ level between the two models considered.
We also estimated the Bayesian evidence 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 9859c 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 9859b has a small median RV semiamplitude 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 semiamplitude “detections” in that their measured values are robustly >0 m s^{−1}. These values are K_{c} = 2.21 ± 0.28 m s^{−1} and K_{d} = 1.67 ± 0.31 m s^{−1} and represent 7.9 and 5.4σ 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 nondetection 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 phasefolded RVs are shown in Fig. 3 along with the MAP L 9859c and d Keplerian orbital solutions and the median L 9859b Keplerian orbit. The periodic modulation from L 9859c and d are clearly discernible. Meanwhile the median value of K_{b} corresponds to a ≲1.5σ detection and is not discernible in the phasefolded RVs.
The planetary masses corresponding to the RV semiamplitudes measured with our data are m_{p,b} < 1.01 M_{⊕} (at 95% confidence), M_{⊕}, and M_{⊕}. K19 measure planet–star radius ratios(r_{p,b}∕R_{s} = 0.0234 ± 0.0009, r_{p,c} ∕R_{s} = 0.0396 ± 0.0010, and r_{p,d} ∕R_{s} = 0.0462 ± 0.0029) from which they derive planetary radii of r_{p,b} = 0.80 ± 0.05 R_{⊕}, r_{p,c} = 1.35 ± 0.07 R_{⊕}, and r_{p,d} = 1.57 ± 0.14 R_{⊕}. The measured planetary masses and radii result in constraints on the planets’ bulk densities of ρ_{p,b} < 12.7 g cm^{−3}, g cm^{−3}, g cm^{−3}. Our mass measurements allow the L 9859 planets to be added to the planetary massradius plane in Fig. 4. In doing so we confirm that the bulk composition of the middle planet L 9859c is consistent with terrestrial and is therefore a bonafide superEarth whose interior appears to be dominated by silicate plus an iron core. The bulk composition of the outermost planet L 9859d is inconsistent with a terrestrial composition. This subNeptune planet must instead contain either a significant fraction of its size in water or have retained a significant gaseous atmosphere. Although the nondetection of K_{b} prevents a precise value of its bulk density from being derived, the close proximity of L 9859b 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 ρ_{p,b} may be realized as we note that the upper limit on ρ_{p,b} < 12.7 g cm^{−3} from the 95% confidence interval of the ρ_{p,b} marginalized posterior exceeds the bulk density of a pure iron ball the size of L 9859b (12.2 g cm^{−3}; Zeng & Sasselov 2013). This implies that the true RV semiamplitude of L 9859b is likely ≲ 1 m s^{−1} and that the detection of K_{b} will require much more stringent RV followup with an instrument whose performance on L 9859 is similar to or better than HARPS, such as ESPRESSO (Pepe et al. 2010).
Model parameters for the L 9859 threeplanet system.
Fig. 2
Time series of each physical RV component in our L 9859 models. Left column: raw RVs (top panel), RV activity (second panel), three L 9859 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 9859c and d planet curves are their MAP Keplerian orbital solutions while the L 9859b curve is its median Keplerian orbit. Right column: the GLSP corresponding to each RV component. 
Fig. 3
Phasefolded RVs for each known L 9859 planet. Each set of RVs has been correctedfor stellar activity and the two planets not depicted in its panel. Only the two outermost planets are detected with semiamplitudes that are inconsistent with 0 m s^{−1}. 
Fig. 4
Planetary mass and radius plane for small planets. The three L 9859 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 massradius curves from two component models of fully differentiated planetary interiors with fractional compositions by mass in water (H_{2} O), silicate (MgSiO_{3}), and/or iron (Fe) (Zeng & Sasselov 2013). Each model’s interior composition is annotated above its respective curve. 
5 Dynamical stability and eccentricity constraints
The presence of three planets in a compact configuration around L 9859 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 longlived but as the initial eccentricities were increased to just 0.1, many of their simulated planetary systems became unstable in ≲ 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} realizations of the L 9859 planetary system and integrating each system forward in time using the WHFast symplectic integrator (Rein & Tamayo 2015) within the opensource REBOUND Nbody package (Rein & Liu 2012). In each realization, the stellar mass is draw from M_{⊙} which in turn prescribes each planet’s initial semimajor 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 . Lastly, the orbital eccentricities of the planets in each realization are treated as free parameters and are drawn from 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:
 2.
Any planet that travels beyond the imposed maximum barycentric distance of 0.2 AU (~4a_{p,d}).
 3.
The integration reaches its stopping time of 10^{6} orbits of the outermost planet; ~2 × 10^{4} 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 9859 (> 1 Gyr; K19), these simulations are not intended to provide a detailed overview of the system’s longterm 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 threeparameter space exhibits a stability fraction that is consistent with zero. In particular, for the two more massive planets in the system (i.e., L 9859c and d), planetary systems for which e_{c} or e_{d} ≳ 0.1 have a stability fraction of < 1%. Conversely, a small fraction of planetary systems with 0.1 ≤ e_{b} ≲ 0.2 can remain stable as the median mass of L 9859b 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_{i} from the PDFs of h_{i} and k_{i} obtained from our MCMC analysis (Fig. A.1). The resulting e_{i} posteriorsare 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 9859 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 9859 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).
Fig. 5
L 9859 stability maps. The stability fraction of the L 9859 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 Nbody 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. 
Fig. 6
Orbital eccentricity marginalized posterior PDFs of each L 9859 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. 
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Γ, lnP_{GP}, s}) plus the star’s systemic velocity ({γ_{0}}) and three Keplerian planet solutions ({P_{i}, T_{0,i}, K_{i}, h_{i}, k_{i}} for i = b, c, d). 
6 Discussion and conclusions
In this study we conducted an intensive HARPS RV followup campaign of the L 9859 multiplanet system to characterize the masses of its three known transiting planets (K19). We measure planet masses of the two outermost planets of M_{⊕} and M_{⊕} 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 superEarth L 9859c is consistent with a bulk terrestrial composition (see Fig. 4) while the bulk density of the outer subNeptune L 9859d 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 9859b 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 9859b will likely require additional RVs with a similar level of precision as our HARPS RVs to measure the semiamplitude of L 9859b at 3σ (Cloutier et al. 2018) given its expected semiamplitude of 0.32 m s^{−1} (Chen & Kipping 2017).
With the precise RV planet masses presented in this study, L 9859c 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 9859c and d are among the smallest TESS planets to have precisely measured masses via RV followup observations.
The nearby L 9859 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 9859 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 cloudfree atmospheres. We find that TSM_{b} > 14.6, TSM_{c} = 23.6 ± 5.2, and TSM_{d} = 212 ± 76. For L 9859c 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 terrestrialsized exoplanet from the preTESS era (Morley et al. 2017). Hence both L 9859c and L 9859d belong to an important set of targets from the TESS mission^{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 dayside temperatures are equivalent to their equilibrium temperatures assuming an Earthlike albedo (i.e., A = 0.3). We find ESM_{b} = 2.2 ± 0.4, ESM_{c} = 3.4 ± 0.6, ESM_{d} = 1.5 ± 0.4 such that the L 9859 planets are somewhat less favorable for emission spectroscopy characterization compared to GJ 1132b (ESM_{GJ 1132b} = 3.6 ± 0.5) although L 9859c still represents a viable target for such observations. The disfavourability of the L 9859 planets compared to GJ 1132b is largely due to the larger stellar radius of L 9859. Nevertheless, the close proximity of L 9859 (d = 10.6 pc) continues tomake the superEarth L 9859c 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 9859 planets in reflected light using nearIR imagers onboard the next generation of extremely large telescopes, the planets’ small angular separations (θ_{b} = 0.0021″, θ_{c} = 0.0030″, θ_{d} = 0.0048″) make them difficult to resolve despite having modest planetstar contrasts (0.6, 1.0, 0.5 × 10^{−6} respectively). Thus despite its close proximity, the orbital architecture of the known planets around L 9859 is likely too compact for any of the planets to be directly imagable within the next decade.
Lastly, we emphasize that L 9859 is slated to be reobserved 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. 337591ExTrA. J.S.J acknowledges support by FONDECYT grant 1161218 and partial support from CONICYT project Basal AFB170002. This work is supported by the French National Research Agency in the framework of the Investissements d’Avenir program (ANR15IDEX02), through the funding of the “Origin of Life” project of the Univ. GrenobleAlpes. Z.B. acknowledges CONICYTFONDECYT/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 Competitividade e Internacionalização by these grants: UID/FIS/04434/2019; PTDC/FISAST/32113/2017 & POCI010145FEDER032113; PTDC/FISAST/28953/2017 & POCI010145FEDER028953. 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 Followup 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 HighEnd 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).
Appendix A Joint and marginalized posterior PDFs from the RV analysis
The joint and marginalized posteriors PDFs from the RV analysis 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^{4} posterior samples in Table A.1.
References
 AstudilloDefru, N., Delfosse, X., Bonfils, X., et al. 2017a, A&A, 600, A13 [Google Scholar]
 AstudilloDefru, N., Forveille, T., Bonfils, X., et al. 2017b, A&A, 602, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AstudilloDefru, N., Díaz, R. F., Bonfils, X., et al. 2017c, A&A, 605, L11 [Google Scholar]
 Ballard, S. 2019, AJ, 157, 113 [NASA ADS] [CrossRef] [Google Scholar]
 Barclay, T., Pepper, J., & Quintana, E. V. 2018, ApJS, 239, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Barnes, R., & Greenberg, R. 2006, ApJ, 647, L163 [Google Scholar]
 Benedict, G. F., Henry, T. J., Franz, O. G., et al. 2016, AJ, 152, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Bonfils, X., Almenara, J. M., Cloutier, R., et al. 2018, A&A, 618, A142 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brahm, R., Espinoza, N., Jordán, A., et al. 2019, AJ, 158, 45 [NASA ADS] [CrossRef] [Google Scholar]
 Cañas, C. I., Stefansson, G., Monson, A. J., et al. 2019, ApJ, 877, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, H., & Rogers, L. A. 2016, ApJ, 831, 180 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, J., & Kipping, D. 2017, ApJ, 834, 17 [Google Scholar]
 Cloutier, R. 2019, AJ, 158, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Cloutier, R., AstudilloDefru, N., Doyon, R., et al. 2017, A&A, 608, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cloutier, R., Doyon, R., Bouchy, F., & Hébrard, G. 2018, AJ, 156, 82 [NASA ADS] [CrossRef] [Google Scholar]
 Cloutier, R., AstudilloDefru, N., Doyon, R., et al. 2019, A&A, 621, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, 2MASS All Sky Catalog of point sources [Google Scholar]
 Cutri, R. M., Wright, E. L., Conrow, T., et al. 2013, Explanatory Supplement to the AllWISE Data Release Products, Tech. rep [Google Scholar]
 Desort, M., Lagrange, A.M., Galland, F., Udry, S., & Mayor, M. 2007, A&A, 473, 983 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dittmann, J. A., Irwin, J. M., Charbonneau, D., et al. 2017, Nature, 544, 333 [NASA ADS] [CrossRef] [Google Scholar]
 Dragomir, D., Teske, J., Günther, M. N., et al. 2019, ApJ, 875, L7 [NASA ADS] [CrossRef] [Google Scholar]
 Dumusque, X., Turner, O., Dorn, C., et al. 2019, A&A, 627, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Duncan, D. K., Vaughan, A. H., Wilson, O. C., et al. 1991, ApJS, 76, 383 [NASA ADS] [CrossRef] [Google Scholar]
 Espinoza, N., Brahm, R., Henning, T., et al. 2019, A&A, 627, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Evans, D. W., Riello, M., De Angeli, F., et al. 2018, A&A, 616, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ, 790, 146 [NASA ADS] [CrossRef] [Google Scholar]
 Ford, E. B. 2006, ApJ, 642, 505 [NASA ADS] [CrossRef] [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [CrossRef] [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2018, A&A, 616, A1 [Google Scholar]
 Gandolfi, D., Barragán, O., Livingston, J. H., et al. 2018, A&A, 619, L10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Giles, H. A. C., Collier Cameron, A., & Haywood, R. D. 2017, MNRAS, 472, 1618 [Google Scholar]
 Gomes da Silva, J., Santos, N. C., Bonfils, X., et al. 2012, A&A, 541, A9 [Google Scholar]
 Hadden, S., & Lithwick, Y. 2014, ApJ, 787, 80 [Google Scholar]
 Henden, A. A., Templeton, M., Terrell, D., et al. 2016, VizieR Online Data Catalog: II/336 [Google Scholar]
 Huang, C. X., Shporer, A., Dragomir, D., et al. 2018a, AJ, submitted [arXiv:1807.11129] [Google Scholar]
 Huang, C. X., Burt, J., Vanderburg, A., et al. 2018b, ApJ, 868, L39 [NASA ADS] [CrossRef] [Google Scholar]
 Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, in Software and Cyberinfrastructure for Astronomy IV, Proc. SPIE, 9913, 99133E [Google Scholar]
 Jin, S., Mordasini, C., Parmentier, V., et al. 2014, ApJ, 795, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Jones, M. I., Brahm, R., Espinoza, N., et al. 2019, A&A, 625, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kempton, E. M. R., Bean, J. L., Louie, D. R., et al. 2018, PASP, 130, 114401 [NASA ADS] [CrossRef] [Google Scholar]
 Kipping, D., Nesvorný, D., Hartman, J., et al. 2019, MNRAS, 486, 4980 [NASA ADS] [Google Scholar]
 Kostov, V. B., Schlieder, J. E., Barclay, T., et al. 2019, AJ, 158, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Li, J., Tenenbaum, P., Twicken, J. D., et al. 2019, PASP, 131, 024506 [NASA ADS] [CrossRef] [Google Scholar]
 Lindegren, L., & Dravins, D. 2003, A&A, 401, 1185 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lindegren, L., Hernández, J., Bombrun, A., et al. 2018, A&A, 616, A2 [Google Scholar]
 Lissauer, J. J., Ragozzine, D., Fabrycky, D. C., et al. 2011, ApJS, 197, 8 [NASA ADS] [CrossRef] [Google Scholar]
 Lopez, E. D., & Fortney, J. J. 2014, ApJ, 792, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Lopez, E. D., & Rice, K. 2018, MNRAS, 479, 4353 [Google Scholar]
 Louie, D. R., Deming, D., Albert, L., et al. 2018, PASP, 130, 044401 [NASA ADS] [CrossRef] [Google Scholar]
 Lovis, C., & Pepe, F. 2007, A&A, 468, 1115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Luque, R., Pallé, E., Kossakowski, D., et al. 2019, A&A, 628, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mann, A. W., Feiden, G. A., Gaidos, E., Boyajian, T., & von Braun, K. 2015, ApJ, 804, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Mann, A. W., Dupuy, T., Kraus, A. L., et al. 2019, ApJ, 871, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20 [NASA ADS] [Google Scholar]
 Ment, K., Dittmann, J. A., AstudilloDefru, N., et al. 2019, AJ, 157, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Morley, C. V., Kreidberg, L., Rustamkulov, Z., Robinson, T., & Fortney, J. J. 2017, ApJ, 850, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Muirhead, P. S., Dressing, C. D., Mann, A. W., et al. 2018, AJ, 155, 180 [NASA ADS] [CrossRef] [Google Scholar]
 Nelson, B. E., Ford, E. B., Buchner, J., et al. 2018, AAS J., submitted [arXiv:1806.04683] [Google Scholar]
 Newton, E. R., Irwin, J., Charbonneau, D., et al. 2016, ApJ, 821, 93 [NASA ADS] [CrossRef] [Google Scholar]
 Nielsen, L. D., Bouchy, F., Turner, O., et al. 2019, A&A, 623, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Owen, J. E., & Wu, Y. 2013, ApJ, 775, 105 [Google Scholar]
 Owen, J. E., & Wu, Y. 2017, ApJ, 847, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Pepe, F. A., Cristiani, S., Rebolo Lopez, R., et al. 2010, in Groundbased and Airborne Instrumentation for Astronomy III, Proc. SPIE, 7735, 77350F [CrossRef] [Google Scholar]
 Perrakis, K., Ntzoufras, I., & Tsionas, E. G. 2014, Comput. Stat. Data Anal., 77, 54 [Google Scholar]
 Queloz, D., Henry, G. W., Sivan, J. P., et al. 2001, A&A, 379, 279 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Quinn, S. N., Becker, J. C., Rodriguez, J. E., et al. 2019, AAS J., submitted [arXiv:1901.09092] [Google Scholar]
 Rein, H., & Liu, S. F. 2012, A&A, 537, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rein, H., & Tamayo, D. 2015, MNRAS, 452, 376 [NASA ADS] [CrossRef] [Google Scholar]
 Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, J. Astron. Telesc. Instrum. Syst., 1, 014003 [NASA ADS] [CrossRef] [Google Scholar]
 Rodriguez, J. E., Quinn, S. N., Huang, C. X., et al. 2019, AJ, 157, 191 [NASA ADS] [CrossRef] [Google Scholar]
 Stassun, K. G., Oelkers, R. J., Pepper, J., et al. 2018, AJ, 156, 102 [NASA ADS] [CrossRef] [Google Scholar]
 Sullivan, P. W., Winn, J. N., BertaThompson, Z. K., et al. 2015, ApJ, 809, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Twicken, J. D., Catanzarite, J. H., Clarke, B. D., et al. 2018, PASP, 130, 064502 [NASA ADS] [CrossRef] [Google Scholar]
 Vanderburg, A., Plavchan, P., Johnson, J. A., et al. 2016, MNRAS, 459, 3565 [NASA ADS] [CrossRef] [Google Scholar]
 Vanderspek, R., Huang, C. X., Vanderburg, A., et al. 2019, ApJ, 871, L24 [NASA ADS] [CrossRef] [Google Scholar]
 Van Eylen, V., & Albrecht, S. 2015, ApJ, 808, 126 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, S., Jones, M., Shporer, A., et al. 2019, AJ, 157, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Winters, J. G., Medina, A. A., Irwin, J. M., et al. 2019, AJ, accepted [arXiv:1906.10147] [Google Scholar]
 Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zeng, L., & Sasselov, D. 2013, PASP, 125, 227 [NASA ADS] [CrossRef] [Google Scholar]
Also of note are LTT 1445Ab (Winters et al. 2019) and GJ 357b (Luque et al. 2019).
All Tables
All Figures
Fig. 1
Generalized LombScargle periodograms of HARPS spectroscopic time series. Left column: GLSPs of our HARPS RV time series, its window function, and Hα, Hβ, Hγ, sodium doublet, Sindex, 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 9859 at P_{rot} ~ 80 days and its first harmonic P_{rot}∕2. Right column: false alarm probabilities computed from bootstrapping with replacement. 

In the text 
Fig. 2
Time series of each physical RV component in our L 9859 models. Left column: raw RVs (top panel), RV activity (second panel), three L 9859 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 9859c and d planet curves are their MAP Keplerian orbital solutions while the L 9859b curve is its median Keplerian orbit. Right column: the GLSP corresponding to each RV component. 

In the text 
Fig. 3
Phasefolded RVs for each known L 9859 planet. Each set of RVs has been correctedfor stellar activity and the two planets not depicted in its panel. Only the two outermost planets are detected with semiamplitudes that are inconsistent with 0 m s^{−1}. 

In the text 
Fig. 4
Planetary mass and radius plane for small planets. The three L 9859 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 massradius curves from two component models of fully differentiated planetary interiors with fractional compositions by mass in water (H_{2} O), silicate (MgSiO_{3}), and/or iron (Fe) (Zeng & Sasselov 2013). Each model’s interior composition is annotated above its respective curve. 

In the text 
Fig. 5
L 9859 stability maps. The stability fraction of the L 9859 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 Nbody 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. 

In the text 
Fig. 6
Orbital eccentricity marginalized posterior PDFs of each L 9859 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. 

In the text 
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Γ, lnP_{GP}, s}) plus the star’s systemic velocity ({γ_{0}}) and three Keplerian planet solutions ({P_{i}, T_{0,i}, K_{i}, h_{i}, k_{i}} for i = b, c, d). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.