The L 98-59 System: Three Transiting, Terrestrial-Sized Planets Orbiting a Nearby M-dwarf

TL;DR

This paper reports the discovery of three terrestrial-sized planets orbiting the nearby M dwarf L 98-59 using TESS data, highlighting their potential for follow-up studies and system dynamical analysis.

Contribution

First detailed characterization of a multi-planet system around a bright M dwarf with planets suitable for follow-up observations.

Findings

Three planets with radii 0.8-1.6 Earth radii

Orbital periods from 2.25 to 7.45 days

System near a 2:1 resonance

Abstract

We report the Transiting Exoplanet Survey Satellite (TESS) discovery of three terrestrial-sized planets transiting L 98-59 (TOI-175, TIC 307210830) -- a bright M dwarf at a distance of 10.6 pc. Using the Gaia-measured distance and broad-band photometry we find that the host star is an M3 dwarf. Combined with the TESS transits from three sectors, the corresponding stellar parameters yield planet radii ranging from 0.8REarth to 1.6REarth. All three planets have short orbital periods, ranging from 2.25 to 7.45 days with the outer pair just wide of a 2:1 period resonance. Diagnostic tests produced by the TESS Data Validation Report and the vetting package DAVE rule out common false positive sources. These analyses, along with dedicated follow-up and the multiplicity of the system, lend confidence that the observed signals are caused by planets transiting L 98-59 and are not associated with other sources in the field. The L 98-59 system is interesting for a number of reasons: the host star is bright (V = 11.7 mag, K = 7.1 mag) and the planets are prime targets for further follow-up observations including precision radial-velocity mass measurements and future transit spectroscopy with the James Webb Space Telescope; the near resonant configuration makes the system a laboratory to study planetary system dynamical evolution; and three planets of relatively similar size in the same system present an opportunity to study terrestrial planets where other variables (age, metallicity, etc.) can be held constant. L 98-59 will be observed in 4 more TESS sectors, which will provide a wealth of information on the three currently known planets and have the potential to reveal additional planets in the system.

Figures

Tables

Table 1. Table 1: Stellar Parameters

Parameter	Value	Notes
Identifying Information
Name	L 98-59
TIC ID	307210830
TOI ID	175
$\alpha$ R.A. (hh:mm:ss)	08:18:07.62	Gaia DR2
$\delta$ Dec. (dd:mm:ss)	-68:18:46.80	Gaia DR2
${\mu }_{\alpha }$ (mas yr-1)	$94.767\pm 0.054$	Gaia DR2
${\mu }_{\delta }$ (mas yr-1)	$-340.470\pm 0.052$	Gaia DR2
Distance (pc)	$10.623\pm 0.003$	Gaia DR2
Photometric Properties
$B$ (mag) ……….	$13.289\pm 0.027$	APASS DR9
$V$ (mag) ……….	$11.685\pm 0.017$	APASS DR9
$G$ (mag) ……….	$10.598\pm 0.001$	Gaia DR2
$g^{\prime }$ (mag) ……….	$12.453\pm 0.019$	APASS DR9
$r^{\prime }$ (mag) ……….	$11.065\pm 0.044$	APASS DR9
$T$ (mag) ……….	9.393	TIC
$J$ (mag) ……….	$7.933\pm 0.027$	2MASS
$H$ (mag) ……….	$7.359\pm 0.049$	2MASS
$K_s$ (mag) ………	$7.101\pm 0.018$	2MASS
$W1$ (mag) ………	$6.935\pm 0.062$	ALLWISE
$W2$ (mag) ………	$6.767\pm 0.021$	ALLWISE
$W3$ (mag) ………	$6.703\pm 0.016$	ALLWISE
$W4$ (mag) ………	$6.578\pm 0.047$	ALLWISE
Stellar Properties
Spectral Type ……….	M3V $\pm$ 1	This Work
$T_{\mathrm{eff}}$ (K) ……….	$3367\pm 150$	This Work
$[$Fe/H$]$ ……….	$-0.5\pm 0.5$	This Work
$M_{star}$ ($M_{\odot }$) ……….	$0.313\pm 0.014$	This Work
$R_{star}$ ($R_{\odot }$) ……….	$0.312\pm 0.014$	This Work
$L_{star}\left(L_{\odot }\right)$ ……….	$0.0113\pm 0.0006$	This Work

Table 2. Table 2: Planet Parameters

Parameter	-1$\sigma$	Median	+1$\sigma$
Model Parameters
Star
$\ln \rho$ [g cm-3]	2.57	2.76	2.91
Limb darkening $q_1$	0.41	0.65	0.83
Limb darkening $q_2$	-0.34	-0.11	0.28
L 98-59 b
$T_0$ (BJD - 2457000)	1366.1694	1366.1701	1366.1707
$\ln Period$ [days]	0.812318	0.812326	0.812334
Impact parameter	0.13	0.36	0.55
$\ln R_p/R_{\ast }$	-3.79	-3.75	-3.72
eccentricity	0.03	0.10	0.27
$\omega$ [radians]	-2.2	0.3	2.4
L 98-59 c
$T_0$ (BJD - 2457000)	1367.2752	1367.2755	1367.2759
$\ln Period$ [days]	1.305791	1.305795	1.305798
Impact parameter	0.09	0.29	0.49
$\ln R_p/R_{\ast }$	-3.25	-3.23	-3.20
eccentricity	0.02	0.09	0.25
$\omega$ [radians]	-2.5	-0.4	2.2
L 98-59 d
$T_0$ (BJD - 2457000)	1362.7367	1362.7375	1362.7382
$\ln Period$ [days]	2.008323	2.008329	2.008334
Impact parameter	0.75	0.89	0.93
$\ln R_p/R_{\ast }$	-3.16	-3.07	-3.01
eccentricity	0.04	0.20	0.52
$\omega$ [radians]	-1.9	0.7	2.3
Derived Parameters
L 98-59 b
Period [days]	2.25312	2.25314	2.25316
$R_p/R_{\ast }$	0.0226	0.0234	0.0243
Radius $\left[R_{\oplus }\right]$	0.75	0.80	0.85
Insolation	19.5	23.9	29.2
$a/R_{\ast }$	15.2	16.2	17.0
$a$ [AU]	0.0216	0.0233	0.0250
Inclination (deg)	88.0	88.7	89.5
Duration (hours)	0.89	1.02	1.19
L 98-59 c
Period [days]	3.690607	3.690621	3.690634
$R_p/R_{\ast }$	0.0388	0.0396	0.0407
Radius $\left[R_{\oplus }\right]$	1.28	1.35	1.43
Insolation	10.1	12.4	15.2
$a/R_{\ast }$	21.1	22.5	, 23.6
$a$ [AU]	0.0300	0.0324	0.0347
Inclination (deg)	88.8	89.3	89.7
Duration (hours)	1.07	1.24	1.36
L 98-59 d
Period [days]	7.45081	7.45086	7.45090
$R_p/R_{\ast }$	0.0426	0.0462	0.0492
Radius $\left[R_{\oplus }\right]$	1.43	1.57	1.71
Insolation	3.96	4.85	5.93
$a/R_{\ast }$	36.2	37.4	38.5
$a$ [AU]	0.048	0.052	0.056
Inclination (deg)	88.0	88.5	88.7
Duration (hours)	0.74	0.91	1.68

Table 3. Table 3: Photocenter Analysis of the three planet candidates

Parameter	Row [pixels]	Column [pixels]
L 98-59 b
Out of Transit Image Centroid	$338.83\pm 0.02$	$664.01\pm 0.01$
Difference Image Centroid	$338.34\pm 1.39$	$664.04\pm 0.55$
Offset	$-0.49\pm 1.39$	$0.03\pm 0.55$
Offset/$\sigma$	0.35	0.05
L 98-59 c
Out of Transit Image Centroid	$338.83\pm 0.02$	$664.01\pm 0.01$
Difference Image Centroid	$338.76\pm 0.31$	$664.11\pm 0.24$
Offset	$-0.07\pm 0.31$	$0.1\pm 0.24$
Offset/$\sigma$	0.23	0.42
L 98-59 d
Out of Transit Image Centroid	$338.93\pm 0.02$	$664.02\pm 0.02$
Difference Image Centroid	$339.27\pm 0.21$	$664.06\pm 0.34$
Offset	$0.34\pm 0.21$	$0.04\pm 0.34$
Offset/$\sigma$	1.62	0.12

