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

TESS observed L 98-59 (TIC 307210830, TOI-175; R.A. = 08:18:07.62, decl. = −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 website⁶² and the Mikulski Archive for Space Telescopes (MAST) TESS alerts page⁶³ 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, respectively,⁶⁴ 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 photometry (SAP) and presearch data conditioned (PDCSAP) light curves (Smith et al. 2012; Stumpe et al. 2014) of L 98-59 are shown in Figure 1.

Zoom In Zoom Out Reset image size

We use the methods appropriate for M dwarfs previously used by Berta-Thompson et al. (2015), Dittmann et al. (2017), and Ment et al. (2019) 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 ± 0.014 M_⊙. We then use single star mass–radius relations (Boyajian et al. 2012) to find a stellar radius of 0.312 ± 0.014 R_⊙. We calculate the bolometric correction in K from Mann et al. (2015, erratum) to be 2.7 ± 0.036 mag, resulting in a bolometric luminosity for L 98-59 of 0.011 ± 0.0004 L_⊙. We calculate the correction in V from Pecaut & Mamajek (2013) to be −2.0 ± 0.03 mag,⁶⁵ resulting in a bolometric luminosity of 0.0115 ± 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 ± 0.0006 L_⊙ (i.e., 4.31e24 W). From the Stephan–Boltzmann law, we find an effective temperature T_eff = 3367 ± 150 K. As a comparison, we also used the relations in Mann et al. (2015) to determine an effective temperature of 3419 ± 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_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_eff and stellar metallicity [Fe/H], and we set the extinction A_V ≡ 0, due to the very close distance of the system. The resulting best fit is shown in Figure 2, with a reduced χ² of 3.8 for eight degrees of freedom. The best-fit parameters are T_eff = 3350 ± 100 K and [Fe/H] = −0.5 ± 0.5. Integrating the SED gives the bolometric flux at Earth as F_bol = 2.99 ± 0.18 × 10⁻⁹ erg s⁻¹ cm⁻². Finally, adopting the Gaia DR2 parallax and the correction of 80 μas from Stassun & Torres (2018), we calculate a stellar radius of 0.305 ± 0.018 R_⊙ using the Stefan–Boltzmann law. These are consistent with the adopted parameters discussed above.

Zoom In Zoom Out Reset image size

We also estimated a prior on the stellar density (ρ_*) by estimating the stellar mass. Here we used the empirical relations of Mann et al. (2015), which provide M_star ≈ 0.32M_⊙ 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 ∼3%; here we conservatively adopt an uncertainty of 10%. Together with the radius determined above from the SED and parallax, we obtain ρ_* = 15.9 ± 3.3 g cm⁻³ (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 Section 2.4) indicate that it is likely an old M dwarf with age >1 Gyr. The stellar parameters of L 98-59 are consistent with a spectral type of M3 ± 1.

We opted to create our own apertures from target pixel file data from Sectors 2, 5, and 8 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)⁶⁶ to extract light curves 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σ; the corresponding mask shape is shown in Figure 3. Thus produced, the resulting light curve 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 light curve by dividing by the median and subtracting one to center the flux about zero. We did the PLD detrending separately for each sector.

Zoom In Zoom Out Reset image size

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. 2019), 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₀) and $\mathrm{log}({\omega }_{0})$, where the power spectrum of the GP is

