# Twenty-five Years of Accretion onto the Classical T Tauri Star TW Hya

Gregory Herczeg, Yuguang Chen, Jean-François Donati, Andrea Dupree, Frederick Walter, Lynne Hillenbrand, Christopher Johns-Krull, Carlo Manara, Hans Moritz Günther, Min Fang, et al.

► To cite this version:

Gregory Herczeg, Yuguang Chen, Jean-François Donati, Andrea Dupree, Frederick Walter, et al. Twenty-five Years of Accretion onto the Classical T Tauri Star TW Hya. The Astrophysical Journal, 2023, 956 (2), pp.102. 10.3847/1538-4357/acf468. hal-04252101

HAL Id: hal-04252101

https://hal.science/hal-04252101v1

Submitted on 26 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

Twenty-five Years of Accretion onto the Classical T Tauri Star TW Hya

Gregory J. Herczeg (沈雷歌)<fnref n="1" />,<fnref n="2" />,<fnref n="3" />, Yuguang Chen (陈昱光)<fnref n="4" />, Jean-François Donati<fnref n="5" />, Andrea K. Dupree<fnref n="6" />, Frederick M. Walter<fnref n="7" />, Lynne A. Hillenbrand<fnref n="8" />, Christopher M. Johns-Krull<fnref n="9" />, Carlo F. Manara<fnref n="10" />, Hans Moritz Günther<fnref n="11" />, Min Fang (房敏)<fnref n="12" />,<fnref n="13" />, P. Christian Schneider<fnref n="14" />, Jeff A. Valenti<fnref n="15" />, Silvia H. P. Alencar<fnref n="16" />, Laura Venuti<fnref n="17" />, Juan Manuel Alcalá<fnref n="18" />, Antonio Frasca<fnref n="19" />, Nicole Arulanantham<fnref n="15" />, Jeffrey L. Linsky<fnref n="20" />, Jerome Bouvier<fnref n="21" />, Nancy S. Brickhouse<fnref n="6" />, Nuria Calvet<fnref n="22" />, Catherine C. Espaillat<fnref n="23" />, Justyn Campbell-White<fnref n="10" />, John M. Carpenter<fnref n="24" />, Seok-Jun Chang<fnref n="25" />, Kelle L. Cruz<fnref n="26" />,<fnref n="27" />,<fnref n="28" />, S. E. Dahm<fnref n="29" />, Jochen Eislöffel<fnref n="30" />, Suzan Edwards<fnref n="31" />, William J. Fischer<fnref n="15" />, Zhen Guo (郭震)<fnref n="32" />,<fnref n="33" />,<fnref n="34" />,<fnref n="35" />, Thomas Henning<fnref n="36" />, Tao Ji (纪涛)<fnref n="1" />,<fnref n="2" />, Jessy Jose<fnref n="37" />, Joel H. Kastner<fnref n="38" />, Ralf Launhardt<fnref n="36" />, David A. Principe<fnref n="11" />, Connor E. Robinson<fnref n="39" />, Javier Serna<fnref n="40" />, Michal Siwak<fnref n="41" />,<fnref n="42" />, Michael F. Sterzik<fnref n="10" />, and Shinsuke Takasao<fnref n="43" />

Abstract

Accretion plays a central role in the physics that governs the evolution and dispersal of protoplanetary disks. The primary goal of this paper is to analyze the stability over time of the mass accretion rate onto TW Hya, the nearest accreting solar-mass young star. We measure veiling across the optical spectrum in 1169 archival high-resolution spectra of TW Hya, obtained from 1998–2022. The veiling is then converted to accretion rate using 26 flux-

<figure><img src="image_3.png" /></figure>

<footnote n="1">¹ Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, People's Republic of China; gherczeg1@gmail.com</footnote>

<footnote n="2">² Department of Astronomy, Peking University, Beijing 100871, People's Republic of China</footnote>

<footnote n="3">³ Visiting astronomer, Department of Astronomy; California Institute of Technology; Pasadena, CA 91125, USA</footnote>

<footnote n="4">⁴ Department of Physics & Astronomy, University of California Davis, 1 Shields Avenue, Davis, CA 95616, USA; yugchen@ucdavis.edu</footnote>

<footnote n="5">⁵ IRAP-UMR 5277, CNRS & Université de Toulouse, 14 Avenue Edouard Belin, F-31400, Toulouse, France</footnote>

<footnote n="6">⁶ Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138 USA</footnote>

<footnote n="7">⁷ Department of Physics & Astronomy, Stony Brook University, Stony Brook NY 11794-3800, USA</footnote>

<footnote n="8">⁸ Department of Astronomy, California Institute of Technology; Pasadena, CA 91125, USA</footnote>

<footnote n="9">⁹ Department of Physics and Astronomy, Rice University, 6100 Main Street, Houston, TX 77005, USA</footnote>

<footnote n="10">¹⁰ European Southern Observatory, Karl-Schwarzschild-Strasse 2, D-85748 Garching bei München, Germany</footnote>

<footnote n="11">¹¹ MIT Kavli Institute for Astrophysics and Space Research, 77 Massachusetts Avenue, Cambridge, MA 02139, USA</footnote>

<footnote n="12">¹² Purple Mountain Observatory, Chinese Academy of Sciences, 10 Yuanhua Road, Nanjing 210023, People's Republic of China</footnote>

<footnote n="13">¹³ University of Science and Technology of China, Hefei 230026, People's Republic of China</footnote>

<footnote n="14">¹⁴ Hamburger Sternwarte, Gojenbergsweg 112, D-21029, Hamburg, Germany</footnote>

<footnote n="15">¹⁵ Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA</footnote>

<footnote n="16">¹⁶ Departamento de Física, ICEx-UFMG, Avenida Antonio Carlos 6627, 30270-901 Belo Horizonte, MG, Brazil</footnote>

<footnote n="17">¹⁷ SETI Institute, 339 Bernardo Avenue, Suite 200, Mountain View, CA 94043, USA</footnote>

<footnote n="18">¹⁸ INAF-Osservatorio Astronomico di Capodimonte, via Moiariello 16, I-80131 Napoli, Italy</footnote>

<footnote n="19">¹⁹ INAF-Osservatorio Astrofísico di Catania, via S. Sofia, 78, I-95123 Catania, Italy</footnote>

<footnote n="20">²⁰ JILA, University of Colorado, Boulder, CO 80309-0440, USA</footnote>

<footnote n="21">²¹ Université Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France</footnote>

<footnote n="22">²² Department of Astronomy, University of Michigan, 1085 South University Avenue, Ann Arbor, MI 48109, USA</footnote>

<footnote n="23">²³ Institute for Astrophysical Research, Department of Astronomy, Boston University, 725 Commonwealth Avenue, Boston, MA 02215, USA</footnote>

<footnote n="24">²⁴ Joint ALMA Observatory, Avenida Alonso de Córdova 3107, Vitacura, Santiago, Chile</footnote>

<footnote n="25">²⁵ Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Strasse 1, D-85748 Garching bei München, Germany</footnote>

<footnote n="26">²⁶ Department of Physics and Astronomy, Hunter College, City University of New York, 365 Fifth Avenue, New York, NY 10016, USA</footnote>

<footnote n="27">²⁷ Physics | Graduate Center of the City University of New York, 365 Fifth Avenue, New York, NY 10016, USA</footnote>

<footnote n="28">²⁸ Department of Astrophysics, American Museum of Natural History, 200 Central Park West, New York, NY 10024, USA</footnote>

<footnote n="29">²⁹ Gemini Observatory/NSF's NOIRLab, 950 N. Cherry Avenue, Tucson, AZ 85719, USA</footnote>

<footnote n="30">³⁰ Thüringer Landessternwarte, Sternwarte 5, D-07778 Tautenburg, Germany</footnote>

<footnote n="31">³¹ Five College Astronomy Department, Smith College, Northampton, MA 01063, USA</footnote>

<footnote n="32">³² Instituto de Física y Astronomía, Universidad de Valparaíso, Avenida Gran Bretaña, 1111, Casilla 5030, Valparaíso, Chile</footnote>

<footnote n="33">³³ Núcleo Milenio de Formación Planetaria (NPF), Avenida Gran Bretaña, 1111, Casilla 5030, Valparaíso, Chile</footnote>

<footnote n="34">³⁴ Centre for Astrophysics Research, University of Hertfordshire, Hatfield AL10 9AB, UK</footnote>

<footnote n="35">³⁵ Departamento de Física, Universidad Tecnicá Federico Santa María, Avenida España 1680, Valparaíso, Chile</footnote>

<footnote n="36">³⁶ Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany</footnote>

<footnote n="37">³⁷ Indian Institute of Science Education and Research (IISER) Tirupati, Rami Reddy Nagar, Karakambadi Road, Mangalam (P.O.), Tirupati 517507, India</footnote>

<footnote n="38">³⁸ Center for Imaging Science, School of Physics and Astronomy, and Laboratory for Multiwavelength Astrophysics, Rochester Institute of Technology, Rochester, NY 14623, USA</footnote>

<footnote n="39">³⁹ Department of Physics & Astronomy, Amherst College, C025 Science Center 25 East Drive, Amherst, MA 01002, USA</footnote>

<footnote n="40">⁴⁰ Instituto de Astronomía, Universidad Nacional Autónoma de México Ensenada, B.C., México</footnote>

<footnote n="41">⁴¹ Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, Eötvös Loránd Research Network (ELKH), Hungarian Academy of Sciences, Konkoly–Thege Miklós út 15–17, 1121 Budapest, Hungary</footnote>

<footnote n="42">⁴² CSFK, MTA Centre of Excellence, Budapest, Konkoly Thege Miklós út 15-17, H-1121, Hungary</footnote>

<footnote n="43">⁴³ Department of Earth and Space Science, Graduate School of Science, Osaka University, Toyonaka, Osaka 560-0043, Japan</footnote>

<footnote>Received 2023 April 23; revised 2023 August 3; accepted 2023 August 14; published 2023 October 12</footnote>
calibrated spectra that cover the Balmer jump. The accretion rate measured from the excess continuum has an average of $2.51 \times 10^{-9} M_{\odot} yr^{-1}$ and a Gaussian distribution with an FWHM of 0.22 dex. This accretion rate may be underestimated by a factor of up to 1.5 because of uncertainty in the bolometric correction and another factor of 1.7 because of excluding the fraction of accretion energy that escapes in lines, especially $Ly\alpha$. The accretion luminosities are well correlated with He line luminosities but poorly correlated with $H\alpha$ and $H\beta$ luminosity. The accretion rate is always flickering over hours but on longer timescales has been stable over 25 years. This level of variability is consistent with previous measurements for most, but not all, accreting young stars.

*Unified Astronomy Thesaurus concepts: Classical T Tauri stars (252); Stellar accretion disks (1579); Protoplanetary disks (1300); High resolution spectroscopy (2096); Variable stars (1761)*

*Supporting material:* data behind figure, machine-readable table

# 1. Introduction

The evolution of protoplanetary disks and the final outcome of any planet formation depends in part on how gas flows through the disk (see review by Manara et al. 2023). Since the flows within the disk are challenging to measure, we often infer the global flow rate by measuring accretion from the disk onto the star. The disk-to-star accretion rate appears to vary on all timescales and with a wide range of amplitude (see reviews by Hartmann et al. 2016; Fischer et al. 2023).

The accretion rate and flow properties can be measured using many emission lines and in continuum emission, both of which are produced as material flows from the disk along the stellar magnetosphere and crashes at the stellar surface. The energy of the infalling gas is initially deposited at the base of the accretion flow, which is shock heated to ~$10^6$ K (e.g., Calvet & Gullbring 1998; Lamzin 1998). The accretion shock itself occurs near the stellar surface and heats the surrounding photosphere to ~$10^4$ K (e.g., Drake 2005; Brickhouse et al. 2012; Bonito et al. 2014). Most of the accretion energy is reprocessed and escapes as hydrogen continuum emission from the heated photosphere (e.g., Calvet & Gullbring 1998), which may then be used to measure accretion rates after scaling by a bolometric correction (e.g., Valenti et al. 1993; Hartigan et al. 1995; Gullbring et al. 1998). The funnel flows and accretion shock also produce line emission (e.g., Muzerolle et al. 2000; Kurosawa et al. 2006; Donati et al. 2014), which may be converted into an accretion luminosity or accretion rate using correlations with the accretion continuum measurements (e.g., Natta et al. 2004; Fang et al. 2009; Alcalá et al. 2017).

Each of the observable diagnostics of accretion rate has its own advantages and disadvantages. Photometry covers large samples on short and long timescales, though changes in accretion are difficult to unambiguously distinguish from other phenomena, such as chromospheric flares, changing spot coverage fractions, and extinction (e.g., Bouvier et al. 1993; Cody et al. 2014; Hillenbrand et al. 2022). Flux-calibrated spectra covering the Balmer jump provide an instantaneous measurement of accretion, but with repeated observations for only a few objects (e.g., Robinson & Espaillat 2019). Multi-epoch measurements of veiling in optical high-resolution spectra provide a consistent set of accretion rates, assuming that the underlying photospheric emission remains constant (e.g., Johns-Krull & Basri 1997; Alencar et al. 2012). Variability may also be inferred from changes in line emission (e.g., Scholz et al. 2005; Costigan et al. 2014), although for one accreting young star, XX Cha, large changes in the accretion continuum did not result in significant changes in line luminosities (Claes et al. 2022). Some recent work has combined spectroscopy with extensive photometric monitoring

to interpret the photometric changes with more precision and across longer periods of time (e.g., Bouvier et al. 2020; Venuti et al. 2021; Fiorellino et al. 2022; Zsidi et al. 2022; Bouvier et al. 2023), with typical changes in accretion rate of ~0.3 dex.

In this paper, we contribute to these efforts to understand accretion variability by analyzing 1169 high-resolution spectra of the classical T Tauri star (CTTS) TW Hya, obtained over 25 years. As one of the closest and brightest CTTSs, TW Hya has been repeatedly observed because it is a remarkable object, a cornerstone for studies of accretion, and a common target of radial velocity searches for young exoplanets. We compile and analyze 1169 high-resolution optical spectra from CFHT/ESPaDOnS, MPG-ESO 2.2 m/FEROS, Magellan/MIKE, La Silla 3.6 m/HARPS, VLT/UVES, VLT/ESPRESSO, Keck/HIRES, McDonald 2.7 m/2coude, and CTIO 4 m/Cassegrain Echelle to conduct the largest analysis of veiling for a single star. The veiling measurements are then converted to accretion luminosity, guided by 26 flux-calibrated spectra, most obtained at low resolution. The paper is organized as follows: Section 2 describes the properties of TW Hya; Section 3 introduces the archival data and data reduction process; Section 4 describes the veiling decomposition technique and how veiling is converted into accretion rate; Section 5 describes our use of low-resolution spectra to convert veiling to accretion rate; Section 6 compares veiling measurements to emission line diagnostics of accretion; and finally Section 7 discusses the accretion variability and searches for timescales.

# 2. Properties of TW Hya

TW Hya was initially discovered as an emission line object by Henize (1976) and has since become the namesake of the young, nearby TW Hya Association (de la Reza et al. 1989; Kastner et al. 1997; Hoff et al. 1998). The disk of TW Hya is one of the few disks in the TW Hya Association (e.g., Weinberger et al. 2013; Luhman 2023), as it has survived for longer than the typical disk dissipation timescale (e.g., Hernández et al. 2008; Mamajek 2009; Fedele et al. 2010). The disk is massive and large, with a series of rings and gaps and an inner hole in micron-sized dust (e.g., Calvet et al. 2002; Bergin et al. 2013; Andrews et al. 2016; van Boekel et al. 2017).

The spectral type of TW Hya in the literature ranges from K6/K7 when measured in blue/optical wavelengths and M2 when measured in the near-IR (e.g., Webb et al. 1999; Yang et al. 2005; Vacca & Sandell 2011; Debes et al. 2013; McClure et al. 2013). This wavelength dependence is likely explained by spots that cover some of the visible surface (e.g., Debes et al. 2013; Gully-Santiago et al. 2017; Gangi et al. 2022). Since the TW Hya disk and magnetosphere are both viewed nearly pole-on (e.g., Qi et al. 2008; Donati et al. 2011), the visible spot coverage
is expected to be roughly steady, without significant rotational modulation (although some photometric and spectroscopic variability is detected and attributed to the spot, see discussion below). We adopt an intermediate spectral type of M0.5 and a corresponding temperature of 3810 K, based on analyses of low-resolution optical spectra (Herczeg & Hillenbrand 2014, see also, e.g., Sokal et al. 2018 and Venuti et al. 2019).

The parameters $T_{\text{eff}} = 3810$ K, $J = 8.217 \pm 0.024$ mag from the Two Micron All Sky Survey (Cutri et al. 2003), J-band bolometric correction from Pecaut & Mamajek (2013), an estimated J-band veiling of 0.1 (Fischer et al. 2011), and $d = 60.14 \pm 0.05$ pc (Brown et al. 2021) yield $\log L_{\text{phot}}/L_{\odot} = -0.54$ and radius $1.23 R_{\odot}$. This temperature and luminosity correspond to a mass of $0.87 M_{\odot}$ and an age of 9.8 Myr for the Somers et al. (2020) evolutionary tracks for stars with 50% spot coverage and $0.59 M_{\odot}$ and an age of 4.3 Myr for 0% spot tracks (the 0% tracks are similar to the results from the Baraffe et al. 2015 models). We adopt the $0.87 M_{\odot}$ mass, based on consistency for a few low-mass stars between mass estimates from the 50% spot models and dynamical masses from the rotation of their disks (Pegues et al. 2021). For the 50% spot models and corresponding bolometric corrections, a temperature of ~4100 K would lead to a mass of $0.96 M_{\odot}$ and an age of 16 Myr; a lower temperature of ~3600 K leads to $0.70 M_{\odot}$ and 8.5 Myr.

The dynamical mass of $0.81 \pm 0.17 M_{\odot}$ measured by Teague et al. (2019) from Atacama Large Millimeter/submillimeter Array (ALMA) gas observations is more consistent with the 50% spot model than the 0% spot model. The adopted mass and radius lead to a surface gravity $\log g$ of 4.2, consistent with a near-IR measurement of surface gravity by Sokal et al. (2018), slightly larger than the 4.02 measured by Mentuch et al. (2008) and lower than 4.46 of Venuti et al. (2019), both from optical spectra. Line emission and absorption in He I and He II lines is detected at ~400–450 km s⁻¹ on the red side of the line (see Section 6), which requires a mass larger than $0.8 M_{\odot}$ for freefall from 5 $R_*$ (or ~0.9 $M_{\odot}$ for freefall from 3.5 $R_*$).