Table 4. Table 4: Observation Log

Planet	Date	Telescope${\dagger}$${\dagger}$Telescopes: LCO-SSO-1: Las Cumbres Observatory - Siding Spring (1.0 m) LCO-CTIO-1: Las Cumbres Observatory - Cerro Tololo Interamerican Observatory (1.0 m) LCO-CTIO-0.4: Las Cumbres Observatory - Cerro Tololo Interamerican Observatory (0.4 m) LCO-SAAO-1: Las Cumbres Observatory - South African Astronomical Observatory (1.0 m) LCO-SAAO-0.4: Las Cumbres Observatory - South African Astronomical Observatory (0.4 m) SSO-Europa: SPECULOOS South Observatory - Europa (1.0 m) SSO-Io: SPECULOOS South Observatory - Io (1.0 m) SSO-T17: Siding Spring Observatory - T17 (0.4 m) MEarth: MEarth-South telescope array (0.4 m $\times$ 5 telescopes)	Filter	ExpT	Exp	Dur.	Transit	Aperture	FWHM
	(UTC)			(sec)	(N)	(min)	coverage	(arcsec)	(arcsec)
b	2018-10-19	LCO-SSO‡	r′	120	41	147	Ingr.+71$\%$	3.89	2.21
	2018-11-11	SSO-Europa	r′	15	764	315	Full	5.25	4.53
	2018-11-11	LCO-CTIO-1	i′	25	120	154	Full	5.83	3.49
	2018-11-18	LCO-SAAO-0.4‡	i′	70	125	175	Full	9.14	7.03
	2018-11-20	LCO-CTIO-1	r′	30	146	178	Full	3.89	2.20
	2018-11-29	LCO-CTIO-1	r′	30	186	223	Full	5.83	3.56
	2018-12-07	SSO-Io	r′	15	998	415	Full	5.60	4.06
	2019-01-26	LCO-SAAO-1‡	r′	12	88	234	Full	3.89	3.08
c	2018-10-16	LCO-CTIO-1	i′	20	63	70	Ingr.+30$\%$	4.27	2.01
	2018-10-22‡	SSO-T17	clear	30	122	86	Ingr.+77$\%$	7.10	2.40
	2018-11-11	SSO-Europa	r′	15	764	315	Full	5.25	4.53
	2018-11-22	MEarth	RG715	45	1682	380	Full	20.16	8.00
	2018-12-25	LCO-SSO-1	i′	22	108	113	Full	5.05	1.89
	2019-01-20	LCO-CTIO-1	g′	100	85	197	Full	6.22	2.75
	2019-01-20	LCO-CTIO-1	zs	30	170	197	Full	9.36	4.43
d	2018-11-07	LCO-CTIO-0.4	i′	70	119	170	Full	6.85	4.48
	2018-11-22	LCO-CTIO-1‡	i′	25	108	113	Full	6.22	5.05
	2018-11-22	MEarth	RG715	45	1682	380	Full	20.16	8.00
	2019-01-13	LCO-CTIO-0.4	i′	14	143	73	Full	5.14	2.32
	2019-01-20	LCO-CTIO-1	r′	30	132	151	Full	6.22	2.89
	2019-01-28	LCO-CTIO-1	g′	50	153	223	Full	7.78	2.35

Table 5. Table 5: Follow-up R p / R ∗ subscript 𝑅 𝑝 subscript 𝑅 {R_{p}/R_{*}} minus TESS R p / R ∗ subscript 𝑅 𝑝 subscript 𝑅 {R_{p}/R_{*}} as a function of wavelength

TOI	Date Obs	Filter (Observatory)	$R_p/R_{\ast }$ minus TESS $R_p/R_{\ast }$
L 98-59 b	UT 2018-11-20	Sloan r’ (LCO)	$-{0.0008}_{-0.0085}^{+0.0054}$
L 98-59 b	UT 2018-11-29	Sloan r’ (LCO)	${0.0018}_{-0.0045}^{+0.0039}$
L 98-59 c	UT 2019-01-20	Sloan g’ (LCO)	$-0.0033\pm 0.003$
L 98-59 c	UT 2018-11-22	RG715 (MEarth)	$-0.0007\pm 0.002$
L 98-59 c	UT 2019-01-20	Sloan z’ (LCO)	$-0.0056\pm 0.0028$
L 98-59 d	UT 2019-01-28	Sloan g’ (LCO)	${0.0072}_{-0.0064}^{+0.0078}$
L 98-59 d	UT 2019-01-20	Sloan r’ (LCO)	$-{0.003}_{-0.012}^{+0.01}$
L 98-59 d	UT 2018-11-07	Sloan i’ (LCO)	${0.0072}_{-0.0072}^{+0.0088}$
L 98-59 d	UT 2018-11-22	RG715 (MEarth)	${0.0104}_{-0.0043}^{+0.0044}$

Equations

$$S(\omega )=\sqrt{\frac{2}{\pi }}\frac{S_0{\omega }_0^4}{({\omega }^2-{\omega }_0^2{)}^2+{\omega }^2{\omega }_0^2/Q^2},$$

Full text

planets and satellites: detection, techniques: photometric, targets: TIC 307210830, TOI-175

††facilities: AAVSO, LCO, TESS††software: AstroImageJ (Collins et al., 2017), astropy (Astropy Collaboration et al., 2013, 2018), celerite (Foreman-Mackey et al., 2017; Foreman-Mackey, 2018), emcee (Foreman-Mackey et al., 2013), exoplanet (Foreman-Mackey, 2018), DAVE (Kostov et al., 2019), forecaster (Chen & Kipping, 2017), IPython (Perez & Granger, 2007), Jupyter (Kluyver et al., 2016), Lightkurve (Lightkurve Collaboration et al., 2018), Matplotlib (Hunter, 2007), NumPy (van der Walt et al., 2011), Pandas (McKinney, 2010), PyMC3 (Salvatier et al., 2016), SciPy (Oliphant, 2007), STARRY (Luger et al., 2018), Tapir (Jensen, 2013), Theano (Theano Development Team, 2016), TTVFast (Deck et al., 2014), TTV2Fast2Furious (Hadden et al., 2018),

1 Introduction

The Transiting Exoplanet Survey Satellite (TESS, Ricker et al., 2015), a near all-sky transit survey that began science operations July 2018, is expected to find thousands of planets. This includes hundreds of small planets with radii ${\rm R<4\,R_{\oplus}}$, around nearby, bright stars (Barclay et al., 2018; Huang et al., 2018a). During the 2-year primary mission, TESS will monitor more than 200,000 pre-selected stars at 2-min cadence and will observe additional targets spread over most of the sky ($\approx$85%) in 30-min cadence Full-Frame-Image (FFI) mode (Ricker et al., 2015). The spacecraft carries four identical wide-field cameras that combine to produce a nearly continuous $24^{\circ}\times 96^{\circ}$ field-of-view (FOV). TESS uses this large FOV to observe thirteen partially overlapping sectors per ecliptic hemisphere, per year and started its survey in the southern ecliptic hemisphere. The spacecraft observes each sector for two consecutive orbits that cover an average time baseline of 27.4 days111The orbital period of TESS is not constant due to 3-body gravitational interactions between TESS, the Earth, and the Moon. This leads to slightly different baselines in each sector.. The increasing overlap of sectors toward the ecliptic poles provides Continuous Viewing Zones (CVZs) surrounding the poles where targets receive $\approx$350 days of coverage. The long observing duration of the TESS CVZs will enable the detection of smaller and longer period planets. It will also overlap with the CVZs of the James Webb Space Telescope (JWST), providing key targets for detailed characterization. In about a hundred days of observations, TESS has already identified more than a hundred planet candidates, provided key observations to confirm several new planets, and provided new data on known transiting systems (Huang et al., 2018b; Gandolfi et al., 2018; Vanderspek et al., 2018; Wang et al., 2018; Nielsen et al., 2018; Shporer et al., 2018; Dragomir et al., 2018; Quinn et al., 2019; Rodriguez et al., 2019).