$$\begin{eqnarray}&&S(\omega )=\sqrt{\displaystyle \frac{2}{\pi }}\displaystyle \frac{{S}_{0}\,{\omega }_{0}^{4}}{{\left({\omega }^{2}-{\omega }_{0}^{2}\right)}^{2}+{\omega }^{2}\,{\omega }_{0}^{2}/{Q}^{2}},\end{eqnarray} \tag{ 1 }$$

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 $\mathrm{log}({S}_{0})$ and $\mathrm{log}({\omega }_{0})$ with the means of the log of the 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 its flexible nature and it smoothly varies (it is once mean square differentiable), enabling 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 $\mathrm{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 of the stellar density, and stellar radius for the three planets. Each individual planet was parameterized in terms of the log of the orbital period, time of first transit, the log of the 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 in addition required to be positive. The log of the mean stellar density, in cgs units, had a Gaussian prior with a mean of $\mathrm{log}15$ and standard deviation of 0.2 dex (as per Section 2.1). The limb darkening followed the Kipping (2013a) parameterization.

The log of the orbital periods, the time of first transits, and the log of the 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 α = 0.867 and β = 3.03 (as suggested by Kipping 2013b), and was bound between zero and one. The periastron angle at transit was sampled in vector space to avoid the sampler seeing a discontinuity at values of π.

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 four 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 the 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 with the 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_⊕. A chain of small terrestrial-size planets is common among M dwarfs (Muirhead et al. 2015), and L 98-59 is reminiscent of other systems such as TRAPPIST-1, Kepler-186, and Kepler-296 (Quintana et al. 2014; Barclay et al. 2015; Gillon et al. 2017). The stellar density obtained from the transit model is fully consistent with that determined from the stellar parameters in Section 2.1 (${15.8}_{-2.7}^{+2.6}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ for the former versus 15.9 ± 3.3 g cm⁻³ for the latter).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We repeated this analysis using the systematics-corrected light curves from the TESS pipeline (PDCSAP; Stumpe et al. 2014; Jenkins et al. 2016) 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 transit depths were lower in the PDCSAP data at the <1σ 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 ignoring 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 ± 0.8, the outer planet is a candidate Venus-zone planet (Kane et al. 2014).

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 minute 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σ 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 light curve 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)–(c) analysis, we used the Modelshift module of DAVE—an automated package designed to emphasize features in the light curve 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 (see e.g., Welsh et al. 2011), Modelshift first convolves the light curve with the transit model of the planet candidate. The module then computes the significance of the primary transits; odd–even differences; and secondary, tertiary, and positive features assuming white noise in the light curve; and compares the ratio between each of these and the systematic red noise F_red to the false alarm thresholds ${\mathrm{FA}}_{1}=\sqrt{2}\mathrm{erfcinv}({T}_{\mathrm{dur}}/({\text{}}P\times {\text{}}N))$ (assuming 20,000 objects evaluated) and ${\mathrm{FA}}_{2}=\sqrt{2}\mathrm{erfcinv}({T}_{\mathrm{dur}}/P)$ (for two events), where T_dur, P, and 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 Sec/F_red > FA₁. The Modelshift results are shown in Figures 6–8, where the panels show the phase-folded light curve (first row), the phase-folded light curve 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, and the most prominent secondary, tertiary, and positive features in the light curve (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 ("FA₁," "FA₂"), 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 nonzero 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–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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 a 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 light curves for nearby field stars. Our analysis revealed that a nearby field star ∼80'' NW of L 98-59 (2MASS 08175808-6817459, TIC 307210817, Tmag = 13.45, i.e., ≈4 mag fainter than L 98-59) is in fact an eclipsing binary (EB), manifesting both primary and secondary eclipses at a period of ≈10.43 days, with T₀ = 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 R.A. = 124.49217149500, decl. = −68.29612321000, ID = 5271055685541797120, and parallax = 0.1888 mas; and source 2 with R.A. = 124.49126472300, decl. = −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.

Zoom In Zoom Out Reset image size

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),⁶⁷ 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 subgroups (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 multiplanet 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'' for TESS versus 4'' for Kepler) and the point-spread function of TESS could be as large as 1', 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 by 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 25 away from L 98-59, we made use of the TFOP SG-1. 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 SG-1 data are the Las Cumbres Observatory (LCO) telescope network (Brown et al. 2013), SPECULOOS South Observatory (SSO; Burdanov et al. 2018; Delrez 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.

Table 4.  Observation Log

L 98-59 Planet	Date	Telescopea	Filter	ExpT	Exp	Dur.	Transit	Aperture	FWHM
	(UTC)			(s)	(N)	(minutes)	Coverage	(arcsec)	(arcsec)
b	2018 Oct 19	LCO-SSOb	r'	120	41	147	Ingr.+71%	3.89	2.21
	2018 Nov 11	SSO-Europa	r'	15	764	315	Full	5.25	4.53
	2018 Nov 11	LCO-CTIO-1	i'	25	120	154	Full	5.83	3.49
	2018 Nov 18	LCO-SAAO-0.4b	i'	70	125	175	Full	9.14	7.03
	2018 Nov 20	LCO-CTIO-1	r'	30	146	178	Full	3.89	2.20
	2018 Nov 29	LCO-CTIO-1	r'	30	186	223	Full	5.83	3.56
	2018 Dec 7	SSO-Io	r'	15	998	415	Full	5.60	4.06
	2019 Jan 26	LCO-SAAO-1b	r'	12	88	234	Full	3.89	3.08
