TESS Hunt for Young and Maturing Exoplanets (THYME): A Planet in the 45 Myr Tucana–Horologium Association

2.1.1. TESS

2.1.2. Spitzer

2.1.3. WASP

2.2.1. Southern Astrophysical Research (SOAR)/Goodman

2.2.2. Archival Data from HARPS, UVES, and FEROS

2.2.3. South African Extremely Large Telescope (SALT)/High Resolution Spectrograph (HRS)

2.2.4. Network of Robotic Echelle Spectrographs (NRES)/Las Cumbres Observatory's (LCO)

We performed H-band integral field spectroscopy of both stars using the Gemini Planet Imager (GPI; Macintosh et al. 2014). As part of the GPI Exoplanet Survey (GPIES), DS Tuc B was observed on 2016 November 18 (program code GS-2015B-Q-500) and DS Tuc A was observed on 2016 October 22 (GS-2015B-Q500) under poor conditions, aborted after nine images, and then observed again under better conditions on 2016 November 18 (GS-2015B-Q-500). A high-order adaptive optics (AO) system compensated for atmospheric turbulence, and an apodized Lyot coronagraph was used to suppress starlight. Using 59.6 s integration times, we obtained 37.78 minutes of data with 149 of parallactic angle rotation for DS Tuc B and 4.97 minutes and 35.79 minutes of data with 50 and 152 of parallactic angle rotation for the two observations of DS Tuc A. All three data sets were reduced using the GPIES automated data reduction pipeline (Wang et al. 2018). Briefly, the data were dark subtracted, a bad-pixel correction was applied, the microspectra positions determined using an Argon arc lamp snapshot taken right before each sequence, 3D spectral data cubes were extracted using wavelength solutions derived from deep Argon arc lamp data, the images were distortion corrected, and fiducial diffraction spots (satellite spots) were used to locate the position of the star in each image. The stellar PSF was then subtracted from each image using both angular differential imaging (Marois et al. 2006) and spectral differential imaging (Sparks & Ford 2002) to disentangle the stellar PSF from any potential companions, and principal component analysis to model the stellar PSF (Soummer et al. 2012; Wang et al. 2015). The resulting image was then used to search for point sources (Section 4.2).

To better characterize the properties of each component we drew resolved photometry and astrometry for DS Tuc A and DS Tuc B from the literature. Specifically, we adopted optical B_T and V_T photometry from the Tycho-2 Survey (Høg et al. 2000), optical G, BP, and RP photometry from the Gaia second data release (DR2; Evans et al. 2018), near-infrared J, H, and K_S photometry from The Two Micron All Sky Survey (2MASS; Skrutskie et al. 2006), and mid-infrared W1, W2, W3, and W4 photometry from the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010). We also adopted proper motions and parallaxes for each component from DR2 (Lindegren et al. 2018), and J2000 positions from Tycho-2.

All photometry and astrometry from the literature used in our analysis is listed in Table 1.

Table 1.  Parameters of DS Tuc

Parameter	DS Tuc A	DS Tuc B	Source
Identifiers
TOI	200.01
Gaia DR2	6387058411482257536	6387058411482257280	Gaia DR2
TIC	410214986	410214984	Stassun et al. (2018)
2MASS	J23393949-6911448	J23393929-6911396	2MASS
HD	222259A	222259B	Cannon & Pickering (1924)
Astrometry
α R.A. (hh:mm:ss J2000)	23:39:39.49	23:39:39.27	Tycho-2
δ decl. (dd:mm:ss J2000)	−69:11:44.88	−69:11:39.51	Tycho-2
μα (mas yr−1)	79.464 ± 0.074	78.022 ± 0.064	Gaia DR2
μδ (mas yr−1)	−67.440 ± 0.045	−65.746 ± 0.037	Gaia DR2
π (mas)	22.666 ± 0.035	22.650 ± 0.030	Gaia DR2
Photometry
BT (mag)	9.320 ± 0.017	10.921 ± 0.060	Tycho-2
VT (mag)	8.548 ± 0.012	9.653 ± 0.030	Tycho-2
G (mag)	8.3193 ± 0.0010	9.3993 ± 0.0014	Gaia DR2
GBP (mag)	8.7044 ± 0.0049	9.9851 ± 0.0059	Gaia DR2
GRP (mag)	7.8137 ± 0.0036	8.7082 ± 0.0044	Gaia DR2
J (mag)	7.122 ± 0.024	7.630 ± 0.058	2MASS
H (mag)	6.759 ± 0.023	7.193 ± 0.034	2MASS
Ks (mag)	6.68 ± 0.03	7.032 ± 0.063	2MASS
W1 (mag)	6.844 ± 0.060	7.049 ± 0.081	WISE
W2 (mag)	6.748 ± 0.030	7.107 ± 0.037	WISE
W3 (mag)	6.777 ± 0.023	7.056 ± 0.029	WISE
W4 (mag)	6.668 ± 0.094	6.958 ± 0.119	WISE
Kinematics
Barycentric RV (km s−1)	8.05 ± 0.06	6.41 ± 0.06	This work
U (km s−1)	−8.71 ± 0.04	−9.27 ± 0.04	This work
V (km s−1)	−21.50 ± 0.04	−20.28 ± 0.04	This work
W (km s−1)	−1.53 ± 0.04	−0.47 ± 0.04	This work
Physical Properties
Spectral type	G6V ± 1	K3V ± 1	Torres et al. (2006)
Rotation period (days)	${2.85}_{-0.05}^{+0.04}$	unknown	This work
Teff (K)	5428 ± 80	4700 ± 90	This work
Fbol (10−8 erg cm−2 s−1)	1.2026 ± 0.017	0.542 ± 0.008	This work
M* (M⊙)	1.01 ± 0.06	0.84 ± 0.06	This work
R* (R⊙)	0.964 ± 0.029	0.864 ± 0.036	This work
L* (L⊙)	0.725 ± 0.013	0.327 ± 0.010	This work
Age (Myr)	45 ± 4	45 ± 4	Bell et al. (2015)
v sin i* (km s−1)	17.8 ± 0.2	14.4 ± 0.3	This work
i* (deg)a	>82°	⋯	This work