X-ray observations reveal the dense shock where the accretion flow heats the star (e.g., Kastner et al. 2002; Stelzer & Schmitt 2004; Brickhouse et al. 2010; Argiroffi et al. 2017). Infall velocities of ~450 km s⁻¹ are seen in He I λ10830 and hot ultraviolet lines (e.g., Herczeg et al. 2002; Dupree et al. 2005; Johns-Krull & Herczeg 2007; Ardila et al. 2013). The X-ray emission is redshifted by only 38 km s⁻¹, consistent with emission in the post-shock region (Argiroffi et al. 2017). Zeeman Doppler Imaging and polarimetry indicate that the accretion morphology is dipole-like, with polar spots where the flow strikes the stellar surface (e.g., Donati et al. 2011; Johns-Krull et al. 2013). Previously measured accretion rates from continuum and line analyses range from $4 \times 10^{-10} - 3 \times 10^{-9} M_{\odot}$ yr⁻¹ (e.g., Muzerolle et al. 2000; Alencar & Batalha 2002; Herczeg et al. 2002; Donati et al. 2011; Robinson & Espaillat 2019), with differences caused by real variability and also by methodological differences.

Monitoring in line and continuum emission reveal changes seen on hours and days timescales, usually attributed to variability in accretion (e.g., Alencar & Batalha 2002; Batalha et al. 2002; Huélamo et al. 2008; Dupree et al. 2012). The stellar rotation period of 3.568 days is measured from sinusoidal changes in the radial velocity and could indicate the presence of a hot Jupiter (Setiawan et al. 2008) but is more likely caused by spot modulation (Huélamo et al. 2008). High time resolution observations with MOST obtained across

several years show frequent bursts on short timescales and no signatures of extinction variability (Siwak et al. 2011, 2018), as can be seen for some accreting systems with disks that are viewed at higher inclinations. Flares are also detected in Si IV and C IV lines (Hinton et al. 2022) and in X-rays (e.g., Brickhouse et al. 2012).

Brγ emission, produced in the accretion flow, extends across at least 3–4 $R_*$, indicating an inner disk truncation radius of at least that size (Garcia Lopez et al. 2020), consistent with the truncation radius inferred from accretion rates and the magnetic field strengths (Johns-Krull 2007; Donati et al. 2011). The disk truncation radius is inside the corotation radius of ~8 $R_*$, which corresponds to the 3.56 day period. In the magnetospheric framework of D’Angelo & Spruit (2010), the accretion should be stable over the long term, while in the simulations of Blinova et al. (2016), the accretion onto TW Hya (and most accreting young stars) should be in the ordered unstable regime, with multiple irregular accretion tongues.

# 3. Observations and Data Reduction

In this section, we provide an overview of the data used in this paper, highlighting aspects of the data and reductions that are most relevant. A brief summary of the spectrographs and representative settings is listed in Table 1.

## 3.1. High-resolution Spectra

The analysis for this survey is built from a foundation of 284 spectra obtained with the Echelle SpectroPolarimetric Device for the Observation of Stars (ESPaDOnS) at the Canada-France-Hawaii Telescope (CFHT). We also downloaded reduced science spectra of TW Hya obtained with the Fiber-fed Extended Range Optical Spectrograph (FEROS; Kaufer et al. 1999) at MPG/ESO 2.2 m telescope, the High Accuracy Radial velocity Planet Searcher (HARPS; Mayor et al. 2003) at the ESO La Silla 3.6 m telescope, in addition to many spectra obtained with other instruments.

The ESPaDOnS spectra were reduced by the automatic data reduction pipeline Libre-ESPRIT (Donati et al. 1997), with some published by Donati et al. (2011). In the ESPaDOnS spectra, an additional post-reduction step eliminated noisy regions near order edges to improve the signal-to-noise ratio (S/N). The individual orders were backed out from the 1D spectra and then recombined using the hrs_merge IDL routine, with final spectra that are oversampled at ~0.01 Å (a factor of a few smaller than the actual pixel scale for ESPaDOnS). Most MIKE spectra were obtained and published by Dupree et al. (2012) and Dupree et al. (2014), with detailed descriptions of the reductions and calibrations. Most HARPS, UVES, and FEROS spectra were automatically reduced and then downloaded from the ESO archive.<fnref n="44" /> Of the FEROS spectra, 54 spectra were reduced and published by Alencar & Batalha (2002) and Batalha et al. (2002). Most Keck/HIRES spectra were obtained from the Keck Observatory Archive, with automated reductions using the MAKEE<fnref n="45" /> pipeline. The Keck spectra obtained on 2008 January 23 and 2008 May 23 were reduced using a custom-written code in IDL. The ESPRESSO data were obtained as part of the PENELLOPE program, with the reduction described in Manara et al. (2021).

<footnote n="44">⁴⁴ https://www.eso.org/sci/facilities/lasilla/instruments/feros/tools/DRS.html</footnote>

<footnote n="45">⁴⁵ see https://sites.astro.caltech.edu/tb/makee/</footnote>
<table><caption>**Table 1**
Summary of Instruments</caption><thead><tr><th>Telescope</th><th>Instrument</th><th>λ Range (Å)</th><th>Resolution</th><th>Aperture</th><th>Sky Subtr.</th><th>#</th><th>Years</th><th>Instr. Reference</th></tr></thead><tbody><tr><td colspan="9" style="text-align:center; font-weight:bold;">High-resolution Spectra</td></tr><tr><td>CFHT</td><td>ESPaDOnS</td><td>3800–10000</td><td>68,000</td><td>1"7 fiber</td><td>No</td><td>284</td><td>2008–2016</td><td>Donati et al. (2006)</td></tr><tr><td>Magellan</td><td>MIKE-red<sup>a</sup></td><td>4850–9100</td><td>35,000</td><td>0"75 slit</td><td>Yes</td><td>467</td><td>2004–2007</td><td>Bernstein et al. (2003)</td></tr><tr><td>ESO-MPG 2.2 m</td><td>FEROS</td><td>3800–9200</td><td>48,000</td><td>2" fiber</td><td>Yes</td><td>241</td><td>1998–2019</td><td>Kaufer et al. (1999)</td></tr><tr><td>ESO 3.6 m</td><td>HARPS</td><td>3800–6900</td><td>115,000</td><td>1" fiber</td><td>No</td><td>34</td><td>2005–2009</td><td>Mayor et al. (2003)</td></tr><tr><td>VLT</td><td>UVES<sup>a</sup></td><td>5000–7000<sup>a</sup></td><td>57,000<sup>a</sup></td><td>0"7<sup>a</sup> slit</td><td>Yes</td><td>47</td><td>2000–2022</td><td>Dekker et al. (2000)</td></tr><tr><td>VLT</td><td>ESPRESSO</td><td>3800–7900</td><td>140,000</td><td>1" fiber</td><td>Yes</td><td>5</td><td>2021</td><td>Pepe et al. (2010)</td></tr><tr><td>SMARTS</td><td>CHIRON</td><td>4100–8900</td><td>30,000</td><td>2"7 fiber</td><td>No</td><td>68</td><td>2021–2022</td><td>Tokovinin et al. (2013)</td></tr><tr><td>Keck I</td><td>HIRES<sup>a</sup></td><td>4800–9200<sup>a</sup></td><td>36,000<sup>a</sup></td><td>1.15<sup>a</sup> slit</td><td>Yes</td><td>18</td><td>2000–2015</td><td>Vogt et al. (1994)</td></tr><tr><td>McDonald 2.7 m</td><td>2D-Coude</td><td>4000–10000</td><td>50,000</td><td>1.15 slit</td><td>Yes</td><td>4</td><td>2006</td><td>Tull et al. (1995)</td></tr><tr><td>CTIO 4 m</td><td>Echelle</td><td>4375–7600</td><td>25,000</td><td>1"5 slit</td><td>Yes</td><td>1</td><td>1998</td><td>...</td></tr><tr><td colspan="9" style="text-align:center; font-weight:bold;">Low- and Medium-resolution Spectra</td></tr><tr><td>VLT</td><td>X-Shooter</td><td>3000–25000</td><td>10,000<sup>a</sup></td><td>1"0 slit<sup>a</sup></td><td>Yes</td><td>6</td><td>2010–2021</td><td>Vernet et al. (2011)</td></tr><tr><td>UH88</td><td>SNIFS</td><td>3100–10000</td><td>1000</td><td>IFU</td><td>Yes</td><td>1</td><td>2014</td><td>Lantz et al. (2004)</td></tr><tr><td>HST</td><td>STIS</td><td>3000–5700</td><td>1000</td><td>0"2 slit</td><td>Yes</td><td>4</td><td>2000–2015</td><td>Woodgate et al. (1998)</td></tr><tr><td>Keck</td><td>LRIS</td><td>3100–9500</td><td>1000</td><td>1"0 slit</td><td>Yes</td><td>1</td><td>2008</td><td>Oke et al. (1995)</td></tr><tr><td>Palomar</td><td>DBSP</td><td>3200–9000</td><td>1000</td><td>2"0 slit</td><td>Yes</td><td>6</td><td>2008</td><td>Oke & Gunn (1982)</td></tr></tbody></table>

**Note.**

We also obtained high-resolution spectra with CHIRON, a bench-mounted, fiber-fed, cross-dispersed echelle spectrograph (Tokovinin et al. 2013) on the 1.5 m telescope at the Cerro Tololo Inter-American Observatory (CTIO) and is part of the Small and Moderate Aperture Research Telescope System (SMARTS). The CHIRON data were taken in *fiber mode*, for a fiber with a diameter of 2"7 on the sky, with 4 × 4 on-chip binning yielding a resolution $\lambda/\delta\lambda \sim 27,800$. Wavelength coverage is complete from 4080–8262 Å in 70 orders, with incomplete coverage to 8900 Å due to interorder gaps between the last five orders. The data were reduced using a pipeline coded in IDL,<sup>46</sup> with cosmic-ray removal using the L. A. Cosmic algorithm (van Dokkum 2001).

No telluric correction is performed on any high-resolution spectrum. All measurements in this paper from the high-resolution spectra are made after normalizing the continuum, including corrections for blaze functions.

This paper uses almost all high-resolution spectra that we could find that cover <6000 Å and had available, high-quality reductions. About 10 high-resolution optical spectra are excluded from this paper because veiling measurements from independent diagnostics have large standard deviations; several of these excluded spectra were obtained with other instruments not included in this paper.

No selection was made for adverse weather conditions. Several spectrographs, including ESPaDOnS (in polarimetry mode) and CHIRON, do not have sky fibers, so sky subtraction is not possible. For bright targets, night-sky emission is generally negligible apart from narrow [O I] and Na D lines and some OH airglow lines at longer wavelengths. However, a bright moon and thin clouds can introduce artificial veiling with a blue color. The CHIRON spectra are excluded from our analysis of [O I] emission because the spectral resolution is low enough that telluric emission blends easily with the narrow line from TW Hya.

Some high-resolution spectra saturate in $H\alpha$, especially affecting UVES and CHIRON spectra. Our $H\alpha$ analysis includes only ESPaDOnS, ESPRESSO, and some HARPS and FEROS observations. Only 10% of HARPS and FEROS spectra are excluded due to saturation, so any introduction of bias is minimal. We confirmed using the 2D images that the $H\alpha$ profiles analyzed in this paper are not affected by saturation or linearity. The only exceptions, where the profile was used but not checked, are old FEROS spectra, which have been previously published and where the line profiles look reasonable by eye. The exclusion of those old FEROS spectra does not affect the results.

## 3.2. Balmer Continuum Spectra

Accretion rates are most accurately measured from flux-calibrated spectra that cover the Balmer jump at ~3646 Å (see review by Hartmann et al. 2016). Our Balmer continuum measurements of TW Hya are listed in Table 1 (see also Table 6). The STIS data were obtained from the archive in a fully reduced and flux-calibrated format. The X-Shooter data were reduced following procedures described in Alcalá et al. (2014). The DBSP, LRIS, and SNIFS data were reduced with custom-build routines in IDL, following Herczeg & Hillenbrand (2014) and Guo et al. (2018). The ground-based spectra have fluxes calibrated to ~10% with a contemporaneous (within ~1 hr) spectrum of the spectrophotometric standard LTT 3864, obtained at similar airmass.

## 3.3. Spectra for Photospheric Templates

Measurements of the excess accretion flux require the subtraction of a photosphere. The primary spectral template used for this paper is TWA 25, a young star in the same association as TW Hya and with a similar spectral type of M0.5 (Herczeg & Hillenbrand 2014). Our default TWA 25 spectrum is obtained by stacking 18 ESPaDOnS spectra from 2016–2017 as part of Programs 16AP18 and 17AF95 (PI Donati). For HARPS

<footnote><sup>a</sup> Typical setting, some spectra obtained with other wavelength settings and slit widths</footnote>

<footnote n="46"><sup>46</sup> http://www.astro.sunysb.edu/fwalter/SMARTS/CHIRON/ch_reduce.pdf</footnote>
spectra of TW Hya, we use as our template HARPS spectra of TWA 25, obtained in 2005–2006 in Programs 076.C-0010 and 074.C-0037. We also use ESPaDOnS spectra of TAP 45 (K6), V819 Tau (K8), and LkCa 7 (M1.2), each selected for spectral type and relatively narrow photospheric lines.

For low-resolution spectra, we measure the accretion luminosity of TW Hya after subtracting a flux-calibrated low-resolution spectrum of TWA 25 that was obtained with Palomar/DBSP (Herczeg & Hillenbrand 2014). A gap from 5500–6300 Å and emission at >8700 Å is filled in with a flux-calibrated Keck/LRIS spectrum of TWA 14.

# 4. Measuring the Veiling of TW Hya

The optical spectrum of TW Hya consists of photospheric emission combined with line and continuum emission produced by the accretion flow and shock. Accretion processes produce strong H and He emission lines, H recombination emission, and an H⁻ continuum, as described and modeled by Calvet & Gullbring (1998). The accretion continuum is blue while the photosphere of a young low-mass star is red, so the ratio of accretion to photospheric flux decreases to longer wavelengths (e.g., Basri & Batalha 1990; Johns & Basri 1995; Dupree et al. 2012). This ratio of accretion flux $F_{\text{acc},\lambda}$ to photospheric flux $F_{\text{phot},\lambda}$ at a wavelength $\lambda$ is defined as the veiling, $r_{\lambda} = F_{\text{acc},\lambda}/F_{\text{phot},\lambda}$. Veiling is measured by comparing the depth of photospheric features in the accreting spectrum to the depth in the spectrum of a photospheric template (e.g., Frasca et al. 2017). In Section 4.3 (see also Figure 6), a properly scaled accretion continuum will be calculated from our measurements of $r_{\lambda}$ and a photospheric template.

Traditionally the veiling is measured by comparing the spectrum to another star, chosen to be a low accretion analog of the target star. Here we use a two-step process, first calculating the relative veiling compared to weak accretion epochs of TW Hya and then calibrating those veiling measurements to a zero-point. We take three different and mostly independent approaches to measure the relative veiling: (1) best fits to spectra over short wavelength segments after adding a flat continuum to weak veiling spectra of TW Hya; (2) equivalent widths of strong lines, and (3) spectral indices obtained by comparing regions with weak fluxes to nearby regions. The analyses are built on the 284 ESPaDOnS spectra, which have high S/N, consistent spectral resolution, and consistent, high-quality reductions, and are then tailored to spectra from other instruments. The final veiling measurements used for subsequent analysis are for the 5000–5100 Å spectral region.

## 4.1. Methodology for Relative Veiling Measurements

In this section, we develop methods to measure veiling, relative to a weak accretion spectrum of TW Hya. We first measure veiling by fitting spectra with a featureless continuum added to the weak accretion spectrum (Section 4.1.1), which works well for fiber-fed spectrographs. We subsequently measure veiling from stacked lines (Section 4.1.2) and spectral indices (Section 4.1.3), using ESPaDOnS to establish relationships to the veiling measurements and feature depths.

The weak accretion spectrum is obtained by summing the 15 ESPaDOnS spectra of TW Hya with the lowest veiling (as measured in Section 4.1.1) and high S/N. We also create a high veiling template from the strongest 11 spectra. Figure 1 shows the comparison between the low and high veiling spectra (for further

details on these spectra, see Appendix A). For other fiber-fed spectra, we also obtain a weak accretion spectrum for the specific instrument and then used stacked lines and spectral indices to scale them to the weak accretion ESPaDOnS spectrum.

### 4.1.1. Veiling Measurements from Spectral Regions

Veiling in the ESPaDOnS, HARPS, FEROS, and CHIRON spectra is measured relative to the weak accretion spectrum of TW Hya, obtained with the same instrument in 25 Å intervals across the full observed wavelength region. This set of values is then converted to a final veiling from 5000–5100 Å, $r_{5050}$, adopted for the remainder of the analysis. This wavelength range is selected because it is blue enough for a wide range of veilings, is covered by and well separated from dichroics in low-resolution spectra, and is red enough to be covered by most high-resolution spectra.

To measure the veiling in each 25 Å interval, the weak accretion template is normalized by the median flux in that segment. A flat featureless emission spectrum is then added to the weak accretion spectrum, and the combined spectrum is renormalized by dividing by $(1+r_{\lambda})$. The best-fit veiling is obtained by minimizing $\chi^2$ in the difference between the science and veiled template spectrum (see also Figure 1).

Veiling measurements from all spectral regions are tightly correlated with the veiling from 5000–5100 Å. Six examples are provided in Figure 2. We obtain final veiling measurements for the 5000–5100 Å region, $r_{5050}$, by combining the highest quality intervals across the full spectrum, as follows. First an average veiling, $\bar{r}_{5050}$, is obtained by averaging the veiling from the four intervals between 5000 and 5100 Å. Relationships between $\bar{r}_{5050}$ and $r_{\lambda}$ (the set of veiling in 25 Å intervals, as in Figure 2) are then fit with a second-order polynomial (not shown in the figure). These relationships are then used to convert the set of $r_{\lambda}$ to a set of veilings $r'_{5050}$, now all on the same scale (veiling between 5000 and 5100 Å).

For each spectrum, the final veiling measurement, $r_{5050}$, is finally obtained by taking the median $r'_{5050}$ from 80 distinct wavelength intervals. These wavelength intervals are selected because the scatter between the measured and the preliminary $r_{5050}$ is less than 0.2. These regions are all located between 4350 and 6650 Å. Regions at longer wavelengths have large scatter, either because of temporal changes in visual spot coverage or temperature or in the shape of the accretion continuum. Comparisons between independent combinations of $r'_{\lambda}$ are consistent with a standard deviation of 0.003, which is adopted as the precision in veiling estimates from ESPaDOnS spectra.

For the individual veiling measurements, many spectral regions are identified as unreliable and are not used. Telluric absorption lines, strong emission lines, and deep, gravity-sensitive lines are straightforward to identify and avoid. The region from 6525–6605 Å is excluded from the fit because that entire region is contaminated by Hα emission; a similar region is avoided around Hβ and Hγ. The region at <4150 Å is ignored because of low S/N and veilings that are often too high to be accurately measured with our automated method. Regions with weak emission lines<fnref n="47" /> are identified and avoided by visually

<footnote n="47">⁴⁷ These emission lines fill in absorption lines and are only detectable after subtracting off a template. These lines are more easily detected in heavily veiled spectra (e.g., Gahm et al. 2008) but the process has been demonstrated with models to also be important for other accretors (Dodin & Lamzin 2012). These lines are not identified here but are likely Fe I and Fe II lines (see, e.g., Hamann & Persson 1992; Beristain et al. 1998; Stempels & Piskunov 2003).</footnote>
Normalized Flux

