TESS Hunt for Young and Maturing Exoplanets (THYME). IV. Three Small Planets Orbiting a 120 Myr Old Star in the Pisces–Eridanus Stream^*

After the initial discovery of the signal in the TESS data (Figure 1), we obtained follow-up transit photometry from the Las Cumbres Observatory (LCO; Brown et al. 2013), the Perth Exoplanet Survey Telescope (PEST), and the Spitzer Space Telescope⁴⁰ (Figure 2). Table 1 lists the photometric data modeled in Section 5.1; this subset of the available data provided the best opportunity for constraining the transit model. These and additional ground-based lightcurves⁴¹ ruled out eclipsing binaries on nearby stars, found the transit to be achromatic, and confirmed TOI 451 as the source of the three candidate event.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

2.1.1. TESS

2.1.2. Spitzer

2.1.3. Investigation into Spitzer Systematics

2.1.4. LCO

2.1.5. PEST

2.1.6. WASP-South

We obtained spectra with the Southern Astrophysical Research (SOAR) telescope/Goodman, South African Extremely Large Telescope (SALT)/High Resolution Spectrograph (HRS), Las Cumbres Observatory (LCO)/Network of Robotic Echelle Spectrographs (NRES), and Small and Moderate Aperture Research Telescope System (SMARTS)/CHIRON. The Goodman spectrum was used to fit the spectral energy distribution (SED; Section 3.1.1), and the SALT, LCO, and CHIRON spectra to measure radial velocities (RVs; Section 3.3). The S/N of the spectra used for RVs are given in Table 3.

Table 2.  Properties of the Host Star TOI 451

Parameter	Value	Source
Identifiers
TOI	451
TIC	257605131
TYC	7577-172-1
2MASS	J04115194-3756232
Gaia DR2	4844691297067063424
Astrometry
α	04 11 51.947	Gaia DR2
δ	−37 56 23.22	Gaia DR2
μα (mas yr−1)	−11.167 ± 0.039	Gaia DR2
μδ (mas yr−1)	12.374 ± 0.054	Gaia DR2
π (mas)	8.0527 ± 0.0250	Gaia DR2
Photometry
GGaia (mag)	10.7498 ± 0.0008	Gaia DR2
BPGaia (mag)	11.1474 ± 0.0027	Gaia DR2
RPGaia (mag)	10.2199 ± 0.0017	Gaia DR2
BT (mag)	11.797 ± 0.074	Tycho-2
VT (mag)	11.018 ± 0.064	Tycho-2
J (mag)	9.636 ± 0.024	2MASS
H (mag)	9.287 ± 0.022	2MASS
KS (mag)	9.190 ± 0.023	2MASS
W1 (mag)	9.137 ± 0.024	ALLWISE
W2 (mag)	9.173 ± 0.020	ALLWISE
W3 (mag)	9.117 ± 0.027	ALLWISE
W4 (mag)	8.632 ± 0.292	ALLWISE
Kinematics and Position
Barycentric RV (km s−1)	19.87 ± 0.12	This paper
Distance (pc)	123.74 ± 0.39	Bailer-Jones et al. (2018)
U (km s−1)	−10.92 ± 0.05	This paper
V (km s−1)	−4.18 ± 0.08	This paper
W (km s−1)	−18.81 ± 0.09	This paper
X (pc)	−41.56 ± 0.14	This paper
Y (pc)	−73.61 ± 0.24	This paper
Z (pc)	−90.43 ± 0.29	This paper
Physical Properties
Rotation period (days)	5.1 ± 0.1 days	This paper
$v\sin {i}_{* }$(km s−1)	7.9 ± 0.5 km s−1	This paper
i* (°)	${69}_{-8}^{11}\circ$	This paper
Fbol(erg cm−2 s−1)	(1.23 ± 0.07) × 10−8	This paper
Teff (K)	5550 ± 56	This paper
M⋆ (M⊙)	0.950 ± 0.020	This paper
R⋆ (R⊙)	0.879 ± 0.032	This paper
L⋆ (L⊙)	0.647 ± 0.032	This paper
ρ⋆ (ρ⊙)	1.4 ± 0.16	This paper
Age (Myr)	125 ± 8	Stauffer et al. (1998)a
	112 ± 5	Dahm (2015)a
	120	Curtis et al. (2019)
	134 ± 6.5	Röser & Schilbach (2020)
E(B − V) (mag)	${0.02}_{-0.01}^{+0.04}$	This paper

Note.

^aAge references denoted are ages for the Pleiades. i_* adopts the convention i_* < 90°.

Download table as:  ASCIITypeset image

2.2.1. SOAR/Goodman

2.2.2. SALT/HRS

We observed TOI 451 twice using NRES (Siverd et al. 2018) on the LCO system. NRES is a set of cross-dispersed echelle spectrographs connected to 1 m telescopes within the Las Cumbres network, providing a resolving power of R = 53,000 over the range 3800–8600 Å. We took both observations at the Cerro Tololo node, the first on 2019 March 27 and the second on 2019 July 31 (UT). The March observation consisted of two back-to-back exposures. The standard NRES pipeline⁴⁵ reduced, extracted, and wavelength-calibrated both observations.

We obtained a single spectrum with the CHIRON spectrograph (Tokovinin et al. 2013) on SMARTS, from which we measured the radial velocity and rotational broadening of the star. CHIRON is an R = 80,000 high-resolution, fiber-bundle-fed spectrograph on the 1.5 m SMARTS telescope, located at located at CTIO, Chile. Data reduction is described in Tokovinin et al. (2013).