Note.

^aWith the convention i < 90.

Download table as:  ASCIITypeset image

Age: DS Tuc was one of the original systems used to define the Tuc-Hor moving group (then called the Tucanae association; Zuckerman & Webb 2000). The group has consistent age estimates based on isochronal fitting (45 ± 4 Myr; Bell et al. 2015) and the lithium-depletion boundary (40 Myr; Kraus et al. 2014). Here we adopt the age estimate from Bell et al. (2015).

Luminosity, effective temperature, and Radius: We first determined the bolometric flux (F_bol), T_eff, and angular diameter of DS Tuc A and DS Tuc B by fitting the resolved spectral energy distributions (SEDs) for each component with unreddened optical and near-infrared template spectra from the cool stars library (Rayner et al. 2009). A demonstration can be seen in Figure 2.

Zoom In Zoom Out Reset image size

Our SED-fitting procedure followed the technique outlined in Mann et al. (2015), which we briefly summarize here. Our comparison assumed zero reddening, as DS Tuc lands within a region near the Sun of low interstellar extinction (the Local Bubble; Sfeir et al. 1999). We simultaneously compared each template spectrum to our optical spectra from SOAR/Goodman (Section 2.2.1) and archival photometry (Section 2.4 and Table 1) using the appropriate system zero-point and filter profile (Cohen et al. 2003; Jarrett et al. 2011; Mann & von Braun 2015; Maíz Apellániz & Weiler 2018). Gaps in each template spectrum are filled with a BT-SETTL atmospheric model (Allard et al. 2012) using the model interpolation and fitting procedure described in Gaidos et al. (2014). This procedure simultaneously provided an estimate of T_eff based on the BT-SETTL model comparison to the observed spectrum. To compute F_bol, we integrated each template/model combination over all wavelengths.

We combined the derived F_bol with the Gaia DR2 distance (d) to determine the total luminosity (L_*) for each component star. We then calculated a stellar radius (R_*) from L_* and T_eff using the Stefan–Boltzmann relation. Errors on each parameter were assigned accounting for both the measurement uncertainties (e.g., in the photometry) as well as the range of possible templates (and their assigned T_eff values) that can fit the data. Final parameters and uncertainties are give in Table 1.

As part of our above procedure, the BT-SETTL model is scaled to match the photometry and template. Assuming perfect models, this multiplicative scale factor is equal to ${R}_{* }^{2}/{d}^{2}$ (Cushing et al. 2008), which provided another estimate of R_* given the Gaia DR2 distance. This technique is similar to the infrared-flux method (Blackwell & Shallis 1977). Radii derived from this scale factor are not totally independent of the above method, as they rely on the same photometry and models, but the latter technique is less sensitive to the assigned T_eff.

The first technique (Stefan–Boltzman) yielded a radius of 0.964 ± 0.029R_⊙, and the scaling (infrared-flux method) yielded a consistent radius of 0.951 ± 0.020R_⊙ for DS Tuc A. We adopt the former value for all analyses.

Mass: We estimated the masses of DS Tuc A and DS Tuc B by interpolating our luminosity estimates onto a modified isochrone grid from the Dartmouth Stellar Evolution Program (Dotter et al. 2008). These grids were adjusted to include the effects of magnetic fields and where the boundary conditions are applied, as described in more detail in Muirhead et al. (2014), Feiden & Chaboyer (2014), and Feiden (2016). We assumed solar metallicity, which is typical within a scatter of ∼0.1 dex for the young stellar populations in the Solar neighborhood (e.g., Spina et al. 2014 and references therein). We used both 40 and 50 Myr grids, using the spread to approximate errors introduced by the age uncertainty for the Tuc-Hor moving group. This interpolation yielded mass estimates of 1.01 ± 0.06M_⊙ for DS Tuc A and 0.84 ± 0.06M_⊙ for DS Tuc B. We considered these errors to be slightly underestimated, as systematic differences between model grids can exceed 10% at this age.