Planets discovered around bright, nearby stars provide ideal targets for mass measurements via Doppler spectroscopy, emission and transmission spectroscopy for atmospheric characterization, and for precise stellar characterization. Multi-planet systems provide an additional layer of information on planet formation and evolution, orbital dynamics, planetary architectures (e.g. Lissauer et al., 2011; Fabrycky et al., 2014), and in some cases mass measurements via transit-timing variations (e.g. Hadden & Lithwick, 2016; Hadden et al., 2018). While NASA’s Kepler and K2 missions successfully discovered thousands of planets around stars in the Kepler field and in the vicinity of the ecliptic plane (e.g. Rowe et al., 2014; Morton et al., 2016; Livingston et al., 2018), TESS will perform a nearly all-sky survey focused on stars in the solar neighborhood and find the touchstone planets that will be prime targets for observations with the Hubble Space Telescope (HST), JWST, and future ground-based observatories (Louie et al., 2018; Kempton et al., 2018).

Here we report the TESS discovery of three small planets transiting the bright (K = 7.1 mag), nearby (10.6 pc) M3-dwarf L 98-59. This paper is organized as follows. In Section 2 we describe the TESS observations and data analysis, as well as our ground-based follow-up efforts. In Section 3 we discuss the properties of the system, and draw our conclusions in Section 4.

2 Observations and Data Analysis

2.1 TESS Observations and Stellar Parameters

TESS observed L 98-59 (TIC 307210830, TOI-175; RA = 08:18:07.62, Dec = -68:18:46.80 (J2000)) in Sectors 2, 5, and 8 with Camera 4. The target was added to the TESS Candidate Target List—a list of targets prioritized for short-cadence observations (Stassun et al., 2018)—as part of the specially curated Cool Dwarf list (Muirhead et al., 2018). The TESS data were processed with the Science Processing Operations Center Pipeline (SPOC; Jenkins et al., 2016) and with the MIT Quick Look Pipeline. The three candidates identified by the SPOC pipeline passed a series of data validation tests (Twicken et al., 2018; Li et al., 2019) summarized below, and were made publicly available on the MIT TESS Data Alerts website222https://tess.mit.edu/alerts/ and the Mikulski Archive for Space Telescopes (MAST) TESS alerts page333https://archive.stsci.edu/prepds/tess-data-alerts/ as TOI-175.01, -175.02 and -175.03. These candidates had periods P = 3.690613 days, 7.451113 days, and 2.253014 days and transit epochs (BTJD) = 1356.203764, 1355.2864 and 1354.906208, respectively444BTJD = BJD - 2457000, and are referred to in the rest of the manuscript as L 98-59 c, L 98-59 d, and L 98-59 b respectively. The SPOC simple aperture photometery (SAP) and pre-search data conditioned (PDCSAP) lightcurves (Smith et al., 2012; Stumpe et al., 2014) of L 98-59 are shown in Figure 1.

We use the methods appropriate for M dwarfs previously used by (Berta-Thompson et al., 2015; Dittmann et al., 2017; Ment et al., 2018) to determine the stellar parameters of the host star, and adopt these parameters throughout our analysis. We estimate the mass of the star using the mass-luminosity relation in the $K-$band from Benedict et al. (2016) to be 0.313$\pm$0.014 M*⊙. We then use single star mass-radius relations (Boyajian et al., 2012) to find a stellar radius of 0.312$\pm$0.014 R⊙. We calculate the bolometric correction in $K$ from Mann et al. (2015, erratum) to be 2.7$\pm$0.036 mag, resulting in a bolometric luminosity for L 98-59 of 0.011$\pm$0.0004 L⊙. We calculate the correction in $V$ from Pecaut & Mamajek (2013) to be -2.0$\pm$0.03 mag555We assume the uncertainty on the bolometric correction in $V$ is that of the ($V-K$) color., resulting in a bolometric luminosity of 0.0115$\pm$0.0005 L⊙. We adopt the mean of the two bolometric luminosities from which we calculate the luminosity of the host star to be 0.0113$\pm$0.0006 L⊙* (i.e. 4.31e24 W). From the Stephan-Boltzmann Law, we find an effective temperature ${\rm T_{eff}=3367\pm 150K}$. As a comparison, we also used the relations in Mann et al. (2015) to determine an effective temperature of 3419$\pm$77 K for L 98-59, in agreement with the $T_{eff}$ derived from the Stefan-Boltzmann Law.

In addition, following the procedures described in Stassun & Torres (2016) and Stassun et al. (2017) to fit a NextGen stellar atmosphere model (Hauschildt et al., 1999) to broadband photometry data from Tycho-2, Winters et al. (2015), Gaia, 2MASS, and WISE, we performed a full fit of the stellar spectral energy distribution (SED) to estimate the stellar $T_{\rm eff}$ and [Fe/H] which, together with the Gaia DR2 parallax, provides an estimate of the stellar radius. The free parameters of the fit were $T_{\rm eff}$ and stellar metallicity [Fe/H], and we set the extinction $A_{V}\equiv 0$ due to the very close distance of the system. The resulting best fit is shown in Figure 2, with a reduced $\chi^{2}$ of 3.8 for 8 degrees of freedom. The best fit parameters are $T_{\rm eff}=3350\pm 100$ K and [Fe/H] $=-0.5\pm 0.5$. Integrating the SED gives the bolometric flux at Earth as $F_{\rm bol}=2.99\pm 0.18\times 10^{-9}$ erg s*-1* cm*-2*. Finally, adopting the Gaia DR2 parallax and the correction of 80 $\mu$as from Stassun & Torres (2018), we calculate a stellar radius of $0.305\pm 0.018$ R⊙ using the Stefan-Boltzmann law. These are consistent with the adopted parameters discussed above.

We also estimated a prior on the stellar density ($\rho_{*}$) by estimating the stellar mass. Here we used the empirical relations of Mann et al. (2015), which provides $M_{star}\approx 0.32M_{\odot}$ from the absolute $K_{S}$ magnitude ($M_{K_{S}}$) determined from the observed 2MASS $K_{S}$ magnitude and the Gaia DR2 parallax (corrected for the offset from Stassun & Torres, 2018). The quoted uncertainty in the Mann et al. empirical relation is $\sim$3%; here we conservatively adopt an uncertainty of 10%. Together with the radius determined above from the SED and parallax, we obtain $\rho_{*}=15.9\pm 3.3$ g cm*-3* (this value is used as a prior in the transit model). We also determined the stellar temperature and radius using empirical relations calibrated using low-mass stars with interferometrically measured radii and precise distances (see Feinstein et al., 2019, and references therein). These alternative stellar parameter estimates were consistent with those determined from the empirical M-dwarf relations and from the SED fitting. An additional set of stellar parameters for L 98-59 were previously derived from a medium-resolution optical spectrum in the CONCH-SHELL survey (Gaidos et al., 2014). This work also provides parameters consistent with our estimates. We compile the stellar parameters used in subsequent analyses and other identifying information for L 98-59 in Table 1. While it is difficult to pin down the ages of old M dwarfs due to their long main-sequence lifetime, the lack of a rapid rotation signal in the TESS SAP light curve and the low activity of L 98-59 (see § 2.4) indicate that it is likely an old M dwarf with an age $>$1 Gyr. The stellar parameters of L 98-59 are consistent with a spectral type of M3 $\pm$ 1.

2.2 Lightcurve Analyses