To rule out unresolved companions that might impact our interpretation of the transit, we obtained speckle imaging using the High-Resolution Camera (HRCam) on the SOAR telescope. We searched for sources near TOI 451 in SOAR speckle imaging obtained on 17 March 2019 UT in the I band, a similar visible bandpass to TESS. Further details of observations from the SOAR TESS survey are available in Ziegler et al. 2020. We detected no nearby stars within 3'' of TOI 451 within the 5σ detection sensitivity of the observation (Figure 3).

Zoom In Zoom Out Reset image size

2.6.1. A Wide-binary Companion

2.6.2. Limits on Additional Companions

We gathered optical and near-infrared (NIR) photometry from the literature for use in our determination of the stellar parameters. Optical photometry comes from Gaia DR2 (Evans et al. 2018; Lindegren et al. 2018), AAVSO All-Sky Photometric Survey (APASS, Henden et al. 2012), and SkyMapper (Wolf et al. 2018). NIR photometry comes from the Two-Micron All-Sky Survey (2MASS, Skrutskie et al. 2006), and the Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010).

We summarize our derived stellar parameters in Table 2.

3.1.1. Luminosity, Effective Temperature, and Radius

To determine the L_*, T_eff, and R_* of TOI 451, we simultaneously fit its SED with the photometry listed in Table 2, our SOAR/Goodman optical spectrum, and Phoenix BT-Settl models (Allard et al. 2011). Significantly more detail of the method can be found in Mann et al. (2015) for nearby unreddened stars, with details on including interstellar extinction in Mann et al. (2016).

We compared the photometry to synthetic magnitudes computed from our SOAR spectrum. We used a Phoenix BT-Settl model (Allard et al. 2011) to cover gaps in the spectra and simultaneously fitting for the best-fitting Phoenix model and a reddening term (because reddening impacts both the spectrum and photometry). The Goodman spectrum is not as precisely flux calibrated as the data used in Mann et al. (2015), so we included two additional free parameters to fit out wavelength-dependent flux variations. The bolometric flux of TOI 451 is the integral of the unreddened spectrum. This flux and the Gaia DR2 distance yield an estimate of L_*.

We show the best-fit result in Figure 4 and adopted stellar parameters in Table 2. Our fitting resulted in two consistent radius estimates: the first from the Stefan–Boltzmann relation (with T_eff from the model grid) and the second from the ${R}_{* }^{2}/$(distance)² scaling (i.e., how much the BT-Settl model needs to be scaled to match the absolutely calibrated spectrum). The latter method is similar to the infrared-flux method (IRFM; Blackwell & Shallis 1977). Both measurements depend on a common parallax and observed spectrum, and hence are not completely independent. However, the good agreement (<1σ) was a useful confirmation of the final fit. Our derived parameters were T_eff = 5550 ± 56 K, L_* =0.647 ± 0.032L_⊙, and R_* = 0.879 ± 0.032R_⊙ from Stefan–Boltzmann and R_* = 0.863 ± 0.024R_⊙ from the IRFM. We adopt the former R_* for all analyses for consistency with previous work.

Zoom In Zoom Out Reset image size

The two reddest bands, W3 (12 μm) and W4 (22 μm), were both brighter than those derived from the best-fit template. We estimated an excess flux of 26% ± 7% at W3, and 70% ± 30% at W4 (3.7 and 2.3σ, respectively, assuming Gaussian errors). This suggests a cool ≲300 K debris disk. The W2 excess was not significant, both because of the large uncertainty in the W4 magnitude (8.63 ± 0.29) and because the SED analysis is sensitive to the choice of template at long wavelengths.

The frequency of infrared excesses decreases with age, declining from tens of percent at ages less than a few hundred megayears to a few percent in the field (Siegler et al. 2007; Meyer et al. 2008; Carpenter et al. 2009). In the similarly aged Pleiades cluster, Spitzer 24 μm excesses are seen in 10% of FGK stars (Gorlova et al. 2006). This excess emission suggests the presence of a debris disk, in which planetesimals are continuously ground into dust (see Hughes et al. 2018, for a review).

3.2.1. Mass

We used high-resolution optical spectra from SALT/HRS, NRES/LCO, and SMARTS/CHIRON to determine stellar radial velocities (RVs). We did not include Gaia because the RV zero-point has not been established in the same manner as our ground-based data.

We computed the spectral line broadening functions (BFs; Rucinski 1992; Tofflemire et al. 2019) through linear inversion of our spectra with a narrow-lined template. For the template, we used a synthetic PHOENIX model with T_eff  = 5400 K and $\mathrm{log}g$  = 4.5 (Husser et al. 2013). The BF accounts for the RV shift and line broadening. We computed the BF for each echelle order and combined them weighted by their S/N. We then fit a Gaussian profile to the combined BF to measure the RV. The RV uncertainty was determined from the standard deviation of the best-fit RV from three independent subsets of the echelle orders.

For HRS epochs, which consisted of three individual exposures, and the first NRES epoch, which consisted of two individual exposures, the RV and its uncertainty were determined from the error-weighted mean and standard error of the three individual spectra.

The resulting RVs are listed in Table 3. The RV zero points were calculated from the spectra obtained in this work and Rizzuto et al. (2020) and are based on telluric features. The zero points are 0.05 ± 0.10 for HRS (28 spectra), 0.32 ± 0.09 for NRES (11 spectra), and −0.05 ± 0.16 for CHIRON (1 spectrum). The S/N was assessed at  ∼6580 Å.