We used high-resolution data from HARPS, UVES, FEROS, SALT/HRS, and NRES/LCO to determine stellar RVs. We measured RVs by computing the spectral line broadening function (BF; Rucinski 1992) between DS Tuc A or B observations and a zero-velocity template. The BF represents the function that, when convolved with the template, returns the observed spectrum, carrying information on RV shifts and line broadening. Throughout the analysis we used the HARPS G2 binary mask as our template (e.g., Pepe et al. 2002). A Gaussian profile was fit to the BF to determine the stellar RV. In each case the BF is single peaked and smooth, indicating a contribution from only one star.

For each echelle order we computed a "first pass" BF, which was used to shift the observed spectrum near zero velocity. Orders that survive a 3σ-clipping algorithm were then stitched into three equal-length wavelength regions where the final BFs were computed. Our geocentric RV measurement and uncertainty were computed from the mean and standard deviation across these three regions. For archival observations that are provided as a single stitched spectrum, we created 150 Å wide initial "orders."

Finally, for each epoch we computed the BF for telluric absorption features using a continuum normalized A0 star as our template. These offsets were applied to our measured RVs. We have measured RVs for all archival data following the above procedure. While the HARPS pipeline provides more precise RVs, we preformed our own measurements to ensure the same zero-point corrections across different instruments. We found a ∼70 m s⁻¹ offset from the HARPS observations, similar to our measurement uncertainty, but recovered the same epoch-to-epoch variability. Our final RVs are corrected for barycentric motion and listed in Table 2. As noted in the introduction, DS Tuc B was previously identified as a binary based on its RV variability and the presence of two spectral components. Our spectra are inconsistent with DS Tuc B having two near-equal spectral type components; for both stars at each epoch, there is only one peak in the BF. While the previous work did not give sufficient information to test the proposed scenario of RV variability, we also do not see evidence for RV variations in excess of reasonable jitter levels for young stars in either star.

We measured the projected rotation velocity (v sin i_*) for DS Tuc A and B by fitting the BF with a rotationally broadened absorption line profile that has been convolved with the instrumental profile (Figure 3). We did not include additional broadening components such as microturbulence, though these factors should have minimal impact given the large v sin i_* values. For DS Tuc A, we find v sin i_* = 17.8 ± 0.2 km s⁻¹ using the HARPS spectra; the value is consistent when using SALT/HRS. From SALT/HRS observations of DS Tuc B, we measure v sin i_* = 14.4 ± 0.3 km s⁻¹.

Zoom In Zoom Out Reset image size

Rotation period: A photometric rotation period of 2.85 days for DS Tuc was previously reported by Kiraga (2012), and is clearly visible in both the TESS and WASP light curves. Using ground-based monitoring with the LCO, we associate this signal with DS Tuc A. We break the WASP light curve into four 200 day observing seasons and measure the rotation period and amplitude of variability in each season. The period is consistently 2.85 days with high variability in the semi-amplitude (2%–2.6%), but the phase shifts. The periodogram shows power at the period and the first harmonic, and no additional signals are seen that could be associated with DS Tuc B.

The TESS light curve of DS Tuc shows consistent rotational modulation with a semi-amplitude of 1%–2%. We modeled the TESS light curve with a Gaussian process (GP) using the celerite package from Foreman-Mackey et al. (2017). We used a kernel composed of a mixture of simple harmonic oscillators and a jitter term. Our GP model has a term to capture the periodic brightness modulation caused by spots on the stellar surface. This kernel is a mixture of two stochastically driven, damped harmonic oscillator models and has two modes in Fourier space: one at the rotation period of the star and one at half the rotation period. We initially included an additional damped harmonic oscillator with a period of 20 days to capture long-term trends in the light curve, but the fitted power of the signal indicated that it was unnecessary.

We used a Lomb–Scargle periodogram to identify the candidate rotation period. We then fit the stellar rotation model using least squares, iterating five times and rejecting 3σ outliers each pass. This served to remove smaller flares. We then started an MCMC fit using the affine-invariant MCMC implemented in the package emcee (Foreman-Mackey et al. 2013), beginning half the chains at the candidate rotation period identified in the periodogram, and a quarter each at half and twice the rotation period. We use 50 walkers and a burn-in of 5000 steps. We end the run when the autocorrelation timescale τ of all chains changes by <0.1 and the length of the chain is >100τ. We measure a rotation period of ${2.85}_{-0.05}^{+0.04}\,\mathrm{days}$.