We opted to create our own apertures from Sectors 2, 5 and 8 target pixel file data to analyze the photometric time-series from TESS, instead of using the pipeline apertures used to first identify the transiting planet candidates. Our primary motivation for performing our own photometry is that we can avoid any attenuation to the transit signals by explicitly masking them during the systematic correction step. We first used the lightkurve package (Lightkurve Collaboration et al., 2018)666https://github.com/KeplerGO/lightkurve to extract lightcurves from each of the three sectors using the threshold method, which selects pixels that (a) are a fixed number of standard deviations above the background, and (b) create a contiguous region with the central pixel in the mask. We used a threshold value of 3$\sigma$; the corresponding mask shape is shown in Figure 3. Thus produced, the resulting lightcurve still contains low-level instrumental systematic signals. To identify and subtract instrumental signals, we used a second order pixel-level de-correlation (PLD), which is a technique based on Spitzer and K2 analysis methods (Deming et al., 2015; Luger et al., 2016). During the PLD step, we masked out transits to avoid attenuating the signals. Finally, we normalized the lightcurve by dividing by the median and subtracting one to center the flux about zero. We did the PLD detrending separately for each sector.

We used the exoplanet toolkit for probabilistic modeling of the exoplanet transits (Foreman-Mackey, 2018). The model we built consisted of four elements: three planet transit components with Keplerian orbits and limb-darkened transits, and a Gaussian Process (GP) component that models residual stellar variability. The planet models were computed with exoplanet using STARRY (Luger et al., 2018), while the GP was computed using celerite (Foreman-Mackey et al., 2017; Foreman-Mackey, 2018). The GP component is described as a stochastically-driven, damped harmonic oscillator with parameters of $\log(S_{0})$ and $\log(\omega_{0})$, where the power spectrum of the GP is

[TABLE]

and a white noise term, with a model parameter of the log variance. We fixed Q to $1/\sqrt{2}$ and put wide Gaussian priors on $\log(S_{0})$ and $\log(\omega_{0})$ with means of the log variance, and one log of one tenth of a cycle, respectively, and a standard deviation on the priors of 10. For each Sector we used separate GP parameters. This form of GP was chosen because of the flexible nature and smoothly varies (it is once mean square differentiable) enables us to use it to model a wide range of low frequency astrophysical and instrumental signals. The white noise term carried the same prior as $\log(S_{0})$. Each sector of data had a separate parameter for the mean flux level.

The planet model was parameterized in terms of consistent limb darkening, log stellar density, and stellar radius for the three planets. Each individual planet was parameterized in terms of log orbital period, time of first transit, log planet-to-star radius ratio, impact parameter, orbital eccentricity and periastron angle at time of transit. The stellar radius had a Gaussian prior with mean 0.312 and 0.014 standard deviation, with solar units, and is additionally required to be positive. The log mean stellar density, in cgs units, had a Gaussian prior with a mean of $\log{15}$ and standard deviation of 0.2 dex (as per Section 2.1). The limb darkening followed the Kipping (2013a) parameterization.

The log orbital periods, time of first transits, and log planet-to-star radius ratio of the three planets had Gaussian priors with means at the values found in the TESS alert data, and standard deviations of 0.1, 0.1, and 1, respectively. The impact parameter had a uniform prior between zero and one plus the planet-to-star radius ratio. Eccentricity had a beta prior with $\alpha=0.867$ and $\beta=3.03$ (as suggested by Kipping, 2013b), and was bounded between zero and one. The periastron angle at transit was sampled in vector space to avoid the sampler seeing a discontinuity at values of $\pi$.

We sampled the posterior distribution of the model parameters using the No U-turn Sampler (NUTS, Hoffman & Gelman, 2014) which is a form of Hamiltonian Monte Carlo, as implemented in PyMC3 (Salvatier et al., 2016). We ran 4 simultaneous chains, with 5000 tuning steps, and 3000 draws in the final sample. The effective number of independent samples of every parameter was above 1000, and most parameters were above 5000. The Gelman–Rubin diagnostic statistic was within 0.005 of 1.000 for each parameter in the model. The impact parameter for the outer planet is relatively high, which caused this parameter along with the orbital eccentricity to be most time consuming to sample independently.

Figure 4 shows the GP model of the low-level variability in the upper panels, and the best fitting transit model in the central panels. The phase-folded transits of the three candidates, along best-fitting transit models, are shown in Figure 5, and the model parameters are provided in Table 2. The transit modeling reveals that the candidate planets have small radii ranging from 0.8 to 1.6 $R_{\oplus}$. A chain of small terrestrial-sized planets is common among M-dwarfs (Muirhead et al., 2015), and L 98-59 is reminiscent of other systems such as the TRAPPIST-1, Kepler-186 and Kepler-296 (Gillon et al., 2017; Barclay et al., 2015; Quintana et al., 2014). The stellar density obtained from the transit model is fully consistent with that determined from the stellar parameters in Section 2.1 (${\rm 15.8^{+2.6}_{-2.7}~{}g~{}cm^{-3}}$ for the former vs ${\rm 15.9\pm 3.3~{}g~{}cm^{-3}}$ for the latter.

We repeated this analysis using the systematics-corrected light curves from the TESS pipeline (PDCSAP Jenkins et al., 2016; Stumpe et al., 2014) rather than using the respective target pixel files. We found consistent results, aside from different GP parameters owing to the different systematics corrections applied. The transits depths were lower in the PDCSAP data at the ${\rm<1\sigma}$ level, which we attribute to masking out transits in the systematics-correction technique we applied. We did not include any flux contamination from nearby stars in our models because there are no bright nearby stars to contaminate our pixel mask – the TIC estimates that the contamination fraction for L 98-59 is 0.002. Even if the TIC contamination is dramatically underestimated, it is highly likely that the stellar radius uncertainty will be the dominant term in the planet radius uncertainty, therefore we feel comfortable neglecting it.

There are significant impact parameter differences between the inner two planets and the outer planet. The outer planet transits close to the limb of the star, although it is not grazing. This manifests in the light curve as a shorter transit duration for the outer planet than the two inner planets and indicates a modest mutual inclination of at least one degree between the inner two and outer planets. All three planets are subject to significantly more flux than Earth receives from the Sun, and are therefore unlikely to be good astrobiology targets. However, with an insolation flux of $4.5\pm 0.8$, the outer planet is a candidate Venus-zone planet (Kane et al., 2014).

2.3 Potential false positive scenarios

The TESS Data Validation Report performs a series of tests designed to rule out various false positive scenarios. The results from these tests are as follows:

i) L 98-59 c and L 98-59 d pass the difference image centroiding test, which employs PSF-based centroiding on the difference images (expected to be more precise and accurate than a brightness-weighted moment on the difference images). While L 98-59 b does not quite pass the difference image centroiding test, its transits are much shallower compared to the other two candidates. Thus it is likely that the centroiding errors are underestimated to some degree due to the variable pointing performance at timescales less than the 2 minutes observation cadence. We expect the analysis of this candidate to improve with new data;

ii) All three candidates pass the odd-even difference tests; iii) Secondary eclipses are ruled out at the 3.6-, 2.6- and 1.8-${\sigma}$ levels;

iv) A bootstrap analysis of the out-of-transit data is used to quantify the probability of false alarms due to stellar variability and residual instrumental systematics. In the case of L 98-59, the light curve is well behaved and the analysis excludes the possibility of a false alarm at the 2.45E-25, 6.6E-62, and 2.2E-25 levels (as extrapolations of the upper tail of the bootstrap distribution to the observed maximum Multiple Event Statistics (MES) that triggered the detections of these candidates in the pipeline (Jenkins et al., 2017));

v) All three candidates pass a ghost diagnostic test, designed to flag instances of scattered light, other instrumental artifacts or background eclipsing binaries.

For completeness, we also applied the vetting pipeline DAVE (Kostov et al., 2019) to the TESS lightcurve of L 98-59. Briefly, DAVE evaluates whether detected transit-like events produced by the candidate are real or false positives by analyzing the data for (a) odd-even differences between consecutive transits; (b) secondary eclipses; (c) stellar variability mimicking a transit; and (d) photocenter shifts during transit;