c	2018 Oct 16	LCO-CTIO-1	i'	20	63	70	Ingr.+30%	4.27	2.01
	2018 Oct 22b	SSO-T17	clear	30	122	86	Ingr.+77%	7.10	2.40
	2018 Nov 11	SSO-Europa	r'	15	764	315	Full	5.25	4.53
	2018 Nov 22	MEarth	RG715	45	1682	380	Full	20.16	8.00
	2018 Dec 25	LCO-SSO-1	i'	22	108	113	Full	5.05	1.89
	2019 Jan 20	LCO-CTIO-1	g'	100	85	197	Full	6.22	2.75
	2019 Jan 20	LCO-CTIO-1	zs	30	170	197	Full	9.36	4.43
d	2018 Nov 7	LCO-CTIO-0.4	i'	70	119	170	Full	6.85	4.48
	2018 Nov 22	LCO-CTIO-1b	i'	25	108	113	Full	6.22	5.05
	2018 Nov 22	MEarth	RG715	45	1682	380	Full	20.16	8.00
	2019 Jan 13	LCO-CTIO-0.4	i'	14	143	73	Full	5.14	2.32
	2019 Jan 20	LCO-CTIO-1	r'	30	132	151	Full	6.22	2.89
	2019 Jan 28	LCO-CTIO-1	g'	50	153	223	Full	7.78	2.35

Notes.

^aTelescopes: LCO-SSO-1: Las Cumbres Observatory—Siding Spring (1.0 m), LCO-CTIO-1: Las Cumbres Observatory—Cerro Tololo Inter-American Observatory (1.0 m), LCO-CTIO-0.4: Las Cumbres Observatory—Cerro Tololo Inter-American 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 × 5 telescopes.) ^bObservations not shown in Figures 13–15 due to intrinsically high scatter in the light curve and/or because they were a deep exposure search for eclipsing binaries in nearby stars.

Download table as:  ASCIITypeset image

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–15.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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., ∼7 mag fainter than TIC 307210830, producing no detectable eclipses in the light curve of the latter), and we used deep exposures to confirm that it is not the origin of the L 98-59 b transits, providing a high level of confidence on 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.

Table 5.  Ground-based Follow-up R_p/R_* Less TESS R_p/R_* as a Function of Wavelength

TOI	Date Obs	Filter (Observatory)	Rp/R* Less TESS Rp/R*
L 98-59 b	UT 2018 Nov 20	Sloan r' (LCO)	$-{0.0008}_{-0.0085}^{+0.0054}$
L 98-59 b	UT 2018 Nov 29	Sloan r' (LCO)	${0.0018}_{-0.0045}^{+0.0039}$
L 98-59 c	UT 2019 Jan 20	Sloan g' (LCO)	−0.0033 ± 0.003
L 98-59 c	UT 2018 Nov 22	RG715 (MEarth)	−0.0007 ± 0.002
L 98-59 c	UT 2019 Jan 20	Sloan z' (LCO)	−0.0056 ± 0.0028
L 98-59 d	UT 2019 Jan 28	Sloan g' (LCO)	${0.0072}_{-0.0064}^{+0.0078}$
L 98-59 d	UT 2019 Jan 20	Sloan r' (LCO)	$-{0.003}_{-0.012}^{+0.01}$
L 98-59 d	UT 2018 Nov 7	Sloan i' (LCO)	${0.0072}_{-0.0072}^{+0.0088}$
L 98-59 d	UT 2018 Nov 22	RG715 (MEarth)	${0.0104}_{-0.0043}^{+0.0044}$

Download table as:  ASCIITypeset image

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 12⁶⁸ and on UT 2018 November 20 using the slicer mode with the CTIO HIgh ResolutiON (CHIRON) spectrograph (Tokovinin et al. 2013; R ≃ 80,000) on the Cerro Tololo Inter-American Observatory (CTIO)/ Small and Moderate Aperture Research Telescope System (SMARTS) 1.5 m telescope. CHIRON has a spectral range of 410–870 nm. We obtained one 7.5 minute exposure for the first epoch and 3 × 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 ± 0.1 km s⁻¹ using the TiO molecular bands at 7065–7165 Å.⁶⁹ Our analyses of the spectra reveal no evidence of double lines. We see negligible rotational broadening (v sin i = 0.0 ± 1.9 km s⁻¹) and do not see Hα 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 ± 52 m s⁻¹. For comparison, a 13 Jupiter-mass companion in a circular, nine-month period would induce a velocity semiamplitude of 863 m s⁻¹ on this star. This semiamplitude 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 ± 0.003 pc from the Gaia second data release is in agreement with the photometric distance estimate of 12.6 ± 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–Russell 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.