<figure><img src="image_4.png" /></figure>

**Figure 1.** Coadded ESPaDOnS spectra of TW Hya during visits with strong accretion ($F_{\text{high veil, blue}}$) and weak accretion ($F_{\text{low veil, red}}$) over five wavelength ranges. The residual spectrum ($F_{\text{residual}}$, black, scaled by a factor of 3 for visualization) is calculated by veiling the weak accretion spectrum and then subtracting it from the strong accretion spectrum, as $F_{\text{residual}} = F_{\text{high veil}} - (F_{\text{low veil}} + r_{\lambda})/(1+r_{\lambda})$, with $r_{\lambda}$ as the veiling at wavelength $\lambda$ and with all spectra normalized to the continuum level within the spectral region (normalization shifted slightly in this Figure for visualization). The residual spectrum reveals emission lines that would be undetectable in any single spectrum, affecting the spectrum by filling in photospheric absorption lines (see examples in Gahm et al. 2008; Dodin & Lamzin 2012). Many of the narrow emission lines in the residual spectrum are Fe I (e.g., Hamann & Persson 1992; Beristain et al. 1998), although line identification is beyond the scope of this paper. Appendix A describes how the strong, weak, and residual spectra are calculated.

(The data used to create this figure are available.)

comparing the strong accretion and weak accretion templates.
The masked regions and the 2% of pixels with the largest
differences between the template and science spectrum are

excluded from the $\chi^2$ calculation. All measurements are made
after applying small wavelength adjustments ($\lesssim 0.5$ km s$^{-1}$ for
ESPaDOnS) from line centroids, as measured in Section 4.1.2.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Left: the relative veiling (veiling relative to the low-veiling spectrum of TW Hya) in six wavelength regions (y-axis) compared to the veiling from 5000–5100 Å (x-axis), as measured from ESPaDOnS spectra. Each veiling measurement plotted here is the average of four to eight different 25 Å spectral regions around the labeled wavelength. The line shows a 1:1 relationship and not a fit. The veiling at bluer wavelengths is generally higher than the veiling at red wavelengths, with exceptions. The veiling at 4700 Å is similar to that at 5050 Å because of similar photospheric fluxes. Right: correlated differences between measured and expected veiling at long wavelengths. The inset shows the correlation between veiling at 7000 Å and 5000–5100 Å. The excess veiling is calculated by fitting a line to the $r_{\lambda-r_{5050}}$ relationship (as can be seen in the left panel) and then subtracting each point from that line. The main plot shows this excess veiling at 7000 Å vs. that at 8700 Å. When the veiling at 7000 Å is higher than expected, the veiling at 8700 Å is also higher than expected. The correlation between excess veiling at 7000 and 8700 Å demonstrates that the scatter in the correlations between veiling measurements is real and not due to S/N or other statistical uncertainties.</figcaption></figure>

### 4.1.2. Line Equivalent Widths as Veiling Measurements

As veiling increases, the equivalent width of photospheric absorption lines decreases. In this subsection, we measure veiling from the depth of coadded photospheric absorption lines,<fnref n="48" /> following the process in the left panels of Figure 3.

We first identify ~260 lines that are between 4500 and 6500 Å and are isolated enough to have single-peak absorption profiles. We then measure equivalent widths for each of the 260 lines from fits with Gaussian profiles. From this analysis, we select 82 lines (see Appendix B) that have equivalent widths in ESPaDOnS spectra greater than 0.08 Å in a median accretion spectrum and that are well correlated with veiling in ESPaDOnS spectra. These selected lines are then separated by wavelength into 5–6 distinct groups, each with eight to 15 lines that are normalized and coadded in velocity space; the exact lines and groupings depend on the wavelength coverage of the instrument. The equivalent widths for each set of coadded lines are measured by fitting Gaussian profiles to the absorption spectrum.

The set of equivalent widths from coadded lines is converted to veiling at 5000–5100 Å by fitting a third-degree polynomial to the relationship between veiling and the inverse of equivalent width, as measured from ESPaDOnS spectra (Figure 3). The final veiling from equivalent widths, $r_{5050}^{EW}$, is then measured from the average of the set of 5–6 coadded equivalent widths (see Appendix B). The number of line sets and the number and categorization of lines into those sets are tailored to the spectral coverage and spectral resolution of the science question. For CHIRON spectra, some lines that are used with ESPaDOnS are excluded because the lower spectral resolution makes them harder to measure. The polynomial fits are recalculated for a set of lines that are easily measured with CHIRON and after convolving the ESPaDOnS spectra to the resolution of CHIRON.

When veiling is measured against a template, the veiling depends on the line, as demonstrated in Figure 1. Rei et al. (2018) found that higher veiling values are measured from stronger photospheric lines. Since our equivalent widths are converted to relative veiling through correlations rather than directly comparing line depths to non-accreting templates, they are robust to effects introduced by line-dependent veiling.

### 4.1.3. Spectral Indices as Veiling Measurements

Across the spectrum, some regions of strong line blends have fluxes that are much lower than nearby regions. The spectral index for these regions (the weak region divided by the strong region, or vice versa) depends on veiling by the accretion continuum. In this subsection, we follow similar procedures as in Section 4.1.2 to measure the flux ratios and then convert the spectral index to veiling, following the right panels of Figure 3.

We measure spectral indices in seven relatively narrow regions with faint emission divided by nearby spectral regions with brighter emission (see Appendix B). In the ESPaDOnS data, seven spectral indices are tightly correlated with veiling. Six additional regions were excluded because the correlations with veiling had large scatter. For several instruments, we discard some of the five indices because either the spectral region or background continuum region is located near the edge of an echelle order. We also avoid TiO bands because they are sensitive to spots and are located at long wavelengths, where weak veiling limits the lever arm and leads to large fractional uncertainties.

The relationship between each spectral index and veiling, $r_{5050}$ for our ESPaDOnS sample, is fit with a fourth-degree polynomial. Each spectral index is then converted into an estimate for veiling, with the median value adopted as the veiling measurement, $r_{5050}^{rat}$. The relative veiling measurements are only applied for the range in ESPaDOnS veiling. The flux ratios and polynomial fits are recalculated for each instrument by degrading either the resolution of either the ESPaDOnS spectra or the science spectrum so that they match.

<footnote n="48">⁴⁸ The procedure and results should be the same by treating each line individually, but results were more robust when first coadding sets of lines.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Top left: the relationship between the equivalent width of coadded lines (one set shown here for lines from 5250–5400 Å) vs. veiling, established from ESPaDOnS spectra and with residuals in ESPaDOnS and FEROS spectra shown below. Top right: similar plots as on the left for the flux ratio for a spectral dip around 5205–5215 Å compared with a continuum region. In both cases, we calculate a best-fit relationship between equivalent width (or flux ratios) and veiling with ESPaDOnS. We then apply those relationships to FEROS spectra. The bottom panels on the left and right show the scatter between the veiling calculated from these relationships and the veiling measured by comparing the FEROS spectrum to a low-veiling FEROS template. Bottom: the final comparison of veiling obtained from the combination of line equivalent widths and flux ratios to the veiling from a low-veiling FEROS template. All FEROS spectra are shifted by ~0.02 to place them on the same scale as the ESPaDOnS veiling measurements. The scatter of 0.013 + 0.045 × r5050 (with a minimum error of 0.02) is applied as an uncertainty to all spectra where veiling is measured from line equivalent widths and flux ratios.</figcaption></figure>

## 4.2. Applying Relative Veiling Methods to Spectra

Figure 3 and Table 2 summarize how the methods and the associated errors are applied to spectra from each instrument. The veiling measurements from ESPaDOnS spectra form the

backbone of our analysis. The FEROS, HARPS, and CHIRON spectra all have veiling initially measured relative to low-veiling spectra from the same instrument. These veilings are then offset to the ESPaDOnS values (calculated only for
<table><caption>**Table 2**
Methods and Errors</caption><thead><tr><th rowspan="2">Instrument</th><th rowspan="2">Direct<br>Scatter</th><th colspan="2">Equivalent Width</th><th colspan="2">Flux Ratios</th><th colspan="3">Adopted</th></tr><tr><th>Offset<sup>a</sup></th><th>Scatter</th><th>Offset<sup>a</sup></th><th>Scatter</th><th>Relative Error</th><th>Random Error</th><th>Method</th></tr></thead><tbody><tr><td>ESPaDOnS</td><td>0.003</td><td>-0.002</td><td>0.021</td><td>0.000</td><td>0.024</td><td>0.02<sup>b</sup></td><td>0.003</td><td>Direct</td></tr><tr><td>FEROS</td><td>0.006</td><td>0.002</td><td>0.043</td><td>0.040</td><td>0.031</td><td>0.02</td><td>0.006</td><td>Direct, shifted</td></tr><tr><td>HARPS</td><td>0.003</td><td>0.068</td><td>0.026</td><td>0.099</td><td>0.019</td><td>0.02</td><td>0.003</td><td>Direct, shifted</td></tr><tr><td>CHIRON</td><td>0.011</td><td>0.128</td><td>0.065</td><td>0.076</td><td>0.073</td><td>0.04</td><td>0.011</td><td>Direct, shifted</td></tr><tr><td>ESPRESSO</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.02</td><td>EW+ratios</td></tr><tr><td>MIKE</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.03</td><td>EW</td></tr><tr><td>UVES</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.02</td><td>EW+ratios</td></tr><tr><td>HIRES</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.02</td><td>EW+ratios</td></tr><tr><td>2coude</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.03</td><td>EW</td></tr></tbody></table>

**Notes.**

spectra with $r_{5050} < 0.5$) by using the correlations of veiling with equivalent widths and flux indices. For the slit-based spectra and for ESPRESSO (with only five spectra), the veiling measurements are obtained by averaging the veiling estimated from the equivalent widths and spectral indices. For most MIKE spectra and one FEROS spectrum, the veiling is measured from the equivalent width of coadded photospheric lines, which is the method most robust to uncertainties in the relative wavelength solution.

The analysis of uncertainties is split into the precision of the veiling measurements relative to each other (Section 4.2.1) and calibrating all of the measurements against a pure photosphere (Section 4.2.2).

### 4.2.1. Relative Errors between Measurements

For ESPaDOnS, the veiling measurements are adopted directly from the set of ~80 veiling values. Random combinations of different sets of veilings indicate a precision of 0.003. Similarly, the precision of FEROS, HARPS, and CHIRON veiling measurements ranges from 0.003–0.006. Those sets of veiling values are scaled to ESPaDOnS, with a relative error of ~0.02, by using results from the equivalent width and spectral index analysis for spectra with $r_{5050} < 0.4$.

For the slit spectra and ESPRESSO, the combined use of equivalent widths and spectral indices leads to uncertainties of ~0.02 for $r_{5050} < 0.3$ and ~0.05 for $r_{5050} > 1$. The scatter versus veiling is well approximated by $0.013 + 0.045 \times r_{5050}$ and is adopted (with a minimum error of 0.02) as the assessed uncertainty. These values are calculated from the dispersion in the individual measurements and from the scatter in the comparison between the FEROS values for the direct veiling measurements and the veilings measured from line equivalent widths and spectral indices.

On eight occasions, two different instruments observed TW Hya within the same 3 hr window (Table 3). The most significant deviation occurred for two observations separated by 62 minutes on MJD 59310, when the veiling was high. This large discrepancy may be caused by a rapid decline of the veiling as an accretion burst faded or by artificial differences due to larger errors when the veiling is high. For other epochs, the average difference in veiling is 0.06.

<table><caption>**Table 3**
Cross-instrument Comparisons of Relative Veilings</caption><thead><tr><th>MJD</th><th>Inst. 1</th><th>r</th><th>Δt (minute)</th><th>Inst. 2</th><th>r</th></tr></thead><tbody><tr><td>54158.382</td><td>FEROS</td><td>0.98</td><td>&lt;20<sup>a</sup></td><td>MIKE</td><td>0.95</td></tr><tr><td>54159.114</td><td>FEROS</td><td>0.46</td><td>&lt;20<sup>b</sup></td><td>MIKE</td><td>0.47</td></tr><tr><td>59280.173</td><td>CHIRON</td><td>0.59</td><td>163</td><td>ESPRESSO</td><td>0.48</td></tr><tr><td>59308.039</td><td>ESPRESSO</td><td>0.61</td><td>131</td><td>CHIRON</td><td>0.61</td></tr><tr><td>59309.135</td><td>CHIRON</td><td>1.24</td><td>11</td><td>ESPRESSO</td><td>1.39</td></tr><tr><td>59310.086</td><td>ESPRESSO</td><td>1.47</td><td>62</td><td>CHIRON</td><td>1.23</td></tr><tr><td>59313.146</td><td>CHIRON</td><td>0.08</td><td>102</td><td>ESPRESSO</td><td>0.06</td></tr><tr><td>59667.139</td><td>CHIRON</td><td>0.71</td><td>127</td><td>UVES</td><td>0.84</td></tr></tbody></table>

**Notes.**

This analysis, though limited, supports our error estimates. A few spectra are likely affected by clouds and the moon and may have uncertainties that are underestimated.

### 4.2.2. Zero-point Measurements for Veiling

The primary shortcoming of using TW Hya as its own template is in calibrating the veiling measurements to a zero-point, $r_{zpt}$. Our analysis described above provides veiling measurements relative to a low accretion spectrum of TW Hya. To convert these relative veiling measurements into an absolute veiling relative to a non-accreting star, the veiling of the low accretion spectrum needs to be calculated through comparisons with non-accreting spectral templates.

We measure a wavelength-dependent zero-point by fitting stacked TW Hya spectra with the spectra of non-accreting stellar templates (Table 4) with spectral types from K6–M1.5 to bracket that of TW Hya. These non-accreting stars are all young,<fnref n="49" /> selected to match the gravity and chromospheric effects of TW Hya (see, e.g., discussions of templates in Ingleby et al. 2013 and Manara et al. 2013).

The spectral fits are similar to those described in Section 4.1.1, with additional free parameters $v \sin i$ and radial velocity alongside the veiling. The rotational broadening kernel

<footnote><sup>a</sup> Difference between veiling from line ratios/equivalent width and veiling measured from a low-resolution template</footnote>

<footnote><sup>b</sup> Overall offset for all data from zero-point analysis</footnote>

<footnote><sup>a</sup> Median of six spectra obtained within 20 minutes of FEROS.</footnote>

<footnote><sup>b</sup> Median of seven spectra obtained within 20 minutes of FEROS.</footnote>

<footnote n="49">⁴⁹ Previous analyses had adopted the K6 dwarf star GJ 1172 as a spectral template for TW Hya (e.g., Alencar & Batalha 2002; Dupree et al. 2012); however, GJ 1172 in particular has much shallower TiO absorption than TW Hya and does not provide a good fit over the full optical wavelength range.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Examples of the correlation between absolute veiling measured from WTTS templates and relative veiling measured from weak-accreting spectra of TW Hya itself. Colored dots indicate the measured $r_{abs}$ and $r_{rel}$ in each sextile with different colors representing different template stars. The solid line and shaded region in the same color as dots denote the best-fit linear relation between $r_{abs}$ and $r_{rel}$ and its $1\sigma$ error for each template. The yellow stars at $r_{rel} = 0$ is the best-fit $r_{zpt}$ determined as the inverse-variance weighted mean of the y-intercepts. The $r_{zpt}$ measured from different templates are consistent with each other.</figcaption></figure>

<table><caption>**Table 4**
Spectral Templates</caption><thead><tr><th>Star</th><th>SpT<sup>a</sup></th><th>v sin i<sup>b</sup></th><th>Reference<sup>c</sup></th></tr></thead><tbody><tr><td>TW Hya</td><td>M0.5</td><td>4</td><td>Donati et al. (2011)</td></tr><tr><td>TAP 45 (V1076 Tau)</td><td>K6</td><td>7.7</td><td>Nguyen et al. (2012)</td></tr><tr><td>V819 Tau</td><td>K8</td><td>9.5</td><td>Donati et al. (2015)</td></tr><tr><td>TWA 25 (V1249 Cen)</td><td>M0.5</td><td>11.9</td><td>Nicholson et al. (2021)</td></tr><tr><td>LkCa 7 (V1070 Tau)</td><td>M1.2</td><td>14.7</td><td>Nguyen et al. (2012)</td></tr></tbody></table>

**Notes.**

<sup>a</sup> From Herczeg & Hillenbrand (2014).

<sup>b</sup> km/s

<sup>c</sup> Reference for $v \sin i$ measurements.

is adopted from Gray (2005) and does not include limb darkening. The TW Hya spectra used in this fit are mean-stacked (ESPaDoNS) spectra in six equal-sized bins of relative veilings. The absolute veiling ($r_{abs}$) measured from this fit is correlated with the mean relative veiling ($r_{rel}$) of each sextile measured from weak-accreting TW Hya spectra (Figure 4).

Figure 5 shows the relation between $r_{zpt}$ and wavelength, which closely resembles the shape of relative veiling as shown in Figure 6. The $r_{zpt}$ is not correlated with the spectral type of the template. The systematic error of $r_{zpt}$, estimated as the median rms of the $r_{zpt}$ determined from different templates, is $\sim 0.02$.

### 4.3. Summary of Veiling Measurements and the Accretion Spectrum

The combined approaches from the analysis above yield veiling measurements for 1169 spectra. The final veiling values (presented in Table 5) for the 5000–5100 Å region range from 0.13–2.13 with a median of 0.67, after excluding observations obtained within 3 hr of another observation.

The uncertainties in the veiling are small enough to be negligible relative to other errors. The relative error in most accretion rate estimates is dominated by uncertainties in how the bolometric correction varies with time and by the change in the photospheric spectrum due to cool photospheric spots. In principle, an increase (or decrease) in accretion spots (or a change in accretion spot coverage on the visible surface) could also decrease (or increase) the observed photospheric emission,

<figure><img src="image_8.png" /><figcaption>**Figure 5.** The zero-point of relative veiling as a function of wavelength. Yellow stars represent the final $r_{zpt}$ determined by averaging the individually measured $r_{zpt}$ from different stellar templates in colored dots in the background. The shape resembles the relation between relative veiling and wavelength.</figcaption></figure>

although usually the accretion spots only cover a few percent of the stellar surface (e.g., Calvet & Gullbring 1998). The zero-point error in the veiling of 0.02 leads to a fractional uncertainty of 10%–15% during periods of the weakest accretion.

The ESPaDoNS spectra yield veiling measurements across the 4150–9000 Å wavelength region, so a veiling at any wavelength is accurately converted to a veiling at any other wavelength, with a relative uncertainty of $\sim 0.01–0.03$. Figure 6 shows that the veiling increases to short wavelengths, as expected for hot emission against a cooler photosphere. Spectral features in the photosphere, including at 5200 Å and the TiO 7140 band, are seen as sharp changes in the veiling because the nearly flat accretion continuum is divided by a photosphere with features.