To perform the (a), (b) and (c) analysis we used the Modelshift module of DAVE—an automated package designed to emphasize features in the lightcurve that resemble the shape, depth and duration of the planetary transit but located at different orbital phase. To identify secondary eclipses and odd-even transit differences, or flares and heartbeat stars (Welsh et al., 2011, see e.g.), Modelshift first convolves the lightcurve with the transit model of the planet candidate. The module then computes the significance of the primary transits, odd-even differences, secondary, tertiary and positive features assuming white noise in the lightcurve, and compares the ratio between each of these and the systematic red noise ${\rm F_{red}}$ to the false alarm thresholds ${\rm FA_{1}=\sqrt{2}erfcinv(T_{dur}/(P\times N))}$ (assuming 20,000 objects evaluated), and ${\rm FA_{2}=\sqrt{2}erfcinv(T_{dur}/P)}$ (for two events), where ${\rm T_{dur},P,N}$ are the duration, period and number of events (see Coughlin et al. 2014 for details). For example, a secondary feature is considered significant if ${\rm Sec/F_{red}>FA_{1}}$. The Modelshift results are shown in Figures 6, 7, and 8 where the panels show the phase-folded lightcurve (first row), the phase-folded lightcurve convolved with the best-fit transit model (second row), as well as the the best-fit to all primary transits, all odd and all even transits, the most prominent secondary, tertiary and positive features in the lightcurve (lower two rows). The tables above the figures list the individual features evaluated by the module: the significance of the primary (“Pri”), secondary (“Sec”), tertiary (“Ter”) and positive (“Pos”) events assuming white noise, along with their corresponding differences (“Pri-Ter”, “Pri-Pos”, “Sec-Ter”, “Sec-Pos”), the significance of the odd-even metric (“Odd-Evn”), the ratio of the individual depths’ median and mean values (“DMM”), the shape metric (“Shape”), the False Alarm thresholds (“${\rm FA_{1}}$”, “${\rm FA_{2}}$”), and the ratio of the red noise to the white noise in the phased light curve at the transit timescale (“Fred”). Our analysis shows that there are no secondary eclipses or odd-even differences for any of the L 98-59 planet candidates. We note that the significant secondary and tertiary eclipses identified by DAVE for L 98-59 d (Figure 7) are due to the transits of L 98-59 c and thus not a source of concern.

To perform the (d) analysis for each candidate we used the photocenter module of DAVE, following the prescription of Bryson et al. (2013). Specifically, for each candidate we: 1) create the mean in-transit and out-of-transit images for each transit (ignoring cadences with non-zero quality flags), where the latter are based on the same number of exposure cadences as the former, split evenly before and after the transit; 2) calculate the overall mean in-transit and out-of-transit images by averaging over all transits; 3) subtract the overall mean out-of-transit image from the overall in-transit image to produce the overall mean difference image; and 4) measure the center-of-light for each difference and out-of-transit image by calculating the corresponding $x$- and $y$-moments of the image. The measured photocenters for the three planet candidates are shown in Figures 9, 10, and 11, and listed in Table 3. We detect no significant photocenter shifts between the respective difference images and out-of-transit images for any of the planet candidates (see Table 3), which confirms that the target star is the source of the transits. We note that some of the individual difference images for L 98-59 b deviate from the expected Gaussian profile, and thus so does the mean difference image.

Overall, our DAVE results rule out false-positive features for all three planet candidates of L 98-59, are consistent with the analysis of the Data Validation Report, and indicate that the detected events are genuine transits associated with the star in question. We also note that while the automated vetting of Osborn et al. (2019) flagged L 98-59 with “a high likelihood of being astrophysical false positives” (their Table 3), their subsequent manual vetting lists the system as planet candidate.

Additionally, to investigate whether one or more of the transits associated with L 98-59 may result from nearby sources (e.g. a background eclipsing binary), we used lightkurve to extract lightcurves for nearby field stars. Our analysis revealed that a nearby field star ${\sim 80\arcsec}$ NW of L 98-59 (2MASS 08175808-6817459, TIC 307210817, Tmag = 13.45, i.e. $\approx 4$ magnitudes fainter than L 98-59) is in fact an eclipsing binary (EB), manifesting both primary and secondary eclipses at a period of ${\approx 10.43}$ days, with ${T_{0}=4.4309}$ (BJD-2,455,000) (see Figure 12). This field star could be associated with one of two sources in the Gaia catalog: Source 1 with RA = 124.49217149500, Dec = -68.29612321000, ID = 5271055685541797120, and parallax=0.1888 mas; and Source 2 with RA = 124.49126472300, Dec = -68.29602748080, ID = 5271055689840223744, and parallax = 0.9977 mas. Given the corresponding approximate distances of 1 and 5 kpc, neither of these targets can be physically associated with L 98-59 as they lie deep in the background. Regardless of which of these sources hosts the detected EB, the faintness of the host compared to L 98-59, the measured EB orbital parameters, and the dilution-corrected eclipse depths are inconsistent with the properties of the candidate planets and effectively rule it out as the potential source of any of these signals.

2.4 Follow-up Observations

We pursued ground-based follow up of the three candidates to rule out potential sources of false positives and strengthen the evidence of their planetary nature. Our L 98-59 follow-up program was organized through the TESS Follow-up Observing Program (TFOP) Working Group (WG)777https://tess.mit.edu/followup/ which facilitates follow-up of TESS candidate systems. The primary goal of the TFOP WG is to provide follow-up observations that will advance the achievement of the TESS Level One Science Requirement to measure masses for 50 transiting planets smaller than 4 Earth radii. A secondary goal of the TFOP WG is to foster communication and coordination for any science coming out of TESS. Our L 98-59 follow-up was conducted by three TFOP Sub Groups (SGs): SG-1, seeing limited photometry; SG-2, reconaissance spectroscopy; and SG-3, high-resolution imaging.

2.4.1 Seeing-Limited Photometry from the TFOP WG

Analysis of multi-planet systems from Kepler has shown that these have a higher probability of being real planets (e.g., Lissauer et al., 2012), lending credibility to the planetary nature of the transit events associated with L 98-59. However, the pixel scale of TESS is larger than Kepler’s (${21\arcsec}$ for TESS vs. ${4\arcsec}$ for Kepler) and the point spread function of TESS could be as large as ${1\arcmin}$, both of which increase the probability of contamination by a nearby eclipsing binary (EB). For example, a deep eclipse in a nearby faint EB might mimic a shallow transit observed on the target star due to dilution. Thus it is critical to explore the potential contamination of relatively distant neighbors in order to confirm transit events detected on a TESS target.

To identify potential false positives due to variable stars such as EBs up to ${2.5\arcmin}$ away from L 98-59, we made use of the TFOP SG1. Specifically, we observed the target with ground-based facilities at the predicted times of the planet transits to search for deep eclipses in nearby stars at higher spatial resolutions. We used the TESS Transit Finder, which is a customized version of the Tapir software package (Jensen, 2013), to schedule photometric time-series follow-up observations. The facilities we used to collect TFOP SG1 data are: Las Cumbres Observatory (LCO) telescope network (Brown et al., 2013); SPECULOOS South Observatory (SSO) (Delrez et al., 2018; Burdanov et al., 2018); MEarth-South telescope array (MEarth) (Irwin et al., 2015); and Siding Spring Observatory T17 (SSO T17). Detailed observation logs are provided in Table 4.

We used the AstroImageJ software package (Collins et al., 2017) for the data reduction and the aperture photometry in most of these follow-up photometric observations. For the observations carried out at SSO, the standard calibration of the images and the extraction of the stellar fluxes were performed using the IRAF/DAOPHOT aperture photometry software as described in Gillon et al. (2013). The results are shown in Figures 13, 14 and 15.

For all three planet candidates, we confirmed that the target star is the source of the transits, and ruled out nearby EBs which could mimic the transits. In addition, observations of L 98-59 c and d in different filters showed no chromatic dependence, which strengthens the hypothesis that the candidates are real planets. Our follow-up did detect a nearby EB at a separation of 54″ (TIC 307210845, Tmag = 16.042, i.e. $\sim 7$ magnitudes fainter than TIC 307210830, producing no detectable eclipses in the lightcurve of the latter), and we used deep exposures to confirm that it is not the origin of the L 98-59 b transits, providing high level of confidence about the planetary nature of this candidate.

The measured transit depths revealed by follow-up transit photometry are consistent with the transit depths measured from TESS. The differences in the follow-up and TESS transit depth measurements (in terms of $Rp/R*$) are listed in Table 5 as a function of wavelength, where we have included only the transits with scatter low enough to reasonably detect the events.