We used Linear Orbits for the Impatient via the python package lofti_gaiaDR2 (LOFTI; Pearce et al. 2020) to constrain the orbit of TOI 451 and TOI 451 B. Briefly, the lofti_gaiaDR2 retrieves observational constraints for the components from the Gaia archive and fits Keplerian orbital elements to the relative motion using the Orbits for the Impatient rejection sampling algorithm (Blunt et al. 2017). Unresolved binaries such as TOI 451 B can pose issues for LOFTI. While TOI 451 B's RUWE suggests binarity, it is still relatively low (1.2) and only on the edge of where astrometric accuracy compromises LOFTI (Pearce et al. 2020). We applied LOFTI to the system but caution that the unresolved binary may influence the results to an unknown, but likely small, degree. Our LOFTI fit constrained the orbit of TOI 451 and TOI 451 B to be close to edge on: we found an orbital inclination of i = 938 ± 116.

3.5.1. Projected Rotation Velocity

3.5.2. Rotation Period

In WASP-South, each of the six seasons of data shows clear modulation at 5.2 days, with variations in phase and the amplitude varying from 0.01 to 0.023 (Figure 5). The mean period from WASP-South is 5.20 days and the standard deviation is 0.02.

Zoom In Zoom Out Reset image size

We measured the stellar rotation period from the TESS data using GPs (e.g., Angus et al. 2018) as implemented in celerite (Foreman-Mackey et al. 2017). We used the same rotation kernel as in Newton et al. (2019), which is composed of a mixture of two stochastically driven, damped harmonic oscillators. The primary signal is an oscillator at the stellar rotation period P_*. The secondary signal is at half the rotation period. We also included a jitter term. The parameters we fit for are described in detail in Section 5.1, where the GP was fit simultaneously with the transits. We used the period from WASP-South to place a wide prior on the GP fit (see Section 5.1). The period we measured from the GP is $\mathrm{ln}{P}_{* }={1.635}_{-0.024}^{+0.027}$ or 5.1 ± 0.1 day.

3.5.3. Stellar Inclination

The original sample from Meingast et al. (2019) required RVs from Gaia for membership. Curtis et al. (2019) extended the sample to two dozen hotter stars by incorporating RVs from the literature. Recent searches have identified candidate Psc–Eri members without RVs. Röser & Schilbach (2020) adapted the convergent point method (van Leeuwen 2009) to the highly elongated structure of the stream. After placing distance and tangential velocity constraints on stars in the vicinity of the Meingast et al. (2019) sample, they identified 1387 probable stream members. Ratzenböck et al. (2020) identified around 2000 new members with a machine-learning classifier, trained on the originally identified sample of stream members. Röser & Schilbach (2020) calculated a bulk Galactic velocity for the Psc–Eri members identified by Meingast et al. (2019) of (U, V, W) = (−8.84, −4.06, −18.33) ± (2.2, 1.3, 1.7) km s⁻¹. This agrees with the value we have calculated for TOI 451 of (−10.92, −4.18, −18.81) km s⁻¹ within 1σ. Based on the space velocities and using the Bayesian membership selection of Rizzuto et al. (2011), we computed a Psc–Eri membership probability of 97% for TOI 451 and 84% for the companion TOI 451 B.

TOI 451 was included (as Gaia DR2 4844691297067063424) in the original membership list from Meingast et al. (2019) and in the subset with rotation periods from TESS data identified by Curtis et al. (2019). It also was listed as a member in Ratzenböck et al. (2020) and Röser & Schilbach (2020).

We used the high-resolution (R ∼ 46,000) spectra obtained with the SALT telescope to derive stellar abundances, as well as T_eff and $\mathrm{log}g$. The spectra cover ∼3700–8900 Å. We median-stacked the spectra to obtain a final spectrum with an S/N around 170 in the continuum at ∼5000 Å. We derived the T_eff, $\mathrm{log}g$, [Fe/H], and the microturbulent velocity (v_micro) using the Brussels Automatic Code for Characterizing High accUracy Spectra (BACCHUS; Masseron et al. 2016) following the method detailed in Hawkins et al. (2020a). To summarize, we set up BACCHUS using the atomic line list from the fifth version of the Gaia-ESO line list (Heiter et al. 2020, submitted) and molecular information for the following species were also included: CH (Masseron et al. 2014); CN, NH, OH, MgH and C₂ (T. Masseron 2020, private communication); and SiH (Kurucz line lists⁴⁷ ). We employed the MARCS model atmosphere grid (Gustafsson et al. 2008) and the TURBOSPECTRUM (Plez 2012) radiative transfer code. BACCHUS uses the standard Fe excitation–ionization balance technique to derive T_eff, $\mathrm{log}g$, and [Fe/H]. We refer the reader to Section 3 of Hawkins et al. (2020b) for a more detailed description of BACCHUS. The stellar atmospheric parameters derived from Fe excitation–ionization balance are T_eff  = 5556 ± 60, $\mathrm{log}g$  = 4.62 ± 0.17 dex, [Fe/H]  = −0.02 ± 0.08 dex. The T_eff and implied $\mathrm{log}g$ derived from a simultaneous fit of the star's SED are T_eff  = 5550 ± 56 and $\mathrm{log}g$  = 4.53 ± 0.04. Encouragingly, these values are, within the uncertainties, consistent with the physical properties outlined in Table 2, which were determined without high-resolution spectra.