The accretion spectrum is obtained by convolving veiling with the flux-calibrated template spectrum of TWA 25. For low and modest veiling values ($r < 0.5$, corresponding to 0.025 $L_\odot$, see Section 5.1), the accretion spectrum is consistent with a constant flux across the optical spectral range. When the veiling is high ($r > 1.2$, or $L_\odot \sim 0.06 L_\odot$), the accretion spectrum becomes stronger at blue wavelengths.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Veiling spectra for strong (pink circles, average relative veiling $r_{5050} = 1.52$), above average (blue diamonds, $r_{5050} = 0.95$), below average (red squares, $r_{5050} = 0.52$), and weak (black triangles, $r_{5050} = 0.21$) veiling epochs measured in ESPaDOnS spectra. The left panel shows the relative veiling measurement, with features in the veiling spectrum caused by the shape of the photospheric template, TWA 25. The right panel shows the veiling spectrum multiplied by TWA 25. The spectra with moderate accretion are flat from 4000–9000 Å, while the strongest accretion spectrum has a bluer slope (increasing flux with decreasing wavelength).</figcaption></figure>

<table><caption>**Table 5**
Veiling and Accretion Measurements</caption><thead><tr><th>MJD</th><th>r<sub>5050</sub></th><th>σ(r<sub>5050</sub>)</th><th>L<sub>☉</sub></th><th>M<sub>☉</sub></th><th>Inst.</th></tr></thead><tbody><tr><td>53842.20872</td><td>1.07</td><td>0.056</td><td>0.057</td><td>3.15e-09</td><td>MIKE</td></tr><tr><td>53842.96101</td><td>0.59</td><td>0.037</td><td>0.032</td><td>1.81e-09</td><td>MIKE</td></tr><tr><td>54225.09375</td><td>0.28</td><td>0.003</td><td>0.017</td><td>9.38e-10</td><td>FEROS</td></tr><tr><td>59241.30469</td><td>0.24</td><td>0.012</td><td>0.015</td><td>8.43e-10</td><td>ESPaDOnS</td></tr><tr><td>59309.13672</td><td>1.41</td><td>0.018</td><td>0.074</td><td>4.10e-09</td><td>ESPaDOnS</td></tr></tbody></table>

Note. Five dates selected randomly. The full table is available online. (This table is available in its entirety in machine-readable form.)

Figure 2 also shows that the excess veiling at red wavelengths is correlated with excess veiling at other red wavelengths. This excess veiling is the difference in veiling between the measured veiling at wavelength $\lambda$ and the veiling expected from the final veiling $r_{5050}$ and $\lambda$. If the scatter in those correlations were a consequence of measurement uncertainties, then the excess veiling at two wavelengths would be uncorrelated. The correlation between the excess veiling at 7000 and 8700 Å demonstrates that the scatter at red wavelengths is real. This scatter may be caused by either spots or by changes in the temperature of the accreting gas.

## 5. Converting Veiling Measurements to Accretion Rate

The primary goal of this paper is to analyze the stability of the accretion rate with time. In this section, we convert the veiling measurements in Section 4 to accretion rates by first measuring accretion luminosities from flux-calibrated spectra of TW Hya, and then find a relationship between the accretion luminosity and the veiling at 5000–5100 Å (Section 5.1). Finally, we use that relationship to calculate accretion rates from the sample of high-resolution spectra (Section 5.2).

The accretion rates are measured from broadband, flux-calibrated spectra. Assuming that all gravitational energy is converted into luminosity and following (Gullbring et al. 1998), the accretion rate is calculated by

$$ \dot{M}_{\text{acc}} = \left(1 - \frac{R_{*}}{R_{\text{in}}}\right)^{-1} \frac{L_{\text{acc}} R_{*}}{GM_{*}} \sim 1.25 \frac{L_{\text{acc}} R_{*}}{GM_{*}}, \quad (1) $$

where $R_*$ and $M_*$ are the stellar radius and mass. For convention and consistency with other estimates, we adopt $R_{\text{in}} \sim 5 R_*$ as the disk truncation radius (e.g., Johns-Krull 2007; Johnstone et al. 2014). The truncation radius of $\sim 3.5 R_*$ measured by Garcia Lopez et al. (2020) would increase the accretion rate by 12%.

The flux-calibrated spectra are fit with the combination of an accretion spectrum and a low-resolution spectrum of the photospheric template TWA 25. The spectra of TW Hya show excess emission at all optical wavelengths, with an increase shortward of the Balmer jump at 3646 Å. The Balmer jumps measured in this paper are the ratio of the excess emission at 3600 Å to the excess emission at 4200 Å, after subtracting the photospheric template. The accretion spectrum across the Balmer jump is modeled with a plane-parallel slab (Valenti et al. 1993, as implemented by Herczeg & Hillenbrand 2014), while the Paschen and $H^-$ continua are assumed to be featureless and constant in flux, when possible. The best-fit combination of a photosphere and accretion continuum is determined by eye to minimize residuals, with a veiling that is generally accurate to $\sim 0.05$. The spectra obtained at medium or high resolution are convolved to low resolution to match the spectrum of TWA 25. The accretion luminosity for each epoch is calculated by summing the total luminosity in the best-fit model.

Figure 7 shows examples of these fits. The residual after subtracting the photosphere should have no significant absorption feature. In periods of strong accretion, the residual emission (the accretion continuum, after subtracting the photosphere) declines in flux to longer wavelengths, while epochs with weak or moderate accretion have an accretion continuum consistent with a constant flux. The accretion continuum on 2019 July 7 is described by $F_{\text{acc}} = F_{\text{acc},5000} \times (1 - 2.5 \times 10^{-4}(\lambda - 5000))$, consistent with the slope in Figure 6. The example epoch with weak accretion, on 2021 April 8, has a slope that is consistent with a constant flux.

The accretion luminosities depend on the assumptions in the broadband accretion spectrum, especially where the accretion emission is unobserved. Table 7 shows the comparison between our accretion rates and literature accretion rates, as measured from the same spectrum. Our accretion luminosities for the X-Shooter spectra are almost exactly the same
<figure><img src="image_10.png" /><figcaption>**Figure 7.** The accretion luminosity measured from flux-calibrated spectra of TW Hya. The left panel shows the spectrum of TW Hya (black), modeled (light blue) as the sum of a stellar photosphere (red) and accretion continuum (blue), with a residual accretion spectrum (pink) calculated by subtracting the scaled photosphere from the TW Hya spectrum. The right panel shows the TW Hya spectrum and the residual accretion spectrum on an epoch of strong accretion (2019 July 7) and weak accretion (2021 April 8). The weak accretion spectra tend to have larger Balmer jumps, while the strong accretion spectra have a bluer slope (increasing flux to shorter wavelengths).</figcaption></figure>

(difference of ~3%) as the values from Manara et al. (2014) and Venuti et al. (2019), both of which applied an independent model for hydrogen slab emission.<fnref n="50" /> However, recent accretion models include multiple accretion shocks with a range in temperatures to better explain the veiling at red wavelengths (e.g., Ingleby et al. 2013; Robinson & Espaillat 2019; Espaillat et al. 2022; Pittman et al. 2022). The introduction of these components significantly enhances the flux at red wavelengths, with accretion luminosities that are 30%–60% higher in Robinson & Espaillat (2019) than the accretion luminosities measured here. The veiling of TW Hya in the near-IR is low and much more stable (Sousa et al. 2023), because either the cooler accretion component or the disk emission varies less than the hotter emission.

## 5.1. Bolometric Corrections for High-resolution Spectra

The fits to the flux-calibrated spectrum of TW Hya combine an accretion spectrum and a photospheric spectrum. In this section, we use those fits to derive a relationship to convert veiling at 5000–5100 Å to accretion luminosity, with results presented in Table 6.

The accretion luminosity is calculated by measuring the accretion spectrum and veiling in 26 flux-calibrated blue spectra of TW Hya (see examples in Figure 7). The veiling measured at high resolution is 20% lower than the veiling measured from low-resolution spectra, in a limited comparison of one Keck/HIRES and two VLT/UVES spectra. The difference in veiling is likely attributed to the definition and normalization of the photospheric continuum level, which should be lower when many lines blend together.

Figure 8 shows our conversions from flux to luminosity. The average bolometric correction is $F_{5050}/F_{acc} = (1.33 \pm 0.12) \times 10^{-4}$, where this uncertainty is calculated from the standard deviation between the data points and the best-fit line. This bolometric correction is applied to the veiling for TW Hya, as

$$L_{acc} = 0.050 \times r_{lowres} = 0.062 \times r_{highres} \quad (2)$$

<table><caption>**Table 6**
Accretion Rates and Veiling in Low-resolution Spectra</caption><thead><tr><th>Instr.</th><th>Date</th><th>r<sub>low</sub><sup>a</sup></th><th>F<sub>acc</sub><sup>b</sup></th><th>BJ</th><th>L<sub>acc</sub>/L<sub>☉</sub><sup>c</sup></th><th>Reference</th></tr></thead><tbody><tr><td>STIS</td><td>2000-05-07</td><td>1.08</td><td>5.99e-14</td><td>2.03</td><td>0.045</td><td>H04</td></tr><tr><td>STIS</td><td>2002-07-19</td><td>0.64</td><td>4.23e-14</td><td>2.82</td><td>0.038</td><td>HH08</td></tr><tr><td>DBSP</td><td>2008-01-18</td><td>0.52</td><td>3.27e-14</td><td>2.17</td><td>0.026</td><td>HH14</td></tr><tr><td>DBSP</td><td>2008-01-19</td><td>0.61</td><td>4.12e-14</td><td>2.12</td><td>0.033</td><td>HH14</td></tr><tr><td>DBSP</td><td>2008-01-20</td><td>0.69</td><td>5.07e-14</td><td>1.82</td><td>0.038</td><td>HH14</td></tr><tr><td>LRIS</td><td>2008-05-28</td><td>0.39</td><td>2.43e-14</td><td>3.36</td><td>0.024</td><td>H09</td></tr><tr><td>DBSP</td><td>2008-12-28</td><td>0.45</td><td>2.92e-14</td><td>2.18</td><td>0.023</td><td>HH14</td></tr><tr><td>DBSP</td><td>2008-12-29</td><td>0.56</td><td>3.48e-14</td><td>2.47</td><td>0.031</td><td>HH14</td></tr><tr><td>DBSP</td><td>2008-12-30</td><td>0.72</td><td>4.20e-14</td><td>1.70</td><td>0.031</td><td>H14</td></tr><tr><td>STIS</td><td>2010-01-28</td><td>1.44</td><td>7.41e-14</td><td>1.58</td><td>0.051</td><td>RE19</td></tr><tr><td>STIS</td><td>2010-02-04</td><td>0.37</td><td>2.31e-14</td><td>2.33</td><td>0.019</td><td>RE19</td></tr><tr><td>XSH</td><td>2010-04-07</td><td>0.75</td><td>3.51e-14</td><td>1.70</td><td>0.026</td><td>V19</td></tr><tr><td>XSH</td><td>2010-05-03</td><td>0.59</td><td>3.85e-14</td><td>2.07</td><td>0.029</td><td>M14</td></tr><tr><td>STIS</td><td>2010-05-28</td><td>0.85</td><td>4.68e-14</td><td>1.48</td><td>0.031</td><td>RE19</td></tr><tr><td>STIS</td><td>2015-04-18</td><td>0.89</td><td>5.31e-14</td><td>1.66</td><td>0.041</td><td>RE19</td></tr><tr><td>HIRES</td><td>2008-05-23</td><td>0.43</td><td>3.01e-14</td><td>2.37</td><td>0.024</td><td>F18</td></tr><tr><td>SNIFS</td><td>2014-11-27</td><td>1.22</td><td>6.10e-14</td><td>1.33</td><td>0.051</td><td>G18</td></tr><tr><td>SNIFS</td><td>2014-11-29</td><td>1.94</td><td>1.21e-13</td><td>1.52</td><td>0.081</td><td>G18</td></tr><tr><td>SNIFS</td><td>2014-12-08</td><td>0.72</td><td>4.47e-14</td><td>2.47</td><td>0.042</td><td>G18</td></tr><tr><td>SNIFS</td><td>2014-12-10</td><td>1.56</td><td>6.43e-14</td><td>1.37</td><td>0.050</td><td>G18</td></tr><tr><td>SNIFS</td><td>2014-12-13</td><td>1.22</td><td>6.79e-14</td><td>1.72</td><td>0.049</td><td>G18</td></tr><tr><td>XSH</td><td>2019-07-06</td><td>1.70</td><td>8.46e-14</td><td>1.43</td><td>0.061</td><td>...</td></tr><tr><td>XSH</td><td>2019-07-07</td><td>1.86</td><td>1.21e-13</td><td>1.47</td><td>0.105</td><td>...</td></tr><tr><td>XSH</td><td>2021-04-02</td><td>1.27</td><td>7.71e-14</td><td>1.60</td><td>0.055</td><td>M21</td></tr><tr><td>XSH</td><td>2021-04-06</td><td>1.70</td><td>1.15e-13</td><td>1.45</td><td>0.092</td><td>M21</td></tr><tr><td>XSH</td><td>2021-04-08</td><td>0.22</td><td>2.47e-14</td><td>2.60</td><td>0.020</td><td>M21</td></tr></tbody></table>

**Notes.**

References: H04: Herczeg et al. (2004) HH08: Herczeg & Hillenbrand (2008) H09: Herczeg et al. (2009) HH14: Herczeg & Hillenbrand (2014) and 2023, in preparation RE19: Robinson & Espaillat (2019) F18: Fang et al. (2018) V19: Venuti et al. (2019) M14: Manara et al. (2014) G18: Guo et al. (2018) M21: Manara et al. (2021)

with an uncertainty of ~10%, which is consistent with typical uncertainties in flux calibration. The low-resolution calculation is obtained directly from the average value on the right panel in Figure 8, while the conversion for high resolution (the right

<footnote n="50">⁵⁰ Our bolometric correction for 5000–5100 Å also matches the bolometric correction of a 9000 K blackbody.</footnote>

<footnote>a Veiling at 5050 Å</footnote>

<footnote>b Accretion continuum flux at 5050 Å.</footnote>

<footnote>c Accretion luminosity in L<sub>☉</sub>.</footnote>
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Left: the size of the Balmer jump (here defined as the flux ratio of 3600–4200 Å) vs. the veiling estimated from low-resolution spectra. The Balmer jump tends to be higher when the veiling is lower. When the Balmer jump is larger, more emission escapes at short wavelengths, so the bolometric correction is larger. Right: the bolometric correction, used to convert the measured accretion flux at 5050 Å to the total accretion luminosity, vs. the veiling, as measured from low-resolution spectra.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Histograms of accretion rates, plotted as a distribution in linear (left) and logarithmic (right) values, with corresponding veiling on the top axis. The average log $M_{\text{acc}}$ is $10^{-8.65} M_{\odot} \text{yr}^{-1}$ with a standard deviation of 0.22 dex. The differences in veiling are much larger than the precision and cross-instrument uncertainties of ~0.02.</figcaption></figure>

<table><caption>**Table 7**
Comparison of Accretion Luminosities</caption><thead><tr><th>Instr.</th><th>Date</th><th>This Work<br>L<sub>acc</sub></th><th>Lit.<br>L<sub>acc</sub></th><th>Reference</th></tr></thead><tbody><tr><td>STIS</td><td>2010-01-28</td><td>0.051</td><td>0.071</td><td>RE19</td></tr><tr><td>STIS</td><td>2010-02-04</td><td>0.019</td><td>0.030</td><td>RE19</td></tr><tr><td>STIS</td><td>2010-05-28</td><td>0.031</td><td>0.048</td><td>RE19</td></tr><tr><td>STIS</td><td>2015-04-18</td><td>0.043</td><td>0.057</td><td>RE19</td></tr><tr><td>XSH</td><td>2010-05-03</td><td>0.029</td><td>0.030</td><td>M14</td></tr><tr><td>XSH</td><td>2010-04-07</td><td>0.026</td><td>0.027</td><td>V19</td></tr></tbody></table>

**Note.** RE19: Robinson & Espaillat (2019) M14: Manara et al. (2014) V19: Venuti et al. (2019), 60 s spectrum

side of the equation) adjusts for the ~20% difference between the veiling measured from the high- and low-resolution spectra (described above). This equation now allows us to convert our set of 1169 veiling measurements to an accretion luminosity.

The bolometric correction depends on the shape of the accretion continuum, which may vary with time and accretion rate. Figure 7 shows that at the highest veiling measurements, the shape of the accretion continuum (after subtracting off the photosphere) becomes bluer. This same result is seen in our fits to the flux-calibrated low-resolution spectra (Figure 2). The size of the Balmer jump in the accretion continuum is also smaller when the accretion rate is higher, a trend that is consistent with higher temperatures or opacities in the accretion shock (e.g., Calvet & Gullbring 1998). At low-modest veilings, when $r_{5050} < 1$, the consistency in the shape of the accretion continuum indicates that the increase in accretion rate is at a constant temperature, so the spot size likely increases. At high veilings, the bluer spectrum indicates that the increase in accretion rate corresponds to a higher temperature.

## 5.2. Final Calculation of Continuum-based Accretion Rates

Figure 9 shows a histogram of accretion rates, as measured from the continuum luminosity (see also Table 5). In our time
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Correlations between accretion luminosity and line luminosities (red circles from ESPaDOnS, black circles from other instruments), with a best-fit (solid pink line) and a comparison to the best-fit correlation (dashed blue line) for the relationships from Alcalá et al. (2017). For Hα 10% width, the comparison is to the correlation calculated by Natta et al. (2004), adjusted to accretion luminosity given the parameters for TW Hya. Some of the lines are poorly correlated with the accretion luminosity.</figcaption></figure>

series from 1998–2022, the average veiling is 0.73 with a standard deviation of 0.36, the average accretion luminosity is 0.045 $L_\odot$ (standard deviation of 0.022 $L_\odot$), and the average accretion rate is $2.51 \times 10^{-9} M_\odot \text{ yr}^{-1}$ (standard deviation $1.25 \times 10^{-9}$); the average $\log M_{\text{acc}}/M_\odot \text{yr}^{-1}$ is $-8.65$ (standard deviation 0.22). The full range of veiling spans from 0.13–2.13, corresponding to accretion rates from $0.46–7.4 \times 10^{-9} M_\odot \text{ yr}^{-1}$. The veiling never goes below zero or even close to the uncertainty in the zero-point, so excess continuum emission due to accretion is detected in every spectrum.

These averages are calculated after reducing the contributions of observations that were obtained contemporaneously. Veiling measurements are randomly selected into a sample. When one point is selected, all other points obtained within 3 hr are excluded. The final histogram and averages are obtained by averaging 10,000 different selections. An average without any weighting would lead to slightly higher measurements because many Magellan/MIKE spectra were obtained on a night when the veiling was high. The average of the logarithm of the accretion rate leads to a slightly smaller accretion rate than the linear average.

The variability of the log of the accretion rate is ~0.22 dex, as found here from the distribution and also in Section 7.2 from a structure function analysis. This variability is consistent with

