TOI 560: Two Transiting Planets Orbiting a K Dwarf Validated with iSHELL, PFS, and HIRES RVs

TOI 560 (TIC 101011575; HD 75383; GAIA EDR3 5746824674801810816) was observed first in TESS Sector 8 from UT 2019 February 2 to UT 2019 February 28, then again in Sector 34 during the TESS extended mission from UT 2021 January 13 to UT 2019 February 9. Stellar properties are listed in Table 1. TOI 560 is in the Hydra constellation and is relatively nearby (31.6 pc) and bright (V = 9.67 mag) making it an ideal candidate for study by TESS. The data collection pipeline was developed by the TESS Science Processing Operations Center (SPOC; Jenkins et al. 2016) and uses a wavelet-based matched filter (Jenkins 2002; Jenkins et al. 2010, 2020). One transit signal was detected during the Sector 8 observations, and fitted with a limb-darkened transit model (Li et al. 2019) with various vetting and diagnostic tests (Twicken et al. 2018) and labeled TOI 560 b (Guerrero et al. 2021), which we hereafter refer to as TOI 560 b. This prompted us to explore the target further with RV measurements (Section 2.4) to verify the TESS candidate and perhaps uncover additional companions in the system. With the release of the Sector 34 light curve, we independently identified by eye a second transiting candidate in the TOI 560 system, which was vetted and adopted as a TOI by the TESS mission as TOI 560 c, which hereafter we refer to as TOI 560 c.

In this work, we specifically analyze the detrended Presearch Data Conditioning Simple Aperture Photometry (PDC-SAP) light curves for Sectors 8 and 34 (Smith et al. 2012; Stumpe et al. 2012, 2014) obtained from the Mikulski Archive for Space Telescopes (MAST). ⁷¹ After gathering this data, we normalize the detrended flux in electrons per second so that the median for each sector is unity. In Figure 1, we show the TESS target pixel files (TPF) around the target star in Sectors 8 and 34, where orange outlines show the pixels used to create the TESS light curve. There are other point sources (notated by circular red points) that are located within the TESS apertures, so these needed to be subtracted from the SAP flux when creating the PDC-SAP flux.

Zoom In Zoom Out Reset image size

TOI 560 b was observed with Spitzer on 2019 August 20, using Director's Discretionary Time (Crossfield et al. 2018). A single transit was observed using the 4.5 μm channel (IRAC2; Fazio et al. 2004) in subarray mode with an integration time of 0.36 s. The transit observation spanned 6 hr 15 min totaling 823 frames with short observations taken before and after transit to check for bad pixels. Peak-up mode was used to place the star as close as possible to the well-characterized "sweet spot" of the detector.

To extract photometry from the Spitzer observations, we use the Photometry for Orbits Eclipses and Transits (POET ⁷² ) package (Cubillos et al. 2013; May & Stevenson 2020). In summary, POET creates a bad pixel mask and discards bad pixels based on the Spitzer Basic Calibrated Data. Outlier pixels are also discarded using sigma-rejection. Then, the center of the point-spread function (PSF) is determined using a 2D Gaussian fitting technique. After the center of the PSF is found, the light curve is extracted using aperture photometry in combination with a BiLinearly Interpolated Subpixel Sensitivity (BLISS) map described in Stevenson et al. (2012). The resulting data are then simultaneously fit with a model that accounts for both the light curve itself and a temporal ramp-like trend attributed to "charge trapping." The posterior distribution is sampled using a Markov Chain Monte Carlo (MCMC) algorithm with chains initialized at the best-fit values.

Aperture photometry was performed with various aperture sizes (ranging from 2–6 pixels in increments of 1 pixel). Visual inspection of the raw data indicated the temporal ramp model was likely either linear or constant. The optimal aperture size was found to be 3 pixels as this size returned the lowest standard deviation of the normalized residuals (SDNRs) for both ramp models. We then tested bin sizes of 0.1, 0.03, 0.01, and 0.003 pixels square for the BLISS map, finding that a bin size of 0.01 minimized the SDNR. Both the linear ramp and constant models gave similar SDNRs for these tests, so we chose the constant model, the simplest of the two.

In 2019 and 2020, seven partial or full transit observations of TOI 560 b were collected with five observatory sites from the Las Cumbres Observatory Global network of Telescopes (LCOGT; Brown et al. 2013). These observations are summarized in Table 2. Data were collected in either the B or Sloan $z^{\prime}$ filters with exposure times of 15 or 30 s to check for transit-depth chromaticity such as would be produced by a false-positive eclipsing binary scenario. For one transit (UT 2019/04/27), data were collected simultaneously in both filters. For a second transit (UT 2020/01/27), data was collected simultaneously in the same filter from two sites. The other five transits were observed in a single band at a single telescope/site. Data were analyzed using AstroImageJ, and the data reduction was done using the BANZAI LCOGT facility pipeline (McCully et al. 2018).

The Perth Exoplanet Survey Telescope (PEST) near Perth, Australia is a 0.3 m telescope equipped with a 1530 × 1020 SBIG ST-8XME camera with an image scale of 12/pixel, resulting in a 31'' × 21'' field of view. Observations of a transit of TOI 560 b were obtained on UT 2020 January 14 in the Rc filter with an exposure time of 30 s. A custom pipeline based on C-Munipack was used to calibrate the images and extract differential photometry.

On UT 2019 December 7, the TOI 560 system was observed by the Next Generation Transit Survey (NGTS; Wheatley et al. 2018) using four telescopes (see Bryant et al. 2020), consisting of 0.2 m diameter robotic telescopes, located at ESO's Paranal Observatory, Chile. A custom NGTS filter in the wavelength range 520–890 nm was used. A transit was detected and confirmed for planet b. The NGTS data were reduced using a custom aperture photometry pipeline version 2, which performs source extraction and photometry using the SEP Python library (Bertin & Arnouts 1996; Barbary 2016) and is detailed in Bryant et al. (2020). The pipeline uses Gaia DR2 (Barbary 2016; Gaia Collaboration et al. 2018) to automatically identify comparison stars that are similar in brightness, color, and CCD position to TOI 560.

The SuperWASP team monitored TOI 560 for four seasons from 2009–2012 included as part of a larger all-sky survey (Pollacco et al. 2006). Data were reduced and light curves were generated using the standard SuperWASP analyses (Maxted et al. 2011).

Reconnaissance spectroscopy for TOI 560 was obtained with the Las Cumbres Observatory (LCOGT) Network of Robotic Echelle Spectrographs (NRES). NRES is a fiber-fed echelle spectrograph at several sites mounted on the respective LCOGT 1 m telescopes with wavelength range 380–860 nm and spectral resolution R ∼ 53,000. Observations were executed at the Sutherland Observatory, South Africa, and the McDonald Observatory, USA, on the nights 2019 November 4, 2019 October 29, and 2019 May 12. We took dark, bias, flat, and arc lamp calibration images at the beginning of each night. The target was centered on one fiber, and the second fiber was used to observe a Th-Ar calibration lamp simultaneously. On each night, we took three consecutive exposures with 480 s on each individual spectrum. The S/Ns for each night were 33, 30, and 33, respectively. We obtained the wavelength-calibrated spectra using the NRES commissioning IDL pipeline and obtained improved and stacked spectra for each night using the Stage2 IDL pipeline with a total integration time of 1440 s and a resulting S/N of ∼35.

The Tillinghast Reflector Echelle Spectrograph (TRES; Fűrész 2008) team obtained two reconnaissance (median S/N = 33.25) spectra of TOI 560 near opposite quadratures of TOI 560 b (phases 0.25 and 0.72 on the initial Sector 8 released ephemeris) at 2019 April 16 04:37 and 2019 April 19 03:49 UT. TRES has a resolving power R ∼ 44,000, and the wavelength range of the spectrograph is 390–910 nm. The Stellar Parameter Classification (SPC) tool specifically uses a ∼310 Å region of the spectrum (∼5050–5360 Å) to derive stellar parameters (Section 3.1.1). Spectra are processed following methods described in Buchhave et al. (2010).