2.4.2 Reconnaissance Spectroscopy

To investigate the magnetic activity and rotation of L 98-59 and rule out spectroscopic binary companions, we obtained two epochs of optical spectra of L 98-59 on UT 2018 February 12888The first epoch of spectroscopic data were obtained as part of the M dwarf spectroscopic program described in Winters et al. (2019), before $TESS$ began observations. and on UT 2018 November 20 using the slicer mode with the CTIO HIgh ResolutiON (CHIRON) spectrograph (Tokovinin et al., 2013) (R $\simeq$ 80,000) on the Cerro Tololo Inter-American Observatory (CTIO)/ Small and Moderate Aperture Research Telescope System (SMARTS) 1.5m telescope. CHIRON has a spectral range of 410-870 nm. We obtained one 7.5-minute exposure for the first epoch and 3 $\times$ 2.5-minute exposures for the second observation, yielding a signal-to-noise ratio of roughly 13 in order 44 for both epochs. As described in Winters et al. (2018), we used an observed template of Barnard’s Star to derive a radial velocity of -5.8$\pm$0.1 km s*-1* using the TiO molecular bands at 7065–7165 Å999We note that the total uncertainty on the systemic velocity should include the 0.5 km s*-1* uncertainty on the Barnard’s Star template velocity.. Our analyses of the spectra reveal no evidence of double lines. We see negligible rotational broadening ($v$sin$i$ = $0.0\pm 1.9$ km/s), and do not see H$\alpha$ in emission, providing evidence that the star is inactive and not host to unresolved, close-in, stellar companions. We are also able to rule out the presence of a brown dwarf companion to the host star. The radial velocity difference between the two observations, separated in time by roughly nine months, is 53$\pm$52 m s*-1*. For comparison, a 13 Jupiter-mass companion in a circular, nine-month period would induce a velocity semi-amplitude of 863 m s*-1* on this star. This semi-amplitude is eight times larger than our velocity uncertainty and would have been readily detectable. Thus, it is highly unlikely that there is a low-mass stellar or brown dwarf companion around L 98-59 at periods shorter than nine months.

We note, as well, that the trigonometric distance of 10.623$\pm$0.003 pc from the Gaia second data release is in agreement with the photometric distance estimate of 12.6$\pm$1.9 pc reported in Winters et al. (2019). Undetected equal-luminosity companions would contribute light to the system, making the system overluminous and resulting in an underestimated photometric distance estimate. Because the two distances are in agreement, this lends further support to the host star being a single star.

We also placed L 98-59 on an observational Hertzsprung-Russel color-magnitude Diagram (blue star; see Figure 16). Because it is not elevated above the main sequence or among the blended photometry binary sequence (red points), even more strength is given to the argument of this star being single.

In addition, we obtained a near-IR spectrum of L 98-59 on 2018 December 22 with the Folded-port InfraRed Echellete (FIRE) spectrograph (Simcoe et al., 2008) on the 6.5 Baade Magellan telescope at Las Campanas observatory. FIRE covers the 0.8-2.5 micron band with a spectral resolution of $R=6000$. The target was observed under favorable conditions, with an average seeing of $\sim 0\farcs 6$. L 98-59 was observed twice in the ABBA nod patterns at 40s integration time for each frame using the $0\farcs 6$ slit. Reductions and telluric corrections, using the nearby A0V standard HIP 41451, were completed with the FIREhose IDL package. We derived stellar parameters following the empirical methods derived by Newton et al. (2015). For L 98-59, we infer: $T_{eff}=3620\pm 74$K, $R_{star}=0.37\pm 0.027R_{\odot}$, and $L=0.021\pm 0.004L_{\odot}$., consistent with the SED analysis.

2.4.3 High Resolution Imaging

Photometric contamination from nearby sources can result in various false positive scenarios (e.g. background eclipsing binaries), and can bias the measured planetary radius from photometric analysis (see e.g. Ciardi et al., 2015; Furlan & Howell, 2017; Ziegler et al., 2018). In this work we use several high resolution images to tightly constrain the possible background sources and companion stars present near L 98-59. Previous speckle observations of the target were collected with Gemini/DSSI on 2018 March 31 as part of the M dwarf speckle program described in Winters et al. (2019). Once the candidate planets in this system had been identified by TESS, we collected additional speckle images with Gemini/DSSI (Horch et al., 2011) on 2018 November 01, and AO images with VLT/NaCo on 2019 January 28. Both epochs of Gemini/DSSI data are collected simultaneously through R and I band filters (692nm and 880nm respectively), while the VLT/NaCo data are collected in the Br$\gamma$ filter. At 0.5” from the host, these data rule out companions 5.0, 6.6 and 5.8 magnitudes fainter than the host star in the R, I and Br$\gamma$ bands respectively. The 5-$\sigma$ contrast curves for each of these observations are presented in Figure 17.

Due to the high proper motion of the target (354 mas/yr), the target undergoes significant motion even over the relatively modest time baseline of these observations. The on-sky position of the target is displaced by 338 mas between the observations on 2018 March 31 and 2019 January 28. This motion is significantly more than the PSF width in the high resolution images, and we are therefore able to rule out the presence of stationary background objects within $\sim$6 magnitudes of the host at any separation: any background objects obscured by the target in the first observation would be clearly visible in the final observation and vice versa. The motion of the target is demonstrated in Figure 17. These data also allow tight constraints to be placed on the presence of co-moving companions beyond $\sim$150 mas (=1.6 AU at the distance of this target). An object with $\Delta$Br$\gamma\sim$6 mag and at the distance of L 98-59 would have a mass of $\sim 75$MJ at an age of 10 Gyr (Baraffe et al., 2003)101010Given the difficulty in estimating M-dwarf ages, we use an age of 10 Gyr so as to calculate a conservative mass limit., and we can therefore rule out any stellar companions to this host with a projected separation greater than 1.6 AU, while stellar companions closer than this could be easily detected in radial velocity data with a sufficient baseline.

Overall, the follow-up efforts demonstrate that there are tentative transit detections of the various candidates from the ground but it is difficult to confirm them since they are all shallow events. More importantly, none of these detections were identified as eclipsing binary scenarios (either on the host star or on nearby stars), supporting the planetary nature of the three transit candidates.

3 Discussion

Our lightcurve and false positive analyses, follow-up observations, and the multiplicity of the system provide strong evidence that the detected transit signals are planetary in nature. We consider the planets L 98-59b, c, and d to be a validated system of terrestrial planets orbiting a very nearby, bright M dwarf. Here we explore additional properties of the system to place constraints on the planet masses, orbital dynamics and evolution, and discuss their potential for future characterization.

3.1 Planet Mass Constraints

In the absence of radial velocity measurements, we placed constraints on the masses of the planets using the forecaster package for probabilistic mass forecasting (Chen & Kipping, 2017). From the mean and standard deviation of each planet’s radius, we generated a grid of 5000 masses within the entire mass range of the conditioned model, which spans dwarf planets to high-mass Jovians, and sampled 50000 times from a truncated normal distribution. For each sampled radius, forecaster computes a vector of probabilities given each element in the mass grid and a randomly chosen set of hyper-parameters from the hyper-posteriors of the model (which include transition points and intrinsic dispersion in the mass-radius relation). From this vector, the package returns the median mass and $\pm 1\sigma$ values. From the calculated radius values of 0.8 [0.05] ${R_{\oplus}}$ (L 98-59 b), 1.35 [0.07] ${R_{\oplus}}$ (L 98-59 c), and 1.57 [0.14] ${R_{\oplus}}$ (L 98-59 d), we determined mass values of 0.5 [+0.3, -0.2] ${M_{\oplus}}$, 2.4 [+1.8, -0.8] ${M_{\oplus}}$, and 3.4 [+2.7, -1.4] ${M_{\oplus}}$, respectively. The large errors on these values suggest that better constrained radii from continued follow-up observations, combined with precise radial velocity measurements, are necessary to constrain the true masses. We note that given the brightness of the host star, the L 98-59 planets should be great targets for mass measurements to establish the M-R relation for M-dwarf planets. Using the forecaster masses, the expected radial velocity semi-amplitude, K, for the three planets are 0.54, 2.22, and 2.48 m/s for L 98-58 b, L 98-58 c, and L 98-58 d respectively—the outer two comparable to the amplitude of the measured radial velocity signal produced by e.g. GJ 581 b (Mayor et al., 2009) and Proxima Centauri b (Anglada-Escudé et al., 2016).