the variability found previously in most single-object analyses and surveys (e.g., Biazzo et al. 2012; Costigan et al. 2014; Fiorellino et al. 2022; Zsidi et al. 2022), though EX Lup-type objects have bursts that are much larger than found in generic samples (see review by Fischer et al. 2023). The factor of 90 change in the accretion rate of XX Cha (Claes et al. 2022) is also larger than the full range of accretion rates measured here for TW Hya, so some accreting young stars seem to have accretion rates that are more variable than that of TW Hya.

## 6. Line Emission and Accretion Rate Variability

The luminosity of emission lines in CTTS spectra is often used to measure accretion through correlations with accretion luminosity. These correlations were developed with single-epoch spectra from a large number of stars and are sometimes applied to spectral monitoring of spectroscopic features of individual stars to infer variability. However, the response of lines to changes in accretion rate for an individual star does not necessarily follow the same relationship as found for the global correlations obtained in large samples of stars.

In this section, we study the reliability of these correlations as a variability indicator for a single star, TW Hya. Correlations between line and accretion luminosities are described in Table 8
<table><caption>**Table 8**
Empirical Parameters of Accretion Rate Indicators<sup>a</sup></caption><thead><tr><th>Line</th><th>a</th><th>σ(a)</th><th>b</th><th>σ(b)</th><th>σ</th><th>p<sup>b</sup></th></tr></thead><tbody><tr><td>Hα λ6563</td><td>9.3</td><td>1.0</td><td>4.7</td><td>0.4</td><td>0.41</td><td>0.41</td></tr><tr><td>Hα 10% width</td><td>-2.82</td><td>0.05</td><td>0.0034</td><td>0.0001</td><td>0.18</td><td>0.60</td></tr><tr><td>Hβ λ4860</td><td>6.72</td><td>0.38</td><td>2.47</td><td>0.11</td><td>0.29</td><td>0.53</td></tr><tr><td>Hγ λ4340</td><td>5.18</td><td>0.24</td><td>1.81</td><td>0.07</td><td>0.19</td><td>0.72</td></tr><tr><td>Paschen-10 λ8598</td><td>2.91</td><td>0.12</td><td>0.85</td><td>0.02</td><td>0.20</td><td>0.58</td></tr><tr><td>He I λ4471</td><td>4.39</td><td>0.17</td><td>1.23</td><td>0.04</td><td>0.17</td><td>0.75</td></tr><tr><td>He I λ5876</td><td>4.45</td><td>0.17</td><td>1.30</td><td>0.04</td><td>0.18</td><td>0.69</td></tr><tr><td>He I λ6678</td><td>5.11</td><td>0.20</td><td>1.36</td><td>0.04</td><td>0.16</td><td>0.77</td></tr><tr><td>He II λ4686</td><td>5.35</td><td>0.21</td><td>1.37</td><td>0.04</td><td>0.14</td><td>0.83</td></tr><tr><td>O I λ6300</td><td>36</td><td>5</td><td>7.5</td><td>1.0</td><td>0.54</td><td>0.35</td></tr><tr><td>O I λ8446</td><td>8.12</td><td>0.40</td><td>2.03</td><td>0.09</td><td>0.30</td><td>0.46</td></tr><tr><td>Ca II λ8498</td><td>5.90</td><td>0.27</td><td>1.52</td><td>0.06</td><td>0.18</td><td>0.71</td></tr></tbody></table>

**Notes.**

<sup>a</sup> log L<sub>acc</sub> = a + b × log L<sub>line</sub> (or log L<sub>acc</sub> = a + b × 10% width).

<sup>b</sup> Pearson correlation coefficient

and Figure 10. The lines are selected to avoid severe blends. We also avoid lines located at wavelengths shorter than 4300 Å, where high veilings lead to uncertain line luminosities. The Ca II λ8498 line is selected as the member of the Ca II infrared triplet that is least blended with Paschen lines. Paschen-10 at 8598 Å is selected as the lowest Paschen line in optical spectra that is not affected by telluric absorption and is not blended with any of the Ca II infrared triplet lines. The line emission is integrated over the entire profile, despite narrow and broad components that may respond differently to changes in accretion rate. We restrict this analysis to ESPaDOnS, FEROS, HARPS, UVES (4900–7000 Å), ESPRESSO, CHIRON, and HIRES spectra, since they have the most reliable measurements and consistent set of reductions.

As with the veiling measurements, this analysis assumes a constant photospheric flux. Most of these lines can also be seen in emission in chromospherically active stars. The equivalent width is measured only after subtracting off a chromospheric template, which provides a natural correction for the chromospheric emission. The chromospheric line emission is much fainter than line emission from accretion and will have a negligible effect on this analysis.

## 6.1. Line Luminosity Measurements

Line luminosities are estimated by first measuring veiling-corrected equivalent widths and subsequently converting them to luminosities from flux-calibrated photospheric templates. The equivalent width of each line is measured by integrating emission in the line profile after subtracting off a high-resolution spectrum of the photospheric template TWA 25. The measured equivalent width is then corrected for veiling by multiplying by $1+r_\lambda$ and converted to a flux using the scaled photospheric template. This correction is necessary here but not in other studies (e.g., Alcalá et al. 2017), because most of the high-resolution spectra presented here are not flux calibrated. The assumption in our calculation is that the photospheric emission is constant.

When subtracting the photosphere, the spectra of TW Hya are first rotationally broadened to match the $v \sin i$ of TWA 25. Then either the TWA 25 or TW Hya spectrum (depending on the instrument resolution) is broadened with a Gaussian profile to account for the minor differences in resolution between, for example, HARPS and ESPaDOnS, so that the line widths

match. The TWA 25 spectra are then shifted to the velocity of TW Hya. The stellar photosphere is then scaled to each line location based on the veiling measurement for that spectrum at 5000–5100 Å. The veiling at the respective wavelength is calculated based on the wavelength dependence developed from the ESPaDOnS spectra in Section 3.

For most lines, the equivalent width is calculated by integrating emission from $\pm 500$ km s$^{-1}$, wide enough to include emission in the line wings. The continuum spans 5–10 Å intervals around each line, selected to avoid strong features. For [O I] λ6300, the equivalent width is measured from a best-fit Gaussian profile to the line.

For Hα, the integration region$^{51}$ covers $\pm 1950$ km s$^{-1}$ from the line center because of the extreme line width (see Section 6.4.1), while the continuum region is measured from $\sim 3000–4000$ km s$^{-1}$. The 10% width for the Hα line is measured by smoothing the line profile by $\sim 15$ km s$^{-1}$ and calculating the location where the continuum-subtracted flux exceeds 10% of the peak value. For analysis of Hα, we exclude all UVES and CHIRON spectra plus 13 FEROS and HARPS spectra with an Hα peak that is at or near saturation.

The [O I] λ6300 line may be contaminated by telluric emission. The analysis excludes CHIRON, which has no sky subtraction and is at a spectral resolution low enough that the [O I] sky emission blends with the source emission.

## 6.2. Relationships between Line and Accretion Luminosities

The correlations between veiling and line luminosities and equivalent widths are presented in Table 8 and Figure 10, with linear fits to log L<sub>acc</sub> = a × b log L<sub>line</sub>. The quality of these relationships is evaluated from the standard deviation of the difference between measured log L<sub>acc</sub> and that estimated from the best-fit line and also from the Pearson correlation coefficient. The standard deviation of all measured log L<sub>acc</sub> is $\sim 0.22$ dex, so any useful relationship would need to predict the log L<sub>acc</sub> with a standard deviation significantly lower than 0.22 dex. As a general rule, the Pearson correlation coefficient is higher than 0.7 for strong correlations and less than 0.4 for poor or no correlation.

Most lines provide some limited predictive power. The He lines are well correlated with accretion luminosity, with best-fit relationships that have significant power above random noise. The Hγ and Paschen-10 line luminosities are also well correlated with accretion and provide some limited predictive power. Presumably H lines from even higher energy levels also scale well with accretion. Hβ is well correlated with accretion luminosity, though the scatter is higher than that for other H lines and the correlation has a different slope than that measured in the sample of accretors by Alcalá et al. (2017). The luminosity of Hα is only poorly correlated with accretion luminosity and offers no predictive power. The 10% width of Hα is well correlated with accretion luminosity but with a relationship that provides only modest power. Neither the [O I] λ8446 line, likely produced by the heated chromosphere, nor the [O I] λ6300 line (see Section 6.4.3) are correlated with accretion variability.

The Ca II line emission correlates with accretion luminosity but is much weaker than expected from the Alcalá et al. (2017) relationship, which may be related to disk structure (Micolta et al. 2023). The accretion flow of TW Hya is deficient in elements

<footnote>$^{51}$ The [N II] λ6548 and λ6583 lines are not detected and do not contribute to these measurements.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Hα, Hβ, and two He I lines from an epoch with strong accretion (solid black line) and an epoch with weak accretion (red dashed line). The spectra shown here have had the photosphere subtracted and are then normalized to the relevant photospheric continuum to place the lines on the same flux scale. The He lines are brighter but the H lines are fainter during the epoch with stronger accretion, illustrating by example that accretion is more tightly correlated with He emission than with H emission. The H lines are formed in an optically thick medium and have emission that may be absorbed by the wind or accretion flow.</figcaption></figure>

such as Fe and O (e.g., Kastner et al. 2002; see also speculation about Si depletion by Herczeg et al. 2002) and may likewise be deficient in Ca. As described in models of Kama et al. (2016) and Booth & Clarke (2018), volatile and refractory species may be preferentially in ices and grains and not in accretion flows. However, models of the accretion flow and shock are needed to accurately measure abundances to determine whether the weak Ca II lines are best explained by an underabundance or by the physical conditions of the emission region.

To illustrate the utility of He lines and futility of H lines (see also phrasing from Fischer et al. 2023) as accretion rate indicators, Figure 11 compares H and He lines from a strong accretion spectrum from 2016 February 27 with a weak accretion spectrum from 2014 February 14. The He lines are all stronger during the high accretion rate epoch, as expected, while the Hα lines are actually stronger during the epoch of the low accretion rate. The difference in the Hα line flux is explained by high line opacity and absorption in the winds and accretion flows. However, any empirical correction in the line profile by, for instance, scaling the weak accretion spectrum to fill in missing flux near the line center, would only lead to the two epochs having similar line fluxes.

For TW Hya, the fractional change in line luminosity is smaller than the fractional change in the accretion luminosity. The correlations from Alcalá et al. (2017) have slopes of ~1, consistent with the idea that an increase in accretion rate corresponds to an equal increase in line luminosity. The Alcalá et al. (2017) correlations were developed from X-Shooter surveys of accreting stars in the Lupus star-forming regions across a wide range of stellar mass. While those correlations are robust across different objects, applying those correlations to the time series observations of TW Hya would overestimate the level of variability.

Some scatter may be introduced in the correlations because the response in some lines may be delayed from the continuum (Dupree et al. 2012), as well as by blueshifted absorption in the wind and redshifted absorption in the accretion column. The H emission line emissivity extends across the ~3.5 R\* magnetospheric cavity (Garcia Lopez et al. 2020), so stronger emission in H lines may precede stronger continuum emission by the freefall

time, ~0.3 day. The H lines may also vary less because the accretion flow and streams are roughly constant, while the accretion shock itself is a near-instantaneous measurement. With multiple accretion streams, some components and indicators will lag others (Espaillat et al. 2021; Robinson et al. 2022).

Many line profiles have a broad component from the accretion flow and shock and a narrow component from the heated photosphere (e.g., Yang et al. 2005; Donati et al. 2014). Variability in line profiles is not analyzed in detail in this paper, despite potential trends of stronger emission in broad profiles during periods of high accretion. However, in Section 6.4 we present some limited results on line profiles.

## 6.3. Total Luminosity of Emission Lines

The accretion rates measured here and in most other studies are calculated from the continuum emission. Even when emission lines are used, they serve as a proxy for the continuum emission. The energy that escapes in lines is usually excluded from these calculations, despite significant potential contributions from the lines (see discussion in Alcalá et al. 2014). For TW Hya, this exclusion leads to an underestimate of accretion rates of ~70%, based on the following calculations.

At optical wavelengths, the H Balmer series dominates the emission line luminosity. Because the Hα line luminosity is poorly correlated with the accretion (continuum) luminosity, the ratio of the Hα luminosity to the accretion luminosity is ~0.5 at the lowest accretion rates to ~0.05 at the highest accretion rates. The average Hα luminosity is 0.15 times the accretion luminosity. Other optical and near-IR lines, including the Balmer series and Ca II H and K lines add another ~0.05, so this leads to an extra 20% in accretion.

The ultraviolet adds even more energy (see tabulation by Herczeg et al. 2004). From 1230–3000 Å, the emission line luminosity is ~8% of the accretion luminosity. However, the Lyα line dominates the UV flux. In five observations, the detected Lyα emission is 10% of the continuum accretion luminosity at high accretion rates and 30% at low accretion rates. Moreover, about half of the Lyα flux is unobserved,
<figure><img src="image_15.png" /><figcaption>**Figure 12.** The Hα profile shows strong Lorentzian wings with emission that extends to ~2000 km s⁻¹ on both sides of the line profile. The spectrum shown here (black), from an epoch with strong accretion, is normalized by the photospheric continuum level. The photospheric template TWA 25 (blue spectrum) is scaled to the appropriate level based on the veiling and is then subtracted from the full spectrum to calculate the accretion spectrum (red).</figcaption></figure>

obscured by circumstellar and interstellar H I. Corrections lead to total Lyα fluxes that are roughly 2 times higher than the observed flux (e.g., Herczeg et al. 2004; Schindhelm et al. 2012; Arulanantham et al. 2023).

The X-ray luminosity, as measured by Kastner et al. (2002) and dominated by accretion (e.g., Brickhouse et al. 2010; Argiroffi et al. 2017), adds only ~1% to the total accretion luminosity.

## 6.4. Behavior of Specific Emission Lines

### 6.4.1. Extremely Broad Wings on H Lines

Some emission in Hα extends to ±2000 km s⁻¹ (see example in Figure 12), far beyond the full width at 10% of the peak flux of ~400 km s⁻¹ (200 km s⁻¹ on both sides of the line center). The weak line wings are consistent with a Voigt profile with a damping parameter of ~40 km s⁻¹. The total flux outside the 10% width is typically ~17.5% (with a standard deviation of 2.6%) of the total flux in the Hα line. These broad wings are seen consistently across all instruments and levels of veiling.

The Lyα emission from TW Hya has weak wings that extend to even larger velocities, as noted by France et al. (2014). The fluxes on the prominent line wings are well reproduced by models with superimposed Gaussian profiles to represent the stellar and accretion emission and interstellar medium and outflow absorption (e.g., Herczeg et al. 2004; Schindhelm et al. 2012) or resonant scattering through an outflowing H I shell (see analysis of DM Tau by Arulanantham et al. 2023), but those models do not capture the extremely high-velocity wings out to 2000 km s⁻¹. The Hβ line wings extend to ~1200 km s⁻¹, while higher H lines do not have such broad wings (see, e.g., Wilson et al. 2022). This pattern is consistent with expectations for pressure broadening in optically thick lines.

Very broad (several thousands of kilometers per second) wings of Hα have also been detected in symbiotic stars (Van Winckel et al. 1993; Ivison et al. 1994; Selvelli & Bonifacio 2000; Skopal 2006) and planetary nebula (Arrieta & Torres-Peimbert 2003; Miranda et al. 2022). The structure of the objects is similar to T Tauri stars because a hot ionized

<figure><img src="image_16.png" /><figcaption>**Figure 13.** He I absorption, seen here to ~300–400 km s⁻¹ in three lines on 2016 February 27, an epoch with strong accretion. The He II λ4686 line (bottom, pink) has emission that extends out to 450 km s⁻¹.</figcaption></figure>

medium and an optically thick neutral medium coexist. Lee (2000) suggested that the broad Hα wings in symbiotic stars follow the Lorentzian profile if the wings are formed by Raman scattering with atomic hydrogen of UV radiation near Lymanβ. Furthermore, Chang et al. (2018) showed that the broad Raman wing near Hα is broader than that near Hβ. Thus, the broad wings near Hα in TW Hya might originate from Raman scattering with atomic hydrogen in the optically thick protoplanetary disk, although further calculations are required to test this hypothesis.

### 6.4.2. Inverse P Cygni Absorption in He I

He lines provide a powerful probe of gas dynamics (e.g., Beristain et al. 2001). Figure 13 shows that the He I λ5876 line exhibits inverse P Cygni profiles during some epochs of strong accretion, with velocities on the red wing that extend to ~350 km s⁻¹. This detection is challenging because the redshifted absorption blends with the broad wings of photospheric Na I D absorption. However, this detection is supported by the detection in other He I lines of redshifted self-absorption at ~300 km s⁻¹. For classical T Tauri stars, self-absorption in the optical He I lines had only been previously identified in a couple of objects (Beristain et al. 2001).

For comparison, He II λ4686 shows emission on the red wing out to ~450 km s⁻¹, similar to the maximum velocities seen in C IV and other high-temperature lines (e.g., Herczeg et al. 2002; Ardila et al. 2013; Dupree et al. 2014). This fast velocity must mean that most of the gas is crashing onto the star along our line of sight to the star. In an accretion geometry with tongues at equatorial or midlatitudes, the gas would flow more along the plane of the sky, so we would not detect gas at such high velocities.
<figure><img src="image_17.png" /><figcaption>**Figure 14.** The [O I] λ6300 line luminosity (blue circles from FEROS, black squares from other instruments) varies with a range of almost a factor of 2 on timescales of months. CHIRON spectra are excluded from this analysis.</figcaption></figure>

Most past work on He lines focused on the He I λ10830, which has a metastable lower level that leads to P Cygni (blueshifted) and inverse P Cygni (redshifted) absorption components (e.g., Dupree et al. 2005; Edwards et al. 2006; Fischer et al. 2008; Erkal et al. 2022). The lower level of the He I λ5876 line is the upper level of the He I λ10830, so the absorption in He I λ10830 could help to populate the level and lead to the absorption in the λ5876 line. A comprehensive analysis of He lines could help to evaluate the excitation and ionization of the accretion flow.

### 6.4.3. Time Variability in [O I] Emission

The [O I] λ6300 line profile is narrow and consistent with the low-velocity component of the wind or innermost disk (e.g., Fang et al. 2018; Pascucci et al. 2020). Figure 14 shows that the [O I] λ6300 emission varies on timescales of months to years. The line equivalent width, corrected for veiling, varies from an epoch-averaged minimum of ~0.65 Å to a maximum of 1.15 Å, corresponding to luminosities of (0.8–1.3) × 10⁻⁵ L⊙.

These temporal changes are likely caused by real changes in the line strength. Fang et al. (2023) evaluate that the centroid of the narrow line core is stable to within ~1 km s⁻¹ but find that emission wings that extend to ~30 km s⁻¹ on both the red and blue sides of the line are sometimes present and sometimes absent.

Some modest dispersion during each epoch may be caused by differences in the spot coverage on the visible surface of the star, since a change in spot coverage will change the flux of the photosphere at 6300 Å but will not change the flux in the line. However, the long-term temporal changes are too large to be explained by differences in spot coverage. The ESPaDOnS spectra include telluric [O I] emission, which in some cases may blend with the stellar emission and in some epochs may be significant. The shape of the veiling continuum does not vary significantly enough to explain these changes.