We searched for close visual companions to the target using high-resolution imaging with both AO and speckle imaging at Gemini. The AO and speckle images are highly complementary: the speckle images reach higher resolutions in the optical, while the AO images reach deeper sensitivities beyond a few hundred milliarcseconds in the NIR, and are therefore more sensitive to more widely separated low-mass stars. We collected high-resolution AO images of TOI 560 with Gemini/NIRI (Hodapp et al. 2003) on 2019 May 26. Given that the star is bright in the K band, we collected images with individual exposure time 0.9 s with the Brγ filter (2.166 μm), to avoid saturating the detector. Our sequence consisted of nine such images, with the telescope dithered in a grid pattern between each frame, and we also collected daytime flats. The science images themselves can be used to construct a sky background frame, by median combining the dithered images.

TOI 560 was observed twice on 2020 March 16 and 2019 May 22 UT using the Zorro speckle instrument on the Gemini South 8 m telescope. ⁷³ Six sets of 1000 by 0.06 s exposures were collected for TOI 560 and processed with Fourier analysis in our standard reduction pipeline (see Howell et al. 2011).

We searched for stellar companions to TOI 560 with speckle imaging on the 4.1 m Southern Astrophysical Research (SOAR) telescope (Tokovinin 2018) on 18 May 2019 UT, observing in the Cousins I band a similar visible bandpass as TESS. This observation was sensitive to a 7.5 mag fainter star at an angular separation of 1'' from the target. More details of the observation, data reduction, and analysis are available in Ziegler et al. (2020).

We have gathered a total of 204 observations of TOI 560 over 30 nights using the iSHELL instrument at NASA InfraRed Telescope Facility (IRTF) atop Maunkea, Hawaii, USA from UT 2020 January 26 to UT 2021 May 29. We observe with iSHELL in KGAS mode covering the wavelengths of 2.18–2.47 μm. Our exposure times were always set at 300 s, and were repeated anywhere from 6–14 times consecutively per night to attempt obtaining a signal-to-noise ratio (S/N) per spectral pixel of ∼120, though our actual results varied from 46–152 due to variable seeing and atmospheric transparency conditions. A methane isotopologue (¹³CH₄) gas cell is used in the instrument (Plavchan et al. 2013; Cale et al. 2019) to constrain the line-spread function (LSF) and provide a common reference for the optical path wavelength. Along with each set of exposures within a night, we also collect a set of 5 × 15 s flat-field images with the gas cell removed for data reduction purposes, and particularly to mitigate flexure-dependent and time-variable fringing present in the spectra. In Cale et al. (2019), the RV pipeline is adapted from the CSHELL RV code described in Gao et al. (2016), and the CSHELL code was rewritten in a Python script pychell ⁷⁴ to adapt to iSHELL's larger spectral grasp with multiple orders.

When the light passes through a gas cell and an echelle spectrograph, the spectrum consists of multiorder spectra of the target with the gas spectrum superimposed, which can be used to subtract off instrument systematic shifts in the wavelength solution and changes in the LSF (Butler et al. 1996). We extract our spectra following the general procedures outlined in Cale et al. (2019). For each of the 29 orders in the regime that we observe between 2.18 and 2.47 μm, we do optimal spectral extraction (weighted summation) so that in each column eventually we get a single value, which gives us a 1D spectrum for each order.

The extracted spectra from pychell were forward-modeled and analyzed using the methods outlined in Cale et al. (2019), and the updated methods described in Reefe et al. (2021, submitted) briefly described herein. We forward model the extracted spectra independently for each order to account for the stellar spectrum, gas cell, telluric absorption, the LSF, and the spectral continuum to construct a complete spectral model. We find that pychell performs better when providing a synthetic stellar template for the stellar model to start from, which is based on properties of the host star.

We assume a solar metallicity and gravity, and we explore synthetic models ±500 K from an initial temperature of 4700 K, based upon the EXO-FOP-TESS values, in increments of 100 K to identify the optimal synthetic stellar template that produces the lowest rms flux residuals in forward modeling our extracted spectra. We use the Spanish Virtual Observatory's BT-Settl models accessible from their web server, ⁷⁵ which we further refine by Doppler broadening the spectrum to the rotational velocity of the star (3.1 km s⁻¹), as found by TRES extracted spectra (Section 3.1.1). We find that the T_eff = 4900K synthetic model minimizes the flux residuals of our observed iSHELL spectra, even though they're hotter than our adopted stellar temperature, and we use this stellar template as our initial stellar spectrum model guess. Barycenter velocities are also generated as an input via the barycorrpy library (Kanodia & Wright 2018). After each order is forward-modeled, we generate one RV measurement per order and per exposure. We optimally coadd RVs across orders and exposures within a night using the same procedures as in Cale et al. (2019).

We have filtered out three individual spectra from UT 2020 March 7 and three more individual spectra from UT 8 March 2020, due to modeled RVs that were in disagreement with other spectra from the same night by >1 km s⁻¹. We suspect this is due to an initially poor seeing on the nights of observation of 13 that was later improved to 10, giving us overall S/Ns of only 54 and 46, respectively, for the nights. Finally, we discard nights with RV uncertainties >20 m s⁻¹.

We obtained 14 RVs of TOI 560 from PFS on the 6.5 m Magellan II Telescope at Las Campanas Observatory in Chile. These spectra were reduced, and RVs were extracted and published in the Magellan-TESS Survey I paper (Teske et al. 2021). The PFS Spectrograph has a resolution of R = 120,000 and wavelength range of 391–734 nm. The observations were carried out between UT 2019 April 18 and May 24.

We include 14 Keck-HIRES (Vogt et al. 1994) observations of TOI 560 in our analyses. Keck-HIRES is located atop Maunkea, Hawaii, USA with spectral resolution R = 67,000. The majority of these observations took place between UT 2019 October 20 and 2020 January 21, with the final night contemporaneous with the start of iSHELL RVs. Exposure times range from 204–500 s, yielding a median S/N ≈ 234 at 550 nm per spectral pixel. HIRES spectra are processed and RVs computed via methods described in Howard et al. (2010).

Minerva-Australis collected two nights of observations of TOI 560 (Addison et al. 2019, 2021a). Minerva-Australis consists of an array of four independently operated 0.7 m CDK700 telescopes situated at the Mount Kent Observatory in Queensland, Australia (Addison et al. 2019). Each telescope simultaneously feeds stellar light via fiber optic cables to a single KiwiSpec R4–100 high-resolution (R = 80,000) spectrograph (Barnes et al. 2012) with wavelength coverage from 480–620 nm. We obtained a total of 10 individual spectra of TOI 560 on UT 2019 May 23 and UT 2019 May 29 using three and two of the four Minerva-Australis telescopes, respectively—twice per night per telescope—and exposure times of 30–60 minutes. Wavelength calibration is achieved using a simultaneous Th-Ar calibration fiber. RVs were derived by cross-correlation, where the template being matched is the mean spectrum. Given that only two nightly RV data points were obtained, we do not include the Minerva-Australis RVs in the remainder of the analysis presented herein.

The reconnaissance spectroscopy observations reveal a typical late K dwarf. Stellar parameters from the NRES spectra were obtained using SpecMatch (Vanderburg et al. 2019), and summarized in Table 11. For the TRES spectra, there is a small velocity difference that is consistent with the photometric ephemeris, given that the data were collected approximately at opposite quadratures. In Table 12, the results for the stellar parameters for the two TRES observations are presented, as well as values from the SPC tool (Buchhave et al. 2012, 2014). Given the absence of any large bulk RV changes or spectroscopic binarity, we do not coadd the multiepoch spectroscopy to improve the cumulative S/N (and consequently reduce the stellar parameter uncertainties), due to current systematic limitations in these approaches (e.g., Duck et al. 2022). We adopt floor errors of 50 K in T_eff, 0.10 in log(g), 0.08 in [m/H], and 0.5 km s⁻¹ in V_rot. Note that V_rot does not include a correction for the contribution of macroturbulence, and the velocities in Table 12 are on the TRES native system and do not correct for the gravitational redshift. The TRES and NRES analyses yield consistent stellar characterization, with the SPC analysis favoring a slightly higher but not statistically significant rotational velocity and metallicity.