Once the stellar atmospheric parameters were determined, we determined the abundance of Li at 6708 Å, A(Li). The presence (or absence) of large amounts of Li is an age indicator. Li fuses at the relatively low temperature of 2.5  × 10⁶ K. As a star ages, Li mixes downward into regions hotter than this temperature, where it is burned into heavier elements. Therefore, the abundance of Li decreases as the star ages. The amount of depletion varies with mass (or T_eff, given a main-sequence population). Therefore, A(Li) at a given T_eff constrains a star's age. This applies to both the Galactic disk (e.g., Ramírez et al. 2012) and open clusters (e.g., Boesgaard et al. 1998; Takeda et al. 2013; Bouvier et al. 2018; Martín et al. 2018).

To measure A(Li), we used the BACCHUS module abund. Using abund, we generated a set of synthetic spectra at 6708 Å with differing atmospheric abundances. We then used χ² minimization to find the synthetic spectrum that best fits the observed spectrum. We determined A(Li)  = 2.80 ± 0.10 dex.

We compare A(Li) and T_eff of TOI 451 to the observed trends in the Galactic disk (Ramírez et al. 2012), the Pleiades (Bouvier et al. 2018), and the Hyades (Takeda et al. 2013; Figure 6). The A(Li) for TOI 451 closely matches the A(Li) for the Psc–Eri stream, indicating that it is likely ∼120 Myr old and a stream member.

Zoom In Zoom Out Reset image size

As described in Section 3.5, we measured a rotation period of 5.1 ± 0.1 days, consistent with the 5.02 days reported in Curtis et al. (2019). As discussed in that work, this places TOI 451 on the slow sequence, the rotation–color sequence to which the initial distribution of periods converges for a uniform-age population. Figure 7 places TOI 451 in the context of other members of Psc–Eri and Pleiades cluster members.

Zoom In Zoom Out Reset image size

We extracted a lightcurve of TOI 451 B from the TESS 30 minute full-frame images (the lightcurve would be from the composite object if TOI 451 B is itself a binary). TOI 451 and its companion(s) are only separated by 37'', or about two TESS pixels, so the images of these two stars overlap substantially on the detector. The lightcurve of the companion TOI 451 B is clearly contaminated by the 14× brighter primary star. We therefore took a nonstandard approach to extracting a lightcurve for the companion. We started with the flux time series of the single pixel closest to the position of TOI 451 B during the TESS observations, clipped out exposures with flags indicating low-quality data points and between times 1419 < BJD − 2,457,000 < 1424 (when a heater on board the spacecraft was activated), and divided by the median flux value to normalize the lightcurve. We removed systematics and contaminating signals from TOI 451 by decorrelating the single-pixel lightcurve of TOI 451 B with mean and standard deviation quaternion time series (see Vanderburg et al. 2019), a fourth-order polynomial, and the flux from the pixel centered on TOI 451 (to model and remove any signals from the primary star). The resulting lightcurve, shown in Figure 8, shows several flares and a clear rotation signal. Before calculating rotation period metrics, we clipped out flares and removed points with times 1449 < BJD − 2,457,000 < 1454, which showed some residual systematic effects.

Zoom In Zoom Out Reset image size

We measured the companion's rotation period by calculating the Lomb–Scargle periodogram and autocorrelation function from the two-sector TESS lightcurve. The ACF and Lomb–Scargle periodogram both showed a clear detection of a 1.64 day rotation period, though we cannot completely rule out the possibility that TOI 451 B's rotation period is actually an integer multiple of this period.

Rebull et al. (2016), in their analysis of the Pleiades, detect periods for 92% of the members and suggest the remaining nondetections are due to nonastrophysical effects. We have suggested TOI 451 B is a binary, which we might expect to manifest as two periodicities in the lightcurve. We only detect one period in our lightcurve; however, a second signal could have been impacted by systematics removal or be present at a smaller amplitude than the 1.64 day signal, and so we do not interpret the lack of a second period further.

At around 100 Myr, stars of this type (early- to mid-M dwarfs) are in the midst of converging to the slow sequence. They may have a range of rotation periods, but are generally rotating with periods of a few days (Rebull et al. 2016). While the Psc–Eri members studied in Curtis et al. (2019) do not extend to stars as low mass as TOI 451 B, we can consider the similarly aged Pleiades members as a proxy. Figure 7 demonstrates that the 1.64 day rotation period of TOI 451 B is typical for stars of its color at 120 Myr.

Chromospheric and coronal activity depend on stellar rotation and are thus also an age indicator (e.g., Skumanich 1972). We use excess UV emission as an age diagnostic for TOI 451 and TOI 451 B, considering the flux ratio F_NUV/F_J as a function of spectral type. The use of this ratio for this purpose was suggested by Shkolnik et al. (2011). More details on this general technique and applications to the Psc–Eri stream are provided in the Appendix.

TOI 451 and TOI 451 B both were detected by GALEX during short NUV exposures taken for the All-Sky Imaging Survey (AIS; Bianchi et al. 2017). Both stars were also observed with longer NUV exposures for the Medium-depth Imaging Survey (MIS), but only the primary's brightness was reported (m_NUV,p = 16.674 ± 0.006 mag). Using the MIS data, we measured the secondary's brightness to be m_NUV,s = 21.33 ± 0.06 mag (see the Appendix for details).

In Figure 9, we plot F_NUV/F_J versus spectral type for TOI 451 and TOI 451 B. Because both spectral type and color should be unaffected by a near-equal-mass unresolved companion, the fact that TOI 451 B is a binary is not expected to impact this analysis. We also show isochronal sequences for the other members of Psc–Eri and the similarly aged Pleiades cluster, using the latter to define the M-dwarf regime at this age. To illustrate the expected fluxes for older stars, we show the sequence for the Hyades cluster. Details on the derivations of these isochronal sequences are given in the Appendix.

Zoom In Zoom Out Reset image size

