Understanding Accretion Variability Through TESS Observations of Taurus

Interpreting the short-timescale variability of the accreting, young, low-mass stars known as Classical T Tauri stars remains an open task. Month-long, continuous light curves from the Transiting Exoplanet Survey Satellite (\textit{TESS}) have become available for hundreds of T Tauri stars. With this vast data set, identifying connections between the variability observed by \TESS and short-timescale accretion variability is valuable for characterizing the accretion process. To this end, we obtained short-cadence \TESS observations of 14 T Tauri stars in the Taurus star-formation region along with simultaneous ground-based, UBVRI-band photometry to be used as accretion diagnostics. In addition, we combine our dataset with previously published simultaneous NUV-NIR \textit{Hubble Space Telescope} spectra for one member of the sample. We find evidence that much of the short-timescale variability observed in the \TESS light curves can be attributed to changes in the accretion rate, but note significant scatter between separate nights and objects. We identify hints of time lags within our dataset that increase at shorter wavelengths which we suggest may be evidence of longitudinal density stratification of the accretion column. Our results highlight that contemporaneous, multi-wavelength observations remain critical for providing context for the observed variability of these stars.


INTRODUCTION
Classical T Tauri stars (CTTS) are young, low-mass stars that are accreting material from the protoplanetary disk that forms and evolves alongside the star. Accretion is a critical part of the star-disk evolution process because it produces much of the high-energy radiation field that permeates the inner regions of the disk. Accretion also provides a method of probing the inner regions of the disk, which traditional observational techniques lack the ability to resolve. For two recent reviews on accretion onto CTTS, see Hartmann et al. (2016) and Schneider et al. (2020).
Nearly all CTTS show some manner of variability. Sources of the variability include rotational modulation from hot/cool spots, changes in the mass accretion rate, stellar activity, and variable extinction from material in the surrounding disk or dust entrained within the stellar magnetosphere (e.g., Herbst et al. 1994;Cody et al. 2014). Simulations of the inner regions of disks have been able to explain some of the observed variability caused by changes in the mass accretion rate on hourand day-timescales via instabilities in the inner regions of the disk and overdensities in the accretion column powered by turbulence (e.g., Kurosawa & Romanova 2013;Robinson et al. 2021).
Despite the above progress, variability on minute-to hour-timescales and the driving forces behind it remain relatively unclear. There are only a handful of studies that exist with long-term, short-cadence (defined here to be 1-to 2-minute) light curves of accreting young stars (e.g., Siwak et al. 2018) . The Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015) has begun to provide near-continuous, high-precision light curves of CTTS which can probe variability that may have been missed by previous studies (e.g., Zsidi et al. 2022). However, it is currently unclear how much of the variability as observed through the TESS bandpass can be attributed to accretion, and how much is driven by other sources of variability. Sets of simultaneous measurements of the accretion rate are thus critical when interpreting TESS light curves for studies of accretion.

arXiv:2207.00058v1 [astro-ph.SR] 30 Jun 2022
Ground-based photometry can be used to measure the accretion rate by relating the excess emission above photospheric levels to the total accretion luminosities using empirical bolometric corrections (e.g., Gullbring et al. 1998). The accretion luminosity can then be transformed into an estimate of the accretion rate using the stellar mass and radius from stellar evolutionary models (e.g., Baraffe et al. 1998), and by assuming that the material is falling at the free fall velocity under the magnetospheric accretion paradigm.
Space-based UV spectra can also be used to measure the accretion rate. In particular, contemporaneous FUV -NIR Hubble Space Telescope (HST ) observations with the Space Telescope Imaging Spectrograph (STIS) have been a powerful tool in previous studies for characterizing the accretion shock (e.g., Ingleby et al. 2011Ingleby et al. , 2013Ingleby et al. , 2015Robinson & Espaillat 2019). The accretion shock models of  predict that the shape of the UV continuum excess is a function of the kinetic energy of the accretion flow. Under the assumption of the flow being in free fall, the kinetic energy flux is proportional to the density of the accretion column. Thus, the shape of the UV continuum can be used to break the degeneracy between accretion column density and surface coverage that is present for optical measurements of the accretion rate (e.g., Ingleby et al. 2013Ingleby et al. , 2015Robinson & Espaillat 2019).
To test the connection between the variability in TTS as observed by TESS to accretion, we present shortcadence TESS observations and simultaneous groundbased UBVRI photometry of 14 TTS in the Taurus star forming region. In addition, we compare our results to five previously published simultaneous HST NUV-NIR STIS spectra and present one new epoch for one member of the survey, GM Aur (see Espaillat et al. 2021). This object is a young K5 star with a transitional disk containing residual optically thin dust in the cavity (Calvet et al. 2005;Hughes et al. 2009;Andrews et al. 2011;Espaillat et al. 2011;Macías et al. 2018) and has been the focus of several other recent HST UV studies of accretion studies of accretion (Ingleby et al. 2015;Robinson & Espaillat 2019). GM Aur is an ideal target for studies of accretion because of its robust accretion rate, similarity to the young Sun, inner cavity, and moderate inclination of 52 • .77 +0.05 −0.04 (Macías et al. 2018), which limits the chance for disk occultation of the star. In §2, information about the three data sets and reduction processes are presented. In §3, we detail our modeling efforts and the steps taken to calculate accretion rates. In §4, we present the results of our analysis and we discuss the significance of these results in §5. We summarize our key findings in §6.
2. OBSERVATIONS 2.1. TESS TESS offers the opportunity to obtain continuous, (aside from a short gap due to downlinking at orbital perigee) month-long light curves for many TTS. The TESS spectral response curve is very broad, and covers roughly between 0.6 µm and 1.1 µm with a steep cutoff at the blue end. While the peak of excess flux from accretion is in the NUV, excess continuum emission extends into the optical and near IR where TESS is sensitive ). Short-cadence (2-minute) light curves of the 14 targets included in our sample were obtained in the TESS GO program 22216 (PI: C. Robinson) as a part of the TESS Sector 19 observing campaign (Nov. 28th -Dec 23rd, 2019). At the time that our objects were observed, two types of data products were available from TESS : Full Frame Images (FFI) and short-cadence observations. FFIs are wide-FOV images that cover a 24 • × 90 • strip of the sky (defined as a TESS data sector) binned to a cadence of 30 minutes aboard the spacecraft. Short-cadence observations are produced for specific individual stars and the data is binned to a 2 minute cadence. Our sample was formed by comparing the population of T Tauri stars from Andrews et al. (2013) with the TESS Sector 19 footprint and selecting from the remaining members those that we could recover at sufficient SNR in the short-cadence mode of TESS . None of these targets have been previously observed by K2, making these the first short-cadence, long-baseline observations of these targets. The TESS magnitudes for our sample range between 12 and 9. SU Aur was originally included in our sample, but because it landed on the edge of the TESS detector, we exclude it from our analysis.
To verify the quality of our 2-minute cadence observations, the 30-minute FFI were inspected using the Python package Lightkurve (Lightkurve Collaboration et al. 2018). Transient (presumably solar-system) objects occasionally cross through the field of view, but from visual inspection, these transient objects do not significantly impact the quality of the light curves. Because of the large size of the TESS pixels (21 ) contamination can be significant in crowded regions. From visual inspection of previous ground-based images of the sample, we note the stars in our sample are significantly brighter than any secondary sources within the extraction region. To further probe the degree of contamination from nearby stars, we applied a variety of test apertures to the FFI to produce light curves for comparison against the short-cadence light curve produced by the Science Processing Operations Center (SPOC) pipeline (Smith et al. 2012;Stumpe et al. 2014;Jenkins et al. 2016). From this analysis, we found only minor contamination from neighboring stars and that light curves from the FFI with these apertures closely match the short-cadence light curve produced by the SPOC pipeline. Given this, we adopt the SPOC pipeline reduction. We also monitored the background flux in the FFI to determine the length of time around the orbit gaps to mask and found the standard SPOC mask is suitable to avoid contamination. Erroneous observations (e.g., single-frame sudden brightening/dimming events) were identified by visual inspection and masked. The light curves can be seen in Figure 1. Comments on the TESS light curves for each object are included in the appendix (Section 8).