We undertake two independent analyses of the broadband SED of the star, and we describe each in turn. The first approach is particularly useful for assessing the evidence for any near-UV (NUV) excess from consideration of the broadband photometry alone. Our second approach is used for our adopted final stellar parameters and is based upon holistic modeling of the SED jointly with the transit light curves and model isochrones.

First, we performed a single-parameter fit of the SED of the star together with the Gaia EDR3 parallax (with no systematic offset applied; see, e.g., Stassun & Torres 2021), in order to determine an empirical measurement of the stellar radius and following the procedures described in Stassun & Torres (2016), Stassun et al. (2017), and Stassun et al. (2018). We use the UBV magnitudes from the compilation of Mermilliod (2006), the JHK_S magnitudes from the Two Micron Ally Sky Survey (2MASS), the W1–W4 magnitudes from the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010), the G_BP G_RP magnitudes from Gaia, and the NUV magnitude from GALEX. Together, the available photometry spans the full stellar SED over the wavelength range 0.2–22 μm (see Figure 2). We used a NextGen stellar atmosphere model, fixing the model effective temperature (T_eff), metallicity ([Fe/H]), and surface gravity ($\mathrm{log}g$) to the values derived from the TRES SPC analysis in Section 3.1.1 and Table 12. The remaining free parameter is the extinction A_V , which we fixed at zero given the proximity of the TOI 560 system to Earth. The resulting fit (Figure 2) has a reduced χ² of 1.7, excluding the GALEX NUV flux, which indicates a moderate level of activity (see below). Integrating the (unreddened) model SED gives the bolometric flux at Earth, F_bol = 6.01 ± 0.14 × 10⁻⁹ erg s⁻¹ cm⁻². Taking the F_bol and T_eff together with the Gaia parallax gives the stellar radius R_⋆ = 0.656 ± 0.016 R_⊙ (consistent with the NRES analysis). In addition, we can estimate the stellar mass from the empirical relations of Torres et al. (2010), giving M_⋆ = 0.69 ± 0.04 M_⊙, which is consistent with the less precise but empirical estimate of M_⋆ = 0.75 ± 0.08 M_⊙ obtained directly from R_⋆ together with the spectroscopic $\mathrm{log}g$.

Zoom In Zoom Out Reset image size

Second, we determine characteristics of the host star, such as effective temperature, gravity, and metallicity, by performing a joint ameba fit followed by MCMC posterior sampling of both stellar properties and planet properties of TOI 560 b and c assuming a two-planet, single-star scenario (from the TESS transit data) simultaneously with EXOFASTv2. EXOFASTv2 is an (IDL) framework for MCMC simulations of exoplanet transits and RVs created by Eastman et al. (2013). Details on the planet transit analysis are given in Section 3.2, but here we present the stellar modeling analysis steps. We start the MCMC with as few assumptions as possible—namely, we place no priors on the spectral type. To ensure we cover an adequate portion of parameter space to find a believable result, we employ parallel tempering with eight parallel threads, following Eastman et al. (2019). We place priors on V-band extinction, parallax, and metallicity summarized in Table 5 along with the priors used in the transit analysis. We simultaneously fit with Mesa Isochrones and Stellar Tracks (MIST) and an SED function. We then use the posteriors of this run as priors for a second iteration MCMC analysis to examine the stability of the solution and adopt our final stellar parameters.

The results of our EXOFASTv2 stellar characterization are shown in and Table 3. The analysis yields a typical K dwarf consistent with the reconnaissance spectroscopy and SED analysis presented in Section 3.1. We provide two separate estimates for R_* and T_eff corresponding to the MIST results and the SED results, respectively. SED results are marked with the subscript SED, while MIST results have no subscript. We adopt the MIST values for the remainder of the analysis, but both sets of values are consistent with one another and with the analysis carried out with the SED and reconnaissance spectroscopy. Finally, we note that the age estimate from EXOFASTv2 shown in Table 3 is drawn from a very flat posterior distribution and is not a reliable estimate; e.g., the stellar age is unconstrained from this particular analysis. Instead, we adopt the final rotation period of 12.2 ± 0.1 days as described in Section 3.1.3.

Table 3. Median Values and 68% Confidence Interval for the Stellar Host, TOI 560, from Our ExoFASTv2 Analysis

Stellar Parameters	Units	Values
M*	Mass (M⊙)	${0.702}_{-0.025}^{+0.026}$
R*	Radius (R⊙)	0.677 ± 0.017
R*,SED	Radius a (R⊙)	${0.6819}_{-0.0094}^{+0.0099}$
L*	Luminosity (L⊙)	${0.1823}_{-0.0060}^{+0.0062}$
FBol	Bolometric Flux ((erg s−1) cm− 2)	5.85 × 10−9 ± 2 × 10−10
ρ*	Density (g cm− 3)	${3.19}_{-0.24}^{+0.26}$
$\mathrm{log}$ g	Surface gravity (cm s− 2)	${4.623}_{-0.024}^{+0.025}$
Teff	Effective Temperature (K)	${4582}_{-62}^{+64}$
Teff,SED	Effective Temperature b (K)	4568 ± 45
[Fe/H]	Metallicity (dex)	−${0.055}_{-0.052}^{+0.044}$
[Fe/H0]	Initial Metallicity c	$-{0.051}_{-0.062}^{+0.057}$
Age	Age d (Gyr)	${6.7}_{-4.5}^{+4.7}$
EEP	Equal Evolutionary Phase e	${329}_{-30}^{+13}$
Av	V-band extinction (mag)	${0.076}_{0.051}^{0.046}$
σSED	SED photometry error scaling	${1.68}_{-44}^{+0.77}$
$\bar{\omega }$	Parallax (mas)	31.683 ± 0.032
d	Distance (pc)	31.562 ± 0.032
Wavelength Parameters	Units	B	R	z'
u1	linear limb-darkening coeff	${0.50}_{-0.32}^{+0.037}$	${0.68}_{-0.039}^{+0.037}$	${0.58}_{-0.38}^{+0.44}$
u2	quadratic limb-darkening coeff	${0.27}_{-0.044}^{+0.038}$	${0.02}_{-0.038}^{+0.045}$	${0.18}_{-0.046}^{+0.043}$
AD	Dilution from neighboring stars	⋯	⋯	⋯
Telescope Parameters	Units	HIRES	PFS
γrel	Relative RV Offset (m s−1)	${0.1}_{-3.6}^{+3.5}$	$-{10.8}_{-2.3}^{+2.2}$	0.5 ± 3.1
σJ	RV Jitter (m s−1)	${12.7}_{-2.4}^{+3.4}$	${8.0}_{-1.7}^{+2.5}$	${10.9}_{-2.5}^{+3.1}$
${\sigma }_{J}^{2}$	RV Jitter Variance	${160}_{-55}^{+97}$	${64}_{-24}^{+47}$	${119}_{-48}^{+79}$

Notes. See Table 3 in Eastman et al. (2019) for a detailed description of all parameters.

^a This value ignores the systematic error and is for reference only. ^b This value ignores the systematic error and is for reference only. ^c The metallicity of the star at birth. ^d This posterior is relatively flat and unconstrained, and does not take into account the stellar-rotation period analysis estimated in Section 3.1.3 and/or Figure 3. ^e Corresponds to static points in a star's evolutionary history. See Section 2 in Dotter (2016).

Download table as:  ASCIITypeset image

Finally, TOI 560 has a GAIA RUWE value of 1.030, which is consistent with a single star (Stassun & Torres 2021). The RUWE statistic is an indication of the goodness of fit of the astrometric solution, and values >1.4 have been shown to be associated with unresolved binaries. A search of EDR3 also shows no stars with similar parallaxes within 75', so TOI 560 likely does not possess any distant, resolved binary companions either.