While solar-type stars only have excess NUV and X-ray fluxes for a short time, the NUV flux of early-M dwarfs remains saturated to ages  ≲ 300 Myr before sharply declining (Shkolnik & Barman 2014). TOI 451 is consistent with the sequences for other G dwarfs in all three clusters: as expected, its UV excess is not a highly discriminating age diagnostic. However, the wide companion TOI 451 B sits above the Hyades sequence and shows an NUV flux excess that is broadly consistent with the Pleiades sequence that is similar in age to Psc–Eri. This supports the membership of the TOI 451 system in the Psc–Eri stream.

We modeled the transit with misttborn (Mann et al. 2016; Johnson et al. 2018). This routine uses batman (Kreidberg 2015) to produce the transit model of Mandel & Agol (2002) and explores the posterior with the MCMC sampler emcee (Foreman-Mackey et al. 2013).

To model the transits of each of the three planets, we fitted for the planet-to-star radius ratio R_P/R_*, impact parameter b, period P, and the epoch of the transit midpoint T₀. We assumed R_P/R_* is the same in all filters. We also fitted for the mean stellar density (ρ_*/ρ_⊙). We used a quadratic limb-darkening law described by g_1,f and g_2,f for each filter f. This is a reasonable choice given that the host star is Sun-like (Espinoza & Jordán 2016). We fitted the limb-darkening parameters using the Kipping (2013) parameterization (q_1,f, q_2,f). For our first fit, we fixed the eccentricity e to 0. For our second fit, we allowed e and the argument of periastron ω to vary, using the parameterization $\sqrt{e}\sin \omega$ and $\sqrt{e}\cos \omega$ (Ford 2006; Eastman et al. 2013).

We used uniform priors for the planetary parameters, and a Gaussian prior on ρ_*/ρ_⊙ centered at 1.34 (following Table 2) with a 1σ width of 0.5. We placed Gaussian priors on the limb-darkening parameters based on theoretical values for a star with the temperature and radius given in Table 2. We assume solar metallicity based on the mean iron abundance of the Psc–Eri members studied by Hawkins et al. (2020a), who found [Fe/H] = − 0.03 dex with a dispersion of 0.04–0.07 dex. For the Gaussian means of the limb-darkening priors, we use the coefficients calculated with the Limb Darkening Toolkit (Parviainen & Aigrain 2015, ; LDTK), using the filter transmission curves available for LCO z_s⁴⁸ and TESS.⁴⁹ The PEST bandpass is similar to that of MEarth, and following Dittmann et al. (2017), we adopted the filter profile from Dittmann et al. (2016). For Spitzer, we used the tabulated results from Claret & Bloemen (2011). Due to systematic uncertainties in limb-darkening parameters (Müller et al. 2013; Espinoza & Jordán 2015), the potential impact of spots (Csizmadia et al. 2013), and differences between LDTK and Claret coefficients, we used 0.1 for the 1σ width of the priors on q.

For stellar rotation, we used the GP model of a mixture of simple harmonic oscillators introduced in Section 3.5, with the parameters sampled in log space (in the following, logarithms are all natural logarithms). The model includes the power at the rotation period P_* (the primary signal) and at P_*/2 (the secondary signal). We fitted for the period of the primary signal $\mathrm{ln}{P}_{* }$, the amplitude of the primary signal $\mathrm{ln}{A}_{1}$, the relative amplitudes of the primary and secondary signals m (where A₁/A₂ = 1 + e^−m), the decay timescale (or "quality factor") of the secondary signal $\mathrm{ln}{Q}_{2}$, and the difference in quality factors $\mathrm{ln}{\rm{\Delta }}Q$ (where $\mathrm{ln}{\rm{\Delta }}Q=\mathrm{ln}\left({Q}_{1}-{Q}_{2}\right)$). We additionally included a photometric jitter term, σ_GP.

We placed a Gaussian prior on $\mathrm{ln}{P}_{* }$ centered at 1.6487, based on the period from the WASP data, with a width of 0.05. We use log-normal priors on the remaining parameters. We require $\mathrm{ln}{\rm{\Delta }}Q\gt 0$ to ensure that the primary signal has a higher quality than the secondary, and $\mathrm{ln}{Q}_{2}\gt \mathrm{ln}0.5$ because we are modeling a signal with periodic behavior (see Equation (23) and Figure 1 in Foreman-Mackey et al. 2017).

The GP model was applied to data from TESS, LCO, and PEST, but not data from Spitzer. The latter was assumed to have no out-of-transit variability remaining after the corrections described in Section 2.

The autocorrelation length of our e = 0 fit was 470 steps and for our variable-eccentricity fit was 1200 steps. We ran the MCMC chain for 100 times the autocorrelation times with 100 walkers, discarding the first half as burn in. Figures 1 and 2 shows the best-fitting models overlain on the data, and Table 4 lists the median and 68% confidence limits of the fitted planetary and stellar parameters.

Table 4.  Transit Fitting Results for the TOI 451 System

Parameter	Planet b	Planet c	Planet d	Planet b	Planet c	Planet d
	e, ω fixed			e, ω free