## 7. Distribution and Variability of Accretion onto TW Hya

### 7.1. Average Accretion Rate

The full veiling lightcurve of TW Hya is presented in Figure 15. The average accretion rate of TW Hya is 2.51 × 10⁻⁹ M⊙ yr⁻¹, as measured from the continuum

emission in Section 5.2. We adopt this measurement for consistency with previous measurements, however, the exclusion of line emission and the treatment of the accretion flow as a slab rather than a multiple-column flow may cause us to underestimate the accretion rate by a factor of 2–3.

The distribution of the log accretion rate is reasonably well described by a Gaussian profile with an average of −8.65, an FWHM of 0.22 dex, and an excess tail at low accretion rates. In our time series, TW Hya never stops accreting. The veiling relative to non-accreting templates is always >0.13 at 5000–5100, significantly larger than the uncertainty in the zero-point. However, visual inspection of the lightcurve reveals that some epochs have low accretion, including 2014 February, and some have high accretion, including 2016 February and 2022 March–April. In addition, the few nights with extensive monitoring show large changes during the night, as originally reported for Magellan/MIKE monitoring (Dupree et al. 2012) and seen here also in UVES monitoring. Much of the nightly variation seen across the lightcurve may be explained by hours-long bursts, seen in a few high-cadence epochs here as well as in photometric monitoring with MOST (e.g., Siwak et al. 2014, 2018).

TW Hya has often been considered a very weak accretor, an interpretation that is not supported by our comparison to unbiased data sets. Figure 16 compares the accretion rate from TW Hya with accretion rates measured for complete samples of stars in Lupus and Cha I (Alcalá et al. 2017; Manara et al. 2017), updated with all stellar properties calculated with the Gaia DR3 distance (inverted parallaxes from Brown et al. 2021) and placed on the mass scale from Somers et al. (2020) tracks with 50% spots. The accretion rate of TW Hya is only 0.15 dex lower than the median accretion rate expected for a ~0.87 M⊙ star in Lupus or Cha I, despite the much older age of TW Hya.

The misconception that TW Hya is a weak accretor has two explanations. First, the accretion models of Muzerolle et al. (2000) fit to the Hα line yielded an accretion rate of 5 × 10⁻¹⁰ M⊙ yr⁻¹, 0.75 dex lower than our median accretion rate. Muzerolle et al. (2000) measured a veiling of ~0.2 at 7000 Å, consistent with our veiling measurement<fnref n="52" /> of 0.49 at 5000–5100 Å and yields an accretion rate of 1.69 × 10⁻⁹ M⊙ yr⁻¹, or 0.17 dex lower than the average accretion rate. The remaining difference of ~0.6 dex between the Muzerolle et al. 2000 and our accretion rates is caused by methodological differences, with our measurements from the accretion luminosity and Muzerolle et al. (2000) from models of the Hα line profile. This difference is consistent with the offset in the correlation between accretion luminosity and the Hα 10% width for TW Hya, in comparison to the relationship developed by Natta et al. (2004), as can be seen in Figure 10. Second, TW Hya may have also been considered a weak accretor because early accretion rates were measured for biased samples. The median accretion rate is 9.6 × 10⁻⁹ M⊙ yr⁻¹ in Gullbring et al. (1998) and 10.3 × 10⁻⁹ M⊙ yr⁻¹ in Valenti et al. (1993), a factor of 3–5 higher than expected for 0.6–1.0 M⊙ stars from the complete surveys in Lupus and Cha I by Alcalá et al. (2017) and Manara et al. (2017).

Our average accretion rate is remarkably consistent with the accretion rate of 1.67 × 10⁻⁹ M⊙ yr⁻¹ derived from X-ray spectroscopy (Brickhouse et al. 2012, increased here by 1.11 to account for the slight change in adopted distance). These X-ray

<footnote n="52">⁵² Measured for lines from ~6200–6500 Å; the spectrum does not cover wavelengths below the 6000.</footnote>
<figure><img src="image_18.png" /><figcaption>**Figure 15.** The full veiling lightcurve (top), with specific segments highlighted in the bottom six panels, with average veiling (horizontal line) and standard deviation (shaded yellow region). In the six panels, the date range in the x-axis differs in each plot.</figcaption></figure>

estimates use hydrodynamic models and account for the H I column density that attenuates the shock, with an interpretation supported by the measurement of redshifted X-ray lines by Argiroffi et al. (2017). The accretion rate measured in the third X-ray integration is at the low end of our distribution but coincides with epochs with low veiling, with the X-ray-based accretion rate a factor of 2 lower than our measurements. This

agreement is remarkable, since the X-ray method depends only on line ratio diagnostics from He-like Ne IX for electron density, electron temperature, and absorption column density and is independent of the X-ray luminosity. Previous estimates yielded lower accretion rates because the H I column density was underestimated (e.g., Stelzer & Schmitt 2004; Günther et al. 2007).
<figure><img src="image_19.png" /><figcaption>**Figure 16.** Accretion rate of TW Hya (red rectangle), compared with accretion rates from Lupus and Cha I star-forming regions (black circles, Alcalá et al. 2017; Manara et al. 2017). The size of the thick rectangle for TW Hya shows the 1σ scatter of the accretion rates while the thin line shows the full range of accretion rates. The average accretion rate of TW Hya is only ~0.1 dex fainter than and consistent with expectations, given the uncertainties, based on a linear fit to the relationship between log M<sub>acc</sub> and log M<sub>*</sub>.</figcaption></figure>

The accretion luminosity is calculated from the continuum
alone and excludes emission lines (see Section 6.3). The
accretion luminosity would be ~70% higher if emission lines,
in particular Lyα and Hα, were included as accretion
luminosity. The H line luminosities are relatively steady
compared to the fluctuations in the continuum accretion
luminosity, which means that this correction is much larger
for low accretion rates than for high accretion rates. Since the H
lines are produced over a larger area than the accretion shock
(Garcia Lopez et al. 2020), they should provide a more steady
source of emission that may not reflect the short-term changes
in accretion rate. The exclusion of lines is even more important
toward very low-mass stars and brown dwarfs, since the optical
line luminosity to UV continuum luminosity approaches unity
(Alcalá et al. 2014; Zhou et al. 2014). The Lyα luminosity for
most of those sources is unmeasured and may be higher than
Hα and other lines in the Balmer series (Arulanantham et al.
2023).

A second significant uncertainty in our average accretion rate
is the bolometric correction, especially in the correction for
low-density accretion flows (e.g., Ingleby et al. 2013; Robinson
& Espaillat 2019; Espaillat et al. 2021). This uncertainty affects
all epochs of TW Hya and all comparison measurements. The
primary uncertainty in the lower bound of accretion rates is the
zero-point, where the ~0.02 absolute veiling uncertainty leads
to a ~20% error in the accretion rate. However, the zero-point
is a negligible (<5%) uncertainty for epochs of average (or
higher) accretion rates. The width of the accretion rate
distribution in Figure 9 is dominated by real variability.

The present-day disk mass within 5 au is ~7 M<sub>J</sub>, as
measured from the wings on CO lines uniquely for TW Hya
(Yoshida et al. 2022), and ~50 M<sub>J</sub> for the entire disk (Bergin
et al. 2013). If the present-day average accretion rate has been
constant for 10 Myr, then the total mass accreted onto the star
to date would be ~26 M<sub>J</sub>; most of the mass was likely accreted
when TW Hya was very young. The initial and even present-
day disk mass of TW Hya is near the upper envelope of disk
masses around a ~0.87 M<sub>☉</sub> star, as measured from complete
samples of dust emission at millimeter wavelengths (see

comparisons to the compilation by Manara et al. 2023). If the
accretion onto the star stays steady and the inner disk is not fed
additional material from the outer disk, then the inner disk has
enough material to survive for another 3.8 Myr. The current
accretion rate is also a factor of ~10 higher than the
photoevaporation rate of 2 − 3 × 10 − 10 M ⊙ yr − 1 (Pascucci
et al. 2012), although the MHD wind may carry significant
mass (e.g., Wang et al. 2019). These values would all be
modified if the line emission increases the accretion rate
by 70%.

7.2. Reset Timescale for Accretion Rate

To evaluate any characteristic timescales of accretion, we
first focus on a robust determination of the nonperiodic
variations using structure functions. The difference in accretion
rates is compared with the time separation for every two data
points in our sample. This so-called “Δm−Δt” analysis is
widely used to study time variations, including for quasars
(e.g., De Vries et al. 2005), for a similar spectroscopic analysis
of a young accreting star, DF Tau, by Johns-Krull & Basri
(1997), and also on larger samples of accretion rates and
photometry by Costigan et al. (2014), Zsidi et al. (2022), and
Sergison et al. (2020). A detailed description of structure
functions and their limitations for young stars is provided by
Findeisen et al. (2015). A monotonic increase of Δm, the
difference in physical property, with respect to Δt, the time
difference, such that the system is further away from the start
point as time passes, indicates that the system is in a varying
mode. However, if the Δm−Δt relation is flat, it suggests that
the system relaxes and the variation occurs on a shorter
timescale.

For our data set, the distribution of differences in accretion
rates, Δ log ˙M, is given as

$$
|\Delta \log \dot{M}| = |\log \dot{M}_1 - \log \dot{M}_2| = |\log \frac{\dot{M}_1}{\dot{M}_2}|, \quad (3)
$$

where ˙M₁ and ˙M₂ are the accretion rates of each pair of data
points as a function of the time difference between those two
data points. Despite the large diversity of timescales in this data
set, we do not have enough high-cadence data at short
timescales to reliably distinguish the difference between
increasing ˙M and decreasing ˙M; hence, the absolute value of
|Δ log ˙M| is used.

Figure 17 shows the |Δ log ˙M|−|Δt| distribution of the 692,076 pairs from the 1169 ˙M measurements. The median |Δ log ˙M|−|Δt| in bins of |Δt| are presented in orange (binned by every 500 pairs) and green (binned by Δlog t = 0.033) points. At short timescales (≈1 day), the median |Δ log ˙M| increases consistently with |Δt|. However, on long timescales (≈1 day), the median value flattens with respect to |Δt| and is nearly flat at ~0.2 dex, despite a large spread in the |Δ log ˙M| values. This analysis leads to the same level of variation as the distribution of accretion rates presented in Figure 9, with the additional conclusion that the accretion rates are (usually) independent from one day to the next. To quantify the turnaround timescale, we fit the median binned |Δ log ˙M|−|Δt| relations, from two binning methods combined, with a
<figure><img src="image_20.png" /><figcaption>**Figure 17.** Top: the so-called “$\Delta m-\Delta t$” diagram of the log $\dot{M}$ variability. Each black point on the background represents a pair of $\dot{M}$ measurements, with x values being the absolute difference in time ($|\Delta t|$), and y values being $|\Delta \log \dot{M}|$. Individual pairs are bins in the order of $|\Delta t|$ by every 500 pairs (orange points) or by every $|\Delta t| = 0.033$ (green points). The purple dashed line denotes the two-part function that fits the binned points best, with a break at ~1.6 days from a power-law increase to a flat line (see Equation (4) and Table 9 for details). Bottom: the periodogram of the log $\dot{M}$ variability. The horizontal dashed and dotted lines represent the 5% and 1% false-alarm rates, respectively.</figcaption></figure>

two-part function,

$$|\Delta \log \dot{M}| = \begin{cases} A_0 |\Delta t|^\gamma, & |\Delta t| < \tau_0; \\ A_1 \log\left(\frac{|\Delta t|}{\tau_0}\right) + A_0 \tau_0^\gamma, & |\Delta t| \ge \tau_0, \end{cases} \quad (4)$$

where $A_0$, $A_1$, $\gamma$, and $\tau_0$ are set as free parameters. For the equation, the unit of $\dot{M}$ is $M_{\odot} \text{yr}^{-1}$ and the units for $\Delta t$ and $\tau_0$ are days. The best-fit results are listed in Table 9.

According to the fitting result, most of the accretion variability occurs on a timescale within $\tau_0 = 1.6 \pm 0.4$ days. On longer timescales, the accretion rate appears randomly distributed around the mean value. However, the $|\Delta \log \dot{M}| - |\Delta t|$ relation is not entirely flat, with a positive slope ($A_1 = 0.009 \pm 0.003$) that suggests the mean accretion rate slowly (either incoherently or coherently) drifts away from the starting point.

Findeisen et al. (2015) found that the timescale in the $\Delta m-\Delta t$ plot is reliable when the time value is ~30 times greater than the sampling interval and ~1/15 of the total time baseline. The sampling of our data has a minimum interval = 0.008 day and a time baseline of ~25 yr. Therefore, the 1.6 day timescale is generally robust.

TW Hya has been previously found to have periodic (or quasiperiodic) behavior with 2–4.5 day periods, attributed to a combination of accretion hot spots and dark starspots on the stellar surface (e.g., Mekkaden 1998; Lawson & Crause 2005;

<table><caption>**Table 9**
Parameters of the Two-part Function</caption><thead><tr><td>Name<sup>a</sup></td><td>Best Fit</td></tr></thead><tbody><tr><td>A<sub>0</sub></td><td>0.15 ± 0.01</td></tr><tr><td>A<sub>1</sub></td><td>0.009 ± 0.003</td></tr><tr><td>γ</td><td>0.50 ± 0.08</td></tr><tr><td>τ<sub>0</sub></td><td>1.6 ± 0.4</td></tr></tbody></table>

**Note.**

Huélamo et al. 2008; Rucinski et al. 2008). Possible periodic signals could occur on short timescales related to oscillations of the accretion flow (Sacco et al. 2010), rotational timescales, turbulent timescales (Robinson et al. 2021), or potentially a viscous timescale, as indicated by a preliminary analysis by Takasao et al. (2022).

We search for these periods in accretion rate by using the generalized Lomb–Scargle periodograms (Lomb 1976; Scargle 1982; Zechmeister & Kürster 2009). The lower panel of Figure 17 shows the periodogram with the false-alarm probability (FAP) estimated from bootstrap. At high frequencies (period $T \lesssim 1$ day), most power is within the background level with low significance, except a peak with FAP <1% at $T \sim 0.004$ days. We examine this signal with phase-folding diagrams and find no convincing periods. Furthermore, the period, which corresponds to ~5 minutes, is roughly the

<footnote><sup>a</sup> See Equation (4) for definitions.</footnote>
cadence of individual exposures for consecutive observations. Therefore, we conclude that this signal is likely due to the sampling window function of the observations. At lower frequencies ($T \gtrsim$ 1 day), the power spectrum becomes noisy with a few prominent peaks close to the observational cadences of 1 day, 29 days, and 365 days. Therefore, despite the power spectrum being consistently above the 1% FAP, the signal is likely due to the combination of these observational cadences.

### 7.3. Interpreting the Timescales for Variable Accretion onto TW Hya

High-frequency photometric monitoring of TW Hya (e.g., Siwak et al. 2018) reveals a lightcurve littered with constant accretion bursts and decays, so many of these timescales may relate to the rise and decay of these bursts. Frequent stochastic bursts are expected for accreting young stars, based on simulations of the magnetospheric geometry (Romanova et al. 2008; Blinova et al. 2016; Takasao et al. 2022). We see only individual points within those stochastic lightcurves, though in a few cases, the veiling lightcurve covers many consecutive hours (e.g., Dupree et al. 2012).

Although accretion fluctuates on any given day, over the long term the average accretion rate appears stable. The structure function indicates that accretion resets on ~1.6 day timescales. Similar analyses from photometry of larger samples indicate that $\tau_0 = 0.25 \times P_*$, where $P_*$ is the rotation period (Sergison et al. 2020; Venuti et al. 2021), which for TW Hya would have suggested a reset timescale of ~0.9 days.

In analysis of high time-resolution photometry, Siwak et al. (2011) found power on a timescale of 1.3 days, which may be a consequence of the timescale for the reset of accretion rates on the 1.6 day timescale found here. This 1.6 day timescale approximates how quickly the accretion rate changes with time, but this change is not periodic. The stellar dipole field is roughly constant, at least on short timescales (e.g., Donati et al. 2011), however, the amount of gas lifted off the disk and accreted onto the star may be constantly changing because of asymmetric structures in the inner disk. The 1.6 day timescale for resetting accretion is similar to the timescales obtained from structure functions in other protoplanetary disks (e.g., Venuti et al. 2021; Zsidi et al. 2022).

## 8. Conclusions

We measure the veiling in 1169 high-resolution spectra of TW Hya obtained over 25 years. These measurements are converted to accretion rates based on scalings obtained from 26 flux-calibrated low-resolution spectra. From analyzing this data set, we find the following results:

1. The veiling at 5000–5100 Å varies from 0.13 to 2.13, equivalent to mass accretion rates from $0.47–9.43 \times 10^{-9} M_{\odot} \text{yr}^{-1}$, with an average of $2.51 \times 10^{-9} M_{\odot} \text{yr}^{-1}$. These accretion rates are only slightly lower than the average accretion rate for accreting stars in nearby young clusters. With its current mass, the disk could survive for another 3.8 Myr. The accretion rate would be 12% higher if we adopted the truncation radius of 3.5 $R_*$ measured by Garcia Lopez et al. (2020).

2. The distribution of accretion rates is well described by $\log \dot{M}_{\text{acc}} = -8.65$ with a standard deviation of 0.22 dex. Accretion never ceases. This variability is less than the scatter seen in accretion rate–stellar mass relationships

(Manara et al. 2023). The variability is consistent with the level of accretion variability of solar-mass stars in many other sources (e.g., Venuti et al. 2021; Zsidi et al. 2022), but with some important exceptions, such as XX Cha (Claes et al. 2022).

3. The uncertainty in our veiling measurements is <0.05 (~0.003 when comparing broadband spectra from the same instrument), which is a minor contribution to the overall error budget. The uncertainty in accretion rates is instead dominated by bolometric corrections and to a lesser extent spots. Comparisons to measurements from Robinson & Espaillat (2019) indicate that the inclusion of multiple accretion streams at lower densities could increase the accretion rate by as much as ~50%. We also find that the accretion spectrum becomes bluer when the accretion rate is higher, which leads to a bolometric correction that depends on the accretion rate. This dependence is not incorporated into our accretion rate measurements.

4. The lines correlate reasonably well with the accretion luminosity. However, the $H\alpha$ luminosity and 10% width are only weakly correlated with accretion luminosity because of higher opacities and absorption in accretion streams and winds and because of possible temporal offsets between continuum and line emission. The use of $H\alpha$ and $H\beta$ as accretion rate indicators for TW Hya is not recommended. These correlations also fail for XX Cha (Claes et al. 2022) and should be used with caution, although at least some sources show strong correlations between accretion and line emission (GM Aur, Bouvier et al. 2023). Long-term variability is found in the [O I] λ6300 line, which traces the inner disk and MHD wind.