Stellar inclination: Following the method detailed in Morton & Winn (2014), we combined the stellar rotation period measured from the TESS light curve, R_*, and v sin i_* measurements from above to estimate of the stellar inclination for DS Tuc A. Although this measurement is not very precise, this method can identify highly misaligned systems (e.g., Hirano et al. 2012) or be used for statistical studies of large planet populations (e.g., Winn et al. 2017). We determine an equatorial velocity of 17.13 ± 0.6 km s⁻¹, which is consistent with our spectroscopic measurement of v sin i_* = 17.8 ± 0.2 km s⁻¹. This corresponds to a 1σ lower limit on the inclination of i > 82° and a 2σ lower limit of i > 70°. We cannot distinguish between i < 90° and i > 90°, and so adopt the convention i < 90°.

We fit orbital parameters to the motion of the binary pair using a modified implementation of the Orbits for the Impatient (OFTI) rejection-sampling methodology described in Blunt et al. (2017). This implementation is publicly available on GitHub³⁴ and described further in Pearce et al. (2019).

Both objects have a well-defined Gaia DR2 astrometric solution, so we used the positions and proper motions of DS Tuc B relative to DS Tuc A in the plane of the sky. We used the RV measurements of Table 2 to interpolate a relative RV at the Gaia observation epoch of 2015.5. Relative separation and position angle (PA) measurements in the Washington Double Star Catalog (WDS) spanning 126 yr provide additional constraints on the stellar orbital motion. We performed a modified OFTI fit constrained by these measurements.

Previous implementations of OFTI have fit orbital parameters to astrometric observations spanning several epochs (e.g., Blunt et al. 2017; Cheetham et al. 2019; Pearce et al. 2019; Ruane et al. 2019). In this system, the precision of the Gaia solution for both objects allowed us to constrain five of the six position vector elements using just this single epoch, and we additionally have the astrometric measurements provided by WDS; only the line of sight position is not sufficiently constrained to contribute to the fit.

Table 3 displays the orbital parameters we determined for the stellar binary orbit. Figure 4 displays the orbital parameter distributions, joint credible intervals, and a selection of orbits plotted in the plane of the sky. The orbital semimajor axis is 157 < a < 174 au, with a closest approach of 59 < r_peri < 93 au (where the ranges are 1σ credible intervals). The stellar binary is constrained to be nearly edge-on (960 < i < 978), which is likely aligned with both the transiting planet's orbit and the primary star's spin axis.

Zoom In Zoom Out Reset image size

To search for companions in high-contrast imaging data from GPI, we forward modeled the PSF template of a hypothetical companion at each pixel in the image using the Forward Model Matched Filter technique (FMMF; Ruffio et al. 2017a). We then ran a matched filter with the template in an attempt to maximize the signal of a planet at that location in the image. The method accounts for the distortion of the signal due to the speckle subtraction step. The detection limits are expressed in terms of the flux ratio between the point source and the star and were calibrated using simulated point-source injection and recovery. The detection limits are set at six times the standard deviation of the noise in the final image, which is calculated in concentric annuli as a function of separation to the star. This detection threshold ensures a false-positive rate of less than one per 20 sequence of observations. The default matched filter reduction used for GPIES assumes a featureless spectrum, corresponding to hot planets, for the estimation of the point-source brightness. However, Ruffio et al. (2017b) showed that it can be used for the detection of stars without loss of sensitivity. We did not detect any candidate companions above our detection threshold in either data set.

We determined completeness to bound substellar companions using the method described in Nielsen et al. (2019). An ensemble of simulated companions were generated with full orbital parameters at a grid of semimajor axis and planet mass. The projected separation in arcseconds was then computed for each simulated companion given the distance to the star, and the contrast was calculated using the BT-Settl models (Baraffe et al. 2015), the age of the star (45 Myr), and the star's H magnitude. Each simulated companion was compared to the measured contrast curve, and companions lying above the curve were considered detectable. The same simulated companions were compared to multiple contrast curves, advanced forward in their orbits when observations are made at different epochs, as is the case for DS Tuc A. Outside a radius of ∼11, not all PAs fall on the detector; to compensate, we reduce the completeness beyond ∼11 using the fractional coverage as a function of radius.

The depth of search plots, giving completeness as a function of semimajor axis and companion mass, are given for DS Tuc A and B in Figure 5, along with the underlying contrast curves. There are two contrast curves at each epoch: a T-type curve assuming heavy methane absorption in the matched filter step (appropriate to companions as hot as ∼1100 K), and an L type contrast curve assuming a flatter spectrum appropriate to hotter brown dwarfs and stars. Overall, wider separation planets and brown dwarfs are ruled out at high confidence between ∼10–80 au, more massive than ∼5 M_Jup, around both A and B.

Zoom In Zoom Out Reset image size