3.2 Dynamical stability and transit-timing variations

3.2.1 Long-Term Stability

To examine the long-term dynamical stability and orbital evolution of the L 98-59 planets, we integrated the system using Rebound (Rein & Spiegel, 2015) for 1 million orbits of the outer planet (${\rm P=7.45}$ days). We used two sets of initial conditions – planets on circular orbits, and on eccentric orbits with $e=0.1$ – and start all integrations with randomly-selected initial arguments of periastron. Given the large uncertainties on the forecaster masses, we also tested the dynamical stability for two sets of planetary masses: best-fit masses (i.e. 0.5 ${M_{\oplus}}$, 2.4 ${M_{\oplus}}$ and 3.4 ${M_{\oplus}}$ for L 98-58 b, L 98-58 c and L 98-58 d respectively), and best-fit+${1\sigma}$ masses (i.e. 0.58 ${M_{\oplus}}$, 4.2 ${M_{\oplus}}$ and 6.1 ${M_{\oplus}}$ for L 98-58 b, L 98-58 c and L 98-58 d respectively). Overall, we performed 1000 numerical simulations for each set of planetary masses, using the IAS15 non-symplectic integrator (Rein & Spiegel, 2015) with a timestep of 0.01 the orbit of the inner planet (i.e. about 30 min).

Our simulations show that for initially circular orbits, the semi-major axes and eccentricities do not exhibit extreme variations, the system does not exhibit chaotic behavior for the duration of the numerical integrations for either set of planet masses, and the orbits remain practically circular (Figure 18). In contrast, for initially eccentric orbits with $e=0.1$, the system becomes unstable in half of our simulations (with randomly-selected initial arguments of periastron), both for the best-fit and the best-fit+${1\sigma}$ planet masses (Figure 19). Thus we consider orbits with non-negligible eccentricity as unlikely. This is consistent with other compact multiplanet systems where the orbital eccentricities are typically on the order of a few percent (e.g. Hadden & Lithwick, 2014), and is in line with the L 98-59 planets being close to but not in resonance (where the orbits may potentially be eccentric, e.g. Charalambous et al., 2018).

Inspired by the closely-spaced multiplanet systems discovered by Kepler (Muirhead et al., 2015) and other surveys (e.g. Gillon et al. 2017), we also explored the possibility of a fourth, non-transiting planet having a dynamically-stable orbit in-between L 98-59 c and L 98-59 d such that the four planets would form a near-resonant chain of 5:8:12:16 period commensurability similar to e.g. TRAPPIST-1 (Gillon et al. 2017). As an example, we tested a planet with a mass of ${\rm 2.5~{}M_{\oplus}}$ and a 5.7-day orbital period (${\rm\approx 1.55}$ and ${\rm\approx 0.77}$ times the period of L 98-59 c and L 98-59 d respectively), again for two cases of (a) initially circular orbits and (b) initially eccentric orbits with $e=0.1$. For simplicity, we only used the best-fit masses for the three planet candidates. The system is dynamically stable in case (a) for the duration of the integrations (Figure 20, upper panel), and becomes unstable within a few thousand orbits of the outer planet for case (b). Thus such a hypothetical planet is potentially possible if on a circular orbit. Overall, while a comprehensive dynamical analysis for the presence of additional planets is beyond the scope of this work, we will continuously monitor the system as data from future TESS sectors become available.

3.2.2 Transit Timing Variations

If detected, deviations in the times of transits from a linear ephemeris can be a powerful method to constrain the masses and orbital eccentricities of planets in multiplanet systems (e.g. Agol et al. 2005). We measured transit times for each individual transit using two different methods. First, we folded the transits on a linear ephemeris, fitting a transit model using the models of Mandel & Agol (2002). Next, we measured the time of each individual transit by, for each transit, sliding this model across a grid of potential transit midpoints with a time resolution of one second, and measuring the likelihood of each transit fit at each grid point. We then found the maximum likelihood transit time and a 68% confidence interval on the same. The ability to measure TTV signals depends sensitively on our ability to measure precise transit times. For L 98-59 b, the scatter in measured transit times, suggestive of the ultimate transit timing precision we measure, is 5.1 minutes. For L 98-59 b, this is 2.1 minutes; and for L 98-59 d, 1.2 minutes. Our analysis showed that a linear ephemeris is sufficient to reproduce the transit times of the three planet candidates detected in Sector 2. We found no evidence for transit timing variations (TTVs), and no further constraints can be placed on the parameters of the system beyond those already provided by dynamical stability considerations. Given that L 98-59 will be observed in 7 of the 13 sectors that comprise the first year of the TESS mission (Mukai & Barclay, 2017) (Sectors 2, 5, 8, 9, 10, 11, and 12), here we examine how continued TESS observations would affect the transit timing analysis of the system.

Specifically, we simulated continued observations following the nominal TESS schedule, assuming a linear ephemeris for future transits and that every planned sector will be observed as scheduled. We then used the TTVFast package (Deck et al., 2014) to calculate predicted transit times for the three planets in various orbital configurations consistent with the current data and examine what can be ruled out by the data by the end of the mission.

The two outer planets (L 98-59 c and d) are close to first-order period commensurability (period ratio of 2.02) whereas the inner planet is not near a first-order resonance with either of the other planets (1.64 period ratio between L 98-59 b and L 98-59 c, and 3.31 period ratio between L 98-59 b and L 98-59 d). Thus the expected TTV signal for the former planet pair is stronger compared to that for the latter planet pair. Indeed, even for eccentricities of $\sim 0.1$, the expected TTV amplitude for the innermost planet is $\sim 90$ seconds, notably smaller than the observed precision on the measured times of transit. To evaluate the potential for measuring TTVs for the outer two planets, we performed numerical simulations of the system for the first year of the TESS mission (using TTVFast), thus covering all sectors it will be observed in. We allowed the planet eccentricities to vary, and assumed the maximum likelihood masses listed in Section 3.1. Adopting a transit-timing precision of 1-2 minutes, we found that significant TTVs could be detected if the eccentricities of the outer two planets were larger than $\sim 0.03$ (as shown in Figure 21).

Overall, as multiplanet systems typically have orbital eccentricities of a few percent (Hadden & Lithwick, 2016), it is unlikely that TESS will reveal timing variations for this system during its primary mission; doing so would suggest either anomalously large eccentricities (which are unlikely based on dynamical stability consideration) or significantly larger planet masses/densities in this system relative to planets with similar radii in other systems.

To explore TTVs using an alternative software framework, we also used the TTV2Fast2Furious package (Hadden et al., 2018) to project the expected TTV signals of the planets through Sector 12, again adopting the forecaster masses and, for simplicity, assuming circular orbits. Similar to the TTVFast analysis described above, our results show that it is unlikely TTVs are measurable for this system during the TESS prime mission — the maximum TTV amplitudes are $0.09$ minutes for L 98-59 b, $0.17$ minutes for L 98-59 c, and $0.56$ minutes for L 98-59 d (see Figure 22).

Following the approach described in Hadden et al. (2018), we also used TTV2Fast2Furious to project the precision of mass constraints derived from future TESS transit timing measurements. Planet mass constraints derived from TTVs depend on the precision of transit time measurements and we adopt the measured scatter in the transit time measurements taken through Sector 2, $\sigma_{t_{L98-59c}}$ = 2.1 minutes, and $\sigma_{t_{L98-59d}}$ = 1.2 minutes. The planets’ masses are expected to be constrained with precisions $\sigma_{m,L98-59c}=13.1\,{\rm M}_{\oplus}$ and $\sigma_{m,L98-59d}=5\,{\rm M}_{\oplus}$ with transit timing data through Sector 12. 111111If the planets are restricted to circular orbits in the TTV model, e.g., under the assumption that eccentricities are damped away by tidal dissipation, the mass-eccentricity degeneracy (Lithwick et al., 2012) is removed and the measured TTV signals therefore place tighter constraints on the planet masses. In particular, if the TTV model is restricted to circular orbits, the mass measurement precisions of $\sigma_{m,L98-59c}=3.4M_{\oplus}$ and $\sigma_{m,L98-59d}=2.1M_{\oplus}$ are projected.