Lomb-Scargle periodograms of the SuperWASP seasonal light curves show a clear and consistent rotation period in each season for TOI 560 of 12 days, with evolution of the light-curve features between seasons (Figure 3). The periodogram peak is statistically significant in all seasons (p-value < 1%) and the average periodogram peak period from all four seasons yields a rotation period of 12.2 ± 0.1 days. We next use the star's NUV excess (Figure 2) to estimate the star's rotation and age via empirical rotation–activity–age relations. The observed NUV excess implies a chromospheric activity of $\mathrm{log}R{{\prime} }_{\mathrm{HK}}\,=$ −4.5 ± 0.1 via the empirical relations of Findeisen et al. (2011), which is fully consistent with the measured value of −4.45 ± 0.05 from Gomes da Silva et al. (2021) and −4.60 from our HIRES spectra. The HIRES and TRES spectra show core Ca II HK flux in emission, consistent with this moderate level of activity. This value of $\mathrm{log}R{{\prime} }_{\mathrm{HK}}$ in turn implies a stellar-rotation period of P_rot = 14 ± 4 d via the empirical relations of Mamajek & Hillenbrand (2008). This rotation period estimate is also consistent with that estimated from the spectroscopic $v\sin i$ and R_⋆, which gives ${P}_{\mathrm{rot}}/\sin i=11\pm 3$ day. Thus, the observed photometric modulation, NUV excess, and the observed rotational velocity all yield consistent rotation periods. This establishes the rotation period of the TOI 560 host star.

Zoom In Zoom Out Reset image size

However, this 12 day rotation period is roughly a factor of 2 slower than that reported by Canto Martins et al. (2020), who reported a "dubious" rotation period 7.17 days from their independent TESS light-curve analysis. We can also estimate rotation periods for TOI 560 by performing Lomb-Scargle periodograms (Lomb 1976; Scargle 1982) of the TESS Sector 8 and 34 light curves of TOI 560. In Figure 4, the TESS light curves contain significant peaks at 12.7 days for PDCSAP and 14.1 days for SAP, both for the Sector 8, respectively, that are not obvious integer fraction multiples of one another nor of 7.17 days; we do find a secondary peak at 7.29 days for the Sector 8 light curve, consistent with Canto Martins et al. (2020). However, the dominant period in the Sector 8 light curve at 12.7 days is reasonably consistent with the WASP, NUV, and rotational velocity inferred rotation periods. A more detailed FF' modeling of the TESS light curves in Section 3.2.3 with a GP and uniform rotation period prior from 2–20 days also favors a rotation period of ${7.13}_{-0.13}^{+0.31}$ days, also consistent with Canto Martins et al. (2020). However, given that the time segments of the TESS light curve in-between data downlinks are themselves <12 days in time baseline, we find that the TESS light curves themselves are insensitive to a 12 day rotation period and are thus not reliable; the 7.1 and 5.6 day power in the TESS light curves may be from a combination of the first harmonic of the stellar-rotation period at P_rot /2, potentially the spot evolution timescale, and the presearch data conditioning (PDC) of the TESS PDC-SAP light curves. Thus, we adopt the rotation period of 12.2 ± 0.1 days from the SuperWASP light curves in the remainder of our analysis.

Zoom In Zoom Out Reset image size

For this rotation period, using the age–rotation relation from Mamajek & Hillenbrand (2008) based upon the power-law model of Barnes (2007), we derive a stellar age estimate of ${\tau }_{\star }={591}_{-97}^{+130}$ Myr. However, Curtis et al. (2020) recently showed that for K dwarfs like TOI 560, there is an age–rotation "pile-up" at rotation periods of ∼10–15 days where the spin-down of K dwarfs stalls between ages ∼0.6–1.4 Gyr, before resuming for ages >1.4 Gyr. We also verified whether TOI 560 may be a plausible member of a nearby moving group or association with the BANYAN Σ (Gagné et al. 2018) tool, which compares the XYZ Galactic coordinates and UVW space velocities of known stellar associations and assigns a membership probability based on Bayes' theorem with the option to marginalize over unknown quantities such as heliocentric RVs. We find a negligible membership probability to all 27 young associations included in BANYAN Σ. We also compared the XYZ and UVW of TOI 560 to those of 1000 moving groups and open clusters identified in the literature so far within 500 pc of the Sun, and we find no plausible association to which TOI 560 may belong to. Inspecting the position of TOI 560 in a Gaia eDR3 (Gaia Collaboration et al. 2018) color–magnitude diagram (Figure 5) reveals a picture consistent with an age of at least 100 Myr. TOI 560 has been detected in GALEX DR5 (Bianchi et al. 2011), and its position in a GALEX-Gaia color–color diagram (Figure 6) shows that it is located on the UV-bright end of field stars and on the UV-faint end of the Pleiades members with similar Gaia eDR3 G − G_RP color (which tracks spectral types), thus indicating that TOI 560 is likely older than the Pleiades, but plausibly on the younger end of the age distribution of field stars, and consistent with the rotation-based age estimate. Additionally, there is no lithium detected, consistent with this star not being <100 Myr (Delgado Mena et al. 2015). We thus conclude from the rotation period of TOI 560 and the K dwarf rotational evolution pile-up that the age of TOI 560 is between 150 Myr and 1.4 Gyr.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the Gemini/NIRI observations, we carried out the data reduction using a custom set of IDL codes, with which we removed bad pixels, sky-subtracted and flat-corrected the frames, and then aligned the stellar position between frames and coadded the images (Hodapp et al. 2003). For Gemini South, we utilize only the observations from 2019 May, as the March observations had poorer seeing and sky conditions, albeit giving similar results to those obtained about a year earlier. Zorro provides simultaneous speckle imaging in two bands (562 and 832 nm) with output data products including a reconstructed image in Figure 7 with robust contrast limits on companion detection (e.g., Howell et al. 2016).

Zoom In Zoom Out Reset image size

For our Gemini/NIRI observation, we did not identify visual companions anywhere in the field of view, which extends to at least 7'' in all directions from the target star, and TOI 560 appears single to the limit of the Gemini/NIRI resolution. The sensitivity as a function of separation is shown in Figure 7, along with an image of the target. The images are sensitive to companions 8 mag fainter than the host star in the background limited regime (beyond ∼1''). For our Gemini South observations, Figure 7 shows our final 5σ contrast curves and the reconstructed speckle images. We find that TOI 560 is single to within the contrast limits achieved by the observations. No companion is identified brighter than 5–8 mag below that of the target star from the diffraction limit (20 mas) out to 1''. At the distance of TOI 560 (d = 31 pc), these angular limits correspond to spatial limits of 0.6–37 au. Finally, for our SOAR observation, Figure 7 shows the 5σ detection sensitivity and speckle autocorrelation functions from the observations. No nearby stars were detected within 3'' of TOI 560 in the SOAR observations.

Before investing detailed resources in the RV analysis of a TESS candidate planetary system, it is useful to rule out false-positive scenarios caused by eclipsing binaries and other systematics, and there are numerous diagnostics that we employ using the TESS light curves alone. We ran two separate vetting analyses on the TOI 560 data gathered by TESS.

The first of our vetting tests was performed with the EDI-Vetter Unplugged tool (Zink et al. 2020). ⁷⁶ EDI-Vetter Unplugged checks for several diagnostics that could be indicative that the target is a false positive. First, EDI-Vetter checks for neighboring flux contamination from nearby stars in the sky, to make sure that the signal is indeed correlated with the expected target; this is necessary due to the large angular size of TESS pixels. Second, it examines the abundance of outliers from the transit model, as well as the validity of each individual transit—e.g., odd/even test, if transits fall into masked data gaps. Third, it also searches for the statistically significant presence of any secondary eclipses, which could indicate an eclipsing binary. Other measures that are considered are the similarities between transit signals, the phase coverage of transit signals, and the relative length of transit duration in comparison to period. The results of our EDI-Vetter analysis are presented in Table 4.

Second, we also scrutinized the TESS light curves using the DAVE software (Kostov et al. 2019) developed at the NASA Goddard Space Flight Center. This tool evaluates similar criteria as the EDI-Vetter tool, such as odd/even transit differences and secondary eclipses, but also searches for any photocenter shift, or the difference between the TIC position and the measured photocenter of the transits during transit, which is often also a diagnostic of a blended eclipsing binary.