Past AO observations of the DS Tuc system have been limited to an outer working angle of ρ ≲ 10'' (e.g., Kasper et al. 2007), leaving open the possibility of a hierarchical architecture with a very wide tertiary companion. The Gaia catalog reveals that there is one co-moving, codistant candidate Tuc-Hor member within <1 pc of the DS Tuc system, 2MASS J23321028-6926537, which was also suggested to be a candidate low-mass (spectral type around M5) member of Tuc-Hor by Gagné et al. (2015). However, given the very wide separation (ρ = 1.12 × 10⁵ au), this source is likely an unbound member of Tuc-Hor and not a bound companion of DS Tuc. There are no other candidate wide companions in Gaia DR2 within ρ < 1 pc and brighter than a limiting magnitude of G ∼ 20.5 mag, corresponding to a mass limit of M > 15M_Jup at τ = 40 Myr (Baraffe et al. 2015).

We tested the detectability of additional planets in the TESS sector 1 light curve of DS Tuc A using the notch-filter detrending and planet search pipeline of Rizzuto et al. (2017). For this process, we used the SAP light curve, which is not corrected for systematics using the cotrending basis vector method. This choice was made based on the presence of artifacts in the PDCSAP light curve, likely introduced by the presence of a strong stellar rotation signal. We first apply a deblending factor based on the TESS magnitudes for DS Tuc A and B and masked the time interval when fine guiding was lost. We then injected a set of model transiting planets synthesized with the BATMAN model of Kreidberg (2015) with orbital and size parameters chosen randomly. We used orbital periods of 1–20 days and planet radii of 1–10 R_⊕, and allowed orbital phase and impact parameter to take values in the interval [0,1]. Eccentricity was fixed to zero for this process, as it does not significantly influence detectability of a transit, but requires two additional variables over which to marginalize. We injected a total of 1000 trial planets for this test.

For each trial planet, we apply the notch-filter detrending pipeline, and then search for periodic signals with the BLS algorithm (Kovács et al. 2002), retaining signals with power-spectrum peaks above 7σ. We then set tolerance windows of 1% in both injected period and orbital phase to flag a trial planet as recovered. Figure 6 shows the completeness map for additional planets in the DS Tuc A system. Our search and the TESS sector 1 data for DS Tuc A are sensitive to ∼4 R_⊕ planets at period <10 days, and ∼3 R_⊕ at periods <6 days. At periods longer than 10 days, the time baseline and gaps due to the masked section significantly decrease sensitivity to transiting planets.

Zoom In Zoom Out Reset image size

The two components of DS Tuc are separated by 5'' and are not resolved by TESS,³⁵ which has a plate scale of 21'' pixel⁻¹ with 50% of light concentrated within one pixel (Ricker et al. 2014). We examined the measured centroid of the in-transit/out-of-transit difference image, which is calculated by the SPOC pipeline and included in the DV report (from the initial TPS run) that accompanied the alert. The DV report indicated that both DS Tuc A and B are contained within the 3σ confusion radius of the centroid (which we note is dominated by the 25 additional error added in quadrature to the propagated uncertainty) and the centroid analysis averages a transit signal and a spurious event. In the second TPS run, not included in the alert, the centroid offset is consistent with DS Tuc A at 2σ. We also analyzed the image centroids measured by the SPOC pipeline. The scatter in the centroid measurements is too large (≃1 millipixel per 4 hr bin) to detect the expected change in centroid position if the planet were to in fact orbit DS Tuc B (0.5 millipixel over a 3 hr transit). In summary, we found that the TESS data alone cannot conclusively identify which star hosts the transit.

Our Spitzer observations definitively show that the planet orbits DS Tuc A. A 4 × 4 pixel aperture placed on DS Tuc A revealed a transit signal that is consistent with that detected in the TESS data. An equal-sized or smaller aperture centered on DS Tuc B yielded no detectable transit signature (Figure 7).

Zoom In Zoom Out Reset image size

We simultaneously fit the TESS and Spitzer photometry using the transit fitting code misttborn.³⁶ misttborn was first used in Mann et al. (2016a) and has been used for a number of more recent works including Johnson et al. (2018). Briefly, we fit each system using emcee, and produced photometric transit models using batman (Kreidberg 2015), which is based on the transit model of Mandel & Agol (2002). In the MCMC we fit for the following planetary parameters: the planet-to-star radius ratio R_P/R_⋆ (assumed to be the same in all filters), impact parameter b, period P, and the epoch of the transit midpoint T₀. We fix eccentricity to zero. We also fit the following stellar parameters: linear and quadratic limb-darkening parameters for each filter (q₁, q₂) using the triangular sampling method of Kipping (2013), and the mean stellar density (ρ_⋆). We use Gaussian priors for the limb-darkening parameters, using the values in Claret & Bloemen (2011) and Claret (2017). We use uniform priors within physically allowed boundaries for the remaining parameters (most notably, we enforced $| b| \lt 1\,+{R}_{P}/{R}_{\star }$ in order to assure that a transit occurs while allowing grazing transits).

DS Tuc is a visual binary with a separation of ρ ∼ 5''. The TESS photometry is de-blended, but the deblending process may introduce errors, while our Spitzer aperture on DS Tuc A includes a small amount of contamination from DS Tuc B. We included as an additional MCMC parameter the contamination of the aperture by flux from other stars. This is implemented as a (fractional) flux added to the transit model to create a diluted model (LC_diluted) of the form