Fitted Transit Parameters
T0 (BJD)	${1410.9900}_{-0.0037}^{+0.0056}$	${1411.7956}_{-0.0026}^{+0.0048}$	${1416.63478}_{-0.00092}^{+0.00088}$	${1410.9909}_{-0.0042}^{+0.0046}$	${1411.7961}_{-0.0030}^{+0.0039}$	${1416.63499}_{-0.00093}^{+0.00097}$
P (days)	${1.858703}_{-3.5\times {10}^{-5}}^{+2.5\times {10}^{-5}}$	${9.192522}_{-10\times {10}^{-5}}^{+6.0\times {10}^{-5}}$	16.364988 ± 4.4 × 10−5	${1.858701}_{-3.3\times {10}^{-5}}^{+2.7\times {10}^{-5}}$	${9.192523}_{-8.4\times {10}^{-5}}^{+6.4\times {10}^{-5}}$	${16.364981}_{-4.9\times {10}^{-5}}^{+4.7\times {10}^{-5}}$
RP/R⋆	${0.0199}_{-0.0011}^{+0.0010}$	${0.03237}_{-0.00070}^{+0.00065}$	0.04246 ± 0.00044	${0.0203}_{-0.0011}^{+0.0014}$	${0.03206}_{-0.00085}^{+0.00090}$	${0.04205}_{-0.00045}^{+0.00050}$
b	${0.22}_{-0.15}^{+0.20}$	${0.139}_{-0.096}^{+0.121}$	${0.387}_{-0.039}^{+0.047}$	${0.42}_{-0.28}^{+0.32}$	${0.52}_{-0.32}^{+0.23}$	${0.23}_{-0.16}^{+0.18}$
ρ⋆ (ρ⊙)	${1.294}_{-0.088}^{+0.061}$			${1.41}_{-0.16}^{+0.15}$
q1,1	${0.344}_{-0.078}^{+0.080}$			${0.385}_{-0.084}^{+0.083}$
q2,1	${0.394}_{-0.082}^{+0.068}$			${0.391}_{-0.088}^{+0.070}$
q1,2	${0.325}_{-0.096}^{+0.095}$			${0.325}_{-0.086}^{+0.099}$
q2,2	${0.365}_{-0.096}^{+0.080}$			${0.364}_{-0.100}^{+0.082}$
q1,3	${0.058}_{-0.045}^{+0.076}$			${0.058}_{-0.045}^{+0.080}$
q2,3	${0.341}_{-0.094}^{+0.085}$			${0.324}_{-0.092}^{+0.089}$
q1,4	${0.036}_{-0.025}^{+0.042}$			${0.027}_{-0.020}^{+0.038}$
q2,4	${0.155}_{-0.086}^{+0.094}$			${0.159}_{-0.090}^{+0.095}$
$\sqrt{e}\sin \omega$	⋯			$-{0.23}_{-0.21}^{+0.19}$	$-{0.26}_{-0.19}^{+0.16}$	${0.02}_{-0.13}^{+0.11}$
$\sqrt{e}\cos \omega$	⋯			$-{0.09}_{-0.38}^{+0.42}$	${0.01}_{-0.43}^{+0.36}$	$-{0.01}_{-0.29}^{+0.32}$
Fitted Gaussian Process Parameters
$\mathrm{log}{P}_{\mathrm{GP}}$ (day)	${1.636}_{-0.024}^{+0.027}$			${1.637}_{-0.024}^{+0.027}$
$\mathrm{log}{A}_{\mathrm{GP}}$ (%2)	$-{12.08}_{-0.26}^{+0.36}$			$-{12.1}_{-0.25}^{+0.32}$
$\mathrm{log}Q{1}_{\mathrm{GP}}$	${1.62}_{-0.94}^{+0.93}$			${1.59}_{-0.95}^{+0.89}$
$\mathrm{log}Q{2}_{\mathrm{GP}}$	${1.25}_{-0.29}^{+0.33}$			${1.25}_{-0.30}^{+0.34}$
Mix Q1,Q2	4.6  ±  3.5			${4.8}_{-3.2}^{+3.6}$
σGP (%)	$-{15.0}_{-3.4}^{+3.1}$			$-{15.0}_{-3.3}^{+3.1}$
Derived Transit Parameters
a/R⋆	${6.93}_{-0.16}^{+0.11}$	${20.12}_{-0.47}^{+0.31}$	${29.56}_{-0.69}^{+0.46}$	${6.63}_{-0.89}^{+0.51}$	${18.9}_{-2.2}^{+1.4}$	${30.69}_{-1.0}^{+0.82}$
i(°)	${88.2}_{-1.7}^{+1.2}$	${89.61}_{-0.36}^{+0.27}$	${89.25}_{-0.1}^{+0.084}$	${86.5}_{-2.9}^{+2.3}$	${88.49}_{-0.67}^{+0.95}$	${89.56}_{-0.35}^{+0.31}$
δ (%)	${0.0396}_{-0.0042}^{+0.0041}$	${0.1048}_{-0.0045}^{+0.0043}$	0.1803 ± 0.0037	${0.0412}_{-0.0043}^{+0.0059}$	${0.1028}_{-0.0054}^{+0.0059}$	${0.1768}_{-0.0038}^{+0.0042}$
T14 (days)	${0.0849}_{-0.0054}^{+0.0024}$	0.1483 ± 0.0016	0.1707 ± 0.001	${0.082}_{-0.017}^{+0.02}$	${0.137}_{-0.031}^{+0.028}$	${0.171}_{-0.017}^{+0.018}$
T23 (days)	${0.0815}_{-0.0058}^{+0.0024}$	${0.1387}_{-0.0016}^{+0.0015}$	0.1543 ± 0.0011	${0.078}_{-0.018}^{+0.019}$	${0.125}_{-0.032}^{+0.026}$	${0.156}_{-0.017}^{+0.016}$
FPP parameter	${0.119}_{-0.094}^{+0.308}$	${0.044}_{-0.038}^{+0.11}$	${0.354}_{-0.072}^{+0.101}$	${0.44}_{-0.38}^{+1.32}$	${0.7}_{-0.61}^{+1.12}$	${0.12}_{-0.11}^{+0.28}$
Tperi (BJD)	${1410.99}_{-0.0037}^{+0.0056}$	${1411.7956}_{-0.0026}^{+0.0048}$	${1416.63478}_{-0.00092}^{+0.00088}$	${1411.26}_{-0.89}^{+0.39}$	${1410.4}_{-2.1}^{+4.6}$	${1416.8}_{-3.8}^{+3.3}$
g1,1	${0.452}_{-0.100}^{+0.093}$			${0.479}_{-0.100}^{+0.094}$
g2,1	${0.122}_{-0.078}^{+0.100}$			${0.132}_{-0.086}^{+0.113}$
g1,2	0.40 ± 0.12			${0.41}_{-0.13}^{+0.12}$
g2,2	${0.148}_{-0.089}^{+0.116}$			${0.154}_{-0.095}^{+0.111}$
g1,3	${0.151}_{-0.080}^{+0.093}$			${0.145}_{-0.077}^{+0.090}$
g2,3	${0.067}_{-0.047}^{+0.082}$			${0.075}_{-0.050}^{+0.083}$
g1,4	${0.052}_{-0.033}^{+0.048}$			${0.046}_{-0.029}^{+0.045}$
g2,4	${0.122}_{-0.061}^{+0.081}$			${0.106}_{-0.057}^{+0.080}$
(RPR⊕)	1.91 ± 0.12	3.1 ± 0.13	4.07 ± 0.15	${1.94}_{-0.13}^{+0.15}$	3.07 ±  0.14	4.03 ± 0.15
a (au)	${0.0283}_{-0.0012}^{+0.0011}$	${0.0823}_{-0.0036}^{+0.0033}$	${0.1208}_{-0.0052}^{+0.0048}$	${0.0271}_{-0.0038}^{+0.0023}$	${0.0771}_{-0.0093}^{+0.0066}$	${0.1255}_{-0.0065}^{+0.0057}$
Teq (K)	${1491}_{-19}^{+23}$	${875}_{-11}^{+13}$	${722}_{-9}^{+11}$	${1524}_{-60}^{+100}$	${903}_{-36}^{+53}$	${708}_{-12}^{+15}$
e	⋯			${0.19}_{-0.14}^{+0.20}$	${0.2}_{-0.14}^{+0.18}$	${0.057}_{-0.040}^{+0.133}$
ω(°)	⋯			${238}_{-48}^{+85}$	266 ± 63	${170}_{-120}^{+170}$