In our EDI-Vetter analysis, two of the criteria were flagged as a potential false positive: the uniqueness of the transit signal and the planet masks (Table 4). However, the latter is due to the system being a multiplanet system, with two different transiting planets, and the former a transit of TOI 560 c falling into the data download gap of Sector 8. In our DAVE analysis, we show full transit data and phased odd/even, secondary, tertiary, and positive transits for each sector in Figures 38 and 39. We see little to no variation between the primary odd and even transits, and we see no statistically significant drop in flux during the expected ephemerides for a secondary eclipse or tertiary transit. However, we note the appearance of an odd–even transit-depth difference for TOI 560 b. Due to the youth of the star, there is significant out-of-transit variability from the star that is not filtered out by the DAVE vetting tool. This in turn impacts its transit-depth determination, resulting in the even–odd discrepancy. Thus, we effectively ignore this metric present in the DAVE tool as a potential sign of a false positive, which is supported by the multiplanet nature of the system.

We also see no statistically significant positive flux variation that could be indicative of lensing between eclipsing binary stars. The evaluation of the photocenter shift during transit helps validate that the transit belongs to the target star, and not a faint companion. Photocenter difference plots for each sector are shown in Figure 40, and overall show no major variations between the expected and observed locations. Thus, this vetting analysis from DAVE and EDI-Vetter, combined with our high-contrast imaging, rules out a blended or background eclipsing binary companion, and confirms that the transit signals originate from TOI 560, and not another nearby point source. Further, when the increasingly low probability (e.g., <1%, Lissauer et al. 2012) that multiple blended or background eclipsing binaries would be spatially coincident with TOI 560 and produce two different eclipsing signals near a 3:1 orbital period coincidence are combined, we can conclude that TOI 560 is the host of the two eclipsing companions. We next turn to presenting our results on constraining the masses of TOI 560 b and c.

With the release of the TESS Sector 34 light curve, we readily identified by eye two transits of a second planet candidate with an orbital period of ∼18.9 days. Several other teams independently identified this second transiting planet during our analysis, which was identified by the TESS mission as TOI 560 c (P_c = 18.8805 days). After masking the transits of TOI 560 b (P_b = 6.3980 days), we computed BLS (Kovacs et al. 2002) and TLS (Hippke & Heller 2019) periodograms and phased light curves (LCs) of the Sector 34 light curve in Figure 8. Due to unfortunate timing, TOI 560 c did not transit in the Sector 8 TESS light curve (Figure 8). We also excluded a period for TOI 560 c of one-half its value, as a third transit of TOI 560 c could have occurred during a downlink data gap in the TESS Sector 34 light curve. However, under that scenario, a transit of TOI 560 c would have occurred in the Sector 8 light curve, which is ruled out.

Zoom In Zoom Out Reset image size

Stars are not uniform disks; they exhibit RV variations from rotationally modulated activity, e.g., star spots and plages on the surface of the star (Dumusque et al. 2011a, 2011b; Plavchan et al. 2015; Fischer et al. 2016, and references therein). Depending on the temperature contrast between magnetically active regions and the surrounding photosphere of the star, these active regions will induce an apparent RV shift that is quasiperiodic as the active regions evolve on timescales distinct from the rotation period of the star. GPs have been widely and successfully employed in modeling the apparent RVs due to stellar activity (Haywood 2015).

The same active regions also give rise to photometric variations from the hemisphere visible to the observer. Thus, we can model the light curves of TOI 560 out of transit with a GP, and use the light-curve hyperparameter posteriors to constrain the GP model hyperparameter priors for our subsequent analysis of the RVs, following the ${{FF}}^{{\prime} }$ technique (Aigrain et al. 2012). This method allows us to generate what a simulated RV curve would look like from the stellar activity:

$$\begin{eqnarray}&&{\rm{\Delta }}{\mathrm{RV}}_{\mathrm{spots}}(t)=-F(t){F}^{{\prime} }(t){R}_{* }/f.\end{eqnarray} \tag{ 1 }$$

Here, F is the photometric flux, ${F}^{{\prime} }$ is the derivative of F with respect to time, f represents the relative flux drop for a spot at the center of the stellar disk, and R^* is the stellar radius.

Our GP kernel describes the functional covariance between any two measurements separated in time. The kernel by definition must be a square, symmetric covariance matrix "K" of length equal to the number of observations. We follow Haywood (2015) who introduced the quasiperiodic (QP) Kernel composed of a decay and periodic term:

$$\begin{eqnarray}&&{K}_{{QP}}({t}_{i},{t}_{j})={\eta }_{\sigma }^{2}\mathop{\overbrace{\exp \left[-\displaystyle \frac{{\rm{\Delta }}{t}^{2}}{2{\eta }_{\tau }^{2}}\right]}}\limits^{\mathrm{Decay}}\mathop{\overbrace{\exp \left[-\displaystyle \frac{1}{2{\eta }_{l}^{2}}{\sin }^{2}\left(\pi \displaystyle \frac{{\rm{\Delta }}t}{{\eta }_{p}}\right)\right]}}\limits^{\mathrm{Periodic}}\end{eqnarray} \tag{ 2 }$$

where Δt = ∣t_i − t_j ∣. Here, η_P typically represents the stellar-rotation period, η_τ is the mean spot lifetime, and η_ℓ is the relative contribution of the periodic term, which may be interpreted as a smoothing parameter (larger is smoother). η_σ is the amplitude of the autocorrelation of the activity signal. The set of {η_σ , η_τ , η_P , η_L } constitutes the set of GP model hyperparameters that are constrained by the data. There exist families of stellar activity models that generate a particular covariance matrix given by a specific set of hyperparameters, and thus GPs are a flexible framework for characterizing the nondeterministic time evolution of stellar activity. The ${{FF}}^{{\prime} }$ analysis does enable us to estimate the expected amplitude of the RV variations at wavelengths comparable to the light-curve observations, η_σ , but we end up not using that amplitude in our RV modeling. Instead, we do use the ${{FF}}^{{\prime} }$ analysis to get estimates of the remaining hyperparameters, particularly η_τ , η_p , and η_ℓ for use in the RV analysis.

We first apply our GP model to the SuperWASP data for the four seasons together (Figure 3). We use wide uniform priors on hyperparameters η_σ ∼ ${ \mathcal U }(0,50)$, η_τ ∼ ${ \mathcal U }(1,100)$, and η_ℓ ∼ ${ \mathcal U }(0.02,1)$. Instead of accounting for intrinsic uncertainties in the light curve, we also fit a jitter term η_LC . For the hyperparameter η_P , we use a prior of ${ \mathcal U }(11.5,13)$ based on periodogram analysis of the WASP light curves in Section 3.1.3. We recover from the MCMC of the WASP data a spot lifetime η_τ of ${57.96}_{-19.39}^{+23.59}$ days, a rotation period of η_P of ${12.03}_{-0.12}^{+0.13}$ days, and a smoothing hyperparameter ${\eta }_{l}={0.44}_{-0.09}^{+0.11}$, which are all consistent with our young star; such long spot lifetimes are seen for other young systems (Figure 9; e.g., Cale et al. 2021). We adopt these three posteriors as priors in constraining the GP model for our RV stellar activity analysis.

We next analyze both sectors of TESS PDC-SAP and SAP separately in order to see if we can recover similar GP kernel hyperparameters {η_σ , η_τ , η_P , η_L } to our SuperWASP analysis. We first median normalize the TESS SAP light curve, and mask out the transits and the edges of the light-curve data, particularly for Sector 8, where there exists a spacecraft systematic producing a "ramp-up" effect at the start of the sector and after the data downlink gap mid-sector. We then fit the remaining light curve via a cubic-spline regression (scipy.interpolate.LSQUnivariateSpline; Virtanen et al. 2020) for each sector individually with knots sampled in units of 1.3 days (excluding any that happen to fall within the TESS data dump regions). The particular value of 1.3 is chosen to be small with respect to the stellar activity timescales, but long with respect to the cadence, so transits are "smoothed" over to produce a "noiseless" and smooth representation of the light curves. Then, in Figure 10 we apply the ${{FF}}^{{\prime} }$ technique (Aigrain et al. 2012) to the spline regression fits.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We used a uniform prior centered on 12.2 days for the stellar-rotation period hyperparameter. As expected due to the relatively short TESS time baselines of two ∼13 day observing sequences per sector, the TESS light curve does not constrain the stellar-rotation period further than the information inferred from the SuperWASP light-curve analysis.