$$\begin{eqnarray}&&{\mathrm{LC}}_{\mathrm{diluted}}=\displaystyle \frac{{\mathrm{LC}}_{\mathrm{undiluted}}+C}{1+C},\end{eqnarray} \tag{ 1 }$$

where LC_undiluted is the model light curve generated from Batman and our GP model. This is comparable to the method used in Johnson et al. (2011) and Gaidos et al. (2016) to correct for flux dilution from a binary using the measured Δm between components. The key difference is that Equation (1) allows for flux to be subtracted from the model (C < 0) in the case of an overcorrection.

We set a Gaussian prior upon C of 0.00 ± 0.02 for TESS and 0.0217 ± 0.0050 for Spitzer. The width of 0.02 for TESS photometry was estimated based on uncertainties in the derived TESS magnitudes from the TIC. Section 2.1.2 describes how C for Spitzer was calculated from a model of the PSF.

The target displays substantial stellar variability in the TESS bandpass. In addition to the transit model described above, we utilized Gaussian process regression to account for stellar variability in the TESS photometry. This enables us to model the variations in the stellar flux occurring during the transit. Our kernel is a mixture of simple harmonic oscillators, the same as described in Section 3. We included the Gaussian process hyperparameters as fit parameters in our MCMC, and placed priors on those parameters based on the results of our stellar rotation modeling. The parameters are the stellar rotation period P_*, the amplitude A_GP of the primary signal at P_*, the relative strength of the secondary signal at P_*/2 (Mix_Q1,Q2), the decay timescales of the primary and secondary signals (Q1_GP, Q2_GP), and a jitter term to account for white noise (σ_GP).³⁷

We ran the MCMC chain with 100 walkers for 30,000 steps and cut off the first 5000 steps of burn-in, producing a total of 2.5 × 10⁶ samples from the posterior distributions of the fit parameters. The resulting fit is shown in Figure 1, and the best-fitting values are listed in Table 4.

Table 4.  Parameters of DS Tuc Ab

Parameter	Value
Measured Parameters
T0 (TJD)a	1332.30997 ± 0.00026
P (days)	8.138268 ± 1.1 × 10−5
RP/R⋆	0.05419 ± 0.00024
b	${0.18}_{-0.12}^{+0.13}$
ρ* (ρ⊙)	${1.7}_{-0.17}^{+0.07}$
q1,1	${0.284}_{-0.053}^{+0.055}$
q2,1	0.284 ± 0.051
q1,2	${0.0266}_{-0.0091}^{+0.0094}$
q2,2	${0.054}_{-0.013}^{+0.014}$
CTESS	${0.015}_{-0.017}^{+0.018}$
CSpitzer	${0.0208}_{-0.005}^{+0.0049}$
$\mathrm{ln}{P}_{* }$ (day)	${1.0606}_{-0.0098}^{+0.0102}$
$\mathrm{ln}{A}_{\mathrm{GP}}$ (%2)	$-{10.87}_{-0.12}^{+0.11}$
$\mathrm{ln}Q{1}_{\mathrm{GP}}$	${2.57}_{-0.37}^{+0.39}$
$\mathrm{ln}Q{2}_{\mathrm{GP}}$	${0.052}_{-0.026}^{+0.027}$
MixQ1, Q2	${0.15}_{-0.11}^{+0.26}$
σGP	−8.682 ± 0.013
Derived Parameters
RP(R⊕)	5.70 ± 0.17
a/R⋆	${20.35}_{-0.69}^{+0.29}$
i(°)	${89.5}_{-0.41}^{+0.34}$
δ (%)	0.2936 ± 0.0026
T14 (days)	${0.13235}_{-0.00039}^{+0.00049}$
T23 (days)	${0.11818}_{-0.00057}^{+0.00039}$
Tperi (TJD)a	1332.30997 ± 0.00026
g1,1	${0.3}_{-0.054}^{+0.055}$
g2,1	${0.228}_{-0.06}^{+0.066}$
g1,2	${0.0172}_{-0.0051}^{+0.0057}$
g2,2	${0.145}_{-0.028}^{+0.024}$

Notes. We report the median and 68% confidence interval for each parameter. Associated probability distributions for key parameters are shown in Figure 1.

^aTJD is TESS Julian Date, which is BJD-2,457,000.0. ^bAlthough we allow b to explore negative values, the absolute value of b is listed because positive and negative values are degenerate. Similarly, we cannot distinguish between i < 90° and i > 90° and adopt the convention i < 90°.

Download table as:  ASCIITypeset image

As we do not have dynamical (RV) confirmation of DS Tuc Ab, we use our other observations to show that the transits are caused by a real transiting planet. We consider and rule out the following false-positive scenarios.