LDT
We obtained UBVRI optical photometry for each object using the Large Monolithic Imager (LMI) on the 4.3m Lowell Discovery Telescope (LDT) in Happy Jack, AZ. A set of UBVRI images were obtained for each target during each of six nights spread over Dec. 2019(specifically 2019-12-02, 2019-12-07, 2019-12-10, 2019-12-13, 2019-12-18, and 2019. These data were taken simultaneously with the TESS observations (with the exception of some of the observations on 2019-12-10 which occurred near the start of the TESS downlink gap). Exposure times for these observations varied depending on object brightness and weather conditions. Exposure times for B, V, R and I were in the range of 0.1 to 10 s and between 10 and 120 s for U. In addition, we obtained six 1-2 hr monitoring series of ∼ 30 s U-band exposures of GM Aur during the same nights. Standard flat fielding and bias subtraction were applied to the images before extraction of photometry. Differential photometry was performed using the ensemble method for inhomogeneous sets of exposures from Honeycutt (1992) using a custom Jupyter widget 1 and the astropy package photutils (Bradley et al. 2020). After our initial source identification process, we identify and exclude variable stars by searching for stars that lie above the locus of points when plotting mean instrumental magnitude against the standard deviation of the instrumental magnitude. Depending on the population of background stars in the field, we used between 3 and 287 comparison stars to perform our differential photom-1 We designed our differential photometry widget to be flexible, interactive, and easy to use, making it a powerful tool with both research and education applications. Living open source code and instructions on its use can be found at https://github.com/ connorrobinson/ensemble. A frozen release can be found here: Robinson (2022).
etry. Stars that only had a few (< 10) comparison stars, which may make their relative photometry slightly less reliable, include FM Tau (3), CY Tau (4), and CW Tau (4). For additional details on this method of inhomogeneous differential ensemble photometry, see Honeycutt (1992).
To transform our instrumental magnitudes to an absolute system, we measured the zero point using the standard field centered around GD 64 (Landolt 2013) on 2019-12-21 UT. That night was chosen in particular because our observations during that time have the smallest variance in atmospheric extinction. Uncertainty estimates for our instrumental magnitudes were calculated following the procedures of Honeycutt (1992), which includes both an estimate of the measurement uncertainty in the flux (arising from Poisson and read noise) and the uncertainty in the atmospheric extinction present in that that exposure. We then add this uncertainty in our instrumental magnitudes in quadrature with an uncertainty in the measurement uncertainty of our zero point (estimated through typical error propagation techniques) to obtain our final estimate of uncertainty in our photometry. A table containing our UBVRI apparent magnitudes for the sample is presented in Table 1. We show the results from our GM Aur U-band monitoring series in Table 2.

HST
We present six NUV -NIR (1700 − 10000Å) low resolution (R ≈ 500 − 1000) spectra of GM Aur using STIS aboard HST, with a 1-3 d separation between epochs. Analysis for five out of six of these observations was first presented in Espaillat et al. (2021). Here we focus on making connections between these spectra and the simultaneous TESS and contemporaneous LDT observations, as well as add a new epoch and an additional sub-exposure analyses. We refer to these observations as Epochs 9 -14 to remain consistent with the eight previous observations presented in Ingleby et al. (2015) and Robinson & Espaillat (2019). Epoch 9 (2019-12-03) is the newly added epoch. This observing strategy was designed to cover a full rotation period of GM Aur (6.1 d; Percy et al. 2010) and to be simultaneous with TESS. The HST visits were originally scheduled to be daily, but had to be rescheduled due to an HST instrument error causing the gap in observations between Epochs 9 and 10. Epochs 9, 10, 11, 12, and 13 were taken simultaneously with the TESS observations. Epoch 14 Figure 1. Short-cadence (2 minute) TESS light curves for the 14 objects in the sample. The photometry shown here was produced using the standard SPOC pipeline. The timestamps of the data are presented as Modified Barycentric Julian Dates (MBJD) in the Barycentric Dynamical Time (TDB) scale. Flares identified by-eye by their characteristic rapid rise and exponential decay morphology have been flagged in magenta, and single-point outliers have been flagged in cyan. The light curves are shown normalized by their median value (this median includes the flagged points). We exclude these flagged points in the time-lag analysis presented in § 4.4.2. The light curve classification for each object is listed in parentheses to the right of the object's name (see § 3.4 and Table 5). Note-Summary of our UBVRI photometry andṀ measurements obtained using the LDT during December 2019.
Observing times are reported in MBJD in the TDB time frame.Ṁ fixed is calculated directly from our assumed stellar parameters. A discussion on photometric uncertainties is included in § 2.2.Ṁ rand and the associated uncertainties were found by re-sampling measurement uncertainties, while those forṀ sys include re-sampling both measurement uncertainties and system parameters. The reported values forṀ rand andṀ sys are the 50 th percentiles, while the uncertainties correspond to the 16 th and 84 th percentiles. A complete machine readable table is available online. Note-The scatter in adjacent values ofṀ, which is typically on the order of 0.01 to 0.03 × 10 −8 M yr −1 , is a better estimate of true random uncertainty since it is set by the internal precision of our differential photometry. Note that some of that scatter may be caused by real variability. A complete machine readable table is available online.
(2019-12-10) fell into the data downlink gap during the perigee, and thus are contemporaneous rather than simultaneous with the TESS dataset. None of our HST observations occurred perfectly simultaneously with our LDT observations. Table 3 summarizes the six HST STIS observations of GM Aur. Figure 2 shows our six epochs of HST STIS spectra overlaid on the range of the eight previous HST STIS observations. These observations were obtained during the HST Director's Discretionary time GO program 16010 (PI: C. Robinson). We took the spectra using the 52" × 2" slit and the G230L grating with the NUV-MAMA detector and the G430L and G750L gratings with the CCD detector. All three spectral orders were taken sequentially within a single orbit, making the observations contemporaneous to within roughly 50 minutes. Moderate fringing occurs at the red end of the NIR spectrum. To correct this, we follow the standard procedure described in Goudfrooij & Christensen (1998) using a contemporaneous fringe flat taken alongside the observations during terrestrial occultation.
The NUV data were taken in the time-tag mode of the NUV-MAMA detector, which records individual photon arrival times, making it possible to break a single long exposure into several shorter sub-exposures. While most of the analysis in this work relies on the full exposure length for the NUV observations, we take advantage of this in §4.4.1 to study minute-timescale changes in the accretion behavior. The initial delay and exposure times for these sub-exposures were chosen such that they are exactly simultaneous with overlapping TESS exposures (see §2.4 for details on matching time frames between HST, TESS, and LDT).
To isolate emission arising from accretion, we require an accurate model of the emission arising from the nondisturbed stellar photosphere. Non-accreting WTTS are preferred as spectral templates due to significant chromospheric emission in young stars in the UV (see Ingleby et al. 2011). A NUV -optical HST STIS spectrum of the K5 (Luhman & Steeghs 2004) WTTS RECX 1 was used as a spectral template in the analysis of GM Aur. This spectrum was taken with the matching 52" x 2" slit and the G230L and G430L gratings using the CCD detector during the HST GO program 11616 (PI: G. Herczeg).

Transforming to a uniform time scale
Because this is a multi-observatory study with importance placed on simultaneity, the timestamps for each observation must be placed into the same time 3. ANALYSIS

Stellar parameters
With the exception of GM Aur, we adopt stellar parameters for all of our objects from Herczeg & Hillenbrand (2014). Those authors include a veiling continuum emission from accretion and extinction during their spectral classification process, which are both critical for accurate measurements of the accretion rate. Spectral types were converted into effective stellar temperatures by interpolating between the grid of WTTS presented in Tab. 4 of Herczeg & Hillenbrand (2014). Following that work, we include K8 as an intermediate spectral type between K7 and M0. We scaled the stellar luminosities from the distances that were assumed in that work to  (2019), which discuss epochs 1-3 and epochs 4-8 in more detail, respectively. Epochs 10-14 were first presented in Espaillat et al. (2021). † While the other exposure start times are in MBJD in the TDB timescale, the exposure time start for RECX 1 is listed in MJD in UTC because simultaneity is not required.
those derived from parallaxes from Gaia EDR3 (Gaia Collaboration et al. 2016 and obtained new stellar radii based on this correction. Stellar masses were obtained by interpolating between the Mesa Isochrones and Stellar Tracks (MIST) model grid (Choi et al. 2016;Dotter 2016). The updated stellar parameters for each object are presented in Table 4. As discussed by Garufi et al. (2019), the parallax from Gaia DR 2 for RY Tau is known to be inconsistent with previous measurements and is assumed to be erroneous. However, this appears to have been corrected with the release of Gaia EDR3 with a value of 138.2 pc, which we adopt. This new distance is roughly consistent with the previous Hipparcos measurement of 133 pc (esa 1997) which has been adopted by previous works. For GM Aur, we adopt the stellar parameters from Manara et al. (2014) to remain consistent with the analysis of the HST STIS spectra presented in Robinson & Espaillat (2019); Espaillat et al. (2019bEspaillat et al. ( ,a, 2021Espaillat et al. ( , 2022.

MeasuringṀ from U-band excess
Excess emission above photospheric levels from accretion in the wavelengths covered by the broadband U filter can be used to measureṀ through empirical relationships (e.g., Gullbring et al. 1998). For each object (except the LDT GM Aur monitoring efforts and HST observations, see Section 2.3), we approximate the photospheric contribution to the observed U-band flux for each night and object by scaling BT Settl photospheric models Allard (2014) using the stellar parameters (see Section 3.1) and convolving the resulting spectrum with the U transmission curve. After subtracting the photospheric emission, the excess U-band accretion luminosity is converted to a bolometric accretion luminosity using the empirical relationship from Robinson & Espaillat (2019). This is then converted into a mass accretion rate by assuming that the material is falling at free-fall velocities from an inner magnetospheric radius of 5R . This value of 5R arises from balancing the ram pressure from the Keplerian motion of the disk and the magnetic pressure for typical CTTS parameters (see Hartmann et al. 2016). Recent interferometric observations of the magnetosphere of TW Hya found similar results, with a measured of radius of 3.5R (Gravity Collaboration et al. 2020). Additionally, because of the weak scaling between accretion luminosity and magnetospheric truncation radius (L acc ∝ 1 − R R M ) due to the flow approaching the free-fall velocity, small differences in R M between objects will not strongly influence our results.
To estimate uncertainties on our reported mass accretion rates, we use a Monte Carlo approach. We assume that the uncertainty distributions for our parameters can be approximated as Gaussian probability distributions. Following Herczeg & Hillenbrand (2014), we adopt uncertainties of σ Av = 0.2 (appropriate for moderately veiled stars) and a 0.1 dex uncertainty for L . For effective temperature, we adopt an uncertainty of σ T eff = 100 which is approximately 1 spectral class. When sampling L and T , we apply a Gaussian prior on the age of the cluster of 2.0 Myr (Kenyon & Hartmann 1995) with a width of 1 Myr. The inclusion of this prior shifts the center of the re-sampled distributions of  (2019) L and T and thus the distribution ofṀ. We found that this effect was particularly noticeable for stars that Herczeg & Hillenbrand (2014) identified as being heavily veiled (e.g., CW Tau, DS Tau, FM Tau). While we have assumed that a single age is representative of the age of the Taurus star-forming region (e.g., Kenyon & Hartmann 1995;Luhman 2018), other authors suggest that it may contain older populations of stars (Kraus et al. 2017;Krolikowski et al. 2021). With this in mind, we explored an alternative age prior of the functional form of 1 2 1−erf( A−A0 √ 2ω ) , which results in roughly equal probability before encountering a drop-off with a characteristic width ω centered at A 0 . With A 0 = 15 Myr and ω = 2 Myr, we find good agreement between the mean of our sampled distribution and the the stellar parameters from Herczeg & Hillenbrand (2014). However, because the shift in sampled stellar parameters primarily occurs for objects that are heavily veiled (which should be an indicator of youth Hartmann et al. 1998), we suggest that the shifts in the L and T distributions that appear under the Gaussian prior should not be used as evidence for an older population, and are instead simply an indication of the difficulty in accurately measuring stellar parameters for strongly accreting sources. We thus chose to adopt the single-age Gaussian prior and also note that the ultimately this choice of prior does not strongly influence our interpretation of our results.
Additionally, we re-sample our U-band photometry, GAIA parallaxes, and the coefficients in the empirical linear relationship between L U and L acc from Robinson & Espaillat (2019) using their associated uncertainties. Values of M for each re-sampling of L and T are found by interpolating within the MIST model grid Choi et al. (2016); Dotter (2016). When sampling stellar masses lower than 0.09 M , we instead interpolate within the Baraffe et al. (2015, hereafter BHAC) evolutionary tracks. From this analysis, we found that typical uncertainties are on the order of 0.2 dex of the reporteḋ M values. We refer to the 50 th percentile of this distribution asṀ sys . We do not report values ofṀ sys for RY Tau because both its previously measured T , L values lie near the upper edge of the MIST model grid, and would require either extrapolation or rejection of points when sampling, each of which may introduce biases.
If we exclude uncertainties on A v , T eff , D, L (i.e., systematic effects that would be eliminated by studying an individual star), then typical uncertainties are instead on the order of 0.09 dex. We refer to these reported values asṀ rand . Note that our uncertainty estimates for this value are still significantly larger than what the scatter in our U-band photometry alone would suggest because we include the uncertainty on the relationship between L acc and L U . Finally, we calculate the mass accretion rate directly from our stellar parameters and data without resampling, which we refer to asṀ fixed . Each are reported in Table 1. In most cases, these values agree with each other. The largest exception is RY Tau in low accretion states, where it becomes difficult to measure the accretion excess against its bright continuum. We exclude negative accretion rates as being nonphysical, which shifts the 50 th percentile to higher values during low accretion states.
Herczeg & Hillenbrand (2014) could not resolve the four known binaries in our sample (UY Aur, Coku Tau 4, DD Tau, and FO Tau), and thus their stellar parameters arise from their combined spectra. While this should not introduce large systematic errors into our L acc measurements, it may affect our resulting measurement ofṀ if M /R varies significantly from that of our single star approximation. Under the assumption that accretion is occurring primarily onto the primary star, the most extreme deviation from our single star approximation will occur for two equally sized stars. Under the MIST/BHAC models, this results in a systematic multiplicative underestimation of the mass accretion rate for these objects of up to about 30%. We do not include this effect in the reported systematic uncertainties foṙ M sys . This effect does not significantly change the interpretation of our results.

MeasuringṀ from HST
We take a different approach for GM Aur because of our additional simultaneous HST observations and ground-based U-band monitoring. To probe the accretion columns of GM Aur, the mass accretion rate and surface coverage of accretion shocks with different densities were measured with the HST STIS spectra using the accretion shock models and fitting methods of Robinson & Espaillat (2019). These models are an updated ver-sion of the models of  and the modeling efforts for five out of six of these epochs of GM Aur were first presented in Espaillat et al. (2021). The model and fitting methods are briefly described here.
The structure and emission arising from the postshock region is solved under a given kinetic energy flux, F i = 1 2 ρu 3 where ρ is density and u is velocity. These models work under the assumption that the accretion flow is in freefall in the pre-shock region, making the kinetic energy flux directly proportional to the density of the column. Half of the emission is radiated toward the star, where it heats the underlying photosphere, and the other half irradiates the pre-shock region. The outgoing, reprocessed emission from these two regions is then summed and scaled by a multiplicative filling factor, f i , which is treated as a free parameter. This is repeated for each value of F i . A non-accreting WTTS, multiplied by a free parameter scaling factor, s, is used as a template for the non-accreting regions of the stellar photosphere. The scaling factor s is defined as the ratio between the brightest observation and the photospheric template as measured through the V filter. The emission from each component is then summed to produce an estimate of the continuum emission.
In this approach, the photosphere is assumed to remain constant between epochs, i.e., the value of s is multiplied uniformly by the photospheric template for each epoch. Treating s as a single free parameter requires that f i for all epochs and s be fit simultaneously. To identify models that fit the observations well in this high-dimensional parameter space, posteriors for f i , s, and model uncertainty w were sampled using a Markov Chain Monte Carlo approach. In particular, we adopt the affine-invariant Goodman & Weare (2010) sampling algorithm via the Python-based ensemble sampler emcee (Foreman-Mackey et al. 2013). The model uncertainty, w, is treated as a nuisance parameter and can be marginalized. Each parameter is evaluated in log-space (limiting possible values to positive quantities) and we impose priors such that s, w, and f i cannot be larger than 1 to remain physical. Additionally, we impose a Gaussian prior on s with a width of 0.1 centered on the value of s derived from the analysis of the Epochs 1-8 from Robinson & Espaillat (2019). It is important to note that the resulting reported posteriors do not include the systematic uncertainty arising from the stellar parameters and extinction. Robinson & Espaillat (2019) found an additional 0.04 dex systematic uncertainty in the fitted parameters is appropriate in most cases. For more discussion on the fitting technique and the model, see that work, and . We also adopt the scaled template when calculating mass accretion rates using our U-band monitoring data for GM Aur. When resampling the photosphere in our Monte Carlo simulation, we assume an uncertainty of 10% based on the results of Robinson & Espaillat (2019).

TESS Symmetry and periodicity
Q and M are statistical metrics that measure periodicity and symmetry around the mean of a light curve, respectively. We measured these metrics for each for object within our sample using their TESS light curves. The Q and M metrics were developed by Cody et al. (2014) and have been used to characterize several starforming regions (e.g., NGC 2264, Ophiuchus, Cody et al. 2014;Cody & Hillenbrand 2018). A recent study also presented Q and M values measured from K2 light curves of the Taurus star-forming region (Cody et al. 2022). We note that there is no overlap between our sample and objects presented in that work. These metrics can be used to separate objects into empirical variability classes. The classifications include burster (B), purely periodic (P), quasi-periodic symmetric (QPS), stochastic (S), quasi-periodic dipper (QPD) and aperiodic dipper (APD). In addition to these classifications, we also note the possibility of several other variability classes, including multi-periodic (MP), which would not be identified by these metrics. We adopt the regions of Q and M parameter space as defined by Cody & Hillenbrand (2018) to determine variability classifications for the objects. More specifics on the algorithms used to measure these metrics can be found in that work and Cody et al. (2014).

Measuring periods from TESS
One of the steps for measuring Q requires determining the stellar rotation period. Following Cody et al. (2014), we measured the period by interpolating the TESS data onto a uniform grid and autocorrelating the data. The first significant peak of the autocorrelation function is then identified. Next, we applied a Lomb-Scargle periodogram (Lomb 1976;Scargle 1982). The resulting periodogram was then weighted by a Gaussian centered at the peak of the autocorrelation function with a width of 30% of the peak period. This was found to more accurately recover rotational periods of models of accreting young stars ) when compared to the box-car scheme of Cody et al. (2014). We took a different approach for CoKu Tau 4, since by-eye inspection of its light curve reveals that it resembles a MP binary. Instead, we identified the two strongest periods directly using a Lomb Scargle periodogram without any additional weighting. We note that our reported periods for sources that exhibit very strong aperiodic variability may not be reliable (BP Tau, DD Tau, RY Tau, and FM Tau).

LDT color magnitude diagrams
We correct for interstellar reddening in our LDT images using the extinction law of Cardelli et al. (1989) (R V = 3.1). Using the distances from GAIA we convert our apparent magnitudes into absolute magnitudes. The extinction-corrected, color-magnitude diagrams for each object are shown in Figure 3. The objects are shown sorted by their bluest M U − M B slope with an offset for clarity. In Figure 3, we also include A v = 1 extinction vectors, which show the expected slope for changes induced by variability in extinction along the line of sight. If the observed variability was caused solely by extinction, we would expect our photometry to lie along those vectors. We find that for most of the objects included in our sample, the M U − M B , M U slope is inconsistent with changes purely in the extinction along the line of sight. This is expected, since T Tauri stars are known to have variable accretion rates, the effect of which is most apparent at shorter wavelengths.

ConnectingṀ from U-band excess to TESS
We leveraged our multi-wavelength coverage to search correlations between the TESS light curves and accretion. To do this, we identify the TESS exposures that are simultaneous with our U and I photometry. We then deredden our I-band photometry, and convert from magnitudes to flux (F I,λ ) by adopting standard zeropoints. We then make the rough approximation of scaling the mean of the simultaneous TESS exposures to the mean of the dereddened, I-band fluxes. We chose I band because both it and the TESS bandpass are centered near 8000Å. We then subtract off the mean scaled TESS flux and correct for the distance and compare the result to the accretion luminosity measured through the U-band excess in Figure 4.
Within Figure 4, we separate objects with variability classes (see Table 5) that are associated with accretion variability (in our sample QPS, B) from those with variability arising primarily from other sources such as occultation and spot modulation (in our sample, QPD, MP, P). However, accretion and these other drivers of variability can and do occur simultaneously, which results in chimeric light curves with features arising from multiple mechanisms. This is particularly relevant for QPS objects, for which previous work has suggested that multiple sub-categories may exist even within the broader QPS category (Cody & Hillenbrand 2018) that are a mix of accretion and non-accretion processes. This introduces some amount of uncertainty in the source of the variability in our TESS light curves and its connection to the inferred mass accretion rate. We do note that none of the regions flagged by-eye as stellar flares by their characteristic rapid rise and exponential decay in the TESS light curves overlapped with our LDT observations, so we can eliminate that specific source of astrophysical noise as a significant contaminant.
To test whether our inferred values ofṀ and the variability in TESS are related, we calculated Pearson correlation coefficients and their associated p-values and did a linear regression for each source. We find positive slopes for our fits for 10 out of 14 of the sources (but note that some of the objects display large p-values). For the objects that have negative slopes, only two have variability classifications that are typically associated with accretion variability (FN Tau, FO Tau). Both of these objects only have relatively small variations in the TESS light curve compared to their mean brightness (∼ 1% and 5%, respectively) and have large p-values. For the objects that do have positive trends, we find significant scatter amongst the fitted slopes. For cases within our sample that have variability driven by accretion and other sources (particularly some QPS sources), we would expect that both the observed TESS flux and the inferred mass accretion rate would be affected in part by some of the non-accretion related events. This could introduce a non-causal correlation between the two quantities. Previous work also identified differences in color slopes between stochastic accretors and dippers (Venuti et al. 2015), which would explain some of the scatter in the measured slopes between our inferred accretion rates and measured TESS fluxes. Thus, we suggest that while TESS is indeed capable of tracing some facets of the accretion process, it can be strongly influenced by other astrophysical sources noise. We discuss this further in Section 5.1.
One caveat to this analysis is that the TESS bandpass is significantly wider than the I bandpass, so differences in the spectral energy distribution between sources and the shape of the transmission curve may introduce systematic offsets. This does not affect the strength of the correlation coefficients for individual objects, but it may make comparing multiple objects more difficult. However, given the wide dispersion in fitted slopes, this effect is not large enough to strongly influence our interpretation. We explored this effect in part by repeating our analysis but instead normalize using our R-band photometry. Similar to I-band, we find large degrees of scatter in the measured slopes between L acc and the scaled TESS flux.

GM Aur
GM Aur was unique in our sample in that we have additional U-band photometry and contemporaneous HST STIS observations that we used to derive mass accretion rates. Here for the first time, we compare to mass accretion rates derived from our contemporaneous HST and U-band data. Figure 5 shows the TESS light curve and the HST and LDT measurements ofṀ. The top panel shows the TESS light curve, while the middle panel shows the derived mass accretion rates from the U-band excess.Ṁ measurements from the LDT U-band excess and those from HST that are close in time are in excellent agreement given the uncertainties in measuring the mass accretion rate (e.g., 2019-12-07 UT and Epochs 10 & 11).
While correlations are present between the mass accretion rate and the observed TESS flux during most of the individual nights of LDT observations, significant scatter is present when all of the observations are considered as a whole (see Figure 6). The HST observations mirror this and also display significant scatter across the multiple epochs of observations. We excluded the points in the TESS light curve that were flagged as single-point outliers in this analysis (see Figure 1 and we did not identify any flares in the light curve of GM Aur from our by-eye inspection).
The contemporaneous NIR -NUV coverage of the STIS spectrograph allows us to break the degeneracy between accretion column density and accretion shock surface coverage. Table 6 contains the fitted mass accretion rates, accretion shock filling factors, and veilings for HST Epochs 9 -14. The model with the maximum posterior probability is overlaid on the observations in Figure 7. Individual contributions from low-, medium-, and high-density accretion shocks are shown separated from the total. The breakdown for Epochs 10, 11, 12, 13, and 14 were presented by Espaillat et al. (2021). Here were present the distribution for Epoch 9 for the first time. We find that the results for this epoch are are in good agreement with previous measurements of GM Aur in a intermediate accretion state, with a shock surface coverage fraction of about 11% and a mass accretion rate of 0.9 × 10 −8 M yr −1 .

Sub-exposure Time-tag Analysis
The STIS NUV-MAMA detector records photon arrival timestamps while in the time-tag mode which can be used to break long, single exposures into multiple, shorter sub-exposures. By breaking our longer NUV GM Aur exposures into 120 s chunks and adding a delay to the first sub-exposure, we obtain exactly simultaneous NUV-TESS observations which we use to search for correlations on minute-timescales. Figure 8 shows the integrated NUV flux between 1700 − 3100Å plotted against the simultaneous TESS observations. The existence of a trend between HST and TESS within individual HST visits is tested by calculating the Pearson correlation coefficient, r and the associated p value for each visit. We find weak, positive trends between TESS and HST within most individual exposures, but note that the reported p values shows that the trends are not strongly statistically significant. We attribute this in part to random uncertainties which are comparable to the measured differences in flux between sub-exposures and the small number of sub-exposures (typically ∼ 10) obtained during each HST visit. Values of p and r for each visit are reported in Figure 8. Ultimately, the observational uncertainty in our measure-   Figure 6). ments prevents us from firmly establishing or ruling out a strong relationship between the emission in the NUV and in the TESS passband on these short timescales, but hints of such a correlation may be present.

Identifying time-lags on hour-to-day timescales
The models of  find that the excess emission produced by higher-density accretion columns tend to peak at shorter wavelengths when compared to that of lower-density columns. Previous work has found that some objects require both high and low density components to reproduce their observed emission (e.g., Ingleby et al. 2013). Using these multicolumn models to fit Epochs 10, 11, 12, 13 and 14 of the HST STIS spectra of GM Aur, Espaillat et al. (2021) identified a one day lag between the UV emission coming from the high-density columns and the optical emission from the low-densities columns (see Figure 5). This was presented as evidence for longitudinal density stratification in the accretion column. The peaks of the filling fac-tors for the low density shocks are also well-aligned with the peaks in the simultaneous TESS light curve, which is reasonable because the TESS bandpass is centered at red/infrared wavelengths where the emission from the low-density accretion columns is higher (e.g., ). This may also explain some of the scatter in the slopes between the excess TESS luminosity and the values ofṀ measured using the excess in the Uband between objects (Figure 4), since the U bandpass is more sensitive to the higher density accretion accretion column energy fluxes than the TESS bandpass. For a visualization of this longitudinal density structure, see Figures 3 and 4 in Espaillat et al. (2021) which present a model of the magnetosphere and a surface map of the accretion column density for a CTTS with such gradients.
We attempt to identify similar lags by crosscorrelating our TESS data with our UBVRI photometry using the Python-based package Stingray (Bachetti Top: Normalized TESS light curve of GM Aur. Points with overlapping HST /LDT measurements ofṀ are highlighted in color. The dashed magenta lines mark the middle of the HST observations across all panels. Middle:Ṁ from the LDT U-band photometry and HST plotted as a function of time. Note that the peak of the measuredṀ is offset from the peak of the TESS light curve. Bottom: The contribution toṀ from the three accretion shock model components with different energy fluxes for each HST epoch. We find that the peak of the accreted mass attributed to the low energy flux model (log F = 10, blue) is roughly co-located with a peak in the TESS observations. This is in contrast with the the totalṀ, which is somewhat misaligned with the peak in the TESS light curve. et al. 2021). We make the assumption that our missing UBVRI data can be considered Missing Completely At Random (MCAR), and impute the mean value for the missing values. Note that the UBVRI data is very sparse (only 6 observations compared to the ∼ 18000 TESS points in each light curve), so the lags identified by this method may not be unique. We exclude regions in the TESS light curve that we identified as having stellar flares or as single-point outliers (see Figure 1), but note that their inclusion does not significantly change our results. We smoothed our cross-correlation results using a 3rd-order Savitzky-Golay filter (Savitzky & Golay 1964). We then identified peaks in the smoothed function by applying a threshold that the minimum peak value must be greater than 0.2 times the maximum peak and separated by least least ∼ 1.5 days. This process of smoothing and thresholding reduces the number of detected peaks to only those that are most significant.
We show the results from this analysis for our U-band photometry in Figure. 9. Histograms of the identified peaks for each object with absolute value closest to 0 d are shown in Figure. 10. The respective median absolute lags for our UBVRI photometry are 0.802, 0.673, 0.443, 0.352, 0.365 d, while the respective rotational-phase lags are 0.145, 0.122, 0.085, 0.060, 0.064. This anticorrelation between wavelength and median lag supports the idea that longitudinal stratification of accretion column density may be common for CTTS. The ∼ 1 d lag between our U-band photometry and the TESS light curve found here for GM Aur is in good agreement with the lag between high-and low-density accretion columns identified by Espaillat et al. (2021). Because of our very coarse sampling, it is not surprising that we do not find a median lag of zero between our TESS observations and our I-band photometry.

TESS as a tracer of accretion
We identify positive correlations between L acc measured through U-band measurements and the TESS flux scaled using our I-band photometry (see Figure 4). Furthermore, we do not see strong positive correlations between objects with variability classifications not associated with accretion within our sample. We also identify Note-Mass accretion rates, filling factors for low (f 1E10 ), medium (f 1E11 ), and high density (f 1E12 ) accretion columns, and veiling at V band, r v , for GM Aur from the analysis of the HST STIS observations with the accretion shock models presented in  and Robinson & Espaillat (2019). The subscripts following the filling factors refer to the kinetic energy flux of the accretion column (which is directly proportional to the column density) in units of erg s −1 cm −2 . The value of each f for each epoch is the fractional coverage of the visible stellar surface by shocks with that energy flux. Note that the uncertainties presented here are the 16 th and 84 th percentiles of the posteriors (roughly analogous to 1σ) recovered from the MCMC analysis, and do not include systematic uncertainties (e.g., extinction and stellar parameters). Robinson & Espaillat (2019) found that systematic uncertainties on the order of ∼ 10% are typically appropriate in derived parameters from this type of analysis. Results for Epochs 10, 11, 12, 13 and 14 were first presented in Espaillat et al. (2021).
correlations between TESS andṀ during the hour-long monitoring sessions of GM Aur. However, when we examine the slopes of the linear regressions, the idea of a single, global relationship between TESS and L acc (analogous to that of L acc and L U ) begins to break down. This is seen both across our entire sample (see the bottom right panel of 4) and between separate nights of observations for GM Aur (see Figure 6). Likewise, on even shorter scales we find different slopes between TESS flux and integrated NUV flux (a proxy for the mass accretion rate) from the time-tag analysis with the HST data (see Figure 8 and §4.4.1, but note that the correlations are very weak in those observations. The results of the decomposition of the accretion columns into individual components presented in Table 6 and Figure 7 as well as the time-lags shown in § 4.4 also demonstrate that the TESS bandpass misses some of the information about accretion encoded in the shorter wavelengths. In short, we suggest that any global relationship between observed TESS flux and the mass accretion rate would be unreliable, but note that TESS flux does trace facets of the accretion behavior for individual objects particularly on shorter timescales. We suggest that some of the discrepancies between the mass accretion rate and the observed TESS flux may arise because of longitudi-nal stratification of accretion column densities (see Section 4.4, also Espaillat et al. 2021). TESS is primarily sensitive to low-density columns. The TESS bandpasses and the inferredṀ values are also sensitive to other sources of astrophysical variability such as disk obscuration or star spots (which can be large for TTS Bradshaw & Hartigan 2014) which could introduce similar effects. We suggest that this non-accretion-related variability may in part lead to the observed differences in relationships between the two quantities for individual sources (analogous to the scatter in color slopes for the sample presented in Venuti et al. 2015).
We searched for trends between the color slopes shown in Figure 3 and our measured values of Q and M, but we did not find obvious trends (likewise with our by-eye classifications). We also did not identify strong correlations between our mean mass accretion rates and Q and M (even when excluding the two dippers in our sample, CY Tau and RY Tau). While our sample size is limited, both results are in agreement with the wide spread of slopes and accretion rates that have been observed in similar studies with larger samples of CTTS (e.g., Venuti et al. 2015;Sousa et al. 2016

Connecting light curve variability metrics to literature data
The nearby, well-studied nature of our sample allows us to make connections that have not been possible in previous space-based monitoring campaigns (e.g., Cody et al. 2014;Cody & Hillenbrand 2018). Here we compiled ancillary literature data when available. Figure 11 plots the 14 objects in Q-M space and indicates their variability classification as well as their disk morphology as measured from IR/mm observations (see Table 4). While we note that our sample is small, we can make a few observations. None of the disks that are classified as B are transitional or debris disks, which may have limited material in the innermost regions compared to full or pre-transitional disks. Of the three transitional disks in our sample, Coku Tau 4 is classified as MP, GM Aur as QPS, and RY Tau as QPD.
Disk inclinations can be measured through resolved mm studies of disks. The right panels of Figure 12 show disk inclinations reported from several studies in the literature (see Table 4) and the Q and M metrics mea-sured from our TESS data. As noted by previous authors (e.g., Cody & Hillenbrand 2018), distinct trends between these variability classifications are difficult to identify. However, we find a weak inverse trend between disk inclination and Q for objects with variability classifications associated with accretion (QPS and B). We find little evidence for a strong correlation between M and i for these objects. Both are in good agreement with the results of Robinson et al. (2021), which showed that for accreting objects, Q varied across the entire range of inclinations while M was only a strong function of i when viewed nearly edge-on (i 70 • ).
One interesting object to note is CY Tau, which we classify as quasi-periodic dipper. The dips in this object are narrow and deep, and occur at a period similar to the apparent rotation period of the star. This may suggest an disk warp induced by the stellar magnetic field in the inner regions of the disk (perhaps analogous to the geometry of V354 Mon presented by Schneider et al. 2018). However, the disk inclination as traced by mm observations is unambiguously close to face-on (i = 32 +1 −1 • , Simon et al. 2017). This suggests that the stellar inclination and/or the inner disk may be strongly misaligned with mm measurements of the disk inclination. This is something that has been observed for other dippers, including the prototypical example AA Tau (Loomis et al. 2017). Misalignment of inner disks has been suggested to be a possible indicator of companions  or planetary mass objects (Zhu 2019).

Comparing BP Tau to simulations
The simulations of Robinson et al. (2021) made predictions about what system parameters are important for setting the morphology of light curves for young stars with rotation-axis-aligned magnetic field topologies. The parameters tested in that work include the stellar mass, the co-rotation radius, the ratio between octupole and dipole magnetic field coefficients, inclination, and the turbulence mach number in the inner disk. Within our sample, BP Tau is the only non-periodic object within our sample that has a published measurement of the magnetic field topology (Donati et al. 2008). Its magnetic field consists of a 1.2 kG dipole and 1.6 kG octupole both tilted slightly (∼ 10 • ) with respect to the stellar rotation axis with a relatively small toroidal component. Coincidentally, the ratio of the octupole to dipole field strengths and the other system parameters of BP Tau are relatively similar to one of simulations of Robinson et al. (2021) (the model in question assumes a turbulent mach number of 0.1 in the inner disk). We find that our measured Q and M values are consistent with predictions from that simulation within their model  Table 6. Filling fractions for this model for each epoch are included in the legend. The light grey bars mark regions that were masked during the fitting procedure due to line emission.
uncertainties. While potentially interesting, magnetic topology measurements exist for only one of our accreting targets and we caution against further conjecture at this point. Magnetic topology observations do exist for V819 Tau Donati et al. (2015), but that source is a WTTS with a magnetic field composed primarily of a dipole that is significantly tilted with respect to the rotation axis (∼ 30 • ), so the simulations of Robinson et al. (2021) are likely less applicable in that case.

Comparisons to previous GM AurṀ measurements
The values of the mass accretion rates of GM Aur measured from the HST data and the LDT are generally consistent with the previously observed range during Epochs 1-8, which spans from 0.6 − 2.0 × 10 −8 M yr −1 (Ingleby et al. 2015;Robinson & Espaillat 2019). The observed range of NUV emission for Epochs 1-8 is shown in Figure 2. When compared to the other previously an-alyzed epochs, the inferred distribution of accretion column densities during Epoch 9 is typical. The range of veiling values estimated from the s parameter for Epochs 9-14 is comparable to that of Epochs 1-8.
While the reported mass accretion rates from this work span roughly a factor of two, the model of azimuthal density structure across the accretion column from Espaillat et al. (2021) suggests that this may not be an accurate measurement of the global mass accretion rate but instead primarily caused by changes in the viewing geometry. In contrast, Robinson & Espaillat (2019) identified an accretion burst during Epoch 7 that remains remarkable. The peak accretion rate during that event was ∼ 20% higher than the peak accretion measured from our LDT photometry, and ∼ 100% higher than the mean accretion rate. Critically, the inferred distribution of accretion columns during this epoch were also distinct when compared with all other 13 epochs, with the largest contribution to the total mass accretion Epoch 9: r = 0.30, p = 4.3e-01 Epoch 10: r = 0.18, p = 6.2e-01 Epoch 11: r = 0.38, p = 3.1e-01 Epoch 12: r = 0.37, p = 2.9e-01 Epoch 13: r = -0.08, p = 8.5e-01 Figure 8. Integrated NUV flux from 120s sub-exposures of the NUV HST STIS observations of GM Aur plotted against the simultaneous TESS observations. The Pearson correlation coefficient r and its p value for each visit are included in the legend. We find hints of weak positive trends during four out of five epochs, but note that these trends are not strongly significant based on their associated p values. We find no correlation between multiple epochs. rate being from high-density columns. This event apppears to be a true increase in the global accretion rate facilitated by a a particularly dense accretion column, which supports the idea that that large density gradients may be present in the inner disk. This accretion burst with the additional context of rotational modulation of complex azimuthal density structures provided by Epochs 9 -14 demonstrates the strength of multi-wavelength monitoring for better understanding the structure and volatility of the inner disk.

SUMMARY
We observed 14 objects using the short-cadence mode on TESS for 26 d. During that time, we also obtained contemporaneous UBVRI photometry with the LDT on six nights for each object and monitored GM Aur over hour-timescales in U. Additionally, six epochs of HST NUV-NIR STIS data were obtained simultaneously for GM Aur during the TESS observations. With the TESS data, we estimate rotational periods and determined empirical variability classifications through the use of previously developed statistical metrics and by-eye checks. With the LDT data, we measured mass accretion rates through the measured U-band excess, searched for correlations between L acc and readily available TESS data, and identified time-lags with our simulatenous TESS data. With the HST data, we present new analysis for one epoch (the other five epochs were presented in Espaillat et al. 2021), compare the measured mass accretion rate to those derived from our U-band photometry, and searched for connections to TESS on minutetimescales using time-tagged observations. Our primary findings are as follows: 1. Some of the variability present in the TESS light curves can be linked to accretion.
2. However, a single global relationship between TESS flux and L acc is not readily apparent. This is demonstrated by the differences in slope between different objects, and even individual objects on different nights.
3. Our measured accretion rates for GM Aur from our U-band photometry agree very well with those derived from detailed accretion shock modeling of contemporaneous NUV-NIR measurements from HST.
4. We found some additional evidence for longitudinal density stratification of accretion columns by identifying time lags within our dataset that increase at shorter wavelengths.
5. We identified CY Tau as a face-on dipper, which may suggest the presence of a misaligned inner disk and/or magnetic field.
We conclude that TESS light curves do trace some of the accretion behavior on shorter timescales, but are not reliable tracers of the total accretion rate, especially over longer timescales. These results highlight the importance of obtaining simultaneous, multi-wavelength observations when interpreting TESS light curves. This will become increasingly important as TESS light curves continue to become available for CTTS.

ACKNOWLEDGEMENTS
This work was supported by HST grant GO-16010, NASA grant 80NSSC19K1712, and NSF Career grant AST-1455042. Support for program #16010 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.
The authors thank the members of Dr. Connor Robinson's Ph.D. defense committee, Dr. Philip Muirhead, Dr. Merav Opher, Dr. Daniel Clemens, and Dr. James Owens, whose comments greatly improved this document during the preparation of his Dissertation. We also   Figure 10. Histograms of the time lag, absolute time lag, phase lag, and absolute phase lag for the peak in the cross-correlation function regarded as significant by our criteria whose absolute time lag is closest to 0 lag (see the gold vertical lines in Figure 9). We found that the median absolute time lag and the median absolute phase lag tend to increase at short wavelengths, which may indicate longitudinal density stratification of the accretion shock. thank Dr. Nuria Calvet, and our anonymous reviewer whose insightful comments improved this work. This paper includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST). Funding for the TESS mission is provided by NASA's Science Mission directorate.
Based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. These observations are associated with program #16010.
These results made use of the Lowell Discovery Telescope (LDT) at Lowell Observatory. Lowell is a private, non-profit institution dedicated to astrophysical research and public appreciation of astronomy and operates the LDT in partnership with Boston University, the University of Maryland, the University of Toledo, Northern Arizona University and Yale University. The Large Monolithic Imager was built by Lowell Observatory using funds provided by the National Science Foundation (AST-1005313).
This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www. cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www. cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.
This research made use of Photutils, an Astropy package for detection and photometry of astronomical sources (Bradley et al. 2020).

APPENDIX: REMARKS ON INDIVIDUAL TESS
LIGHT CURVES Here, we discuss individual objects with regards to their observed variability classifications and stellar activity. Derived statistical metrics Q and M, periods, and empirical variability classifications can be found in Table 5. This section refers to events occurring in terms of days since MBJD 58815 in the TDB time scale (see Figure 1).
BP Tau Using Q and M, BP Tau is classified as a QPS object. Periodic behavior in the light curve is not immediately obvious by-eye. Some features that resemble both dips and accretion bursts are present. A bright flare was observed near 3.5d. Typical changes in brightness are on the order of 10%.
Coku Tau 4: We distinctly see periodic behavior with two frequencies and thus classify CoKu Tau 4 as an MP object rather than the QPS that Q and M would suggest. The amplitudes of observed changes for this object are small, on the order of 1% and are likely from rotational modulation of stable hot/cold spots. We do not see obvious evidence for flaring in this object. A Lomb-Scargle Periodogram (Lomb 1976;Scargle 1982) was used to identify both periods of the signals present in the light curve. CW Tau: This object exhibits very large changes in flux (up to 50%) with limited obvious periodicity, hence its classification as a QPS. Most of the changes resemble accretion events, but some events may be dips. We do not see evidence of flaring, but note that flares may be masked by the strong variability.
CY Tau: This object exhibits periodicity and shows evidence for both accretion bursts and dips, resulting in a classification of QPD. As discussed in § 5.2, this object has an inclination of (i = 32 +1 −1 • , Simon et al. 2017), which may suggest the presence of a misaligned inner disk and/or stellar inclination. A single flare was identified near 11.5d.
DD Tau: Some degree of periodicity is visible for this object. Most of the variability resembles accretion bursts, resulting in a B classification. Two flares were identified near 1.2d and 24.8d.
DE Tau: This object shows a periodic signal with accretion bursts and was classified as a QPS object. Two flares are present near 2.0d and 24.0d.
DS Tau: A possible period can be identified for this object, but it is somewhat masked by changes in the accretion rate. This object displays large changes in flux (up to 50%) and was assigned a B classification. FM Tau: Clear periodicity is difficult to identify for this object. The light curve is symmetric with both dimming and brightening events, resulting in a QPS classification. Two flares were observed, near 6.8d and 21.0d. Relative changes for this object are moderate (on the order of 25%).
FN Tau: While the changes in this object are quite small (on the order of 1%), a periodic signal is identifiable. We found that the TESS flux uncertainties from the SPOC pipeline were slightly overestimated compared to the observed scatter, causing the algorithm used to measure Q to fail. To account for this, we report an upper limit on Q for this object and note that both P and QPS are possible variability classifications. Two potential flares were identified near 6.0d and 21.7d.
FO Tau: A period is identifiable and accretion bursts are present in the TESS light curve, resulting in a B classification via Q and M. A single flare was identified near 25.0d. No quasi-periodicity is detected above a significance of 6 standard deviations. GM Aur GM Aur is unique within the sample because of the NUV -NIR HST STIS data obtained simultaneously with the TESS observations and the additional U-band monitoring sessions. More analysis on GM Aur is presented in Espaillat et al. (2021). From the TESS light curve, a period can be identified. On top of the periodic behavior, we see moderate increases, likely associated with changes in the accretion rate, resulting in a classification of QPS. No obvious flares were identified for this object.
RY Tau: Deep dips are visible during the first half of the TESS observations which contributes to the classification of QPD. During the second half, almost no period is visible. RY Tau has been previously observed to show dipping behavior thought to be associated with variations in the upper layers of the inner disk (Davies et al. 2020) which matches the observed behavior here. No flares were identified for this object.
UY Aur: A clear periodic signal can be identified, punctuated by flares and accretion events, resulting in a classification of B. We identified 5 events that resemble flares in this object, one of which appears to be at the beginning of a major accretion event.
V819 Tau: V819 Tau is a very periodic source with many flares. We identified 10 events that resemble the characteristic rapid rise and exponential decay of flares. We see very limited variability that resembles dipping/accretion.