5. The line emission is excluded from the calculated accretion rates and is an important correction for TW Hya. The exclusion of lines, and in particular, $Ly\alpha$ emission, from the accretion luminosity would increase the average accretion rate by 70% to $4.3 \times 10^{-9} M_{\odot} \text{yr}^{-1}$. The exclusion of lines such as $Ly\alpha$ emission may lead to systematic underestimates of all accretion rates. Including the lines in instantaneous accretion rate estimates would decrease the scatter in accretion rates.

6. The $\Delta m-\Delta t$ diagram indicates that the reset timescale for the aperiodic accretion variability is ~1.6 days. Large increases in accretion are likely caused by short bursts, as seen in photometry (Siwak et al. 2018). We did not find convincing evidence of periodic signals associated with the bursts. Despite the 1.6 day timescale for resetting the accretion rate, some epochs appear to have different average accretion rates. Accretion onto TW Hya was weaker than average in 2014 February and stronger than average in 2022 March–April.

The accretion variability of TW Hya is likely common and unremarkable in the context of other accretion disks, but this lengthy time series provides new insights into unremarkable variability. The steady accretion likely occurs along dipole-like field lines that land near the pole (e.g., Donati et al. 2011), with high-velocity emission and absorption in He lines that indicate mass loading onto these field lines during bursts. While the accretion rate appears unstable on any given day, the overall accretion rate and fluctuations in that rate have been stable over the past 25 years.
## Acknowledgments