• 1.  
  The transits are caused by instrumental artifacts or residuals from stellar variability: Though there are only two transits in the TESS data set with amplitudes much lower than the amplitude of starspot variability, we confirm the transits with Spitzer, conclusively ruling out an instrumental origin for the signal. The Spitzer detection of the transits in the near-infrared, at the predicted time and with the same depth as in TESS rules out stellar variability as an origin, which should be significantly lower in the Spitzer bandpass and should not produce periodic transit-like signals.
• 2.  
  DS Tuc A is an eclipsing binary: Our RV observations showed no variations large enough to be caused by a stellar companion. To test this, we generated 100,000 binaries with random (uniform) mass ratios, argument of periastron, phase, inclination, and eccentricty. The period was fixed at 8.138 days, and inclination was restricted ensure the companion eclipses (≳70°). We then compared each synthetic binary's predicted velocities to the observed velocities assuming an extra jitter term in the velocities of 100 m s⁻¹ (from stellar variability). All generated binaries down to 20M_J in mass were rejected at >5σ, and >99% were rejected down to 5M_J.
• 3.  
  Light from a physically unassociated eclipsing binary star or transiting planet system is blended with light from DS Tuc: Spitzer confirms that the transit signal detected toward DS Tuc A must originate from within a few arcseconds of the star. We detected no stars nearby DS Tuc in our GPI AO imaging, and other groups have previously detected no nearby stars in their own AO observations (Kasper et al. 2007; Vogt et al. 2015). Crucially, due to its proper motion, DS Tuc has moved over half an arcsecond with respect to stationary background sources between the different AO imaging epochs over the last decade, so we are able to definitively rule out background stars too close to DS Tuc A for GPI to resolve.
• 4.  
  Light from a physically associated eclipsing binary or planet-hosting companion is blended with light from DS Tuc A: For this to be true, DS Tuc A must have a binary companion close enough to escape detection by GPI (inside about 8 au) and bright enough to cause the transit signal that we see. The magnitude difference Δm between DS Tuc A and the faintest companion, which could contribute the transit signal, is given by

  $$\begin{eqnarray}&&{\rm{\Delta }}m\lesssim 2.5{\mathrm{log}}_{10}\left(\displaystyle \frac{{t}_{12}^{2}}{{t}_{13}^{2}\delta }\right)\end{eqnarray} \tag{ 2 }$$

  where t₁₂ is the duration of transit ingress/egress, t₁₃ is the transit duration from first contact (beginning of ingress) to third contact (beginning of egress), and δ is the observed transit depth (Vanderburg et al. 2019). Fitting the TESS light curve with MCMC, but without any constraints from the stellar parameters, yields Δm ≲ 2.4 (95% confidence). From a 45 Myr MIST isochrone (Choi et al. 2016; Dotter 2016) at solar metallically (provided in the TESS bandpass), this magnitude difference corresponds to a companion star with a mass >0.63 M_⊙.To place a dynamical upper limit on the mass of a companion, we perform a Monte Carlo simulation of companion orbits to DS Tuc A with randomly drawn isotropic inclinations, masses below 1M_⊙, and semimajor axes below 8 au (holding the eccentricity to zero). For orbits that produce semimajor amplitudes less than half the range of our RV observations (0.6 km s⁻¹), we find that we can exclude companion masses above 0.28 M_⊙ at 95% confidence. The large discrepancy between these mass limits excludes this scenario at high confidence.

Our observational constraints confidently rule out these false-positive scenarios, so DS Tuc Ab is almost certainly a genuine exoplanet.

With an age of τ ∼ 45 Myr, DS Tuc Ab is one of the few transiting planets with ages τ < 100 Myr, joining the planets K2-33b (David et al. 2016b; Mann et al. 2016b), V1298 Tau b (David et al. 2019) and AU Mic b (P. Plavchan et al. 2019, in preparation). At V = 8.5, DS Tuc A is the brightest of these transiting planet host stars, closely followed by AU Mic at V = 8.6.

Using photometry from TESS and Spitzer, we determined that DS Tuc Ab has a radius of 5.70 ± 0.17 R_⊕, placing it in the sparsely populated realm of super-Neptunes and sub-Saturns. The planet is young enough that it likely still contracting due to internal cooling and may also be losing mass; models from Bodenheimer et al. (2018) suggest that its radius will shrink by 5%–10% over the next few 100 Myr.

DS Tuc is a visual binary, and we find no evidence for additional massive companions in the system. While DS Tuc B has previously been suggested to be a spectroscopic binary, we do not see two components in the spectrum of DS Tuc B at any observed epoch, a visual companion in high-contrast imaging data, or periodic RV variations at the precision of our data (200 m s⁻¹). The detection of planetary or substellar companions orbiting DS Tuc A exterior to DS Tuc Ab could indicate that dynamical interactions played a role in the present orbit of DS Tuc A; however, our high-contrast imaging data from GPI shows no companions with masses more than about 5M_Jup between 10 and 80 au.