Note.

T_eq assumes an albedo of 0 (the planets reflect no light) with no uncertainty. The planetary radii listed above are from the joint fit of TESS, ground-based, and Spitzer data. For the TESS-only fit, R_P/R_⋆ for planets b, c, and d, respectively are: 0.0194 ± 0.0011, 0.0286 ± 0.0016, and ${0.0418}_{-0.0013}^{+0.0012}$. The largest difference is for c. For the Spitzer-only fit of c, R_P/R_⋆ = 0.034 ± 0.001. We suggest adopting the e = 0 fit because the eccentricities from the variable-eccentricity fit are consistent with 0.

Download table as:  ASCIITypeset image

We performed fits of the TESS data using the EXOFASTv2⁵⁰ (Eastman et al. 2019) software package, which simultaneously fits the transits with the stellar parameters. We first removed the stellar variability using the GP model described in Section 3.5, masking the transit prior to fitting. The stellar parameters were constrained by the SED and MESA Isochrones and Stellar Tracks (MIST) stellar evolution models (Dotter 2016; Choi et al. 2016). We enforced Gaussian priors on T_eff, [Fe/H], and $\mathrm{log}g$ following Section 3, and on the Gaia DR2 parallax corrected for a systematic offset (Stassun & Torres 2018). We used a Gaussian prior on the age of 125 ± 15 Myr and restricted extinction to  <0.039 (Schlafly & Finkbeiner 2011). The stellar parameters from EXOFASTv2 are T_eff = 5563 ± 44 K, M_* = 0.94 ± 0.025, and R_* =0.834 ± 0.01. We generally found excellent agreement with the best-fitting planetary parameters from misttborn.

We also fitted the TESS and Spitzer data sets independently. There is a 2σ discrepancy between the Spitzer transit depths of TOI 451 c (R_P/R_⋆ = 0.034 ±  0.001) and the TESS transit depth (R_P/R_⋆ = 0.029 ±  0.002). No significant differences are seen in the transit depths of TOI 451 d, which we might expect if there was dilution in the TESS data. The mildly larger planetary radius for TOI 451 c measured at 4.5 μm relative to that measured at TESS's red-optical bandpass could be due to occulted starspots,⁵¹ unocculted plages, or atmospheric features. Given that the transits of c and d were observed at similar times, we might also expect differences in the depths of TOI 451 d if starspots caused the discrepancy for c; however, the nature of starspots (size, temperature, location/active latitudes, longevity) is not well understood. Alternatively, the difference in transit depths for TOI 451 c could derive from its atmosphere. It is likely, however, that the difference arises from systematics in the Spitzer data, as discussed in Section 2.1.2.

We considered four different false-positive scenarios for each of the candidate planets: (1) there is uncorrected stellar variability, (2) TOI 451 is an eclipsing binary, (3) TOI 451 is a hierarchical eclipsing system, and (4) there is a background or foreground eclipsing system. Relevant to the blend scenarios, transits are visible in the Spitzer data even when shrinking or shifting the aperture, indicating that the signal lands within 2 Spitzer pixels (24) of TOI 451. Speckle imaging and Gaia RUWE rule out companions brighter than these limits (i.e., bright enough to cause the transits) down to 02, so a very close blend is required.