We thank the referee for a careful read and comments, which saved us from publishing a significant error and misinterpretation of possible periods in the veiling. We thank the many people who have contributed time and effort into obtaining this data at the telescope, reducing the data, and developing archives to make the data available. This list includes Ilaria Pascucci, Feng Long, and Claudio Melo. We also thank Kevin France for discussions of Lyα emission. This work benefited from discussions with the ODYSSEUS team (HST AR-16129; Espaillat et al. 2022, https://sites.bu.edu/odysseus/).

G.J.H. and J.T. are supported by the National Key R&D program of China 2022YFA1603102 from the Ministry of Science and Technology (MOST) of China and by general grant 12173003 from the National Natural Science Foundation of China. Z.G. acknowledges support from FONDECYT Postdoctoral 3220029. Z.G. acknowledges support by ANID Millennium Science Initiative Program NCN19_171. Support for this HMG was provided by the National Aeronautics and Space Administration through Chandra Award Number GO1-22007X issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of the National Aeronautics Space Administration under contract NAS8-03060. A.F. and J.A. acknowledge support by the PRIN-INAF 2019 STRADE (Spectroscopically TRACING the Disk dispersal Evolution) and by the Large Grant INAF YODA (YSOs Outflow, Disks and Accretion). This work has been funded by the European Union under the European Union Horizon Europe Research & Innovation Program 101039452 (WANDA). This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No 716155 (SACCRED). J.F.D. acknowledges funding from the European Research Council (ERC) under the H2020 research and innovation program (grant 740651, NewWorlds). SHPA acknowledges financial support from CNPq, CAPES, and Fapemig. J.H.K. is supported by NASA XRP grant 80NSSC19K0292 and NASA ADAP grant 80NSSC22K0625 to RIT.

This work is based in part on observations obtained at the Canada-France-Hawaii Telescope (CFHT), which is operated by the National Research Council of Canada, the Institut National des Sciences de l'Univers of the Center National de la Recherche Scientifique of France, and the University of Hawaii. The CFHT observations were obtained in programs 08AF11, 10AP11, 12Ap12, 14AP18, and 16AP18 (PI Donati) and 15BE97 (PI Malo).

This work is also based in part on observations collected at the European Southern Observatory under ESO programs 106.20Z8 (PI Manara), 074.A-9021 (PI Setiawan), 089-A-9007, 090.A-9013, 092.A-9007 (PI Mohler), 093.A-9029 (PI Gredel), 099.A-9010 (PI Sarkis), 099.A-9008 (PI Mueller), 0101.A-9012 (PI Launhardt), 075.C-0202 (PI Gunther), 081.C-0778, 082.C-0390 (PI Weise), 082.C-0427 (PI Doellinger), 082.C-0218 (PI Melo), 089.C-0299 (PI Pascucci), 085.C-0238 (PI Alcala), 085.C-0764 (PI Güenther), 103.200T (PI Günther), and 60.A-9036 and 60.A-9022 (engineering runs, no PI listed).

This research is based on observations made with the NASA/ESA Hubble Space Telescope obtained from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under

NASA contract NAS 5-26555. These observations are associated with programs GO-8041 (PI Linsky), GO-9093 (PI Johns-Krull), 11608 (PI Calvet), and 13775 (PI Espaillat) and can be accessed via DOI:10.17909/0v8a-vq31.

Some of the data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. This research has made use of the Keck Observatory Archive (KOA), which is operated by the W. M. Keck Observatory and the NASA Exoplanet Science Institute (NExScI), under contract with the National Aeronautics and Space Administration. The Keck data were obtained in programs C199LA (PI Herczeg), C247Hr (PI Carpenter), C252Hr (PI Hillenbrand), C186Hr (PI Hillenbrand), C199Hb (PI Herczeg), and C269Hr (PI Dahm).

Observing time with SMARTS/Chiron obtained by PI Walter was made possible by a Research Support grant from Stony Brook University. We thank Wei-Chun Jao, Leonardo Paredes, and Todd Henry for managing the Chiron spectrograph and their prompt scheduling of the requested observations.

The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.

Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

## Appendix A

### High and Low Veiling Accretion Spectra

In data behind the figures, we provide high S/N spectra by coadding 10–15 spectra of similar veiling, including the high and low veiling spectra and the residual spectrum in Figure 11. These high and low veiling spectra were identified after initial veiling measurements, and then veiling was remeasured to produce the residual spectrum. For each of these spectra, the number of coadds, veiling range, and veiling average are provided in Table 10. The spectra are coadded over every pixel by normalizing the flux over ~30 Å around each pixel.

The residual spectrum is calculated by first measuring the veiling from 5000–5100, $r_{5050}$ in the high accretion spectrum, $f_0$, and the low accretion spectrum, $f_5$. We then apply the correlations between the veiling at other wavelengths and $r_{5050}$.

<table><caption>**Table 10**
Coadded Spectraª</caption><thead><tr><th>Spectrum</th><th>#</th><th>r<sub>5050</sub> Range</th><th>Avg (r<sub>5050</sub>)</th></tr></thead><tbody><tr><td>f<sub>0</sub></td><td>14</td><td>1.1–1.52</td><td>1.25</td></tr><tr><td>f<sub>1</sub></td><td>17</td><td>0.80–0.85</td><td>0.83</td></tr><tr><td>f<sub>2</sub></td><td>12</td><td>0.60–0.65</td><td>0.62</td></tr><tr><td>f<sub>3</sub></td><td>11</td><td>0.4–0.45</td><td>0.43</td></tr><tr><td>f<sub>4</sub></td><td>11</td><td>0.2–0.25</td><td>0.23</td></tr><tr><td>f<sub>5</sub></td><td>17</td><td>-0.01–-0.03</td><td>0.01</td></tr></tbody></table>

Note. Relative veiling values listed.
(as shown in Figure 2), to obtain a veiling spectrum. The veiling spectrum is then converted to flux (following methods in Figure 6), so that excluded wavelengths could be estimated by extrapolating from nearby wavelengths. We then reconvert this flux spectrum into a veiling spectrum, tailored to the photospheric flux level. Finally, we subtract the low veiling spectrum from the high veiling spectrum, after scaling the low veiling spectrum to its estimated contribution to the high veiling spectrum.

The residual spectrum is what remains after this process. Photospheric lines are no longer present in the spectrum. Instead, the spectrum is mostly the flat accretion continuum. However, many narrow emission lines are also present across this spectrum. For many of these lines, the emission fills in the core of the photospheric lines and would not be detectable in a single spectrum, but the excess emission is apparent after subtracting off a high-quality photospheric template, in this case, TW Hya itself during weak accretion epochs.

## Appendix B
Relationships between Veiling Diagnostics and Final Veiling Measurements

This section describes how equivalent widths and spectral indices are used to calculate veiling, with measurements from ESPaDonS spectra. Table 11 lists 82 photospheric absorption lines that are coadded for equivalent width measurements. Each line is normalized by the nearby continuum before being added to other nearby lines. Table 12 describes the lines that are coadded and the relationship between equivalent width and veiling. The lines that are used and coadded differ, depending on spectral coverage and resolution.

Veiling is also measured using the spectral indices in Table 13. Each spectral index is calculated by dividing the median flux in a low-flux region by the median flux of a nearby region with stronger emission. The relationship between each spectral index and veiling is then described with a 4th-order polynomial fit to ESPaDonS data.

<table><caption>Table 11
Equivalent Widths of Photospheric Lines Used for Veiling</caption><thead><tr><th>λ<sub>obs</sub> (Å)</th><th>EW (Å)</th><th>λ<sub>obs</sub> (Å)</th><th>EW (Å)</th><th>λ<sub>obs</sub> (Å)</th><th>EW (Å)</th></tr></thead><tbody><tr><td>4591.587</td><td>0.116</td><td>5426.490</td><td>0.128</td><td>5737.296</td><td>0.093</td></tr><tr><td>4592.824</td><td>0.122</td><td>5434.755</td><td>0.229</td><td>5857.699</td><td>0.197</td></tr><tr><td>4594.290</td><td>0.138</td><td>5455.876</td><td>0.249</td><td>5866.697</td><td>0.116</td></tr><tr><td>4611.438</td><td>0.122</td><td>5490.386</td><td>0.082</td><td>5942.003</td><td>0.100</td></tr><tr><td>4613.517</td><td>0.139</td><td>5497.736</td><td>0.242</td><td>5956.953</td><td>0.090</td></tr><tr><td>4617.468</td><td>0.106</td><td>5501.700</td><td>0.130</td><td>6013.738</td><td>0.097</td></tr><tr><td>4651.480</td><td>0.128</td><td>5506.979</td><td>0.180</td><td>6016.900</td><td>0.113</td></tr><tr><td>4881.779</td><td>0.117</td><td>5535.707</td><td>0.100</td><td>6020.359</td><td>0.081</td></tr><tr><td>4903.494</td><td>0.086</td><td>5569.851</td><td>0.092</td><td>6022.045</td><td>0.114</td></tr><tr><td>4934.267</td><td>0.142</td><td>5573.104</td><td>0.165</td><td>6065.740</td><td>0.119</td></tr><tr><td>4965.066</td><td>0.116</td><td>5582.192</td><td>0.146</td><td>6081.698</td><td>0.102</td></tr><tr><td>4991.293</td><td>0.192</td><td>5588.984</td><td>0.174</td><td>6085.493</td><td>0.108</td></tr><tr><td>5007.420</td><td>0.096</td><td>5590.333</td><td>0.153</td><td>6090.467</td><td>0.103</td></tr><tr><td>5189.038</td><td>0.109</td><td>5594.700</td><td>0.220</td><td>6111.902</td><td>0.103</td></tr><tr><td>5219.918</td><td>0.107</td><td>5601.519</td><td>0.137</td><td>6141.977</td><td>0.126</td></tr><tr><td>5252.306</td><td>0.094</td><td>5603.094</td><td>0.184</td><td>6150.414</td><td>0.099</td></tr><tr><td>5261.925</td><td>0.116</td><td>5671.083</td><td>0.096</td><td>6154.478</td><td>0.092</td></tr><tr><td>5262.455</td><td>0.110</td><td>5672.044</td><td>0.086</td><td>6166.693</td><td>0.102</td></tr><tr><td>5282.017</td><td>0.098</td><td>5682.867</td><td>0.154</td><td>6199.438</td><td>0.108</td></tr><tr><td>5341.218</td><td>0.173</td><td>5688.446</td><td>0.156</td><td>6210.933</td><td>0.112</td></tr><tr><td>5346.033</td><td>0.247</td><td>5698.711</td><td>0.165</td><td>6216.628</td><td>0.114</td></tr><tr><td>5348.536</td><td>0.187</td><td>5703.818</td><td>0.089</td><td>6219.547</td><td>0.087</td></tr><tr><td>5349.724</td><td>0.190</td><td>5707.230</td><td>0.094</td><td>6231.034</td><td>0.148</td></tr><tr><td>5353.646</td><td>0.080</td><td>5709.677</td><td>0.118</td><td>6246.586</td><td>0.081</td></tr><tr><td>5394.902</td><td>0.153</td><td>5727.291</td><td>0.121</td><td>6252.094</td><td>0.111</td></tr><tr><td>5415.437</td><td>0.114</td><td>5727.900</td><td>0.091</td><td>6261.393</td><td>0.126</td></tr><tr><td>5420.589</td><td>0.124</td><td>5731.489</td><td>0.108</td><td>6274.934</td><td>0.093</td></tr><tr><td>6285.420</td><td>0.092</td><td></td><td></td><td></td><td></td></tr></tbody></table>

<table><caption>Table 12
Relationships to Convert Equivalent Widths to Veiling</caption><thead><tr><th rowspan="2">Line Range (Å)</th><th colspan="4">Σa<sub>i</sub> × EW<sub>i</sub></th></tr><tr><th>a<sub>0</sub></th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>a<sub>3</sub></th></tr></thead><tbody><tr><td>4500–5200</td><td>-0.758</td><td>-5.920</td><td>28.355</td><td>62.400</td></tr><tr><td>5200–5450</td><td>-1.783</td><td>-21.250</td><td>-68.016</td><td>-144.046</td></tr><tr><td>5450–5600</td><td>-1.451</td><td>-20.997</td><td>-81.493</td><td>-280.907</td></tr><tr><td>5600–5800</td><td>-1.922</td><td>-19.083</td><td>-54.677</td><td>-132.374</td></tr><tr><td>5800–6100</td><td>-2.540</td><td>-26.538</td><td>-95.976</td><td>-198.139</td></tr><tr><td>6100–6300</td><td>-3.970</td><td>-48.499</td><td>-208.793</td><td>-379.800</td></tr></tbody></table>
<table><caption>**Table 13**
Spectral Indices and Veiling</caption><thead><tr><th rowspan="2">Region of<br>Low Flux (Å)</th><th rowspan="2">Region of<br>Continuum Flux (Å)</th><th colspan="5">∑a<sub>i</sub> × (F<sub>low</sub>/F<sub>cont</sub>)<sup>i</sup></th></tr><tr><th>a<sub>0</sub></th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>a<sub>3</sub></th><th>a<sub>4</sub></th></tr></thead><tbody><tr><td rowspan="2">5166-5185<sup>a</sup></td><td>5186-5188.5, 5190-5191.5</td><td>331.36</td><td>-1848.60</td><td>3855.37</td><td>-3568.97</td><td>1240.63</td></tr><tr><td>5175-5185.5, 5160.0-5162.5,5164-5165</td><td>-20.04</td><td>106.60</td><td>-208.50</td><td>170.49</td><td>-43.19</td></tr><tr><td rowspan="2">5204-5211.5</td><td>5196.7-5200, 5212-5219.5</td><td>-20.04</td><td>106.60</td><td>-208.50</td><td>170.49</td><td>-43.19</td></tr><tr><td>5220.5-5222, 5240-5242</td><td>-20.04</td><td>106.60</td><td>-208.50</td><td>170.49</td><td>-43.19</td></tr><tr><td rowspan="2">5262-5276.5</td><td>5285.5-5294, 5256.5-5260</td><td>1153.54</td><td>-5924.35</td><td>11409.70</td><td>-9775.75</td><td>3147.71</td></tr><tr><td>5277.5-5279.5,5245-5246.5</td><td>1153.54</td><td>-5924.35</td><td>11409.70</td><td>-9775.75</td><td>3147.71</td></tr><tr><td>5323-5330.5<sup>b</sup></td><td>5312-5316,5335-5338</td><td>-200.01</td><td>837.71</td><td>-1245.03</td><td>739.73</td><td>-126.41</td></tr><tr><td>5369-5373</td><td>5354-5361,5374.5-5375.5,5378-5382</td><td>301.27</td><td>-1700.42</td><td>3572.08</td><td>-3324.61</td><td>1160.96</td></tr><tr><td>5396.5-5411.3<sup>c</sup></td><td>5388-5392.5,5401-5403.5,5416.5-542</td><td>543.57</td><td>-2928.97</td><td>5911.19</td><td>-5305.77</td><td>1790.89</td></tr><tr><td>6160.5-6171</td><td>6144-6146,6171.5-6173,6178-6179.5</td><td>577.36</td><td>-3480.02</td><td>7849.18</td><td>-7864.79</td><td>2959.92</td></tr></tbody></table>

**Notes.**

a Excluding 5169-5170.5

b Excluding 5325-5327

c Excluding 5398-5404, 5407-5409

**ORCID iDs**

Gregory J. Herczeg (沈雷歌) https://orcid.org/0000-0002-7154-6065

Yuguang Chen (陈昱光) https://orcid.org/0000-0003-4520-5395

Jean-Francois Donati https://orcid.org/0000-0001-5541-2887

Andrea K. Dupree https://orcid.org/0000-0002-8985-8489

Frederick M. Walter https://orcid.org/0000-0001-7796-1756

Christopher M. Johns-Krull https://orcid.org/0000-0002-8828-6386

Carlo F. Manara https://orcid.org/0000-0003-3562-262X

Hans Moritz Günther https://orcid.org/0000-0003-4243-2840

Min Fang (房敏) https://orcid.org/0000-0001-8060-1321

P. Christian Schneider https://orcid.org/0000-0002-5094-2245

Jeff A. Valenti https://orcid.org/0000-0003-3305-6281

Silvia H. P. Alencar https://orcid.org/0000-0002-5171-8376

Laura Venuti https://orcid.org/0000-0002-4115-0318

Juan Manuel Alcalá https://orcid.org/0000-0001-8657-095X

Antonio Frasca https://orcid.org/0000-0002-0474-0896

Nicole Arulanantham https://orcid.org/0000-0003-2631-5265

Jeffrey L. Linsky https://orcid.org/0000-0003-4446-3181

Nancy S. Brickhouse https://orcid.org/0000-0002-8704-4473

Nuria Calvet https://orcid.org/0000-0002-3950-5386

Catherine C. Espaillat https://orcid.org/0000-0001-9227-5949

Justyn Campbell-White https://orcid.org/0000-0002-3913-3746

John M. Carpenter https://orcid.org/0000-0003-2251-0602

Seok-Jun Chang https://orcid.org/0000-0002-0112-5900

Kelle L. Cruz https://orcid.org/0000-0002-1821-0650

S. E. Dahm https://orcid.org/0000-0002-2968-2418

Jochen Eislöffel https://orcid.org/0000-0001-6496-0252

Suzan Edwards https://orcid.org/0000-0002-3232-665X

William J. Fischer https://orcid.org/0000-0002-3747-2496

Zhen Guo (郭震) https://orcid.org/0000-0003-0292-4832

Thomas Henning https://orcid.org/0000-0002-1493-300X

Jessy Jose https://orcid.org/0000-0003-4908-4404

Joel H. Kastner https://orcid.org/0000-0002-3138-8250

Ralf Launhardt https://orcid.org/0000-0002-8298-2663

David A. Principe https://orcid.org/0000-0002-7939-377X

Connor E. Robinson https://orcid.org/0000-0003-1639-510X

Javier Serna https://orcid.org/0000-0001-7351-6540

Michal Siwak https://orcid.org/0000-0002-0786-7307

Michael F. Sterzik https://orcid.org/0000-0002-5784-4437

Shinsuke Takasao https://orcid.org/0000-0003-3882-3945

**References**

Alcalá, J., Natta, A., Manara, C., et al. 2014, A&A, 561, A2

Alcalá, J. M., Manara, C. F., Natta, A., et al. 2017, A&A, 600, A20

Alencar, S. H., & Batalha, C. 2002, ApJ, 571, 378

Alencar, S. H. P., Bouvier, J., Walter, F. M., et al. 2012, A&A, 541, A116

Andrews, S. M., Wilner, D. J., Zhu, Z., et al. 2016, ApJL, 820, L40

Ardila, D. R., Herczeg, G. J., Gregory, S. G., et al. 2013, ApJS, 207, 1

Argiroffi, C., Drake, J. J., Bonito, R., et al. 2017, A&A, 607, A14

Arrieta, A., & Torres-Peimbert, S. 2003, ApJS, 147, 97

Arulanantham, N., Gronke, M., Fiorellino, E., et al. 2023, ApJ, 944, 185

Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42

Basri, G., & Batalha, C. 1990, ApJ, 363, 654

Batalha, C., Batalha, N., Alencar, S., Lopes, D., & Duarte, E. 2002, ApJ, 580, 343

Bergin, E. A., Cleeves, L. I., Gorti, U., et al. 2013, Natur, 493, 644

Beristain, G., Edwards, S., & Kwan, J. 1998, ApJ, 499, 828

Beristain, G., Edwards, S., & Kwan, J. 2001, ApJ, 551, 1037

Bernstein, R., Shectman, S. A., Gunnels, S. M., Mochnacki, S., & Athey, A. E. 2003, Proc. SPIE, 4841, 1694

Biazzo, K., Alcalá, J. M., Covino, E., et al. 2012, A&A, 547, A104

Blinova, A., Romanova, M., & Lovelace, R. 2016, MNRAS, 459, 2354

Bonito, R., Orlando, S., Argiroffi, C., et al. 2014, ApJL, 795, L34

Booth, R. A., & Clarke, C. J. 2018, MNRAS, 473, 757

Bouvier, J., Alecian, E., Alencar, S. H. P., et al. 2020, A&A, 643, A99

Bouvier, J., Sousa, A., Pouilly, K., et al. 2023, A&A, 672, A5

Bouvier, J., Cabrit, S., Fernandez, M., Martin, E. L., & Matthews, J. M. 1993, A&AS, 101, 485

Brickhouse, N., Cranmer, S., Dupree, A., et al. 2012, ApJL, 760, L21
Brickhouse, N. S., Cranmer, S. R., Dupree, A. K., Luna, G. J. M., & Wolk, S. 2010, ApJ, 710, 1835

Brown, A. G. A., Vallenari, A., (Gaia Collaboration), et al. 2021, A&A, 649, A1

Calvet, N., D'Alessio, P., Hartmann, L., et al. 2002, ApJ, 568, 1008

Calvet, N., & Gullbring, E. 1998, ApJ, 509, 802

Chang, S.-J., Lee, H.-W., Lee, H.-G., et al. 2018, ApJ, 866, 129

Claes, R. A. B., Manara, C. F., Garcia-Lopez, R., et al. 2022, A&A, 664, L7

Cody, A. M., Stauffer, J., Baglin, A., et al. 2014, AJ, 147, 82

Costigan, G., Vink, J. S., Scholz, A., Ray, T., & Testi, L. 2014, MNRAS, 440, 3444

Cutri, R., Skrutskie, M., Van Dyk, S., et al. 2003, VizieR Online Data Catalog: 2MASS All-Sky Catalog of Point Sources (Cutri+ 2003)

D'Angelo, C. R., & Spruit, H. C. 2010, MNRAS, 406, 1208

de la Reza, R., Torres, C. A. O., Quast, G., Castilho, B. V., & Vieira, G. L. 1989, ApJL, 343, L61

De Vries, W., Becker, R., White, R., & Loomis, C. 2005, AJ, 129, 615

Debes, H. H., Jang-Condell, H., Weinberger, A. J., Roberge, A., & Schneider, G. 2013, ApJ, 771, 45

Dekker, H., D'Odorico, S., Kaufer, A., Delabre, B., & Kotzlowski, H. 2000, Proc. SPIE, 4008, 534

Dodin, A. V., & Lamzin, S. A. 2012, AstL, 38, 649

Donati, J. F., Catala, C., Landstreet, J. D., & Petit, P. 2006, in ASP Conf. Ser. 358, Solar Polarization 4, ed. R. Casini & B. W. Lites (San Francisco, CA: ASP), 362

Donati, J.-F., Gregory, S., Alencar, S., et al. 2011, MNRAS, 417, 472

Donati, J.-F., Hébrard, E., Hussain, G., et al. 2014, MNRAS, 444, 3220

Donati, J. F., Hébrard, E., Hussain, G. A. J., et al. 2015, MNRAS, 453, 3706

Donati, J.-F., Semel, M., Carter, B., Rees, D., & Cameron, A. C. 1997, MNRAS, 291, 658

Drake, J. 2005, in 13th Cambridge Workshop on Cool Stars, Stellar Systems and the Sun, Vol. 560, ed. F. Favata, G. A. J. Hussain, & B. Battrick (Paris: European Space Agency), 519

Dupree, A., Brickhouse, N., Cranmer, S., et al. 2012, ApJ, 750, 73

Dupree, A. K., Brickhouse, N. S., Cranmer, S. R., et al. 2014, ApJ, 789, 27

Dupree, A. K., Brickhouse, N. S., Smith, G. H., & Strader, J. 2005, ApJL, 625, L131

Edwards, S., Fischer, W., Hillenbrand, L., & Kwan, J. 2006, ApJ, 646, 319

Erkal, J., Manara, C. F., Schneider, P. C., et al. 2022, A&A, 666, A188

Espaillat, C. C., Robinson, C. E., Romanova, M. M., et al. 2021, Natur, 597, 41

Espaillat, C. C., Herczeg, G. J., Thanathibodee, T., et al. 2022, AJ, 163, 114

Fang, M., Pascucci, I., Edwards, S., et al. 2018, ApJ, 868, 28

Fang, M., Van Boekel, R., Wang, W., et al. 2009, A&A, 504, 461

Fang, M., Wang, L., Herczeg, G. J., et al. 2023, NatAs, 7, 905

Fedele, D., van den Ancker, M., Henning, T., Jayawardhana, R., & Oliveira, J. 2010, A&A, 510, A72

Findeisen, K., Cody, A. M., & Hillenbrand, L. 2015, ApJ, 798, 89

Fiorellino, E., Zsidi, G., Kóspál, Á., et al. 2022, ApJ, 938, 93

Fischer, W., Edwards, S., Hillenbrand, L., & Kwan, J. 2011, ApJ, 730, 73

Fischer, W. J., Hillenbrand, L. A., Herczeg, G. J., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S.-i. Inutsuka et al. (San Francisco, CA: ASP), 355

Fischer, W., Kwan, J., Edwards, S., & Hillenbrand, L. 2008, ApJ, 687, 1117

France, K., Schindhelm, R., Bergin, E. A., Roueff, E., & Abgrall, H. 2014, ApJ, 784, 127

Frasca, A., Biazzo, K., Alcalá, J. M., et al. 2017, A&A, 602, A33

Gahm, G. F., Walter, F. M., Stempels, H. C., Petrov, P. P., & Herczeg, G. J. 2008, A&A, 482, L35

Gangi, M., Antonucci, S., Biazzo, K., et al. 2022, A&A, 667, A124

Garcia Lopez, R., Natta, A., (Gravity Collaboration), et al. 2020, Natur, 584, 547

Gray, D. F. 2005, The Observation and Analysis of Stellar Photospheres (Cambridge: Cambridge Univ. Press)

Gullbring, E., Hartmann, L., Briceno, C., & Calvet, N. 1998, ApJ, 492, 323

Gully-Santiago, M. A., Herczeg, G. J., Czekala, I., et al. 2017, ApJ, 836, 200

Günther, H. M., Schmitt, J. H. M. M., Robrade, J., & Liefke, C. 2007, A&A, 466, 1111

Guo, Z., Herczeg, G. J., Jose, J., et al. 2018, ApJ, 852, 56

Hamann, F., & Persson, S. E. 1992, ApJ, 394, 628

Hartigan, P., Edwards, S., & Ghandour, L. 1995, ApJ, 452, 736

Hartmann, L., Herczeg, G., & Calvet, N. 2016, ARA&A, 54, 135

Henize, K. G. 1976, ApJS, 30, 491

Herczeg, G. J., Cruz, K. L., & Hillenbrand, L. A. 2009, ApJ, 696, 1589

Herczeg, G. J., & Hillenbrand, L. A. 2008, ApJ, 681, 594

Herczeg, G. J., & Hillenbrand, L. A. 2014, ApJ, 786, 97

Herczeg, G. J., Linsky, J. L., Valenti, J. A., Johns-Krull, C. M., & Wood, B. E. 2002, ApJ, 572, 310

Herczeg, G. J., Wood, B. E., Linsky, J. L., Valenti, J. A., & Johns-Krull, C. M. 2004, ApJ, 607, 369

Hernández, J., Hartmann, L., Calvet, N., et al. 2008, ApJ, 686, 1195

Hillenbrand, L. A., Kiker, T. J., Gee, M., et al. 2022, AJ, 163, 263

Hinton, P. C., France, K., Batista, M. G., et al. 2022, ApJ, 939, 82

Hoff, W., Henning, T., & Pfau, W. 1998, A&A, 336, 242

Huélamo, N., Figueira, P., Bonfils, X., et al. 2008, A&A, 489, L9

Ingleby, L., Calvet, N., Herczeg, G., et al. 2013, ApJ, 767, 112

Ivison, R. J., Bode, M. F., & Meaburn, J. 1994, A&AS, 103, 201

Johns, C. M., & Basri, G. 1995, ApJ, 449, 341

Johns-Krull, C. M. 2007, ApJ, 664, 975

Johns-Krull, C. M., & Basri, G. 1997, ApJ, 474, 433

Johns-Krull, C. M., Chen, W., Valenti, J. A., et al. 2013, ApJ, 765, 11

Johns-Krull, C. M., & Herczeg, G. J. 2007, ApJ, 655, 345

Johnstone, C. P., Jardine, M., Gregory, S. G., Donati, J.-F., & Hussain, G. 2014, MNRAS, 437, 3202

Kama, M., Bruderer, S., van Dishoeck, E. F., et al. 2016, A&A, 592, A83

Kastner, J. H., Huenemoerder, D. P., Schulz, N. S., Canizares, C. R., & Weintraub, D. A. 2002, ApJ, 567, 434

Kastner, J. H., Zuckerman, B., Weintraub, D. A., & Forveille, T. 1997, Sci, 277, 67

Kaufer, A., Stahl, O., Tubbesing, S., et al. 1999, Msngr, 95, 8

Kurosawa, R., Harries, T. J., & Symington, N. H. 2006, MNRAS, 370, 580

Lamzin, S. A. 1998, ARep, 42, 322

Lantz, B., Aldering, G., Antilogus, P., et al. 2004, Proc. SPIE, 5249, 146

Lawson, W. A., & Crause, L. A. 2005, MNRAS, 357, 1399

Lee, H.-W. 2000, ApJL, 541, L25

Lomb, N. R. 1976, Ap&SS, 39, 447

Luhman, K. L. 2023, AJ, 165, 269

Mamajek, E. E. 2009, in AIP Conf. Ser. 1158, Exoplanets and Disks: Their Formation and Diversity, ed. T. Usuda, M. Tamura, & M. Ishii (Melville, NY: AIP), 3

Manara, C. F., Ansdell, M., Rosotti, G. P., et al. 2023, in ASP Conf. Ser. 534, Protostars and Planets VII, ed. S.-i. Inutsuka (San Francisco, CA: ASP), 539

Manara, C. F., Frasca, A., Venuti, L., et al. 2021, A&A, 650, A196

Manara, C. F., Testi, L., Herczeg, G. J., et al. 2017, A&A, 604, A127

Manara, C. F., Testi, L., Natta, A., et al. 2014, A&A, 568, A18

Manara, C. F., Testi, L., Rigliaco, E., et al. 2013, A&A, 551, A107

Mayor, M., Pepe, F., Queloz, D., et al. 2003, Msngr, 114, 20

McClure, M., Calvet, N., Espaillat, C., et al. 2013, ApJ, 769, 73

Mekkaden, M. 1998, A&A, 340, 135

Mentuch, E., Brandeker, A., van Kerkwijk, M. H., Jayawardhana, R., & Hauschildt, P. H. 2008, ApJ, 689, 1127

Micolta, M., Calvet, N., Thanathibodee, T., et al. 2023, ApJ, 953, 177

Miranda, L. F., Torrelles, J. M., & Lillo-Box, J. 2022, A&A, 657, L9

Natta, A., Testi, L., Muzerolle, J., et al. 2004, A&A, 424, 603

Muzerolle, J., Calvet, N., Briceno, C., Hartmann, L., & Hillenbrand, L. 2000, ApJL, 535, L47

Nguyen, D. C., Brandeker, A., van Kerkwijk, M. H., & Jayawardhana, R. 2012, ApJ, 745, 119

Nicholson, B. A., Hussain, G., Donati, J. F., et al. 2021, MNRAS, 504, 2461

Oke, J. B., Cohen, J. G., Carr, M., et al. 1995, PASP, 107, 375

Oke, J. B., & Gunn, J. E. 1982, PASP, 94, 586

Pascucci, I., Gorti, U., & Hollenbach, D. 2012, ApJL, 751, L42

Pascucci, I., Banzatti, A., Gorti, U., et al. 2020, ApJ, 903, 78

Pecaut, M. J., & Mamajek, E. E. 2013, ApJS, 208, 9

Pegues, J., Czekala, I., Andrews, S. M., et al. 2021, ApJ, 908, 42

Pepe, F. A., Cristiani, S., Rebolo Lopez, R., et al. 2010, Proc. SPIE, 7735, 77350F

Pittman, C. V., Espaillat, C. C., Robinson, C. E., et al. 2022, AJ, 164, 201

Qi, C., Wilner, D. J., Aikawa, Y., Blake, G. A., & Hogerheijde, M. R. 2008, ApJ, 681, 1396

Rei, A. C. S., Petrov, P. P., & Gameiro, J. F. 2018, A&A, 610, A40

Robinson, C. E., & Espaillat, C. C. 2019, ApJ, 874, 129

Robinson, C. E., Espaillat, C. C., & Owen, J. E. 2021, ApJ, 908, 16

Robinson, C. E., Espaillat, C. C., & Rodriguez, J. E. 2022, ApJ, 935, 54

Romanova, M. M., Kulkarni, A. K., & Lovelace, R. V. E. 2008, ApJL, 673, L171

Rucinski, S. M., Matthews, J. M., Kuschnig, R., et al. 2008, MNRAS, 391, 1913

Sacco, G. G., Orlando, S., Argiroffi, C., et al. 2010, A&A, 522, A55

Scargle, J. D. 1982, ApJ, 263, 835

Schindhelm, R., France, K., Herczeg, G. J., et al. 2012, ApJL, 756, L23
Scholz, A., Jayawardhana, R., & Brandeker, A. 2005, ApJL, 629, L41

Selvelli, P. L., & Bonifacio, P. 2000, A&A, 364, L1

Sergison, D. J., Naylor, T., Littlefair, S. P., Bell, C. P. M., & Williams, C. D. H. 2020, MNRAS, 491, 5035

Setiawan, J., Henning, T., Launhardt, R., et al. 2008, Natur, 451, 38

Siwak, M., Ogloza, W., Moffat, A. F. J., et al. 2018, MNRAS, 478, 758

Siwak, M., Rucinski, S. M., Matthews, J. M., et al. 2011, MNRAS, 410, 2725

Siwak, M., Rucinski, S. M., Matthews, J. M., et al. 2014, MNRAS, 444, 327

Skopal, A. 2006, A&A, 457, 1003

Sokal, K. R., Deen, C. P., Mace, G. N., et al. 2018, ApJ, 853, 120

Somers, G., Cao, L., & Pinsonneault, M. H. 2020, ApJ, 891, 29

Sousa, A. P., Bouvier, J., Alencar, S. H. P., et al. 2023, A&A, 670, A142

Stelzer, B., & Schmitt, J. H. M. M. 2004, A&A, 418, 687

Stempels, H. C., & Piskunov, N. 2003, A&A, 408, 693

Takasao, S., Tomida, K., Iwasaki, K., & Suzuki, T. K. 2022, ApJ, 941, 73

Teague, R., Bae, J., Huang, J., & Bergin, E. A. 2019, ApJL, 884, L56

Tokovinin, A., Fischer, D. A., Bonati, M., et al. 2013, PASP, 125, 1336

Tull, R. G., MacQueen, P. J., Sneden, C., & Lambert, D. L. 1995, PASP, 107, 251

Vacca, W. D., & Sandell, G. 2011, ApJ, 732, 8

Valenti, J. A., Basri, G., & Johns, C. M. 1993, AJ, 106, 2024

van Boekel, R., Henning, T., Menu, J., et al. 2017, ApJ, 837, 132

van Dokkum, P. G. 2001, PASP, 113, 1420

Van Winckel, H., Duerbeck, H. W., & Schwarz, H. E. 1993, A&AS, 102, 401

Venuti, L., Cody, A. M., Rebull, L. M., et al. 2021, AJ, 162, 101

Venuti, L., Stelzer, B., Alcalá, J. M., et al. 2019, A&A, 632, A46

Vernet, J., Dekker, H., D'Odorico, S., et al. 2011, A&A, 536, A105

Vogt, S. S., Allen, S. L., Bigelow, B. C., et al. 1994, Proc. SPIE, 2198, 362

Wang, L., Bai, X.-N., & Goodman, J. 2019, ApJ, 874, 90

Webb, R. A., Zuckerman, B., Platais, I., et al. 1999, ApJL, 512, L63

Weinberger, A. J., Anglada-Escudé, G., & Boss, A. P. 2013, ApJ, 762, 118

Wilson, T. J. G., Matt, S., Harries, T. J., & Herczeg, G. J. 2022, MNRAS, 514, 2162

Woodgate, B. E., Kimble, R. A., Bowers, C. W., et al. 1998, PASP, 110, 1183

Yang, H., Johns-Krull, C. M., & Valenti, J. A. 2005, ApJ, 635, 466

Yoshida, T. C., Nomura, H., Tsukagoshi, T., Furuya, K., & Ueda, T. 2022, ApJL, 937, L14

Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577

Zhou, Y., Herczeg, G. J., Kraus, A. L., Metchev, S., & Cruz, K. L. 2014, ApJL, 783, L17

Zsidi, G., Manara, C. F., Kóspál, Á., et al. 2022, A&A, 660, A108