For the TESS SAP Sector 8 (34; 8+34) light curve, we recover from the MCMC a spot lifetime η_τ of ${12.94}_{-9.03}^{+23.55}$ (10.10${}_{-6.62}^{+23.49}$ days; ${15.97}_{-11.74}^{+22.55}$) days (Figure 11). In other words, the TESS light-curve analysis leads to a significantly shorter spot lifetime analysis than inferred from the SuperWASP analysis as well, even for the SAP fluxes, but we discard these results as a consequence of the shorter TESS time baseline compared to the SuperWASP light curve, even though the former is at higher photometric precision. For the TESS SAP Sector 8 (34; 8+34) light curve, we recover from the MCMC a smoothing hyperparameter η_l = 0.25 ± 0.03 (η_l = 0.25 ± 0.03; η_l = 0.25 ± 0.03), and these are all somewhat smaller than the value inferred from the SuperWASP analysis, and this could potentially be due to the higher precision of the TESS light curves and improved photometric sampling of the rotation period. Finally, from an MCMC of the TESS PDC-SAP light curves, we recover an even shorter spot lifetime η_τ , indicating that the PDC can lead to introducing systematics that lead to inaccurate characterization of the stellar activity (Figure 12).

Zoom In Zoom Out Reset image size

In this section, we perform an EXOFASTv2 analysis of the candidate planet's transit light curves. After normalizing the TESS PDC-SAP data as described in Section 2, we cut out a region of at least 8 hr around each individual transit and create separate data files for each transit and use these as input data for EXOFASTv2.

We jointly model the TESS, Spitzer, and ground-based light curves, and all RVs with EXOFASTv2. The EXOFASTv2 RV results are presented in Appendix H in Figure 41. However, while ExoFAST can jointly model the light curves and RVs, the EXOFASTv2 RV model does not account for the stellar activity manifested in the RV measurements. Thus, we perform an independent modeling of the RVs with a stellar activity model in Section 3.3.

Our minimal set of priors are detailed in Table 5, including the period P, time of conjunction, and planetary radius R_P /R_* for TOI 560 b and c (the EXOFASTv2 full results are in Table 7). We allow other model parameters, i.e., eccentricity e, inclination i, and argument of periastron ω to vary with no imposed priors, starting with circular e and edge-on i values. The posterior values for this initial MCMC run are then used as the initial values for a second iteration run, though we keep the same uniform and Gaussian priors as the initial run. We allow this second run to run longer, and we confirm the second MCMC converges on the same results within 1σ to check for the robustness of the MCMC posteriors. Note that the model is parameterized in the EXOFASTv2 standard basis for transit only fits. This model parameterization is $\{{T}_{C},\mathrm{log}P,{R}_{P}/{R}_{* },\cos i,\sqrt{e}\sin \omega ,\sqrt{e}\cos \omega \}$ (Eastman et al. 2019).

Table 5. Planetary and Stellar Prior Probability Distributions for Our EXOFASTv2 MCMC Simulations

Parameter	Initial	Priors
[Units]	Value
	[P0]
Pb (days)	6.3981	${ \mathcal U }({P}_{0}\pm 10 \% )$
TCb (days)	2458517.69007	${ \mathcal U }({P}_{0}\pm P/3)$
Pc (days)	18.881	${ \mathcal U }({P}_{0}\pm 10 \% )$
TCc (days)	2458533.59092	${ \mathcal U }({P}_{0}\pm P/3)$
Rp,b /R*	0.0388	⋯
Rp,c /R*	0.0382	⋯
Mp,b /M⊙	0.000030	⋯
Mp,c /M⊙	0.000030	⋯
R*	0.674	⋯
R*,SED	0.679	⋯
Teff (K)	4599.94	⋯
Teff,SED (K)	4591.85	⋯
feh	−0.031	${ \mathcal N }(-0.210,0.080)$
initfeh	−0.043	⋯
eep	319.8	⋯
logM*	−0.145	⋯
AV	0.097	${ \mathcal U }(0,0.143)$
errscale	1.22	⋯
distance (pc)	31.567	⋯

Note. ${ \mathcal U }({\ell },r)$ signifies a uniform prior with left bound ℓ and right bound r. A_V is the galactic extinction in the V band, "feh" is metallicity, "initfeh" is the initial metallicity at the time of the zero-age main sequence, "eep" stands for equivalent evolutionary point, logM_* is the log of the stellar mass, and "errscale" is the error scale. The eccentricity e and the argument of periastron ω are assumed to be circular, with no priors.

Download table as:  ASCIITypeset image

To jointly model the RVs from iSHELL, PFS, and HIRES, we again use pychell, this time in combination with the co-dependent package optimize, which is a general-purpose Bayesian analysis tool that pychell expands upon with RV-specific MCMC tools. In Section 5.2 we show that no statistically significant periodic signals were identified in the raw RVs, and so we proceed herein with the assumption of a two-planet model in our Bayesian analysis with the ephemerides from TESS. Each planet in our model consists of five parameters, which comprise a complete orbit basis: the period P, time of conjunction TC, eccentricity e, argument of periastron ω, and RV semiamplitude K, with subscripts denoting the planet that said parameter is associated with (in this case, b and c). We next include a GP model for the stellar activity, which is described further in the next subsection. We also include an absolute RV offset term (γ) to account for the average overall recessional velocity of the star, and a jitter term (σ) that quantifies the variational RV amplitude not accounted for by any modeled planets, stellar activity, or instrumental noise that is nonfrequency dependent on our timescales of interest (e.g., white noise). Analyses of the PDC-SAP TESS transit data (Section 4.1.1, Figure 13 and Table 5) shows a period and time of conjunction consistent with those listed on the ExoFOP-TESS database in the TESS Object of Interest (TESS Project section) ⁷⁷ (Akeson et al. 2013), which are P_b = 6.398069 ± 0.000015 days, P_c = 18.879744 ± 0.000162, TC_b = 2458517.690108 ± 0.000747 days, and TC_c = 2458533.620329 ± 0.005422, all of which have uncertainties that are orders of magnitude finer than we can hope to resolve with our RV cadence. We thus decide to lock both parameters at their nominal values for all analyses in this paper. For all MCMC results quoted, the median value is defined as the 50th percentile, while the lower and upper uncertainty bounds are the 15.9th and 84.1th percentile posterior confidence intervals, respectively. We sample the posterior distributions using the emcee package (Foreman-Mackey et al. 2013) for a subset of models to determine parameter uncertainties, always starting from the MAP-derived parameters (i.e., affine invariant is our actual sampler). We perform a burn-in phase of 1000 steps followed by an MCMC analysis for approximately 50 times the median autocorrelation time (steps) of all chains. The prior, posterior distributions, and results of this model are given in Section 4. Finally, we also perform a model comparison test (Table 9) by calculating the $\mathrm{ln}{ \mathcal L }$, small-sample Akaike Information Criterion (AIC; Akaike 1974; Burnham & Anderson 2002), and Bayesian information criterion (BIC) for each combination of planets in our model. These results are shown in Section 4.2.

Zoom In Zoom Out Reset image size

Rather than applying an independent GP model to each individual RV data set for PFS, HIRES, and iSHELL, we extend the kernel in Equation (1) to utilize a single global (joint) GP model and covariance kernel across multiple spectrographs. We follow Cale et al. (2021) where they already implemented our desired framework in two Python packages. As in Cale et al. (2021), we first reparameterize the GP amplitude through a linear kernel, hereafter referred to as the J₁ kernel as in Cale et al. (2021), as follows:

$$\begin{eqnarray}&&{K}_{{J}_{1}}({t}_{i},{t}_{j})={\eta }_{\sigma ,s(i)}{\eta }_{\sigma ,s(j)}\times \exp [...].\end{eqnarray} \tag{ 3 }$$

Here, η_σ,s(i) and η_σ,s(j) are the effective stellar activity amplitudes for each spectrograph s at times t_i and t_j , respectively, where s(i) represents an indexing set between the observations at time t_i . ⁷⁸ In the J₁ kernel, each instrument is given its own independent amplitude hyperparameter, η_σ,s(i), but the other three hyperparameters are shared between all instruments. Also, the covariance kernel is still square, but now has dimensions corresponding to the total number of observations across all spectrographs, and thus represents a "joint" covariance matrix.