We ruled out stellar variability for all three planets from the Spitzer transits. For any spot contrast, stellar variation due to rotation and spots/plages will always be weaker at Spitzer wavelengths compared to TESS. For TOI 451, out-of-transit data taken by Spitzer show a factor of ≃4 lower variability than the equivalent baseline in TESS data. Thus, if any transit was due purely to uncorrected stellar variation, the shape, duration, and depth would be significantly different or the entire transit would not be present in the Spitzer photometry.

We used the source brightness parameter from Vanderburg et al. (2019) to constrain the magnitude of a putative blended source (bound or otherwise). The parameter, Δm, relates the ingress or egress duration to transit duration and reflects the true radius ratio, independent of whether there is contaminating flux: ${\rm{\Delta }}m\leqslant 2.5{\mathrm{log}}_{10}({T}_{12}^{2}/{T}_{13}^{2}/\delta )$. Here, δ is the transit depth, T₁₂ is the ingress/egress duration, and T₁₃ the time between the first and third contacts. We calculated Δm for the posterior samples for our variable-eccentricity transit fit and took the 99.7% confidence limit. We find Δm < 4.6, 2.0, 1.0 for TOI 451 b, c, and d, respectively.

Using the brightness constraints from our imaging data and the Δm parameter, we statistically rejected the scenario where any of the transits are due to an unassociated field star, either an eclipsing binary or an unassociated transiting planetary system (Morton & Johnson 2011). We drew information on every star within 1° of TOI 451 from Gaia DR2 satisfying the most conservative brightness limits (ΔT < 4.6). This yielded a source density of  ≃3500 stars per square degree, suggesting a negligible  ≃7 × 10⁻⁵ field stars that were missed by our speckle data and still bright enough to produce the candidate transit associated with b. Constraints are stronger for the other planets.

To investigate scenarios including a bound companion blended with the source, we first considered the constraint placed by the multiwavelength transit depths. As explained in Désert et al. (2015), if the transit signals were associated with another star in the aperture, the transit depth observed by Spitzer would be deeper than that observed by TESS, owing to the decreased contrast ratio between the target and the blended star (which must be fainter as is assumed to be cooler) at 4.5 μm compared to 0.75 μm. Thus, the ratio of the Spitzer-to-TESS transit depths (δS/δT) provides a range of possible TESS–Spitzer colors of the putative companion (C_TS,comp) in terms of the combined (unresolved) color (C_TS,combined). Following Tofflemire et al. (2021), we used the 95th percentile range for δS/δT for each transit to derive the putative companion color:

$$\begin{eqnarray}&&{C}_{\mathrm{TS},\mathrm{comp}}\lt {C}_{\mathrm{TS},\mathrm{combined}}+2.5{\mathrm{log}}_{10}\left(\tfrac{\delta S}{\delta T}\right).\end{eqnarray} \tag{ 1 }$$

Adopting the weakest constraints from the three planets gives a C_TS,comp between 1.09 and 1.30 mag, indicating that the putative bound companion must be similar in color to the target star.

We then simulated binary systems following M. Wood et al. (2020, in preparation), which we compared to the observational data and constraints. In short, we generated 5 million binaries following the period and mass ratio from Raghavan et al. (2010) and the eccentricity distribution from Price-Whelan et al. (2020). For each binary, we calculated the expected radial velocity curve, magnitude, and projected separation at the epoch of the speckle data. We then compared the generated models for a given simulated binary to our radial velocities, speckle data, and Gaia imaging and astrometry (Section 2.6.2). To account for stellar variability, we added a 50 m s⁻¹ error to the velocity measurements. Binaries were then rejected based on the probability of the observational constraints being consistent with the binary star parameters by chance (color, source brightness, radial velocities, and all imaging/photometric data). We ran two versions of this simulation, one where a single companion was forced to eclipse the primary (for ruling out eclipsing binaries), and one where the binary's orbital inclination was unrestricted (for hierarchical systems). The former set ruled out all nonplanetary signals at periods matching any of the three planets (Figure 10). In the latter simulation, less than 0.01% of the simulated blends passed all observational constraints for the signals associated with any of the three planets.

Zoom In Zoom Out Reset image size

We now take into consideration the probability that any given star happens to be an eclipsing binary (about 1%; Kirk et al. 2016), compared to the probability of a transiting planet (also about 1%; Thompson et al. 2018). The chance alignment of unassociated field stars is therefore 7 × 10⁻⁵ times as likely as the transiting planet scenario, and the hierarchical triple 1 × 10⁻⁴ times as likely. The probability of the host star itself being an eclipsing binary is negligible. The transiting planet hypothesis is significantly more likely than the nontransiting planet hypothesis, with a false-positive probability of about 2 × 10⁻⁴. The probability of a star hosting a planet discovered by TESS is likely lower (about 0.05% based on current search statistics; Guerrero et al. 2021) than assumed here. On the other hand, studies of the Kepler multiplanet systems indicate that the false-positive rate for systems with two or more transiting planets is low (Lissauer et al. 2012, 2014), about a factor of 15 less for TESS planets (Guerrero et al. 2021), which would more than compensate for our assumption. The false-positive probability of 2 × 10⁻⁴ is therefore an overestimate.

We have ruled out instrumental variability, chance alignment of unassociated field stars (including both eclipsing binaries and planetary systems), and bound eclipsing binaries or hierarchical triples, and therefore reject the false positives identified at the beginning of this section. We consider all three planets validated at high confidence.