With three sectors of data available for this system at the time of writing (Sectors 2, 5, and 8), there is at present no evidence for TTVs. In line with the predictions of the previous paragraph, there is no significant mass constraint beyond what is given from plausible compositions of terrestrial planets. We use TTVFast to compare maximum-likelihood dynamical models of the orbits of the planets, assuming the masses listed in Section 3.1. We also repeat this procedure holding the masses fixed at four times this nominal value, which would imply a density of approximately 20 g/cc. In both scenarios, there is a dynamical model which fits the observed transits.

3.3 Potential for Atmospheric Characterization

Owing to the small, bright host star the three planets of the L 98-59 system are promising targets for follow-up atmosphere characterization. The planets’ small radii suggest that it is unlikely that they retain hydrogen-rich atmospheres (Rogers, 2015; Fulton et al., 2017), but secondary atmospheres could form from volcanic outgassing and/or delivery of volatiles from comets.

To investigate the feasibility of atmosphere studies for the L 98-59 planets, we compared the expected signal-to-noise of atmospheric features to that of GJ 1132b, another small planet around a nearby M-dwarf (Berta-Thompson et al., 2015). Morley et al. (2017) found that a CO2-dominated atmosphere could be detected for GJ 1132 with a modest number of JWST transits or eclipses (11 transits with NIRSpec/G235M or 2 eclipses with MIRI/LRS). To scale these estimates for L 98-59, we used the transmission and emission spectroscopy metrics from Kempton et al. (2018), which calculate expected signal-to-noise for atmospheric features based on planet and star properties. We scaled the signal relative to each planet’s transit/eclipse duration, estimated the brightness of the star based on its K-band magnitude, and assumed zero noise floor. We found that L 98-59 b, L 98-59 c and L 98-59 d have Transmission Spectroscopy Metric (TSM) values of 0.8, 1.4, and 1.0 that of GJ 1132b and Emission Spectroscopy Metric (ESM) values of 0.3, 0.4, and 0.7 respectively (Kempton et al., 2018). This implies that features in the transmission spectrum could be detected with 16, 6, or 11 transits, or 24, 13, or 4 eclipses for L 98-59 b, L 98-59 c and L 98-59 d. Provided that JWST observations reach the photon limit for stars as bright as L 98-59, this system is an exciting opportunity for studying comparative planetology of terrestrial exoplanet atmospheres.

3.4 Planets in the Venus Zone

It is worth noting the possibility that the planets in the system are analogs to Venus in terms of their atmospheric evolution. Venus shares several characteristics with Earth including its relative composition, size, and mass. Although Venus may have previously had temperate surface conditions (Way et al., 2016), Venus eventually diverged significantly from the habitable pathway of Earth and transitioned into a runaway greenhouse state. The planet now has a high pressure, high temperature, and a carbon dioxide dominated atmosphere. In our study of exoplanets and the search for life, it is vitally important that we understand why Earth is habitable and Venus is not (Kane et al., 2014). There is a need to discover planets that may have evolved into a post-runaway greenhouse state so that we can target their atmospheres for characterization with future facilities, such as JWST (Ehrenreich et al., 2012). However, most of the potential Venus analog candidates hitherto discovered orbit relatively faint stars (Barclay et al., 2013; Kane et al., 2013; Angelo et al., 2017; Kane et al., 2018).

The L 98-59 planets receive significantly more energy than the Earth receives from the Sun (a factor of between 4–22 more than Earth’s insolation) and fall into the region that Kane et al. (2014) dubbed the Venus Zone. This is a region where the atmosphere of a planet like Earth would likely have been forced into a runaway greenhouse, producing conditions similar to those found on Venus. The range of incident fluxes within the Venus Zone corresponds to insolations of between 1–25 times that received by the Earth. Planets in the Venus Zone that can be spectroscopically characterized will become increasingly important in the realm of comparative planetology that aims to characterize the conditions for planetary habitability. In that respect, and considering the potential for atmospheric characterization discussed in Section 3.3, L 98-59 could become a benchmark system.

4 Conclusions

We presented the discovery of a system of three transiting, terrestrial-size planets orbiting L 98-59 (TESS Object of Interest TOI-175). The host star is a bright M3 dwarf (K = 7.1) at a distance of 10.6 pc, with ${\rm M_{*}=0.313\pm 0.014M_{\odot}}$, ${\rm R_{*}=0.312\pm 0.014R_{\odot}}$, and ${\rm T_{eff}=3367\pm 150K}$. TFOP-led follow-up observations found no evidence of binarity or significant stellar activity. To thoroughly vet the transit signals detected in the TESS data, we used the software package DAVE. Our analysis ruled out significant secondary eclipses, odd-even differences or photocenter shifts during transits, verifying their planetary nature. Using lightkurve, we also discovered that the nearby field star 2MASS 08181825-6818430, inside the TESS aperture of L 98-59, is an eclipsing binary system with an orbital period of ${\sim 10.43}$ days, manifesting both primary and secondary eclipses. Utilizing Gaia data we confirmed that the eclipsing binary is a background object (likely a red giant) not associated with L 98-59. This battery of tests highlights the importance of comprehensive analysis of all sources inside the TESS aperture.

The planets range in size from slightly smaller to slightly bigger than Earth, with radii of ${\rm 0.8\pm 0.05R_{\oplus}}$, ${\rm 1.35\pm 0.07R_{\oplus}}$, and ${\rm 1.59\pm 0.23R_{\oplus}}$ from inner to outer respectively. The planetary system is quite compact, with orbital periods of ${\rm 2.25~{}days}$, ${\rm 3.69~{}days}$, and ${\rm 7.45~{}days}$ respectively. We estimated their masses using the forecast package for probabilistic mass forecasting, confirmed the dynamical stability of the system for circular orbits, and showed that there are no significant transit-timing variations.

TESS will continue observing the system in upcoming sectors (9, 10, 11, 12), and it is also likely that the system will be observed during a TESS Extended Mission. These observations will allow for refinement of the known planet parameters, searches for additional planets, further investigations of the dynamics of the system, as well as long-term monitoring of the host star activity.

We thank the referee for the insightful comments that helped us improve this manuscript. This manuscript includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST). Funding for the TESS mission is provided by NASA’s Science Mission directorate. We acknowledge the use of TESS Alerts data, as provided by the TESS Science Office. 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. JGW is supported by a grant from the John Templeton Foundation. Acquisition of the CHIRON data and the first epoch of DSSI data was made possible by a grant from the John Templeton Foundation. The opinions expressed in this publication are those of the authors and do not necessarily reflect the views of the John Templeton Foundation. We thank Leonardo Paredes, Hodari James, Rodrigo Hinojosa, and Todd Henry for their work in gathering and processing the CHIRON data, as well as for the management of CTIO / SMARTS 1.5m telescope. We are also grateful to the observer support staff at CTIO, at ESO/VLT (for program number 0102.C-0503(A)), and Gemini (for program number GS-2018B-LP-101). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research made use of observations from the LCOGT network, and the AAVSO Photometric All-Sky Survey (APASS), funded by the Robert Martin Ayers Sciences Fund and NSF AST-1412587. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013) ERC Grant Agreement n$\arcdeg$ 336480, and from the ARC grant for Concerted Research Actions, financed by the Wallonia-Brussels Federation. M.G. and E.J. are FNRS Senior Research Associates. Work by B.T.M. was performed in part under contract with the Jet Propulsion Laboratory (JPL) funded by NASA through the Sagan Fellowship Program executed by the NASA Exoplanet Science Institute. This work is partly supported by JSPS KAKENHI Grant Numbers JP18H01265 and 18H05439, and JST PRESTO Grant Number JPMJPR1775. KH acknowledges support from STFC grant ST/R000824/1.