Second, to first order, we expect the amplitude from stellar activity to be linearly proportional to frequency (or inversely proportional to wavelength; Cale et al. 2021). This approximation is a direct result of the spot-contrast scaling with the photon frequency (or inversely with wavelength) from the ratio of two blackbody functions with different effective temperatures (Reiners et al. 2010). Thus, we also consider a variation of this kernel that further enforces the expected inverse relationship between the amplitude with wavelength. As in Cale et al. (2021), we parameterize the kernel to become:

$$\begin{eqnarray}&&{K}_{{J}_{2}}({t}_{i},{t}_{j},{\lambda }_{i},{\lambda }_{j})={\eta }_{\sigma ,0}^{2}{\left(\displaystyle \frac{{\lambda }_{0}}{\sqrt{{\lambda }_{i}{\lambda }_{j}}}\right)}^{2{\eta }_{\lambda }}\times \exp [...],\end{eqnarray} \tag{ 4 }$$

which hereafter we refer to as the J₂ kernel, as in Cale et al. (2021). Here, η_σ,0 is the effective amplitude at λ = λ₀, and η_λ is an additional power-law scaling parameter with wavelength to allow for a more flexible nonlinear (with frequency) relation. λ_i and λ_j are the effective wavelengths for observations at times t_i and t_j , respectively. Note for this kernel that we are ignoring the chromatic effects of limb-darkening and convective blueshifts on the RVs, which will impact the covariance matrix; however, given the flexibility of the GPs, Cale et al. (2021) found this chromatic kernel effectively recovers the wavelength dependence of stellar activity induced RV variations. For both Equations (3) and (4), the expression within square brackets is identical to that in Equation (2).

We apply Kernels J₁ and J₂ to our RV data, using the priors in Table 6, the priors from our FF${}^{{\prime} }$ analysis in Section 3.2.3 for the model hyperparameters, as well as a set of disjoint quasiperiodic kernels as in Equation (1), one GP for each spectrograph akin to RadVel (Fulton et al. 2018). In the analysis presented herein, we fix η_l and η_τ to the median values from the SuperWASP light-curve analysis posteriors; letting these model hyperparameters vary in our RV analysis—with prior minimums of 0.2 and 20 days, respectively—yields quantitatively identical results for the median posterior values for K_b and K_c , albeit with slightly larger confidence intervals. We also analyze the RVs using no GP, effectively assuming the stellar activity is not present.

Table 6. Model Parameters and Prior Distributions Used in Our RV Model That Considers the Transiting b and c Planets, as Well as the Recovered MAP Fit and MCMC Posteriors for the J2 Kernel

Parameter [units]	Initial Value (P0)	Priors	MAP Value J2	MCMC Posterior J2
Pb (days)	6.3980661	locked a	⋯	⋯
TCb (days)	2458517.68971	locked a	⋯	⋯
eb	0.294	${ \mathcal U }(0,1)$; ${ \mathcal N }({P}_{0},0.13)$	0.0.33	${0.29}_{-0.09}^{+0.09}$
ωb	130π/180	${ \mathcal U }({P}_{0}-\pi ,{P}_{0}+\pi )$; ${ \mathcal N }({P}_{0},45\pi /180)$	4.11	${3.94}_{-0.49}^{+0.26}$
Kb (m s−1)	10	${ \mathcal U }(0,\infty )$	6.57	${6.98}_{-1.82}^{+1.76}$
Pc (days)	18.8805	locked a	⋯	⋯
TCc (days)	2458533.593	locked a	⋯	⋯
ec	0.093	${ \mathcal U }(0,1);$ ${ \mathcal N }({P}_{0},0.13)$	0.19	${0.19}_{-0.10}^{+0.10}$
ωc	− 190π/180	${ \mathcal U }({P}_{0}-\pi ,{P}_{0}+\pi );$ ${ \mathcal N }({P}_{0},45\pi /180)$	−3.87	$-{3.53}_{-0.063}^{+0.62}$
Kc (m s−1)	10	${ \mathcal U }(0,\infty )$	6.80	${6.64}_{-1.42}^{+1.36}$
γiSHELL (m s−1)	MEDIAN(RViSHELL) + π/100 b	${ \mathcal N }({P}_{0},100)$	4.22	${3.52}_{-4.62}^{+4.63}$
γPFS (m s−1)	MEDIAN(RVPFS) + π/100 b	${ \mathcal N }({P}_{0},100)$	−13.28	−${13.30}_{-17.63}^{+17.87}$
${\gamma }_{\mathrm{HIRES}}$ (m s−1)	$\mathrm{MEDIAN}({\mathrm{RV}}_{\mathrm{HIRES}})+\pi /100$ b	${ \mathcal N }({P}_{0},100)$	4.59	${7.07}_{-13.48}^{+13.68}$
ηP	12.03	${ \mathcal N }({P}_{0},0.13)$	12.02	${12.00}_{-0.0}^{+0.09}$
ηℓ	0.44	locked a	⋯	⋯
ητ	57.96	locked a	⋯	⋯
ησ,0	1	${ \mathcal J }(0.67,50)$	23.42	${27.87}_{-4.30}^{+5.58}$
ηλ	1.17	${ \mathcal U }(-1,2)$	0.61	${0.66}_{-0.30}^{+0.30}$

Notes. Initial values for orbital period, timing conjunction, eccentricity, and planet radius for both planets as well as stellar mass come from our ExoFAST analysis in Section 3.2.4, and the stellar activity model hyperparameters and prior distributions come from our SuperWASP light-curve analysis in Section 3.2.3.

^a Locked indicates the parameter is fixed. Gaussian priors are denoted by ${ \mathcal N }(\mu ,\sigma )$, uniform priors by ${ \mathcal U }$(lower bound, upper bound), and Jeffrey's priors by ${ \mathcal J }$(lower bound, upper bound). The +1 on the initial gamma values is in case the RVs are already median-subtracted, and to pseudo-randomize this initial parameter. ^b We want the initial value to be the median of the RVs for that spectrograph; the +1 is used in case the median is already zero, as Nelder–Mead solvers cannot start at zero.

Download table as:  ASCIITypeset image

The TOI 560 TESS light curves for Sectors 8 and 34 are shown in Figure 13, after subtracting off a cubic-spline regression fit for the out-of-transit stellar activity from the PDC-SAP light curves shown in Figure 10. The ExoFASTv2 analysis reveals clear transits of TOI 560 b and c with expected depths and transit times consistent with the TESS mission TOI 560 b and 560.02 candidates. Table 7 shows the median and confidence intervals for the model and derived planetary parameters. Figure 14 shows the phased TESS transits. Of particular note, we confirm a nonzero eccentricity for TOI 560 b at 4.7σ; the eccentricity posterior for TOI 560 c is consistent with zero (Figure 15). The statistical significance of the nonzero eccentricity detection for TOI 560 b is primarily driven by the photoeccentric effect using the Spitzer data (Dawson & Johnson 2012), as excluding this particular data set decreases the statistical significance on e_b > 0 to 1.6σ, but still with a similar nonzero median eccentricity. Finally, the Spitzer light curve is presented in Section 4.1.2, and ground-based PEST, NGTS, and LCO light curves are presented in Figures 26–28 in Appendix C.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 7.  ExoFASTv2 Planetary Parameters: Median Values and 68% Confidence Interval Created for TOI 560 b and c, ExoFASTv2 Commit Number 101011575.2p.2