The orbit of the stellar binary is likely to be closely but not perfectly aligned with both the orbit of the transiting planet and the spin axis of the planet-hosting star. We found a binary orbit inclination of 969 ± 09, a planetary inclination of ${89.5}_{-0.41}^{{+0.34}^{\circ}}$, and a stellar inclination of i > 82° (1σ limit). The latter two quantities use the convention of i < 90; however, i > 90 is equally likely. Although the PAs are presently unconstrained, the chance of all three having the similar inclinations by chance is small, suggesting that the three axes are in fact close to aligned. This is similar to the five-planet Kepler-444ABC system (Campante et al. 2015). Dupuy et al. (2016) found that the orbit of Kepler-444BC and the orbits of the planets around Kepler-444A have the same inclination angle, and suggested that the planets formed in situ in close orbits around Kepler-444A.

The stellar density that we determine from the transit fit differs from that which we calculate from the stellar parameters by 3σ. The most likely reason is either errors in the model-derived stellar mass, or a mild eccentricity (0.05 ≲ e ≲ 0.1). While our mass estimate has formal errors of ≃6%, predictions from different model grids can vary by ≃10%. Moderate eccentricities have been found for some other young planets, including two in the Hyades (Quinn et al. 2014; P. C. Thao et al. 2019, in preparation).

Due to the brightness of DS Tuc A, this system offers an exciting opportunity for detailed characterization of a young planet. Measuring the planetary mass would allow one to compare the planet's density to that of older planets. A distinct possibility is that mass estimates based on field-age planets represent an overestimate for DS Tuc Ab, given that the planet could still retain heat from its formation and might undergo future radius evolution as its atmosphere is sculpted by photoevaporative ultraviolet flux. While these processes would impact the planetary radius, they are not expected to have a substantial impact on the planetary mass.

The Chen & Kipping (2017) mass–radius relation, which is based on field-aged planetary systems, predicts a planetary mass of ${28}_{-13}^{+35}$ M_⊕. The expected RV semi-amplitude produced by DS Tuc Ab would then be ${9}_{-4}^{+11}$ m s⁻¹. As evidenced by the large error bars on the inferred planet mass, there are relatively few planets with sizes between Neptune and Saturn with measured masses; and the planetary mass–radius relation is poorly constrained for planets of this size.

Measuring the Rossiter–McLaughlin effect would determine the sky-projected angle between the stellar rotational and planetary orbital angular momentum vectors, and test our hypothesis that the stellar spin and planetary orbital axes are aligned. We estimate the RV amplitude due to the Rossiter-McLaughin effect using the relation ΔRV ≃ 0.65 v sin i_* ${\left(\tfrac{{R}_{P}}{{R}_{* }}\right)}^{2}\sqrt{1-{b}^{2}}$ (Gaudi & Winn 2007), finding a predicted amplitude of 32 m s⁻¹. Combining a spin–orbit misalignment measurement from Doppler Tomography (e.g., Johnson et al. 2017) or the Rossiter–McLaughlin effect (e.g., Narita et al. 2010) with our measurement of i_* from the rotation period and v sin i_*, one could measure full 3D spin–orbit misalignment ψ. DS Tuc Ab joins a small number of planets where such measurements are possible.

Measuring RV signals on the scales noted above would be well within reach of current high precision RV instruments, but stellar activity poses a major challenge (e.g., Saar & Donahue 1997; Paulson et al. 2004). DS Tuc A is a very magnetically active star, with $\mathrm{log}{R}_{{HK}}^{{\prime} }=-4.09$ (Henry et al. 1996). For stars like DS Tuc A, the stellar activity signal on many-day timescales (i.e., over many stellar rotation periods) is expected to be 100–200 m s⁻¹ based on the sample of active stars monitored with Keck by Hillenbrand et al. (2015). While a jitter of this level would seem to preclude RV measurements of the planetary signal, stellar activity signals can be mitigated by simultaneously modeling the activity and planetary signals using, e.g., Gaussian processes, a process which would be aided by our knowledge of the star's photometric variability (e.g., Haywood et al. 2014; Rajpaul et al. 2015; López-Morales et al. 2016). It is not clear how well the activity signal can be modeled and removed in an intensive RV campaign to measure a planet's mass or Rossiter–McLaughlin effect.

We investigate prospects for atmospheric characterization with the James Webb Space Telescope (JWST) by computing its transmission spectroscopy metric using Equation (1) of Kempton et al. (2018). We assume zero albedo and full day-night heat redistribution to estimate an equilibrium temperature for the planet of 850 K. We find a transmission spectroscopy metric is 264, which can be interpreted as the S/N with which its transmission spectrum is expected to be measured (assuming a cloud-free atmosphere) with a 10 hr observing program with the NIRISS instrument. This makes DS Tuc Ab an excellent target for observations with JWST. Finally, we note that it may be possible to detect the planetary exosphere, e.g., using He 10830 Å transit observations (Oklopčić & Hirata 2018; Spake et al. 2018).