Zoom In Zoom Out Reset image size

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 μm band with a spectral resolution of R = 6000. The target was observed under favorable conditions, with an average seeing of ∼06. L 98-59 was observed twice in the ABBA nod patterns at 40 s integration time for each frame using the 06 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 ± 74 K, R_star = 0.37 ± 0.027R_⊙, and L = 0.021 ± 0.004L_⊙, 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 1 and AO images with VLT/NaCo on 2019 January 28. Both epochs of Gemini/DSSI data are collected simultaneously through the R- and I-band filters (692 nm and 880 nm respectively), while the VLT/NaCo data are collected in the Brγ filter. At 05 from the host, these data rule out companions 5.0, 6.6, and 5.8 mag fainter than the host star in the R, I, and Brγ bands, respectively. The 5σ contrast curves for each of these observations are presented in Figure 17.

Zoom In Zoom Out Reset image size

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 ∼6 mag 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 comoving companions beyond ∼150 mas (=1.6 au at the distance of this target). An object with ΔBrγ ∼ 6 mag and at the distance of L 98-59 would have a mass of ∼75M_J at an age of 10 Gyr (Baraffe et al. 2003),⁷⁰ 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 because 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.

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 50,000 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 hyperparameters from the hyperposteriors 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 ±1σ values. From the calculated radius values of 0.8 [0.05] R_⊕ (L 98-59 b), 1.35 [0.07] R_⊕ (L 98-59 c), and 1.57 [0.14] R_⊕ (L 98-59 d), we determined mass values of 0.5 [+0.3, −0.2] M_⊕, 2.4 [+1.8, −0.8] M_⊕, and 3.4 [+2.7, −1.4] M_⊕, 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 semiamplitude, 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.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 (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_⊕, 2.4 M_⊕, and 3.4 M_⊕ for L 98-58 b, L 98-58 c, and L 98-58 d, respectively), and best-fit+1σ masses (i.e., 0.58 M_⊕, 4.2 M_⊕, and 6.1 M_⊕ 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 nonsymplectic integrator (Rein & Spiegel 2015) with a time step of 0.01 times the orbit of the inner planet (i.e., about 30 minutes).

Our simulations show that for initially circular orbits, the semimajor 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σ planet masses (Figure 19). Thus, we consider orbits with nonnegligible 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).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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, nontransiting 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 2.5 M_⊕ and a 5.7 day orbital period (≈1.55 and ≈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.

Zoom In Zoom Out Reset image size

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 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 ∼0.1, the expected TTV amplitude for the innermost planet is ∼90 s, 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 ∼0.03 (as shown in Figure 21).

Zoom In Zoom Out Reset image size

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).

Zoom In Zoom Out Reset image size

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}_{{\rm{L}}98-59{\rm{c}}}}=2.1$ minutes and ${\sigma }_{{t}_{{\rm{L}}98-59{\rm{d}}}}\,=1.2$ minutes. The planets' masses are expected to be constrained with precisions of ${\sigma }_{m,{\rm{L}}98-59{\rm{c}}}=13.1\,{M}_{\oplus }$ and ${\sigma }_{m,{\rm{L}}98-59{\rm{d}}}=5\,{M}_{\oplus }$ with transit timing data through Sector 12.⁷¹

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 ∼20 g cm⁻³. In both scenarios, there is a dynamical model that fits the observed transits equally well. This is partly due to the large gaps in the data at present: the best opportunity for future TTV detection with TESS will be in upcoming data releases, when the system is observed for five consecutive sectors.

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 ratio 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 CO₂-dominated atmosphere could be detected for GJ 1132 with a modest number of JWST transits or eclipses (11 transits with NIRSpec/G235M or two 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), who calculated the expected signal-to-noise ratio 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 0.8, 1.4, and 1.0 times 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.

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 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 postrunaway 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, 2018; Angelo et al. 2017).

The L 98-59 planets receive significantly more energy than the Earth receives from the Sun (a factor of between 4 and 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 and 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.