Planetary Parameter	Units	Planet b	Planet c
P	Period (days)	${6.3980661}_{-0.0000097}^{+0.0000095}$	${18.8805}_{-0.0011}^{+0.0024}$
RP	Radius (RJ )	${0.253}_{-0.011}^{+0.014}$	${0.2433}_{-0.0096}^{+0.011}$
RP	Radius (RNep)	${0.74}_{-0.03}^{+0.04}$	${0.71}_{-0.03}^{+0.03}$
RP	Radius (R⊕)	${2.84}_{-0.12}^{+0.16}$	${2.73}_{-0.11}^{+0.12}$
Tc	Time of conjunction a (BJDTDB )	${2458517.68971}_{-0.00054}^{+0.00053}$	${2458533.593}_{-0.091}^{+0.041}$
TT	Time of minimum projected separation b (BJDTDB )	${2458517.68973}_{-0.00054}^{+0.00053}$	${2458533.593}_{-0.091}^{+0.041}$
T0	Optimal conjunction Time c (BJDTDB )	${2458703.23362}_{-0.00045}^{+0.00044}$	${2459251.0510}_{-0.0014}^{+0.0016}$
a	Semimajor axis (AU)	${0.05995}_{-0.00072}^{+0.00074}$	0.1233 ± 0.0015
i	Inclination (degrees)	${89.45}_{-0.50}^{+0.38}$	${89.61}_{-0.24}^{+0.26}$
e	Eccentricity	${0.294}_{-0.062}^{+0.13}$	${0.093}_{-0.066}^{+0.13}$
ω*	Argument of Periastron (degrees)	${130}_{-46}^{+30}$	−190 ± 130
Teq	Equilibrium temperature d (K)	${742.7}_{-7.0}^{+7.2}$	${517.8}_{-4.9}^{+5}$
τcir	Tidal circularization timescale (Gyr)	${44}_{-36}^{+69}$	${9500}_{-7100}^{+15000}$
K	RV semiamplitude (m s−1) e	${2.9}_{-1.9}^{+2.7}$	${1.5}_{-1.1}^{+1.9}$
RP /R*	Radius of planet in stellar radii	${0.03803}_{-0.00064}^{+0.00063}$	0.0379 ± 0.0011
a/R*	Semimajor axis in stellar radii	${19.03}_{-0.48}^{+0.51}$	${39.15}_{-0.99}^{+1.0}$
δ	${({R}_{P}/{R}_{* })}^{2}$	0.001446 ± 0.000048	$.{001433}_{-0.000082}^{+0.000086}$
δB	Transit depth in B (fraction)	${0.00192}_{-0.00032}^{+0.00059}$	${0.00186}_{-0.00029}^{+0.00056}$
δR	Transit depth in R (fraction)	${0.00216}_{-0.00047}^{+0.00080}$	${0.00208}_{-0.00043}^{+0.00076}$
${\delta }_{z^{\prime} }$	Transit depth in z' (fraction)	${0.00202}_{-0.00041}^{+0.00085}$	${0.00195}_{-0.00036}^{+0.00081}$
δ4.5μ m	Transit depth in 4.5 μm (fraction)	${0.0046}_{-0.0012}^{+0.0023}$	${0.0042}_{-0.0012}^{+0.0023}$
δTESS	Transit depth in TESS (fraction)	${0.00163}_{-0.00013}^{+0.00021}$	${0.0042}_{-0.0012}^{+0.0023}$
τ	Ingress/egress transit duration (days)	${0.003232}_{-0.000078}^{+0.00020}$	${0.00588}_{-0.00040}^{+0.00098}$
T14	Total transit duration (days)	0.0866 ± 0.0014	${0.1510}_{-0.0032}^{+0.0039}$
TFWHM	FWHM transit duration (days)	${0.0833}_{-0.0013}^{+0.0014}$	${0.1448}_{-0.0031}^{+0.0038}$
b	Transit impact parameter	${0.134}_{-0.093}^{+0.13}$	${0.26}_{-0.17}^{+0.18}$
bs	Eclipse impact parameter	${0.19}_{-0.13}^{+0.18}$	${0.26}_{-0.17}^{+0.13}$
Ts	Ingress/egress eclipse duration (days)	${0.00500}_{-0.00066}^{+0.00049}$	${0.00597}_{-0.00046}^{+0.00047}$
Ts,14	Total eclipse duration (days)	${0.127}_{-0.017}^{+0.012}$	${0.152}_{-0.015}^{+0.013}$
Ts,FWHM	FWHM eclipse duration (days)	${0.122}_{-0.016}^{+0.012}$	${0.147}_{-0.015}^{+0.012}$
Ts,2.5μ m	Blackbody eclipse depth at 2.5 μ m (ppm)	${1.57}_{-0.11}^{+0.12}$	${0.0537}_{-0.0059}^{+0.0065}$
Ts,5.0μ m	Blackbody eclipse depth at 5.0 μ m (ppm)	${26.8}_{-1.2}^{+1.3}$	${4.85}_{-0.36}^{+0.39}$
Ts,7.5μ m	Blackbody eclipse depth at 7.5 μ m	${61.4}_{-2.5}^{+2.7}$	${18.8}_{-1.2}^{+1.3}$
〈F〉	Incident flux (109 erg s−1 m−2)	${0.0629}_{-0.0054}^{+0.0034}$	${0.01600}_{-0.00078}^{+0.00071}$
TP	Time of periastron (BJDTDB )	${2458511.67}_{-0.43}^{+0.15}$	${2458517.0}_{-5.8}^{+4.0}$
TS	Time of eclipse (BJDTDB )	${2458520.15}_{-0.85}^{+0.83}$	${2458524.1}_{-1.7}^{+1.3}$
TA	Time of ascending node (BJDTDB )	${2458516.27}_{-0.49}^{+0.32}$	${2458528.85}_{-1.0}^{+0.70}$
TD	Time of descending node (BJDTDB )	${2458512.20}_{-0.25}^{+0.27}$	2458538.22 ± 0.80
Vc /Ve	Equivalent circular to measured eccentric velocity ratio	${0.791}_{-0.027}^{+0.029}$	${0.991}_{-0.048}^{+0.061}$
$e\cos {\omega }_{* }$	Eccentricity times cosine of the periastron angle	$-{0.18}_{-0.22}^{+0.20}$	$-{0.004}_{-0.14}^{+0.11}$
$e\sin {\omega }_{* }$	Eccentricity times sine of the periastron angle	${0.199}_{-0.069}^{+0.043}$	$-{0.003}_{-0.066}^{+0.047}$
d/R*	Separation at mid-transit	${14.2}_{-1.2}^{+1.1}$	38.6 ± 2.7
PT	A priori nongrazing transit prob	${0.0676}_{-0.0048}^{+0.0061}$	${0.0249}_{-0.0016}^{+0.0018}$
PT,G	A priori transit prob	${0.0730}_{-0.0051}^{+0.0066}$	${0.0269}_{-0.0017}^{+0.0020}$
PS	A priori nongrazing eclipse prob	${0.0437}_{-0.0024}^{+0.093}$	${0.0247}_{-0.0011}^{+0.0025}$
PS,G	A priori eclipse prob	${0.0472}_{-0.0026}^{+0.010}$	${0.0267}_{-0.0012}^{+0.0027}$

Notes. See Table 3 in Eastman et al. (2019) for a detailed description of all parameters.

^a Time of conjunction is commonly reported as the "transit time." ^b Time of minimum projected separation is a more correct "transit time." ^c Optimal time of conjunction minimizes the covariance between TC and period. ^d Assumes no albedo and perfect redistribution. ^e The recovered semiamplitudes for the planets in the ExoFASTv2 analysis are significantly smaller than the those recovered in Table 6 with our RV analysis that includes a stellar activity model, and thus would appear at first to be contradictory. However, they are not statistically significant detections in this ExoFASTv2 analysis, as the posteriors are peaked at 0 m s⁻¹, and the reported confidence intervals herein are consequently misleading. The corresponding 5σ upper limits are ∼15 and 10 m s⁻¹ for K_b and K_c , respectively.

Download table as:  ASCIITypeset image

The Spitzer light curve (Figure 16) shows a transit that is consistent with the expected timing, duration, and depth from TESS for TOI 560 b. This helps rule out false-positive eclipsing binary scenarios due to the lack of any observed chromaticity in the depth and shape of the observed transit. Additionally, the Spitzer light curve is at higher photometric precision and temporal sampling with respect to TESS given the larger aperture of Spitzer; there is less limb-darkening compared to visible wavelengths, and less photometric variations due to stellar activity. Therefore, these observations help further constrain the exoplanet parameters of TOI 560 b and refine the orbital ephemerides. In particular, the Spitzer light curve enables us to detect the eccentricity of TOI 506 b via the photoeccentric effect as deviating from zero at 4.7σ. This eccentricity is marginally detected in the TESS transits, but is not statistically significant without the Spitzer observations.

Zoom In Zoom Out Reset image size