# Ly$\alpha$ Scattering Models Trace Accretion and Outflow Kinematics in T Tauri Systems*

Nicole Arulanantham¹, Max Gronke², Eleonora Fiorellino³,⁴, Jorge Filipe Gameiro⁵,⁶, Antonio Frasca⁷, Joel Green¹, Seok-Jun Chang², Rik A. B. Claes⁸, Catherine C. Espaillat⁹, Kevin France¹⁰, Gregory J. Herczeg¹¹,¹², Carlo F. Manara⁸, Laura Venuti¹³, Péter Ábrahám⁴,¹⁴,¹⁵, Richard Alexander¹⁶, Jerome Bouvier¹⁷, Justyn Campbell-White⁸, Jochen Eislöffel¹⁸, William J. Fischer¹, Ágnes Kóspál⁴,¹⁴,¹⁵,¹⁹, and Miguel Vioque²⁰,²¹

¹Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA; narulanantham@stsci.edu

²Max Planck Institut fur Astrophysik, Karl-Schwarzschild-Straße 1, D-85748 Garching bei München, Germany

³INAF-Osservatorio Astronomico di Capodimonte, via Moiariello 16, I-80131 Napoli, Italy

⁴Konkoly Observatory, Research Centre for Astronomy and Earth Sciences, Eötvös Loránd Research Network (ELKH), Konkoly-Thege Miklós út 15-17, 1121 Budapest, Hungary

⁵Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150-762 Porto, Portugal

⁶Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, rua do Campo Alegre 687, 4169-007 Porto. Portugal
⁷INAF—Osservatorio Astrofísico di Catania, via S. Sofia, 78, I-95123 Catania, Italy

⁸European Southern Observatory, Karl-Schwarzschild-Strasse 2, D-85748 Garching bei München, Germany

⁹Institute for Astrophysical Research, Department of Astronomy, Boston University, 725 Commonwealth Avenue, Boston, MA 02215, USA
¹⁰Laboratory for Atmospheric and Space Physics, University of Colorado Boulder, Boulder, CO 80303, USA

¹¹Department of Astronomy, Peking University, Yiheyuan Road 5, Haidian District, Beijing 100871, People's Republic of China

¹²Kavli Institute for Astronomy and Astrophysics, Peking University, Yiheyuan Road 5, Haidian District, Beijing 100871, People's Republic of China
¹³SETI Institute, 339 Bernardo Ave, Suite 200, Mountain View, CA 94043, USA

¹⁴ELTE Eötvös Loránd University, Institute of Physics, Pázmány Péter Sétány 1/A, 1117 Budapest, Hungary

¹⁵CSFK, MTA Centre of Excellence, Konkoly-Thege Miklós út 15-17, 1121 Budapest, Hungary

¹⁶School of Physics and Astronomy, University of Leicester, Leicester, LE1 7RH, UK
¹⁷Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

¹⁸Thüringer Landessternwarte, Sternwarte 5, D-07778 Tautenburg, Germany
¹⁹Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg, Germany

²⁰Joint ALMA Observatory, Alonso de Córdova 3107, Vitacura, Santiago 763-0355, Chile

²¹National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA

Received 2022 September 20; revised 2022 December 22; accepted 2022 December 22; published 2023 February 23

## Abstract

T Tauri stars produce broad Ly$\alpha$ emission lines that contribute ~88% of the total UV flux incident on the inner circumstellar disks. Ly$\alpha$ photons are generated at the accretion shocks and in the protostellar chromospheres and must travel through accretion flows, winds, and jets, the protoplanetary disks, and the interstellar medium before reaching the observer. This trajectory produces asymmetric, double-peaked features that carry kinematic and opacity signatures of the disk environments. To understand the link between the evolution of Ly$\alpha$ emission lines and the disks themselves, we model HST-COS spectra from targets included in Data Release 3 of the Hubble UV Legacy Library of Young Stars as Essential Standards program. We find that resonant scattering in a simple spherical expanding shell is able to reproduce the high-velocity emission line wings, providing estimates of the average velocities within the bulk intervening H I. The model velocities are significantly correlated with the K-band veiling, indicating a turnover from Ly$\alpha$ profiles absorbed by outflowing winds to emission lines suppressed by accretion flows as the hot inner disk is depleted. Just 30% of targets in our sample have profiles with redshifted absorption from accretion flows, many of which have resolved dust gaps. At this stage, Ly$\alpha$ photons may no longer intersect with disk winds along the path to the observer. Our results point to a significant evolution of Ly$\alpha$ irradiation within the gas disks over time, which may lead to chemical differences that are observable with ALMA and JWST.

Unified Astronomy Thesaurus concepts: Classical T Tauri stars (252); Protoplanetary disks (1300); Hubble Space Telescope (761); Stellar accretion (1578); Stellar winds (1636)

Supporting material: figure set

## 1. Introduction

Classical T Tauri stars (CTTSs) are readily identified by the infrared excess associated with their circumstellar disks and the strong excess UV and optical emission generated at their

accretion shocks, which become visible as the natal clouds disperse (see, e.g., Hartmann et al. 2016; Schneider et al. 2020 for reviews). UV photometry from a few thousand of these targets has been measured with the Galaxy Evolution Explorer (GALEX) satellite, showing that both the FUV (1400–1700 Å) and NUV (1800–2750 Å) broadband fluxes decrease with stellar age, as circumstellar material disperses and accretion slows (see, e.g., Findeisen et al. 2011; Gomez de Castro et al. 2015). The disappearance of accretion signatures (and therefore the protoplanetary disks) roughly sets the upper limit on timescales for giant planet formation ($t < 10$ Myr; see, e.g., Hartmann et al. 1998; Lubow & D’Angelo 2006; Fedele et al. 2010;

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

<footnote>* Based on observations collected at the European Southern Observatory under ESO program 106.20Z8.</footnote>

<footnote>Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.</footnote>
Manara et al. 2019). Residual optically thick, gas-rich inner disks may linger, however, as evidenced by large accretion rates derived from Hα emission, *U* band excess (Najita et al. 2007; Espaillat et al. 2012), and UV-near-infrared continuum fits (Manara et al. 2014) for many transitional disk systems with disk dust gaps or cavities that may be associated with protoplanet formation. Such reservoirs have subsequently been resolved around TW Hya (Ercolano & Pascucci 2017), SR 24S and CIDA 1 (Pinilla et al. 2019, 2021), DM Tau (Hashimoto et al. 2021; Francis et al. 2022), and LkCa 15 (Blakely et al. 2022).

Physical-chemical models of circumstellar disks demonstrate that the UV fluxes from accreting CTTSs play a critical role in regulating the abundances and spatial distributions of volatile-bearing molecules (e.g., C₂H, C₂H₂, HCN, and CN; van Zadelhoff et al. 2003; Walsh et al. 2015; Cazzoletti et al. 2018). For example, larger input UV fluxes are expected to produce brighter, more radially extended rings of CN emission than weaker radiation fields (Cazzoletti et al. 2018), and species like H₂O, CO₂, and HCOOH are expected to have enhanced abundances near water snowlines, due to larger desorption rates at these radii (Ruaud & Gorti 2019). This work is consistent with submm observations of C₂H, which can only be reproduced in models with high UV radiation and C/O ratios greater than unity (Bergin et al. 2016; Bergner et al. 2019; Miotello et al. 2019). However, the observed ratios of submillimeter CN and HCN emission unexpectedly do not show any clear correlation with total UV luminosity across targets spanning two orders of magnitude in irradiation (Bergner et al. 2021), despite models predicting increased photodissociation and lower HCN abundances in strong UV fields (Bergin et al. 2003).

One possible explanation for this discrepancy between model predictions and observations is that the UV radiation fields produced by accreting CTTSs are dominated by the Lyα transition of H I (λ = 1215.67 Å), which spans photodissociation and excitation cross sections for a large number of molecules (including HCN and H₂O; see, e.g., van Dishoeck et al. 2006) and contributes ~88% of the total UV flux on average (France et al. 2014). Many of the bright disks observed at high spatial resolution with ALMA do not have accompanying Lyα spectra (see, e.g., Bergner et al. 2021), meaning that disk models cannot fully account for the Lyα properties of individual targets and may not reveal the impact of Lyα photodissociation on abundance ratios like CN/HCN, for example. The feature is challenging to observe, as the line centers where the fluxes peak are always masked by circumstellar and interstellar absorption. However, broad emission line wings are detected at velocities up to ~1000 km s⁻¹ and carry signatures of the radiation fields incident on the disks (Schindhelm et al. 2012b). Characterizing these observed features is critical for interpreting infrared and submm observations of UV-sensitive molecular gas in protoplanetary disks.

The shapes of Lyα emission lines observed from accreting T Tauri stars fall into three general categories: (a) P Cygni-like, where blueshifted absorption from outflowing material has removed photons from wavelengths at λ < 1215.67 Å; (b) inverse P Cygni-like, consistent with redshifted absorption at λ > 1215.67 Å by gas being accreted onto the central stars; and (c) intermediate, symmetric profiles absorbed primarily near line center (Arulanantham et al. 2021). Absorption below the continuum is not observed at these velocities (v > 300 km s⁻¹).

The ratios of flux from the red (λ > 1215.67Å) and blue (λ < 1215.67Å) Lyα emission line wings are strongly correlated with UV-fluorescent emission from electronic transitions of Lyα-pumped H₂ and CO, demonstrating that the radiation is received by molecular gas in surface layers of the inner disks (Schindhelm et al. 2012a; France et al. 2012; Hoadley et al. 2015). The Lyα flux ratios are also significantly correlated with the infrared SED slopes of circumstellar disks with dust rings, gaps, and cavities (Arulanantham et al. 2021), measured between λ = 13μm and λ = 31μm (Furlan et al. 2009). Because the dust continuum emission at these wavelengths is determined by the temperature of the optically thick disk, which maps to the size of the emitting area, the infrared slopes trace the radii where dust grains have been removed from the disks and are not contributing to the total flux (see, e.g., Espaillat et al. 2010). This implies that Lyα irradiation of the gas disks changes as the dust disks evolve, potentially altering the systems' chemical evolution over the timescales of planet formation.

In order to further explore this connection between Lyα radiation and the circumstellar disks, we model HST spectra from 42 targets included in the Hubble UV Legacy Library of Young Stars as Essential Standards Director's Discretionary program (ULLYSES; Roman-Duval et al. 2020). Lyα is the *n* = 2–1 resonant line of H I, which implies that Lyα photons typically scatter thousands to millions of times before reaching the observer because of the high optical depth in the ground state. Each scattering event leads to a slight frequency shift—mostly simply due to the Doppler effect. In other words, Lyα escape through astrophysical media can be thought of as a double diffusion process through space and frequency. The frequency shift depends on the H I properties of the medium, in particular the kinematics and column density (Adams 1972; Bonilha et al. 1979; Neufeld 1990)—but also less intuitive properties such as the clumping factor, which describes the number of neutral gas clouds along the line of sight (Neufeld 1991; Hansen & Oh 2006; Gronke & Oh 2018). We compare the observed spectra to models including resonant scattering for the first time, using a simplified shell geometry to represent the accretion shocks, disk winds, accretion flows, and disks. Optical and infrared data from the PENELLOPE Large Program on the ESO Very Large Telescope (Manara et al. 2021), which provides spectra acquired contemporaneously with the ULLYSES observations on HST, are included in our analysis when available, to assemble a comprehensive multi-wavelength view of interactions between the CTTSs and disks. This work was done as part of the Outflows and Disks around Young Stars: Synergies for the Exploration of Ulyses Spectra (ODYSSEUS; Espaillat et al. 2022) and PENELLOPE collaborations (Manara et al. 2021).

## 2. Targets and Observations

Our sample consists of 42 low-mass targets (M* = 0.14–1.71 M⊙) with either archival HST spectra or new observations from Data Release 3 of the ULLYSES program (Roman-Duval et al. 2020). This group includes objects in Chamaeleon I (N = 10), Lupus (N = 13), Orion OB1 (N = 7), Taurus-Auriga (N = 9), and σ Orionis (N = 1), TW Hya, and V4046 Sgr, covering a total protostellar age range of ~1–12 Myr. All targets were observed with the G130M λ1291 mode on the Cosmic Origins Spectrograph (HST-COS), which provides spectra from λ = 1132–1436 Å at R ~ 12, 000–16, 000. Data reduction for
the ULLYSES targets was done by the ULLYSES core implementation team at STScI (Roman-Duval et al. 2020), which provides publicly available high-level science products (HLSPs) for the entire sample.<fnref n="22" /> All of the ULLYSES data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute, via 10.17909/t9-jzeh-xy14.

The ULLYSES spectra include the Lyα line (Schindhelm et al. 2012a; McJunkin et al. 2014) and fluorescent emission from H₂ and CO in the protoplanetary disks (Herczeg et al. 2002; France et al. 2011). Because HST-COS is effectively a slitless spectrograph (Green et al. 2012), Lyα photons from the sky background readily enter the aperture. This produces a narrow (FWHM ~ 250 km s⁻¹) geocoronal emission line with a peak flux density of ~10⁻¹² erg s⁻¹ cm⁻¹ Å⁻¹ at line center (λ ~ 1216Å).<fnref n="23" /> The Lyα emission from T Tauri stars is heavily absorbed by the ISM and much weaker than the geocoronal emission at line center. However, the line profiles from our targets are sufficiently broadened (FWHM ~ 550–850 km s⁻¹; Schindhelm et al. 2012b), such that Lyα photons from high velocities where geocoronal emission is not produced are readily detectable. These features appear double-peaked, where separations between line wings are much larger than both the instrumental resolution (Δv ~ 18 km s⁻¹) and the extent of the geocoronal emission (FWHM ~ 250 km s⁻¹). For this reason, we do not fit for any radial velocity shifts from line center but rather measure the properties of the redder wing relative to the bluer wing (see Figure 1). The continuum emission is effectively zero, so we do not fit it, for the sake of model simplicity (see also McJunkin et al. 2014).

Seventeen targets in our sample were also observed through the ESO PENELLOPE program (Manara et al. 2021), which uses the ESPRESSO, UVES, and X-Shooter instruments on the Very Large Telescope (VLT) to acquire contemporaneous UV, optical, and near-IR data for all ULLYSES targets. The spectra include Hα emission, which also originates within the CTTS magnetospheres and the base of the disk winds (see, e.g., Lima et al. 2010). Additionally, the near-IR X-Shooter K-band observations can be used to trace hot inner disk material that is readily exposed to Lyα photons (see, e.g., McClure et al. 2013). A detailed description of the data reduction procedures for the PENELLOPE program is available in (Manara et al. 2021).

Table 1 includes relevant published disk properties for all targets in our sample, including WISE photometry (Wright et al. 2010) from data set 10.26131/IRSA142, Herschel 70 µm flux densities, and outer disk inclinations. Mass accretion rates based on Gaia DR3 distances for the systems in Orion OB1 are from Pittman et al. (2022), a range for XX Cha based on its Gaia DR3 distance from Claes et al. (2022), and the remaining values taken from the ULLYSES documentation (based on stellar positions from Gaia DR2). We note that upcoming work will extract updated measurements from all new spectra (R. A. B. Claes et al. 2023, C. V. Pittman et al. 2023, J. Wendeborn et al. 2023, in preparation) All stellar coordinates, spectral types, masses, and extinction corrections adopted by the ULLYSES implementation team are provided at https://

ullyses.stsci.edu and will be modified as the program is completed.

## 3. Measured Properties of Lyα Emission Lines from CTTS

The observed Lyα emission lines were parameterized by measuring six properties from the spectra in our sample: the integrated fluxes at v > 0 (F<sub>Lyα,red</sub>) and v < 0 (F<sub>Lyα,blue</sub>), and the peak velocities and FWHMs of the red (v<sub>red</sub>, FWHM<sub>red</sub>) and blue (v<sub>blue</sub>, FWHM<sub>blue</sub>) emission line wings (see Table 2). These parameters are marked in Figure 1, which shows the observed Lyα emission lines from DF Tau and DM Tau and the rough dust disk morphologies detected for targets with similar Lyα profiles (see Section 5.1 for a more detailed discussion). Figure 2 shows the measured ratios of red to blue fluxes for the full sample, which indicate whether blueshifted outflows or redshifted accretion flows are the dominant absorbers of Lyα photons. The measured flux ratios span three orders of magnitude, with 0.4 < F<sub>Lyα,red</sub>/F<sub>Lyα,blue</sub> < 93. Approximately 70% of targets have F<sub>Lyα,red</sub>/F<sub>Lyα,blue</sub> > 1, corresponding to P Cygni-like profiles with strong absorption from outflowing H I (e.g., DF Tau). The remaining ~30% have more flux on the blue sides of the line profiles, with red emission suppressed by infalling material (inverse P Cygni-like; e.g., DM Tau). A small subset of targets have F<sub>Lyα,red</sub>/F<sub>Lyα,blue</sub> ~ 1, corresponding to roughly symmetric Lyα emission lines (V4046 Sgr, CVSO 109A, Sz 77, IN Cha, CVSO 146, LkCa 15, CHX18N). These systems are among the weakest accretors in the sample (with the exception of the known binary CVSO 109; Pittman et al. 2022) and have lower integrated Lyα fluxes than the targets with larger accretion rates. The physical mechanisms responsible for producing intermediate Lyα flux ratios are discussed further in Section 5.1.

Figure 3 shows the measured velocity separations between the peaks of the red and blue Lyα emission line wings (v<sub>red</sub> - v<sub>blue</sub>), for all targets with distinct double-peaked line profiles (N = 30/42). Of the remaining 12 systems, nine have P Cygni-like profiles with very little Lyα flux above the continuum at v < 0, two have inverse P Cygni-like profiles that are almost entirely suppressed at v > 0, and one has a symmetric profile with peaks that are hidden by the geocoronal emission. Models of Lyα emission in the extragalactic context show that column densities of neutral, intervening H I are correlated with v<sub>red</sub>, v<sub>blue</sub>, and v<sub>red</sub> - v<sub>blue</sub> measured from double-peaked profiles Verhamme et al. (2015, 2017), which appear similar to the wings observed from T Tauri stars. The lower limit of the grid explored in the extragalactic work, below which the H I medium becomes optically thin, corresponds to N<sub>H I</sub> = 10<sup>17</sup> cm<sup>-2</sup> when v<sub>red</sub> - v<sub>blue</sub> < 300 km s<sup>-1</sup>. Our sample ranges from 721 < v<sub>red</sub> to v<sub>blue</sub> < 2070 km s<sup>-1</sup>, indicative of much larger column densities (N(H I) > 10<sup>20</sup> cm<sup>-2</sup>). This lower limit is roughly consistent with the results from McJunkin et al. (2014), who measured H I column densities between log<sub>10</sub>N(H I) ~ 19.6–21.1 by fitting the red wings of the Lyα emission line profiles from a sample of T Tauri and Herbig Ae/Be stars.

## 4. Description of Models

In order to explore the physical mechanisms responsible for the range in Lyα emission line shapes observed in the ULLYSES sample, we fit the spectra with simple models of

<footnote n="22">²² See https://ullyses.stsci.edu/ullyses-data-description.html for a full description of the available ULLYSES data products.</footnote>

<footnote n="23">²³ See https://www.stsci.edu/hst/instrumentation/cos/calibration/airglow for HST-COS calibration spectra of geocoronal Lyα emission.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Characteristic Lyα emission lines observed in HST-COS spectra from ULLYSES targets. We parameterize the observed Lyα profiles using the integrated fluxes at $v > 0$ ($F_{Lyα,red}$) and $v < 0$ ($F_{Lyα,blue}$) and the peak velocities and FWHMs of the red ($v_{red}$, FWHM$_{red}$) and blue ($v_{blue}$, FWHM$_{blue}$) emission line wings. DF Tau and DM Tau have the largest and smallest ratios of $F_{Lyα,red}/F_{Lyα,blue}$, respectively, within our sample. DF Tau shows P Cygni-like emission, where outflowing H I contained in jets and inner disk winds has suppressed the blue side of the profile. This scenario is illustrated in the schematic on the left, which also shows that the targets have smooth, full dust disks at this stage (Arulanantham et al. 2021). DM Tau shows opposite behavior, where infalling H I has absorbed photons on the red side of the line profile. The schematic on the right illustrates this stage, showing that large dust gaps consistent with models of planet–disk interactions have been cleared in the dust distributions (see, e.g., DM Tau; Hashimoto et al. 2021; Francis et al. 2022). Low-velocity components of [O I] 6300 Å and [Ne II] 12.81 µm emission lines are still observed, consistent with origins in MHD winds (see, e.g., Banzatti et al. 2019; Pascucci et al. 2020) or photoevaporative flows (Pascucci & Sterzik 2009). The green boxes in each panel mark velocities between −250 $< v < +250$ km s$^{-1}$, where the Lyα line profiles are dominated by geocoronal emission. We set these fluxes to zero throughout the analysis presented here and calculate velocities relative to line center at 1215.67 Å.</figcaption></figure>

Lyα radiative transfer through H I in protoplanetary disk systems. The photons originate in the protostellar chromospheres and at the accretion shocks, then propagate through the accretion flows, inner disk winds and jets, circumstellar disks, outer disk winds, and ISM (see, e.g., McJunkin et al. 2014). Although this medium includes multiple distinct populations of gas with a superposition of infalling, outflowing, and rotating velocity signatures, the goal of this work is to constrain the bulk properties of intervening gas along the line of sight to these sources. The two modeling frameworks used to reproduce the line profiles are described in the following sections.

## 4.1. Shell Models with Lyα Emission and Absorption

Previous work on the observed Lyα emission lines from T Tauri stars prioritized constraining column densities of H I along the line of sight (McJunkin et al. 2014; Espaillat et al. 2022), which could be accurately measured from the red wings of the line profiles. By fitting the data with models of superimposed broad and narrow Gaussian emission components and a single Voigt absorption feature, McJunkin et al. (2014) determined much lower ISM extinction values than measured from X-ray, optical, and IR observations. No correlations were detected between the best-fit H I column

densities and disk inclinations, implying that the bulk of the absorbing gas is contained in the ISM instead of within the circumstellar disks themselves. Because those model fits were optimized for the red sides of the observed Lyα profiles ($\lambda > 1216$ Å), the data at $\lambda < 1216$ Å carrying signatures of stellar and disk winds were not explored (see also Herczeg et al. 2004).

Other efforts to model the Lyα emission line wings used fluorescent emission from electronically excited CO as an additional constraint on the low-velocity Lyα flux that reaches the protoplanetary disk surface but is hidden from HST by geocoronal emission and interstellar absorption (Schindhelm et al. 2012a). Because Lyα photons from both the blue and red sides of the line profiles are required to excite the electronic transitions of CO, these Lyα models included a second Voigt absorption to represent outflowing H I in disk winds, in addition to the broad and narrow emission components and Voigt absorption from the ISM used in later work (McJunkin et al. 2014). The best-fit Lyα models constrain the total flux reaching the CO molecules, in order to reproduce the fluorescent CO emission lines. However, targets like DF Tau still required some additional absorption component to fully reproduce the observed Lyα profile, which is almost fully suppressed at $v < 0$ (Schindhelm et al. 2012a).
<table><caption>**Table 1**
Sample of Targets</caption><thead><tr><th>Target Name</th><th>log M<sub>acc</sub><sup>a</sup><br>(M<sub>☉</sub> yr<sup>-1</sup>)</th><th>W<sub>3</sub><sup>b</sup><br>(mag)</th><th>W<sub>4</sub><sup>c</sup><br>(mag)</th><th>F<sub>70 µm</sub><sup>d</sup><br>(mJy)</th><th>i<sub>disk</sub><sup>e</sup><br>(degrees)</th><th>ALMA/SMA<sup>e</sup><br>Dust Substructure</th></tr></thead><tbody><tr><td>2MASS J11432669-7804454</td><td>-8.71</td><td>8.676</td><td>7.052</td><td>...</td><td>...</td><td>...</td></tr><tr><td>2MASS J16083070-3828268</td><td>-8.96</td><td>7.842</td><td>3.07</td><td>3400 ± 27</td><td>74</td><td>Cavity at r = 75 au</td></tr><tr><td>AA Tau</td><td>-7.82</td><td>4.645</td><td>2.505</td><td>1300 ± 300</td><td>59.1 ± 0.3</td><td>Rings at r = 49, 95, 143 au</td></tr><tr><td>CHX18N</td><td>-8.09</td><td>...</td><td>3.584</td><td>300 ± 100</td><td>...</td><td>...</td></tr><tr><td>CS Cha</td><td>-8.29</td><td>7.095</td><td>2.758</td><td>3200 ± 600</td><td>8</td><td>Cavity at r = 37 au</td></tr><tr><td>CVSO 104</td><td>-7.94<sup>f</sup></td><td>7.56</td><td>5.571</td><td>77.41 ± 2.35</td><td>43<sup>+8</sup><sub>-4</sub></td><td>...</td></tr><tr><td>CVSO 107</td><td>-7.32<sup>f</sup></td><td>7.419</td><td>5.128</td><td>105.34 ± 3.32</td><td>...</td><td>...</td></tr><tr><td>CVSO 109A</td><td>-6.76<sup>f</sup></td><td>6.523</td><td>4.527</td><td>69.38 ± 2.69</td><td>...</td><td>...</td></tr><tr><td>CVSO 146</td><td>-8.28<sup>f</sup></td><td>7.072</td><td>4.597</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CVSO 165A</td><td>-8.15<sup>f</sup></td><td>6.571</td><td>4.626</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CVSO 165B</td><td>-7.78<sup>f</sup></td><td>6.571</td><td>4.626</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CVSO 176</td><td>-7.38<sup>f</sup></td><td>7.896</td><td>6.221</td><td>...</td><td>...</td><td>...</td></tr><tr><td>CVSO 90</td><td>-7.90<sup>f</sup></td><td>7.622</td><td>5.345</td><td>...</td><td>...</td><td>...</td></tr><tr><td>DE Tau</td><td>-7.55</td><td>4.895</td><td>2.725</td><td>1300 ± 300</td><td>34.1 ± 1.0</td><td>...</td></tr><tr><td>DF Tau</td><td>-6.75</td><td>3.829</td><td>2.27</td><td>700 ± 100</td><td>24 ± 5</td><td>...</td></tr><tr><td>DK Tau</td><td>-7.47</td><td>3.599</td><td>...</td><td>1170 ± 120</td><td>12.8<sup>+2.5</sup><sub>-2.8</sub></td><td>Smooth disk</td></tr><tr><td>DM Tau</td><td>-8.54</td><td>7.085</td><td>3.571</td><td>780 ± 80</td><td>35.2 ± 0.7</td><td>Gap between r = 4 and 20 au</td></tr><tr><td>DN Tau</td><td>-8.0</td><td>5.166</td><td>3.039</td><td>700 ± 100</td><td>35.2<sup>+0.5</sup><sub>-0.6</sub></td><td>Gap (radius not reported)</td></tr><tr><td>HT Lup</td><td>-8.24</td><td>2.896</td><td>0.856</td><td>3500 ± 49</td><td>48.1 ± 4.5</td><td>Spiral arms</td></tr><tr><td>IN Cha</td><td>-9.34</td><td>...</td><td>5.316</td><td>&lt;287.3</td><td>...</td><td>...</td></tr><tr><td>LkCa 15</td><td>-8.51</td><td>5.696</td><td>3.565</td><td>1230 ± 120</td><td>55</td><td>Cavity at r = 76 au</td></tr><tr><td>MY Lup</td><td>-9.67</td><td>5.164</td><td>2.833</td><td>1100 ± 22</td><td>73</td><td>Rings at r = 8, 20, 30, 40 au</td></tr><tr><td>RECX-11</td><td>-9.7</td><td>5.388</td><td>3.718</td><td>212 ± 21</td><td>70</td><td>...</td></tr><tr><td>RECX-15</td><td>-9.0</td><td>...</td><td>3.409</td><td>204 ± 6</td><td>60</td><td>...</td></tr><tr><td>RU Lup</td><td>-7.3</td><td>2.82</td><td>0.658</td><td>...</td><td>18.5 ± 2</td><td>Eight rings between r = 14 and 50 au</td></tr><tr><td>RW Aur</td><td>-7.50</td><td>3.375</td><td>1.448</td><td>2470 ± 150</td><td>55.51 ± 0.13</td><td>Tidal streams</td></tr><tr><td>RY Lup</td><td>-8.19</td><td>3.647</td><td>1.404</td><td>...</td><td>68</td><td>Cavity at r = 69 au</td></tr><tr><td>Sz 10</td><td>-8.7</td><td>7.816</td><td>5.847</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Sz 111</td><td>-9.12</td><td>8.922</td><td>5.731</td><td>1200 ± 24</td><td>...</td><td>...</td></tr><tr><td>Sz 130</td><td>-9.19</td><td>6.53</td><td>4.514</td><td>...</td><td>53</td><td>Smooth disk</td></tr><tr><td>Sz 45</td><td>-8.09</td><td>6.948</td><td>4.156</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Sz 69</td><td>-9.51</td><td>...</td><td>...</td><td>3500 ± 49</td><td>69</td><td>Smooth disk</td></tr><tr><td>Sz 71</td><td>-9.06</td><td>5.716</td><td>3.671</td><td>240 ± 34</td><td>40</td><td>Smooth disk</td></tr><tr><td>Sz 72</td><td>-8.65</td><td>5.987</td><td>3.911</td><td>...</td><td>53</td><td>Smooth disk</td></tr><tr><td>Sz 75</td><td>-7.67</td><td>4.342</td><td>2.379</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Sz 76</td><td>-9.26</td><td>7.366</td><td>5.274</td><td>...</td><td>...</td><td>...</td></tr><tr><td>Sz 77</td><td>-8.79</td><td>5.362</td><td>3.437</td><td>...</td><td>...</td><td>...</td></tr><tr><td>T Cha</td><td>-8.40</td><td>4.663</td><td>2.567</td><td>...</td><td>73</td><td>Gap between r = 18 and 28 au</td></tr><tr><td>TW Hya</td><td>-8.7</td><td>4.54</td><td>1.516</td><td>...</td><td>7</td><td>Ten rings between r = 1 and 52 au</td></tr><tr><td>UX Tau</td><td>-8.0</td><td>5.71</td><td>2.176</td><td>3300 ± 700</td><td>40</td><td>Cavity at r = 31 au</td></tr><tr><td>V4046 Sgr</td><td>-7.89</td><td>4.966</td><td>1.93</td><td>...</td><td>34</td><td>Cavity at r = 31 au</td></tr><tr><td>V510 Ori</td><td>-8.38</td><td>5.598</td><td>3.258</td><td>...</td><td>...</td><td>...</td></tr><tr><td>XX Cha<sup>g</sup></td><td>-9.06--7.69</td><td>...</td><td>3.569</td><td>190 ± 40</td><td>...</td><td>...</td></tr></tbody></table>

**Notes.**

While the superposition of emission and absorption components used in previous work to represent the stellar-disk-wind environment provides good fits to the Lyα profiles with high S/N, the large number of model parameters makes it

difficult to fit similar models to the noisier targets from the ULLYSES sample. To explore the simplest possible framework for all targets in our sample, we first compare the data to models with a single Gaussian emission line and a Voigt

<footnote><sup>a</sup> Mass accretion rates are taken from Alcala et al. (2014), Alcala et al. (2017; Lupus), Manara et al. (2015), Manara et al. (2017; ρ Ophiuchus, Chamaeleon I), France et al. (2017), and Xu et al. (2021).</footnote>

<footnote><sup>b</sup> WISE photometry, 12 µm (Wright et al. 2010).</footnote>

<footnote><sup>c</sup> WISE photometry, 22 µm.</footnote>

<footnote><sup>d</sup> Herschel 70 µm flux densities from Mauco et al. (2016).</footnote>

<footnote><sup>e</sup> Submm disk inclinations and dust substructure (unless otherwise indicated) from Marel et al. (2018; 2MASS J16083070-3828268), Loomis et al. (2017; AA Tau), Francis & van der Marel (2020; CS Cha, LkCa 15, RY Lup, UX Tau, V4046 Sgr), Frasca et al. (2021; CVSO 104; Hα modeling), Simon et al. (2019; DE Tau), Shakhovskoj et al. (2006; DF Tau), Long et al. (2019; DK Tau, DN Tau), Kudo et al. (2018; DM Tau), Kurtovic et al. (2018; HT Lup), Huang et al. (2018; MY Lup, RU Lup, TW Hya), Lawson et al. (2004; RECX-11, RECX-15; Hα modeling), Rodriguez et al. (2018; RW Aur), Ansdell et al. (2016; Sz 130, Sz 69, Sz 71, Sz 72), and Hendler et al. (2018; TCha).</footnote>

<footnote><sup>f</sup> Mass accretion rates from Pittman et al. (2022).</footnote>

<footnote><sup>g</sup> Range of mass accretion rates from Claes et al. (2022).</footnote>
<table><caption>**Table 2**
Measured Properties of Lyα Emission Lines</caption><thead><tr><th>Target Name</th><th>F<sub>Lyα,red</sub><br>(erg s<sup>-1</sup> cm<sup>-2</sup>)</th><th>F<sub>Lyα,blue</sub><br>(erg s<sup>-1</sup> cm<sup>-2</sup>)</th><th>F&#773;<sub>Lyα,red</sub><br>F&#773;<sub>Lyα,blue</sub></th><th>v<sub>red</sub><br>(km s<sup>-1</sup>)</th><th>v<sub>blue</sub><br>(km s<sup>-1</sup>)</th><th>v<sub>red</sub> - v<sub>blue</sub><br>(km s<sup>-1</sup>)</th><th>FWHM<sub>red</sub><br>(km s<sup>-1</sup>)</th><th>FWHM<sub>blue</sub><br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>J11432669-7804454</td><td>9.3 × 10<sup>-14</sup></td><td>2.1 × 10<sup>-14</sup></td><td>4.4</td><td>389</td><td>-330</td><td>719</td><td>198</td><td>232</td></tr><tr><td>J16083070-3828268</td><td>3.9 × 10<sup>-15</sup></td><td>1.94 × 10<sup>-15</sup></td><td>2.0</td><td>542</td><td>-359</td><td>901</td><td>423</td><td>260</td></tr><tr><td>AA Tau</td><td>2.6 × 10<sup>-13</sup></td><td>4.1 × 10<sup>-14</sup></td><td>6.3</td><td>683</td><td>-971</td><td>1654</td><td>250</td><td>1230</td></tr><tr><td>CHX18N</td><td>2.8 × 10<sup>-14</sup></td><td>2.4 × 10<sup>-14</sup></td><td>1.2</td><td>566</td><td>-439</td><td>1005</td><td>515</td><td>567</td></tr><tr><td>CS Cha</td><td>3.7 × 10<sup>-12</sup></td><td>5.3 × 10<sup>-12</sup></td><td>0.7</td><td>618</td><td>-566</td><td>1184</td><td>308</td><td>371</td></tr><tr><td>CVSO 104</td><td>4.3 × 10<sup>-15</sup></td><td>2.1 × 10<sup>-15</sup></td><td>2.0</td><td>602</td><td>-575</td><td>1177</td><td>303</td><td>355</td></tr><tr><td>CVSO 107</td><td>2.0 × 10<sup>-14</sup></td><td>9.6 × 10<sup>-15</sup></td><td>2.1</td><td>541</td><td>-1160</td><td>1700</td><td>258</td><td>1408</td></tr><tr><td>CVSO 109</td><td>6.7 × 10<sup>-15</sup></td><td>6.9 × 10<sup>-15</sup></td><td>0.96</td><td>599</td><td>-719</td><td>1318</td><td>302</td><td>241</td></tr><tr><td>CVSO 146</td><td>4.9 × 10<sup>-15</sup></td><td>4.8 × 10<sup>-15</sup></td><td>1.0</td><td>644</td><td>...</td><td>...</td><td>589</td><td>...</td></tr><tr><td>CVSO 165</td><td>5.9 × 10<sup>-15</sup></td><td>2.2 × 10<sup>-15</sup></td><td>2.7</td><td>686</td><td>...</td><td>...</td><td>323</td><td>...</td></tr><tr><td>CVSO 176</td><td>7.9 × 10<sup>-15</sup></td><td>1.7 × 10<sup>-15</sup></td><td>4.7</td><td>627</td><td>...</td><td>...</td><td>322</td><td>...</td></tr><tr><td>CVSO 90</td><td>1.9 × 10<sup>-14</sup></td><td>2.0 × 10<sup>-15</sup></td><td>9.5</td><td>663</td><td>...</td><td>...</td><td>353</td><td>...</td></tr><tr><td>DE Tau</td><td>1.2 × 10<sup>-13</sup></td><td>7.0 × 10<sup>-14</sup></td><td>1.7</td><td>654</td><td>-1192</td><td>1846</td><td>351</td><td>1450</td></tr><tr><td>DF Tau</td><td>5.1 × 10<sup>-12</sup></td><td>5.5 × 10<sup>-14</sup></td><td>93</td><td>876</td><td>...</td><td>...</td><td>400</td><td>...</td></tr><tr><td>DK Tau</td><td>1.4 × 10<sup>-13</sup></td><td>5.7 × 10<sup>-14</sup></td><td>2.4</td><td>585</td><td>...</td><td>...</td><td>474</td><td>...</td></tr><tr><td>DM Tau</td><td>3.4 × 10<sup>-13</sup></td><td>8.5 × 10<sup>-13</sup></td><td>0.4</td><td>1215</td><td>-849</td><td>2064</td><td>1385</td><td>389</td></tr><tr><td>DN Tau</td><td>1.2 × 10<sup>-13</sup></td><td>1.7 × 10<sup>-13</sup></td><td>0.69</td><td>724</td><td>-1146</td><td>1870</td><td>...</td><td>...</td></tr><tr><td>HT Lup</td><td>3.2 × 10<sup>-15</sup></td><td>1.4 × 10<sup>-15</sup></td><td>2.2</td><td>644</td><td>...</td><td>...</td><td>156</td><td>...</td></tr><tr><td>IN Cha</td><td>2.5 × 10<sup>-15</sup></td><td>2.4 × 10<sup>-15</sup></td><td>1.0</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>LkCa 15</td><td>9.6 × 10<sup>-14</sup></td><td>9.1 × 10<sup>-14</sup></td><td>1.1</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MY Lup</td><td>6.8 × 10<sup>-15</sup></td><td>3.6 × 10<sup>-15</sup></td><td>1.9</td><td>501</td><td>-365</td><td>865</td><td>318</td><td>409</td></tr><tr><td>RECX-11</td><td>2.9 × 10<sup>-12</sup></td><td>2.1 × 10<sup>-12</sup></td><td>1.4</td><td>397</td><td>-381</td><td>777</td><td>178</td><td>275</td></tr><tr><td>RECX-15</td><td>1.9 × 10<sup>-11</sup></td><td>8.6 × 10<sup>-13</sup></td><td>22</td><td>410</td><td>-345</td><td>755</td><td>243</td><td>182</td></tr><tr><td>RU Lup</td><td>2.8 × 10<sup>-13</sup></td><td>2.1 × 10<sup>-14</sup></td><td>14</td><td>428</td><td>-481</td><td>909</td><td>246</td><td>264</td></tr><tr><td>RW Aur</td><td>2.7 × 10<sup>-13</sup></td><td>3.1 × 10<sup>-14</sup></td><td>8.5</td><td>620</td><td>-386</td><td>1006</td><td>407</td><td>231</td></tr><tr><td>RY Lup</td><td>2.7 × 10<sup>-13</sup></td><td>1.5 × 10<sup>-13</sup></td><td>1.8</td><td>463</td><td>-496</td><td>958</td><td>281</td><td>391</td></tr><tr><td>Sz 10</td><td>4.4 × 10<sup>-14</sup></td><td>5.5 × 10<sup>-14</sup></td><td>0.79</td><td>555</td><td>-472</td><td>1027</td><td>282</td><td>195</td></tr><tr><td>Sz 111</td><td>8.4 × 10<sup>-15</sup></td><td>1.3 × 10<sup>-14</sup></td><td>0.63</td><td>670</td><td>-734</td><td>1404</td><td>246</td><td>408</td></tr><tr><td>Sz 130</td><td>4.4 × 10<sup>-15</sup></td><td>1.0 × 10<sup>-14</sup></td><td>0.44</td><td>531</td><td>-593</td><td>1124</td><td>358</td><td>589</td></tr><tr><td>Sz 45</td><td>4.7 × 10<sup>-14</sup></td><td>8.8 × 10<sup>-14</sup></td><td>0.54</td><td>602</td><td>-677</td><td>1278</td><td>80</td><td>346</td></tr><tr><td>Sz 69</td><td>6.3 × 10<sup>-15</sup></td><td>1.1 × 10<sup>-15</sup></td><td>5.6</td><td>720</td><td>-400</td><td>1119</td><td>434</td><td>340</td></tr><tr><td>Sz 71</td><td>1.5 × 10<sup>-14</sup></td><td>5.4 × 10<sup>-15</sup></td><td>2.8</td><td>529</td><td>-594</td><td>1123</td><td>309</td><td>276</td></tr><tr><td>Sz 72</td><td>2.4 × 10<sup>-14</sup></td><td>1.1 × 10<sup>-14</sup></td><td>2.1</td><td>625</td><td>...</td><td>...</td><td>213</td><td>...</td></tr><tr><td>Sz 75</td><td>5.0 × 10<sup>-13</sup></td><td>1.3 × 10<sup>-13</sup></td><td>3.8</td><td>649</td><td>-1113</td><td>1762</td><td>321</td><td>1407</td></tr><tr><td>Sz 76</td><td>4.2 × 10<sup>-15</sup></td><td>5.6 × 10<sup>-15</sup></td><td>0.74</td><td>430</td><td>-764</td><td>1194</td><td>239</td><td>951</td></tr><tr><td>Sz 77</td><td>2.1 × 10<sup>-14</sup></td><td>2.2 × 10<sup>-14</sup></td><td>0.98</td><td>590</td><td>-714</td><td>1304</td><td>386</td><td>615</td></tr><tr><td>T Cha</td><td>2.9 × 10<sup>-14</sup></td><td>1.8 × 10<sup>-14</sup></td><td>1.6</td><td>...</td><td>-286</td><td>...</td><td>...</td><td>164</td></tr><tr><td>TW Hya</td><td>5.2 × 10<sup>-11</sup></td><td>9.2 × 10<sup>-12</sup></td><td>5.6</td><td>388</td><td>-529</td><td>917</td><td>366</td><td>397</td></tr><tr><td>UX Tau</td><td>1.0 × 10<sup>-13</sup></td><td>2.3 × 10<sup>-13</sup></td><td>0.43</td><td>792</td><td>-754</td><td>1546</td><td>344</td><td>464</td></tr><tr><td>V4046 Sgr</td><td>5.2 × 10<sup>-12</sup></td><td>3.9 × 10<sup>-12</sup></td><td>1.3</td><td>421</td><td>-468</td><td>890</td><td>262</td><td>347</td></tr><tr><td>V510 Ori</td><td>1.3 × 10<sup>-14</sup></td><td>4.3 × 10<sup>-15</sup></td><td>3.0</td><td>710</td><td>...</td><td>...</td><td>306</td><td>...</td></tr><tr><td>XX Cha</td><td>1.04 × 10<sup>-13</sup></td><td>2.4 × 10<sup>-14</sup></td><td>4.3</td><td>740</td><td>-489</td><td>1229</td><td>356</td><td>323</td></tr></tbody></table>

Note. Uncertainties in $v_{red}$, $v_{blue}$, FWHM$_{red}$, and FWHM$_{blue}$ are dominated by the instrument spectral resolution ($\Delta v \sim 17$ km s$^{-1}$). Flux uncertainties are determined by the data calibration pipeline and are on the order of ~5%–10%.

absorption component that includes all intervening H I from the ISM, disk wind, and circumstellar disk. This simplest case model already includes five free parameters: the amplitude ($I_0$), FWHM, and central wavelength ($\lambda_0$) of the intrinsic Gaussian emission line, and the column density ($N_{out}$) and velocity ($v_{out}$) of absorbing H I. The full model profile is calculated as

$$I_{Ly\alpha}(\lambda_i) = I_0 \times \exp\left[\frac{-(\lambda_i - \lambda_0)^2}{2 \times (\text{FWHM}/2\sqrt{2*\ln 2})^2}\right] \times \exp[-\tau_{\text{out}}(N_{\text{out}}, v_{\text{out}})], \quad (1)$$

where $\tau_{out} = \sigma_{Ly\alpha}N_{out}$ is the optical depth of the Voigt absorption and $\sigma_{Ly\alpha}$ is the absorption cross section. Because the low-velocity regions of the Lyα spectra ($|v| < 250$ km s$^{-1}$) are hidden by geocoronal emission (and likely fully absorbed

by the ISM), the central portions of the line profiles were masked by setting the fluxes to zero before comparing the data and models.

The best-fit model parameters were determined by minimizing the mean square error (MSE):

$$\text{MSE} = \frac{1}{N} \sum_{i=0}^{N} (y_i - I_{Ly\alpha,i})^2, \quad (2)$$

where $y_i$ and $I_{Ly\alpha,i}$ represent the observed and model fluxes, respectively, at each wavelength along the Lyα profiles. We chose the MSE test statistic because it does not include the flux measurement uncertainties, which are derived from the HST-COS data reduction pipeline and are anomalously small for low-S/N spectra (see, e.g., Arulanantham et al. 2021). This biases
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Measured ratios of integrated Lyα flux from $\lambda \ge 1215.67$ ($F_{Ly\alpha}$ (red)) and $\lambda < 1215.67$ ($F_{Ly\alpha}$ (blue)) for all targets in our sample, which span three orders of magnitude across the parameter space. Roughly 70% of targets have $F_{Ly\alpha}$(red)/$F_{Ly\alpha}$(blue) > 1, corresponding to outflow absorbed, P Cygni-like emission line profiles. The remaining 30% have $F_{Ly\alpha}$(red)/$F_{Ly\alpha}$(blue) < 1, characteristic of infalling H I suppressing the red sides of the features (inverse P Cygni-like), although a small subset have symmetric profiles primarily absorbed near line center.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Measured separations between peaks of blue and red Lyα emission line wings for all targets in our sample with distinct double-peaked line profiles ($N = 30/42$, or ~70%). All targets sit well above the threshold of $v_{red} - v_{blue} \sim 300$ km s⁻¹, which roughly corresponds to $N_{H\ i} = 10^{17}$ cm⁻² (Verhamme et al. 2015, 2017).</figcaption></figure>

test statistics that include the flux errors (e.g., $\chi^2$) toward models that fit the continuum well but not the high-S/N emission line wings. An MCMC sampler (Foreman-Mackey et al. 2013) with a Gaussian log-likelihood function was then used to determine the uncertainties on the model parameters that produced the smallest MSE, by exploring a set of uniform priors for each variable with 300 walkers over 1500 steps and a correction factor for the underpredicted flux uncertainties (see Figure A2). We assessed convergence of the chain by computing the acceptance fraction and reinitialized the sampler around the parameters that maximized the negative log-likelihood function if fewer than ~30% of samples were accepted.

The most robustly constrained prior is on the H I column densities, based on the results from McJunkin et al. (2014). The smallest value measured in that sample was $\log N(H\ i) = 19.2$ for the weak accretor TWA 3A, and the largest was $\log N(H\ i) = 21.05$ for IP Tau. Here, we explored values between $\log N(H\ i) = 18.5$ and 22.0 in order to capture the

full posterior probability distributions. The Gaussian amplitudes were allowed to vary between $10^{-16} < I_0 < 10^{-11}$ erg s⁻¹ cm⁻² Å⁻¹ and the FWHMs between $250 < \text{FWHM} < 2000$ km s⁻¹. Because the low-velocity regions of the observed Lyα profiles are fully absorbed, we restricted the model velocities of absorbing gas to $|v_{\text{out}}| < 250$ km s⁻¹. These model parameters and priors for the MCMC sampler are summarized in Table A1.

## 4.2. Resonant Scattering in a Spherical Shell

In this work, we resorted to the highly simplified “shell model” in order to include resonant scattering in fits to the observed Lyα spectra (Ahn et al. 2001; Verhamme et al. 2015). This framework consists of a Lyα (and adjacent shell) continuum emitting source surrounded by an isotropic shell of neutral hydrogen and dust. The source can be parameterized by the intrinsic Lyα equivalent width $EW_i$ and the intrinsic width of the Lyα line $\sigma_i$. The shell is characterized by its neutral hydrogen column density $N_{H\ i}$, the (absorbing) dust optical depth $\tau_d$, the (effective) temperature $T$ (which includes thermal as well as turbulent motion), and the outflowing (or infalling if $<0$) velocity $v_{\text{exp}}$.

Despite its simplicity, the shell model can fit a range of observed spectra (in the extragalactic context, e.g., Verhamme et al. 2008; Gronke 2017)—including the ones presented here. However, it is clear that a real scattering geometry is neither isotropic nor does it have constant velocity and density profiles as assumed in the shell model. There is, therefore, an ongoing debate in the literature regarding the usefulness of shell model fitting as well as the physical interpretation of the derived parameters (e.g., Orlitova et al. 2018; Li & Gronke 2022). What is clear from this debate is that shell model parameters cannot be interpreted literally, if at all, but instead represent some weighted quantity probed by Lyα photons. In this case, we can interpret the model parameters as representative of the bulk intervening H I along the dominant photon trajectory, contained in accretion flows, jets, and winds, protoplanetary disks, and the ISM.

In order to recover the physical parameters from the Lyα spectra observed with HST-COS, one can compare the data to a suite of synthetic spectra that stem from radiative transfer calculations through simplified geometries. In this work, we used the improved fitting pipeline originally described in Gronke et al. (2015), which employs an affine invariant Monte Carlo algorithm to map out the likelihood landscape (Foreman-Mackey et al. 2013). First, the pipeline searches for the global maximum of the likelihood function. This is nontrivial due to the multimodal and non-Gaussian nature of the likelihood function, but we use a “basinhopping” algorithm in which we randomly sample the discrete parameters ($N_{H\ i}$, $v_{\text{exp}}$, and $T$) and at each point use a COBYLA optimization algorithm to maximize efficiently (using gradients) over the smooth algorithm (Powell 1994). This is repeated until the global maximum is not changed for 300 steps. Afterward, the MC walkers are initialized around that maximum. If any of the walkers find a point with a greater likelihood, we restart the MC sampling and use that new global maximum as starting point.

Prior to feeding the ULLYSES data into the model fitting pipeline, we masked the observed spectra at velocities where geocoronal emission dominates the observed profiles ($|v| < 250$ km s⁻¹) by setting the fluxes to zero and increasing the associated uncertainties by ~5 orders of magnitude. This effectively removes the geocoronal emission from the likelihood calculation without modifying the pipeline
itself. However, this leaves the model fluxes at line center unconstrained, adding larger uncertainties to the model fits for the observed emission line wings. To correct for this, we added an overall normalization A to the integrated $Ly\alpha$ flux as a free parameter. Figure A3 shows an example corner plot from the scattering model for Sz 111, illustrating the marginalized probability distributions of the parameters over the likelihood space sampled by the MCMC routines and showing intrinsic degeneracies between the model parameters. These degeneracies are further discussed in Section 6.3, and the model parameters and priors for the MCMC sampler are summarized in Table A2.

Double-peaked $Ly\alpha$ emission lines similar to those observed from T Tauri stars can also be reproduced by models where photons scatter through a multiphase outflow (Li & Gronke 2022), with a velocity profile that decreases as a function of radial distance from the $Ly\alpha$ source (Gronke & Dijkstra 2016). Such models are highly dependent on the covering factor and column densities of intervening H I clumps, along with the power-law index describing the wind velocity structure. These models consist typically of many more parameters than the shell model, and a full discussion is beyond the scope of this work. However, the models can be roughly divided into two subgroups, with a covering factor below or above a critical value.<fnref n="24" /> The resulting spectra also predict excess $Ly\alpha$ emission at line center, unfortunately, where geocoronal emission dominates any remaining signal from our sources. Because observational constraints on these parameters are difficult to extract with HST-COS, spatially resolved spectra from, e.g., HST-STIS are required to trace the propagation of $Ly\alpha$ photons through the individual components of T Tauri environments using multiphase models instead of the shell framework. Uncontaminated data from $Ly\alpha$ line center will also enable measurements of how much flux scatters into the molecular gas disk, providing a critical input for physical-chemical models (see, e.g., Bergin et al. 2003; Bethell & Bergin 2011).

# 5. Results

## 5.1. $Ly\alpha$ Properties Derived from Shell Models

We fit two different simple shell models to the $Ly\alpha$ emission lines from HST-COS spectra: one with a single Gaussian emission profile and Voigt absorption, and one with resonant scattering effects included. The observed double-peaked emission line wings are well-reproduced with the resonant scattering models for 27/42 targets. Four of the remaining 15 systems had best-fit models without two distinct peaks, despite showing double-peaked profiles in the observed spectra (AA Tau, HT Lup, J16083070-3828268, Sz 77), and the rest had low-S/N emission line wings for which the scattering models did not converge.

Figure 4 compares the best-fit models from the frameworks with and without resonant scattering for a target with P Cygni-like $Ly\alpha$ emission (DF Tau) and one with inverse P Cygni-like line wings (DM Tau). In the case of DF Tau, the best-fit FWHM from the model without scattering ($FWHM = 1200_{-90}^{+140}$ km s⁻¹) is more than double the value from the resonant scattering model ($FWHM = 526_{-2}^{+2}$ km s⁻¹). The best-fit column density in the model without scattering is correspondingly smaller ($N(H\ I) = 20.72_{-0.2}^{+0.01}$

versus $N(H\ I) = 21.39 \pm 0.06$), but the absorber in this model still removes too much flux from the red wing while leaving too much behind in the blue wing. The scattering model, which allows the photons to be redistributed in frequency space in addition to being removed entirely via absorption, is much better able to dilute the flux in the blue wing while preserving the intensity and shape of the red peak that dominates other P Cygni-like line profiles. A similar effect is seen in the best-fit model parameters for DM Tau, although we note that the narrow FWHM predicted by the scattering model ($FWHM = 42_{-0.6}^{+0.7}$ km s⁻¹) is smaller than expected from non-accreting systems. Again, the model including resonant scattering was better able to redistribute photons to both the red and blue high-velocity emission line wings.

We also consider two intermediate cases, V4046 Sgr and CVSO 109A, which have $F_{Ly\alpha}(\text{red})/F_{Ly\alpha}(\text{blue}) = 1.3$ and 0.96, respectively (see Figure 5). Although the two targets have similar observed flux ratios, the integrated line wing fluxes from CVSO 109 are ~3 orders of magnitude smaller than those measured from V4046 Sgr, and we note that CVSO 109A appears to generate less $Ly\alpha$ emission than the average CTTS (Espaillat et al. 2022). The remaining targets with $Ly\alpha$ flux ratios between ~0.8 and 1.4 are all too noisy to fit with the $Ly\alpha$ scattering models (Sz 77, IN Cha, CVSO 146, LkCa 15, CHX18N) but could be reproduced with the models without scattering included. No correlations are detected between the measured $Ly\alpha$ properties (see Table 2) and mass accretion rates, or the model H I velocities, column densities, and FWHMs for this subset of targets. Further observations are required to fully characterize the intermediate sources, three of which are known binary systems: V4046 Sgr (de La Reza et al. 1986), CVSO 109 (Proffitt et al. 2021), and Sz 77 (Zagaria et al. 2022).

Table 3 lists the H I velocities, column densities, and initial $Ly\alpha$ FWHMs from both sets of models for the full sample, which are compared to the ratios of red to blue $Ly\alpha$ flux ($F_{Ly\alpha}(\text{red})/F_{Ly\alpha}(\text{blue})$) in Figure 6. The velocities of absorbing H I are the only parameters that are statistically significantly correlated with the observed emission line profile shapes, and values extracted from models when scattering is and is not included are almost always consistent within the $1\sigma$ uncertainties. One has to bear in mind that degeneracies between the parameters exist, e.g., between $N(H\ I)$, the intrinsic emission line width, and the effective temperature of the medium (Li & Gronke 2022). It is likely then that the best-fit scattering models with the smallest H I column densities ($N(H\ I) = 16.3, 18.2$ for T Cha and CVSO 146, respectively) diverged from the solutions with larger column densities and narrower FWHMs. We report the solutions favored by the MCMC sampler in Table 3 for consistency across the sample, as the geocoronal emission hides any observational constraints on the intrinsic $Ly\alpha$ FWHMs for this sample. Similarly, scattering models returning the smallest FWHMs ($FWHM = 42, 94, 100$ km s⁻¹ for DM Tau, Sz 69, V510 Ori, respectively) could not be distinguished from those with larger FWHMs and smaller H I column densities. We explore alternate constraints on the intrinsic $Ly\alpha$ FWHMs in Section 6.3.

## 5.2. Measured Disk and Accretion Tracers from PENELLOPE Spectra

The near-IR VLT/X-Shooter data from the PENELLOPE program cover the K bandpass, which traces hot dust at small radii close to the YSOs ($T_{\text{dust}} \sim 2000$ K, $r < 0.1$ au; see, e.g.,

<footnote n="24">²⁴ This critical number of clumps per line of sight depends on other parameters and can be analytically computed (Gronke et al. 2017), with the former and latter predicting high or low flux at line center, respectively.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Comparison of simple shell model fits from a target with P Cygni-like, outflow absorbed Lyα emission (DF Tau; left) and a system with an inverse P Cygni-like, infall absorbed profile (DM Tau; right). Observed fluxes between ±250 km s⁻¹ are masked by geocoronal emission, so these velocities were deprioritized in the model fitting routines. The scattering models are better able to reproduce the observed spectra, because they allow the photons to be redistributed in frequency (velocity) space in addition to being removed from the physical system via absorption.</figcaption></figure>

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

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparison of simple shell model fits from targets with intermediate values of $F_{Lyα}$(red)/$F_{Lyα}$(blue) = 1.3 (V4046 Sgr, left) and 0.96 (CVSO 109A, right). The other five sources with Lyα flux ratios between ~0.8 and 1.4 do not have high enough S/N in the emission line wings to fit with the scattering models. Follow-up observations are required to further characterize these intermediate objects, three of which are known binaries (V4046 Sgr, CVSO 109, and Sz 77).</figcaption></figure>

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

McClure et al. 2013). This component of the inner disk produces excess emission that fills in, or “veils,” absorption lines from the stellar photosphere, in addition to the excess emission generated by accretion (see, e.g., Johns-Krull & Valenti 2001; Fiorellino et al. 2021). To explore this hot dust that is likely spatially adjacent to the Lyα-generating accretion shocks, we computed the veiling in the K band ($r_K$) for all 17 targets with available PENELLOPE data and high S/N in the Lyα line wings. This process selects K-band data where telluric contamination was perfectly removed and emission lines were not present, and it compares the spectra to synthetic stellar templates. BT-Settl synthetic spectra from Allard et al. (2012) were used, with a compatible effective temperature to the targets (see Table 4), solar metallicity, and assuming a $logg = 4$, a common value adopted for low-mass YSOs. Each synthetic template was degraded in resolution and rebinned according to the target spectra. Archival values of veiling ($r_K$) are available for an additional eight targets in our sample, for a total of 25 measurements. Table 4 lists the means and standard deviations of the veiling values, along with integrated Brγ fluxes and Hα FWHMs for the PENELLOPE targets. The Hα

line widths were measured using the STAR-MELT Python package for fitting spectral features (Campbell-White et al. 2021).

# 6. Discussion

## 6.1. Measured Lyα Properties Trace Dust Disk Evolution

Previous work reported a strong correlation between the ratios of red to blue Lyα fluxes and infrared spectral indices (Arulanantham et al. 2021), which trace the presence of dust gaps in transition disks via the slopes of the SEDs between 13 and 31 µm (Furlan et al. 2009; Espaillat et al. 2010). This result indicates that the turnover from P Cygni-like to inverse P Cygni-like Lyα emission line profiles is a physical transition that occurs as the dust disks evolve, rather than an effect of the disk inclinations relative to the line of sight. However, the small sample size of N = 12 targets made it difficult to explore the causation behind the trend, prior to the release of data from the ULLYSES program.

Few additional ULLYSES targets have Spitzer spectra from which the infrared spectral index between 13 and 31 µm can be
<table><caption>**Table 3**
Shell Model Properties of Lyα Emission Lines</caption><thead><tr><th rowspan="2">Target Name</th><th colspan="2">N(H I)</th><th colspan="2">FWHM</th><th colspan="2">H I Velocity</th></tr><tr><th>NS</th><th>S</th><th>NS</th><th>S</th><th>NS</th><th>S</th></tr></thead><tbody><tr><td>J11432669-7804454</td><td>20.2<sup>+0.1</sup><sub>-0.3</sub></td><td>20<sup>+0.075</sup><sub>-0.066</sub></td><td>430<sup>+110</sup><sub>-50</sub></td><td>593<sup>+11</sup><sub>-9</sub></td><td>-75<sup>+40</sup><sub>-100</sub></td><td>-28<sup>+4</sup><sub>-2</sub></td></tr><tr><td>J16083070-3828268</td><td>20.0<sup>+0.9</sup><sub>-1.1</sub></td><td>...</td><td>700<sup>+400</sup><sub>-110</sub></td><td>...</td><td>-3<sup>+4</sup><sub>-160</sub></td><td>...</td></tr><tr><td>AA Tau</td><td>20.8 ± 0.1</td><td>...</td><td>610<sup>+170</sup><sub>-10</sub></td><td>...</td><td>-100<sup>+50</sup><sub>-90</sub></td><td>...</td></tr><tr><td>CHX18N</td><td>19.2<sup>+1.5</sup><sub>-0.9</sub></td><td>...</td><td>1200 ± 600</td><td>...</td><td>60<sup>+150</sup><sub>-170</sub></td><td>...</td></tr><tr><td>CS Cha</td><td>20.3<sup>+0.2</sup><sub>-0.2</sub></td><td>20.60 ± 0.07</td><td>760<sup>+230</sup><sub>-60</sub></td><td>710 ± 4</td><td>30<sup>+30</sup><sub>-50</sub></td><td>20 ± 3</td></tr><tr><td>CVSO 104</td><td>20.0<sup>+1.0</sup><sub>-0.2</sub></td><td>20.04<sup>+0.43</sup><sub>-0.76</sub></td><td>1400<sup>+0</sup><sub>-1000</sub></td><td>976<sup>+546</sup><sub>-255</sub></td><td>-210<sup>+190</sup><sub>-30</sub></td><td>-88<sup>+90</sup><sub>-45</sub></td></tr><tr><td>CVSO 107</td><td>20.7 ± 0.2</td><td>20.59<sup>+0.08</sup><sub>-0.07</sub></td><td>520<sup>+90</sup><sub>-80</sub></td><td>536<sup>+15</sup><sub>-18</sub></td><td>-100<sup>+50</sup><sub>-90</sub></td><td>-81 ± 4</td></tr><tr><td>CVSO 109</td><td>20.8<sup>+0.6</sup><sub>-1.5</sub></td><td>21.01 ± 0.08</td><td>950<sup>+720</sup><sub>-460</sub></td><td>566 ± 29</td><td>-80<sup>+140</sup><sub>-10</sub></td><td>-2<sup>+1</sup><sub>-3</sub></td></tr><tr><td>CVSO 146</td><td>20.5<sup>+0.9</sup><sub>-0.6</sub></td><td>18.2<sup>+0.4</sup><sub>-0.8</sub></td><td>1200<sup>+700</sup><sub>-600</sub></td><td>1648<sup>+178</sup><sub>-357</sub></td><td>50<sup>+170</sup><sub>-100</sub></td><td>-332<sup>+80</sup><sub>-40</sub></td></tr><tr><td>CVSO 165</td><td>20.7<sup>+0.3</sup><sub>-0.4</sub></td><td>20.99 ± 0.08</td><td>780<sup>+740</sup><sub>-270</sub></td><td>508<sup>+68</sup><sub>-92</sub></td><td>-160<sup>+120</sup><sub>-10</sub></td><td>-60<sup>+8</sup><sub>-6</sub></td></tr><tr><td>CVSO 176</td><td>20.8<sup>+0.2</sup><sub>-0.1</sub></td><td>21.20 ± 0.08</td><td>540<sup>+110</sup><sub>-70</sub></td><td>350<sup>+69</sup><sub>-83</sub></td><td>-130<sup>+80</sup><sub>-60</sub></td><td>-30 ± 4</td></tr><tr><td>CVSO 90</td><td>20.9<sup>+0.5</sup><sub>-0.3</sub></td><td>20.79 ± 0.07</td><td>610<sup>+120</sup><sub>-60</sub></td><td>639<sup>+56</sup><sub>-76</sub></td><td>-190<sup>+130</sup><sub>-5</sub></td><td>-179<sup>+23</sup><sub>-32</sub></td></tr><tr><td>DE Tau</td><td>20.7<sup>+0.8</sup><sub>-0.9</sub></td><td>...</td><td>1900<sup>+1100</sup><sub>-100</sub></td><td>...</td><td>-40<sup>+20</sup><sub>-90</sub></td><td>...</td></tr><tr><td>DF Tau</td><td>20.72<sup>+0.01</sup><sub>-0.2</sub></td><td>21.39 ± 0.06</td><td>1200<sup>+140</sup><sub>-90</sub></td><td>526 ± 2</td><td>-200<sup>+80</sup><sub>-20</sub></td><td>-74<sup>+3</sup><sub>-3</sub></td></tr><tr><td>DK Tau</td><td>20.1<sup>+0.6</sup><sub>-1.2</sub></td><td>...</td><td>950 ± 400</td><td>...</td><td>-140<sup>+60</sup><sub>-70</sub></td><td>...</td></tr><tr><td>DM Tau</td><td>20.83<sup>+0.2</sup><sub>-0.04</sub></td><td>21.58<sup>+0.08</sup><sub>-0.07</sub></td><td>1100<sup>+80</sup><sub>-140</sub></td><td>42<sup>+0.7</sup><sub>-0.6</sub></td><td>80<sup>+100</sup><sub>-40</sub></td><td>22<sup>+3</sup><sub>-2</sub></td></tr><tr><td>DN Tau</td><td>20.1<sup>+0.8</sup><sub>-0.9</sub></td><td>20.75<sup>+0.05</sup><sub>-0.04</sub></td><td>1900<sup>+100</sup><sub>-900</sub></td><td>1883<sup>+1</sup><sub>-3</sub></td><td>29<sup>+90</sup><sub>-20</sub></td><td>49 ± 1</td></tr><tr><td>HT Lup</td><td>20.80<sup>+0.3</sup><sub>-0.02</sub></td><td>...</td><td>630<sup>+40</sup><sub>-120</sub></td><td>...</td><td>-90<sup>+60</sup><sub>-70</sub></td><td>...</td></tr><tr><td>IN Cha</td><td>20.0 ± 0.7</td><td>...</td><td>240<sup>+100</sup><sub>-1</sub></td><td>...</td><td>10 ± 140</td><td>...</td></tr><tr><td>LkCa 15</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>MY Lup</td><td>20.4<sup>+0.4</sup><sub>-1.5</sub></td><td>20.60<sup>+0.07</sup><sub>-0.09</sub></td><td>520<sup>+150</sup><sub>-10</sub></td><td>465<sup>+16</sup><sub>-33</sub></td><td>-80<sup>+90</sup><sub>-70</sub></td><td>-28<sup>+5</sup><sub>-2</sub></td></tr><tr><td>RECX-11</td><td>20.1<sup>+0.3</sup><sub>-0.2</sub></td><td>19.80 ± 0.07</td><td>400<sup>+170</sup><sub>-60</sub></td><td>627<sup>+16</sup><sub>-19</sub></td><td>-20<sup>+110</sup><sub>-120</sub></td><td>-40 ± 4</td></tr><tr><td>RECX-15</td><td>20.3<sup>+0.2</sup><sub>-0.4</sub></td><td>19.81 ± 0.08</td><td>560<sup>+160</sup><sub>-10</sub></td><td>517<sup>+5</sup><sub>-6</sub></td><td>-190<sup>+110</sup><sub>-1</sub></td><td>-170 ± 4</td></tr><tr><td>RU Lup</td><td>20.3<sup>+0.1</sup><sub>-0.3</sub></td><td>20.20 ± 0.07</td><td>490 ± 180</td><td>544 ± 3</td><td>-110<sup>+40</sup><sub>-70</sub></td><td>-150<sup>+4</sup><sub>-3</sub></td></tr><tr><td>RW Aur</td><td>20.6<sup>+0.1</sup><sub>-0.2</sub></td><td>20.40 ± 0.07</td><td>860<sup>+10</sup><sub>-190</sub></td><td>716<sup>+17</sup><sub>-16</sub></td><td>-220<sup>+140</sup><sub>-20</sub></td><td>-180 ± 4</td></tr><tr><td>RY Lup</td><td>20.2<sup>+0.3</sup><sub>-0.2</sub></td><td>20.20 ± 0.07</td><td>510<sup>+180</sup><sub>-50</sub></td><td>739 ± 11</td><td>-50<sup>+30</sup><sub>-120</sub></td><td>-50<sup>+4</sup><sub>-3</sub></td></tr><tr><td>Sz 10</td><td>19.7 ± 0.9</td><td>...</td><td>1300<sup>+500</sup><sub>-800</sub></td><td>...</td><td>60<sup>+10</sup><sub>-90</sub></td><td>...</td></tr><tr><td>Sz 111</td><td>20.6<sup>+0.4</sup><sub>-0.2</sub></td><td>20.79<sup>+0.08</sup><sub>-0.07</sub></td><td>820<sup>+200</sup><sub>-670</sub></td><td>691<sup>+25</sup><sub>-24</sub></td><td>7<sup>+140</sup><sub>-10</sub></td><td>13 ± 5</td></tr><tr><td>Sz 130</td><td>20.1<sup>+0.9</sup><sub>-0.6</sub></td><td>...</td><td>930<sup>+320</sup><sub>-450</sub></td><td>...</td><td>70<sup>+110</sup><sub>-50</sub></td><td>...</td></tr><tr><td>Sz 45</td><td>20.3<sup>+0.7</sup><sub>-1.6</sub></td><td>...</td><td>1500<sup>+400</sup><sub>-800</sub></td><td>...</td><td>120<sup>+30</sup><sub>-100</sub></td><td>...</td></tr><tr><td>Sz 69</td><td>21.0<sup>+0.05</sup><sub>-0.2</sub></td><td>21.38<sup>+0.09</sup><sub>-0.2</sub></td><td>620 ± 80</td><td>94<sup>+42</sup><sub>-26</sub></td><td>-170<sup>+120</sup><sub>-20</sub></td><td>-22<sup>+6</sup><sub>-4</sub></td></tr><tr><td>Sz 71</td><td>20.0<sup>+1.7</sup><sub>-1.6</sub></td><td>...</td><td>1300<sup>+500</sup><sub>-800</sub></td><td>...</td><td>-110<sup>+80</sup><sub>-70</sub></td><td>...</td></tr><tr><td>Sz 72</td><td>19.7<sup>+1.3</sup><sub>-0.2</sub></td><td>...</td><td>1800<sup>+5</sup><sub>-300</sub></td><td>...</td><td>-200<sup>+160</sup><sub>-40</sub></td><td>...</td></tr><tr><td>Sz 75</td><td>20.5<sup>+0.1</sup><sub>-0.2</sub></td><td>20.9<sup>+0.4</sup><sub>-0.1</sub></td><td>800<sup>+150</sup><sub>-30</sub></td><td>638 ± 3</td><td>-120<sup>+40</sup><sub>-70</sub></td><td>-98 ± 2</td></tr><tr><td>Sz 76</td><td>20.0<sup>+1.2</sup><sub>-1.3</sub></td><td>...</td><td>1000<sup>+400</sup><sub>-500</sub></td><td>...</td><td>-30<sup>+190</sup><sub>-40</sub></td><td>...</td></tr><tr><td>Sz 77</td><td>20.4<sup>+0.4</sup><sub>-1.9</sub></td><td>...</td><td>720<sup>+30</sup><sub>-190</sub></td><td>...</td><td>-20 ± 50</td><td>...</td></tr><tr><td>T Cha</td><td>20.1<sup>+0.6</sup><sub>-1.4</sub></td><td>16.3<sup>+1</sup><sub>-0.4</sub></td><td>200<sup>+150</sup><sub>-60</sub></td><td>653<sup>+34</sup><sub>-39</sub></td><td>-20<sup>+40</sup><sub>-160</sub></td><td>-101<sup>+30</sup><sub>-28</sub></td></tr><tr><td>TW Hya</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>UX Tau</td><td>20.7<sup>+0.3</sup><sub>-0.6</sub></td><td>19.7<sup>+1.2</sup><sub>-0.5</sub></td><td>830<sup>+140</sup><sub>-90</sub></td><td>108<sup>+10</sup><sub>-19</sub></td><td>50<sup>+20</sup><sub>-130</sub></td><td>25<sup>+12</sup><sub>-6</sub></td></tr><tr><td>V4046 Sgr</td><td>19.65<sup>+0.007</sup><sub>-0.04</sub></td><td>20.21<sup>+0.07</sup><sub>-0.1</sub></td><td>820<sup>+20</sup><sub>-10</sub></td><td>724 ± 2</td><td>-40 ± 10</td><td>-23<sup>+2</sup><sub>-3</sub></td></tr><tr><td>V510 Ori</td><td>20.7 ± 0.4</td><td>20.81 ± 0.08</td><td>870<sup>+610</sup><sub>-300</sub></td><td>100<sup>+24</sup><sub>-17</sub></td><td>-180<sup>+130</sup><sub>-10</sub></td><td>-160<sup>+4</sup><sub>-5</sub></td></tr><tr><td>XX Cha</td><td>20.9<sup>+0.02</sup><sub>-0.3</sub></td><td>20.93<sup>+0.04</sup><sub>-0.1</sub></td><td>760<sup>+150</sup><sub>-50</sub></td><td>724<sup>+1</sup><sub>-2</sub></td><td>-150<sup>+70</sup><sub>-50</sub></td><td>-43<sup>+1</sup><sub>-2</sub></td></tr></tbody></table>

Note. Best-fit Lyα shell parameters when considering a model without resonant scattering (NS) and with resonant scattering (S).

calculated, but 36/42 of the circumstellar disks from our sample were detected at infrared wavelengths with WISE (see, e.g., Koenig & Leisawitz 2014). The remaining six were observed, but the photometry was flagged for contamination and is excluded from the following analysis. Figure 7 compares the ratios of red to blue Lyα flux to the [W3 – W4] colors, which roughly capture the SED slopes between 12 and 22 µm. Large [W3 – W4] can also indicate contributions from circumstellar envelopes (Koenig & Leisawitz 2014), but the range of [W1 – W2] colors from our sample (0.04 < [W1 – W2] < 1.0) place all targets in the Class II/ transition disk regime defined in that work.

We detect a statistically significant negative correlation between [W3 – W4] and the ratios of red to blue Lyα flux ($\rho = -0.3, p = 0.04$), such that larger positive values of [W3 – W4] correspond to inverse P Cygni-like Lyα emission lines. However, we note that the W3 filter is broad enough to capture both optically thin 10 µm silicate emission and the optically thick dust continuum. This trend then implies that targets with more blue than red Lyα flux have less flared disks and/or gaps in the dust distributions that place them in the transition disk category, as illustrated in Figure 1.

We also compare the ratios of red to blue Lyα flux to 70 µm flux densities from the 23 targets that were also observed with
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Ratios of red to blue integrated Ly$\alpha$ fluxes vs. model H I velocities (left), column densities (middle), and emission line FWHMs (right). The velocities of absorbing H I are statistically significantly correlated with the observed emission line profile shapes and generally consistent between models when scattering is and is not included. Outlier values in the upper right and middle panels (N(H I) = 16.3, 18.2 for T Cha and CVSO 146 and FWHM = 42, 94, 100 km s$^{-1}$ for DM Tau, Sz 69, and V510 Ori) likely correspond to models that diverged from the more physical solutions with larger column densities and narrower FWHMs. Spearman $\rho$ and $p$ values are shown between the Ly$\alpha$ flux ratios and parameters from models with scattering included.</figcaption></figure>

<table><caption>**Table 4**
Measured Properties from PENELLOPE Spectra</caption><thead><tr><th>Target Name</th><th>WTTS Template/SpT</th><th>r<sub>K</sub></th><th>Brγ Flux<br>(10<sup>-15</sup> erg s<sup>-1</sup> cm<sup>-2</sup> nm<sup>-1</sup>)</th><th>Hα FWHM<br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>AA Tau</td><td>...</td><td>0.8 ± 0.4<sup>a</sup></td><td>...</td><td>...</td></tr><tr><td>CHX18N</td><td>RX J1540.7-3756/K6</td><td>0.66 ± 0.2</td><td>&lt;11</td><td>220</td></tr><tr><td>CVSO 104</td><td>...</td><td>...</td><td>4 ± 1</td><td>235</td></tr><tr><td>CVSO 107</td><td>RX J1540.7-3756/K6</td><td>0.48 ± 0.25</td><td>5.8 ± 0.6</td><td>200</td></tr><tr><td>CVSO 109A</td><td>...</td><td>0.65 ± 0.2</td><td>13.6 ± 0.8</td><td>168</td></tr><tr><td>CVSO 146</td><td>RX J1540.7-3756/K6</td><td>0.7 ± 0.3</td><td>10.9 ± 0.6</td><td>206</td></tr><tr><td>CVSO 165</td><td>RX J1540.7-3756/K6</td><td>0.7 ± 0.2</td><td>8.9 ± 0.6</td><td>244</td></tr><tr><td>CVSO 176</td><td>TWA 7/M2</td><td>1.3 ± 0.3</td><td>&lt;2.8</td><td>211</td></tr><tr><td>CVSO 90</td><td>SO 879/K7</td><td>2.2 ± 0.4</td><td>28.1 ± 0.3</td><td>244</td></tr><tr><td>DF Tau</td><td>...</td><td>0.9 ± 0.3<sup>a</sup></td><td>...</td><td>...</td></tr><tr><td>DK Tau</td><td>...</td><td>1.8 ± 0.6<sup>a</sup></td><td>...</td><td>...</td></tr><tr><td>DM Tau</td><td>...</td><td>0<sup>a</sup></td><td>...</td><td>...</td></tr><tr><td>DN Tau</td><td>...</td><td>0.1 ± 0.1<sup>a</sup></td><td>...</td><td>...</td></tr><tr><td>IN Cha</td><td>Par-Lup3-2/M6</td><td>0.1 ± 0.1</td><td>&lt;2.8</td><td>154</td></tr><tr><td>RW Aur</td><td>...</td><td>&gt;1.5<sup>a</sup></td><td>...</td><td>...</td></tr><tr><td>Sz 10</td><td>Par-Lup3-2/M6</td><td>0.6 ± 0.4</td><td>5.1 ± 0.4</td><td>148</td></tr><tr><td>Sz 45</td><td>TWA 25/M0.5</td><td>0.25 ± 0.2</td><td>5 ± 1</td><td>150</td></tr><tr><td>Sz 69</td><td>SO 797/M4.5</td><td>1.7 ± 0.4</td><td>10.8 ± 0.7</td><td>173</td></tr><tr><td>Sz 71</td><td>TWA 13B/M1</td><td>0.36 ± 0.09</td><td>10 ± 2</td><td>135</td></tr><tr><td>Sz 72</td><td>TWA 13B/M1</td><td>0.8 ± 0.2</td><td>81 ± 1</td><td>295</td></tr><tr><td>Sz 75</td><td>RX J1540.7-3756/K6</td><td>2.4 ± 0.5</td><td>147 ± 13</td><td>213</td></tr><tr><td>Sz 77</td><td>RX J1540.7-3756/K6</td><td>0.7 ± 0.1</td><td>&lt;7.3</td><td>237</td></tr><tr><td>TW Hya</td><td>...</td><td>0.1 ± 0.1<sup>a</sup></td><td>...</td><td>...</td></tr><tr><td>UX Tau</td><td>...</td><td>0.4<sup>a</sup></td><td>...</td><td>...</td></tr><tr><td>XX Cha</td><td>TWA 9B/M3</td><td>0.3 ± 0.3</td><td>5 ± 1</td><td>136</td></tr></tbody></table>

**Note.**

Herschel (Table 1; see, e.g., Mauco et al. 2016). Disk models fit to the 70 µm emission return best-fit temperatures at $r = 10$ au ranging between $T = 24$ and 236 K, which are strongly correlated with the observed flux densities (Ribas et al. 2017). This relationship makes the Herschel data important probes of large dust grains that continue to produce optically thick continuum emission even after near-IR excess from smaller grains disappears (see, e.g., Rebollido et al. 2015). No significant correlation is detected between the Herschel flux densities and Ly$\alpha$ flux ratios, perhaps indicating that the 70 µm emission remains relatively constant over the timescales when accretion broadened Ly$\alpha$ emission line wings are readily detectable. Although the observed Ly$\alpha$ flux ratios do not

provide constraints on outer disk dust evolution, we conclude that the shapes of the emission line profiles are linked to the inner disk flaring angles and dust gaps within $r < 10$ au.

Figure 2 shows that ~30% of targets have accretion absorbed, inverse P Cygni-like Ly$\alpha$ emission line wings— and therefore larger [W3 – W4] colors in Figure 7. This percentage is slightly higher than the fraction of transition disks with large submm dust cavities (26% with $r > 15$ au; see, e.g., Andrews et al. 2011) and much larger than the fraction predicted from infrared SED analyses (10%; see, e.g., Espaillat et al. 2014), perhaps indicating that some targets in our sample have dust gaps opened between the inner and outer disks at

<footnote><sup>a</sup> Values of r<sub>K</sub> taken from the literature: Folha & Emerson (1999; AA Tau, DF Tau, DK Tau, DN Tau, RW Aur, TW Hya); Espaillat et al. (2010; DM Tau, UX Tau).</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Ratios of red to blue integrated Lyα fluxes vs. WISE colors between 12 (W3) and 22 (W4) µm (left) and mass accretion rates (right). No correlations are detected between the flux ratios and accretion rates, implying that targets with P Cygni-like Lyα profiles are not necessarily undergoing more rapid accretion. Larger values of [W3 − W4] correspond to less flux from 10 µm silicate emission and optically thick dust continuum, both of which are captured in the W3 filter. These colors may trace the disk flaring in smooth disks and the locations of dust gaps in structured disks.</figcaption></figure>

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

radii smaller than 15 au. Although not all ULLYSES targets have been observed at submm wavelengths with sufficient spatial resolution to detect dust substructure, four disks with inverse P Cygni-like Lyα profiles have dust cavities imaged at submm wavelengths: DM Tau (r = 19 au; Andrews et al. 2011), UX Tau (r = 30 au; Andrews et al. 2011), Sz 111 (r = 80 au; Ansdell et al. 2016), and CS Cha (r = 37 au; Francis & van der Marel 2020). Those same sources have a higher detection frequency for the 1600 Å H₂O dissociation bump (France et al. 2017), implying that Lyα photons are more readily propagating through and photodissociating the molecular gas in these disks.

Signatures of fast inner disk winds and jets are no longer detected at the stage of disk evolution when inverse P Cygni-like Lyα emission is observed, although low-velocity emission from [O I] and [Ne II] in slow photoevaporative winds is still detected (see, e.g., Banzatti et al. 2019; Pascucci et al. 2020). Lyα photons must then only travel through the accretion flows, protoplanetary disks, and ISM before reaching the observer. The blue sides of the emission line profiles are no longer absorbed by fast outflowing H I, and redshifted absorption within accretion flows becomes the dominant factor in shaping the observed Lyα features. However, we note that not all targets with dust substructure in our sample have inverse P Cygni-like Lyα profiles (see Table 1). Furthermore, models that produce dust substructure via photoevaporative winds or planet-disk interactions do not make predictions for the impact of these physical mechanisms on Lyα emission (see, e.g., Gárate et al. 2021), making it difficult to distinguish between scenarios in this data set. Submm observations at high spatial resolution for the remaining ULLYSES targets are required in order to further explore this relationship between Lyα flux ratios and dust disk evolution.

Observations of the C II λ1335 doublet in HST-COS spectra of T Tauri stars also show absorption signatures from fast, collimated jets and slow winds originating from the inner disks (Xu et al. 2021). Figure 8 compares the wind velocities derived from the C II features (Xu et al. 2021) to the model H I velocities for targets with P Cygni-like, outflow-dominated Lyα profiles. Of the seven targets included in both samples, only three have C II and H I velocities that are consistent within

the 1σ uncertainties. All four of the remaining systems have smaller H I velocities than the C II measurements, corresponding to slower winds. The lack of correlation between H I and C II wind velocities (ρ = −0.2, p = 0.6) is also consistent with observations of Lyα emission lines from compact, star-forming green pea galaxies, which are always well-reproduced by models with lower H I velocities than indicated by absorption against other low-ionization features (Orlitova et al. 2018).

The C II absorption velocities that Xu et al. (2021) associate with fast winds are inversely correlated with the disk inclinations, as expected for material originating in collimated jets with speeds that increase with distance from the source. Photons traveling from edge-on sources only pass through the slow bases of the jets but must propagate through faster regions to escape more face-on disks. Meanwhile, absorption from slow winds is only detected in disks with intermediate to high inclinations, but no significant correlation is detected between the wind velocities and disk inclinations. To explore whether the Lyα-absorbing H I is more consistent with collimated jets or slow disk winds, Figure 8 compares the disk inclinations to the H I velocities for targets with P Cygni-like, outflow-dominated Lyα emission, which are not statistically significantly correlated (ρ = 0.43, p = 0.07). The targets in our sample that have been observed at submm wavelengths have outer disk inclinations ranging from $i_d = 7^\circ$ to $i_d = 74^\circ$ (see Table 1). However, we note that submm observations are not available for most targets with H I velocities much larger or smaller than the median value of the sample $v \sim -96$ km s$^{-1}$. Furthermore, misalignments between the inner and outer disk inclination angles are frequently detected (see, e.g., Davies 2019; Bohn et al. 2022). The high model H I velocities reported here are consistent with absorbing gas located very close to the stars (see, e.g., Pascucci et al. 2020) and may be more correlated with the inner disk inclinations that are much harder to constrain.

Of the five systems with submm observations and inverse P Cygni-like, infall-dominated Lyα profiles, four have disk inclinations between 35° and 37° that may be in alignment with the launching angle of MHD disk winds (Banzatti et al. 2019). This inclination angle places them at the edge of the population
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Absorption velocities of outflowing H I from best-fit Lyα models vs. C II wind velocities for the targets included in Xu et al. (2021; left) and circumstellar disk inclinations measured from submm observations (right). The median H I velocity (−96 km s⁻¹) is marked as a black, dashed line in the right panel. No statistically significant correlation is detected between the velocities and disk inclinations, although we note that submm observations are not available for most targets with H I velocities much larger or smaller than the median value. Uncertainties on the C II slow wind velocities are <100 km s⁻¹ (Xu et al. 2021).</figcaption></figure>

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

with slow winds absorbing the C II measured by Xu et al. (2021). Targets with disk inclinations ~35° also show the largest blueshifts in both the broad and narrow low-velocity components of [O I] 6300 Å emission line profiles, as expected for coinciding disk inclinations and wind launching angles (Banzatti et al. 2019). It is possible then that Lyα photons from disks at intermediate inclinations do not scatter through the winds before escaping the system, leaving only the signatures of the accretion flows imprinted on the observed emission line profiles. Naturally, the simple, isotropic shell model does not capture these subtleties. In the future, we can model individual systems with more complex frameworks to study this point further. We note that three of the four targets with intermediate disk inclinations and P Cygni-like Lyα emission are known binaries (V4046 Sgr, HT Lup, and Sz 69), which likely complicates the scattering geometry.

## 6.2. Model H I Velocities Probe Kinematics of Accretion and Outflows

Lyα photons are generated at the YSO accretion shocks (see, e.g., Hartmann et al. 2016; Schneider et al. 2020), and previous work reports significant correlations between reconstructed Lyα emission line fluxes at the disk surfaces and mass accretion rates measured from the C IV $\lambda_0 = 1548, 1550$Å resonance doublet (France et al. 2014). We compare the Lyα flux ratios to measured accretion rates for our sample, to explore whether the turnover in flux ratios occurs as accretion slows (see Figure 7). However, no clear relationships are detected, demonstrating that targets with P Cygni-like Lyα profiles are not necessarily undergoing more rapid accretion. We note that there is no statistically significant correlation between mass accretion rates and stellar masses for this sample ($\rho = 0.18, p = 0.3$). The measured accretion rates for the 13 disks with inverse P Cygni-like Lyα emission range from $5 \times 10^{-10} < \dot{M} < 3 \times 10^{-8} M_{\odot} \text{ yr}^{-1}$ (see Table 1), with one outlier and known binary at the high end ($\dot{M} = 1.72 \times 10^{-8} M_{\odot} \text{ yr}^{-1}$ for CVSO 109A; Pittman et al. 2022), confirming that accretion rates can remain high even as material is dissipated from the outer disks. This is consistent with submm observations, which still show significant gas abundances in disks with high accretion rates and large dust

cavities or gaps (see, e.g., Bruderer 2013; Espaillat et al. 2014; Kudo et al. 2018; Francis et al. 2022).

Previous modeling work found that the velocities of expanding (outflowing) or collapsing (accreting) shells are well-constrained for emission lines with distinct, asymmetric double peaks (Li & Gronke 2022). Figure 6 compares the observed Lyα flux ratios to the H I velocities derived from the shell models in this work that do and do not account for resonant scattering. The H I velocities from the two modeling frameworks are nearly always consistent with each other within the $1\sigma$ uncertainties, and they are strongly correlated with the ratios of red to blue Lyα flux (Spearman $\rho = -0.95, p = 3 \times 10^{-13}$ without scattering; Spearman $\rho = -0.59, p = 2 \times 10^{-3}$ when scattering is included). Because almost all targets in our sample have double-peaked Lyα emission lines, we are able to interpret the model values as average outflow (or infall) velocities.

We compare the model H I velocities and Lyα flux ratios to veiling measured from K-band spectra acquired as part of the PENELLOPE program ($r_K$; see Figure 9), which quantifies how much the NIR spectral features are filled in, or veiled, by excess emission from hot gas that is accreting onto the central protostars and inner disk dust (see, e.g., McClure et al. 2013; Fiorellino et al. 2021). The Lyα flux ratios and veiling measurements are statistically significantly positively correlated (Spearman $\rho = 0.60, p = 0.006$), showing that increased veiling generally corresponds to P Cygni-like Lyα emission lines. We also find a strong negative correlation between the H I velocities and the K-band veiling measurements (Spearman $\rho = -0.66, p = 0.0006$), consistent with the disappearance of $T \sim 2000$ K dust that fills in the intrinsic protostellar spectrum (see, e.g., Kidder et al. 2021). This may indicate a depletion of inner disk material as outflowing winds slow by shifting to larger radii and infalling accretion flows become the dominant absorber of Lyα photons, although $r_K$ may also be influenced by the inner disk inclination or scale height that we cannot measure from these observations (see, e.g., Johns-Krull & Valenti 2001). As an additional tracer of accretion from spectra acquired simultaneously—and therefore not impacted by variability—we also compare the K-band veiling to Brγ emission measured from the PENELLOPE spectra in
<figure><img src="image_14.png" /><figcaption>**Figure 9.** K-band veiling ($r_K$) measured from X-Shooter spectra acquired through the PENELLOPE program vs. (a) measured ratios of red to blue Ly$\alpha$ flux, (b) model absorption velocities of outflowing and infalling H I, and (c) Br$\gamma$ fluxes measured from PENELLOPE spectra. Statistically significant correlations are detected between $r_K$ and $F_{Ly\alpha}(\text{red})/F_{Ly\alpha}(\text{blue})$ ($\rho = 0.60, p = 0.0006$), the H I velocities ($\rho = -0.66, p = 0.0006$), and Br$\gamma$ emission ($\rho = 0.62, p = 0.009$), indicating that the turnover in shape of Ly$\alpha$ profiles occurs for targets with spectra that are less veiled by hot dust ($T_{\text{dust}} \sim 2000\text{K}$).</figcaption></figure>

Figure 9. A statistically significant positive correlation is detected (Spearman $\rho = 0.623, p = 0.009$), although the Br$\gamma$ emission is not significantly correlated with the Ly$\alpha$ flux ratios (Spearman $\rho = 0.42, p = 0.16$). These trends will likely become more clear once simultaneous accretion rates have been measured for all ULLYSES targets (C. V. Pittman et al., J. Wendeborn et al., R. A. B. Claes et al. 2023, in preparation). Despite the statistical significance of the correlations between Ly$\alpha$ flux ratios, model H I velocities, and $r_K$, we note that the small *p* values reported here are dependent on the choice of stellar template used to measure the veiling, and they should not be interpreted as absolute.

### 6.3. Constraints on Turbulence within Accretion Shocks

Figure 6 demonstrates that no strong correlations are detected between the Ly$\alpha$ flux ratios and H I column densities (Spearman $\rho = -0.1, p = 0.6$ when scattering is not included; Spearman $\rho = 0.08, p = 0.7$ with scattering) or FWHMs of the intrinsic emission lines (Spearman $\rho = -0.2, p = 0.5$ without scattering; Spearman $\rho = -0.3, p = 0.1$ when scattering is included). Li & Gronke (2022) find that $N(H\,I)$, emission line widths, and effective temperatures within the scattering medium can be degenerate parameters, which may be responsible for the lack of correlations with the observed flux ratios (see corner plot for Sz 111 in Appendix A). The intrinsic line widths are challenging to constrain, because of both the geocoronal emission masking velocities <250 km s$^{-1}$ and turbulence within the accretion shocks themselves. The Ly$\alpha$ photons must propagate through this highly turbulent region before reaching the protoplanetary disks, winds/jets, and ISM, which requires the intrinsic model emission line to be broader than the velocity dispersions of randomly moving clumps within the scattering medium (Li & Gronke 2022). When the clumps are accounted for in multiphase models, the total H I column densities are expected to decrease, as the wider emission lines allow flux to be distributed to higher velocities with fewer scattering events.

Because H$\alpha$ scattering requires a large population of atomic hydrogen in the *n* = 2 level, the process occurs at a much lower rate than Ly$\alpha$ scattering and leads to emission lines with less turbulent broadening than the intrinsic Ly$\alpha$ profiles. To explore the impact of turbulence on the Ly$\alpha$ observations, Figure 10 compares the Ly$\alpha$ FWHMs from the models presented in this work to H$\alpha$ emission line widths measured from the PENELLOPE spectra. The Ly$\alpha$ FWHMs are ~3 times larger than the H$\alpha$

<figure><img src="image_15.png" /><figcaption>**Figure 10.** H$\alpha$ FWHMs measured from PENELLOPE data vs. Ly$\alpha$ FWHMs extracted from models that do and do not include resonant scattering. The Ly$\alpha$ FWHMs are always larger than the measured H$\alpha$ values, with median ratios of FWHM(Ly$\alpha$)/FWHM(H$\alpha$) ~3 under both modeling frameworks. The discrepancy can be explained if the Ly$\alpha$ photons are propagating through a clumpy medium (Li & Gronke 2022), as expected within the accretion shocks surrounding TTSs.</figcaption></figure>

line widths, regardless of whether or not scattering is included in the models. This is consistent with the results of Orlitova et al. (2018), whose models of Ly$\alpha$ emission lines from star-forming green pea galaxies had to have FWHMs that were ~3 times broader than the Balmer lines from the same systems. This discrepancy requires the Ly$\alpha$ photons to propagate through a clumpy medium (Gronke et al. 2017; Li & Gronke 2022), as expected within the turbulent accretion shocks surrounding TTSs.

Alternatively, we explore whether the difference in emission line FWHMs is caused by degeneracies in model parameters by fitting a second set of shell models to the targets with PENELLOPE data and holding the Ly$\alpha$ line widths fixed at the measured H$\alpha$ FWHMs. In the case that H$\alpha$ scattering is minimal, this represents the “true” intrinsic Ly$\alpha$ emission line width prior to passing through the turbulent accretion shock. Figure 11 compares the best-fit H I column densities and shell velocities from the two modeling frameworks. Despite decreasing the Ly$\alpha$ FWHMs by a factor of ~3 on average to match the H$\alpha$ line widths, the best-fit H I column densities and velocities generally do not change significantly between models. This implies that the more significantly unconstrained pair of parameters is likely the effective temperature within the
<figure><img src="image_16.png" /><figcaption>**Figure 11.** Comparison of best-fit H I column densities (left) and shell velocities (right) in models when the Lyα line widths are held fixed at the measured Hα FWHMs (y axes) or left as variable parameters (x axes). Despite changing the Lyα FWHMs by a factor of ~3 on average between models, the column densities and velocities are roughly consistent, implying that the temperature of the scattering medium is the more strongly degenerate parameter. One target, Sz 69, has an H I velocity that is ~17 times larger when the Lyα FWHM is held fixed. The value corresponds to the stronger of two peaks in the posterior distribution extracted from the MCMC sampler. However, the larger value is nonphysical, as it exceeds the expected escape velocity by a factor of ~3.</figcaption></figure>

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

scattering medium and the H I column densities, rather than the H I column densities and intrinsic emission line widths (see also corner plot for Sz 111 in Appendix A). If a significant fraction of the intervening H I is located in the ISM, as reported by McJunkin et al. (2014), the simple shell models used here will not accurately capture the temperature gradient of the scattering medium. However, the dominant velocity signatures setting the observed Lyα flux ratios are still consistent with outflowing and accreting H I within the protostellar/disk/wind environments.

## 6.4. Implications for Chemistry in Protoplanetary Disks

Physical-chemical models that account for Lyα irradiation of protoplanetary gas disks identify several volatile-bearing species that are particularly sensitive to its effects. The large photodissociation cross sections near 1216 Å for H₂O, HCN, and C₂H₂ can result in decreased abundances of all three species (Bergin et al. 2003; Bethell & Bergin 2011), depending on the fraction of the UV radiation fields comprised of Lyα emission (Walsh et al. 2015). Disks around M dwarfs are expected to receive less UV radiation than T Tauri or Herbig Ae/Be disks, resulting in relatively larger abundances of UV-sensitive molecules (Walsh et al. 2015).

We do not find statistically significant correlations between stellar masses and total observed Lyα fluxes ($ρ = 0.0486, p = 0.76$) or Lyα flux ratios ($ρ = -0.0797, p = 0.611$) from our sample of primarily K and M type stars, although the relationships between reconstructed Lyα fluxes at the disk surfaces and stellar properties will be explored in future work (K. France et al. 2023, in preparation). The ratios of Lyα/FUV continuum emission are strongly correlated with the mass accretion rates derived from C IV emission (France et al. 2014), and there is also no correlation between stellar mass and accretion rates for the sample presented here. However, the correlations between Lyα flux ratios, [W3 − W4] colors, and 1600 Å H₂O dissociation “bump” detection rates (France et al. 2017) reported in this work are consistent with models predicting increased Lyα propagation with dust settling (Bethell & Bergin 2011) and observations of CN/HCN abundance ratios peaking within UV-exposed dust gaps (Bergner et al. 2021).

The ALMA Large Program “Molecules with ALMA at Planet-forming Scales” (MAPS) unexpectedly found that the observed CN/HCN ratios from AS 209, GM Aur, IM Lup, HD 163296, and MWC 480 did not show any clear trends with stellar UV luminosity, despite the five targets spanning two orders of magnitude in integrated flux between 910 and 2000 Å (Bergner et al. 2021). However, just two of these targets have Lyα spectra from HST: HD 163296 (P Cygni-like) and GM Aur (intermediate). AS 209 and MWC 480 were observed at low spectral resolution with IUE, and IM Lup does not have short wavelength UV coverage as of ULLYSES DR3 (see, e.g., Dionatos et al. 2019). Furthermore, physical-chemical models have thus far incorporated approximations of Lyα emission lines from different types of host stars (Walsh et al. 2015), as direct observations were not available. The true nature of Lyα as a driver of protoplanetary disk chemistry has not yet been studied from an observational or modeling perspective, but it can now be fully explored using the ULLYSES data.

## 7. Summary and Conclusions

We have used two different modeling frameworks to reproduce observed Lyα emission lines in HST-COS spectra of 42 CTTSs from the ULLYSES program, in an effort to explore the connection between the line profile shapes and outflow and accretion kinematics. Our results demonstrate that:

1. ~70% of targets in our sample have P Cygni-like Lyα emission lines, with blue velocity emission suppressed by outflowing H I. The remaining ~30% have inverse P Cygni-like profiles, a percentage that is roughly consistent with the observed fraction of submm transition disks with dust cavities carved by planet-disk interactions or photoevaporative or MHD winds (see, e.g., Andrews et al. 2011). The subset of our sample with inverse P Cygni-like profiles and spatially resolved submm observations show dust gaps or cavities, and a statistically significant negative correlation is observed between the ratios of red to blue Lyα flux and [W3 − W4] ($ρ = -0.33, p = 0.04$).
2. The observed emission line wings for 27/42 targets are well-reproduced by the models with resonant scattering in a simple shell. Although the geometry of T Tauri systems is more complex than the simple shell model used here, we interpret the best-fit parameters as average properties of the bulk intervening H I within inner disk winds, accretion flows, protoplanetary disks, and the ISM.

3. The model H I velocities are strongly correlated with the K-band veiling when $r_K$ is measured against synthetic spectra ($\rho = -0.66$; $p = 0.0006$), tentatively indicating the Ly$\alpha$ line profiles change shape as the hot inner disk is depleted.

4. No statistically significant correlation is detected between the model H I velocities and outer disk inclinations, indicating that the trends reported here are not solely caused by the viewing angles. We note that (a) submm observations are not available for targets with H I velocities much larger or smaller than the median value $v \sim -94 \text{ km s}^{-1}$, and (b) the Ly$\alpha$ emission lines are more strongly associated with the inner disks, which are frequently misaligned (see, e.g., Davies 2019; Bohn et al. 2022).

The resonant scattering models that were applied in this work use a simple shell geometry that may not be fully representative of the complex environment surrounding CTTSs. However, the trends reported here reveal that the ratios of red to blue Ly$\alpha$ fluxes observed with HST-COS are closely linked to inner disk outflow and accretion kinematics, which go along with the evolution of the dust disks. The Ly$\alpha$ flux ratios across our sample have a range of three orders of magnitude, a spread in emission line shape that may cause significant variations in the photodissociation rates of Ly$\alpha$-sensitive molecules that will be detected with JWST and ALMA (see, e.g., Bergin et al. 2003). Spatially resolved UV spectra (e.g., from HST-STIS) are required to fully explore multiphase Ly$\alpha$ resonant scattering through the protoplanetary disks, accretion flows, outflowing winds, and ISM.

Based on observations obtained with the NASA/ESA Hubble Space Telescope, retrieved from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science Institute (STScI). STScI is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555. Based on data obtained with ESO programmes 106.20Z8.002 and 106.20Z8.009. This work benefited from discussions with the ODYSSEUS team (HST AR-16129, Espaillat et al. 2022), https://sites.bu.edu/odysseus/. This research made use of Astropy,<fnref n="25" /> a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013, 2018). 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). Funded by the European Union under the European Unions Horizon Europe Research and Innovation Programme 101039452 (WANDA). 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. J.F. Gameiro was supported by the Fundação para a Ciência e Tecnologia (FCT) through the research grants UIDB/04434/2020 and UIDP/04434/2020. N.A. acknowledges funding from the Giacconi Fellowship at the Space Telescope Science Institute. R.A. acknowledges funding from the Science & Technology Facilities Council (STFC) through Consolidated Grant ST/W000857/1. This work has been supported by the project PRIN-INAF 2019 "Spectroscopically Tracing the Disk Dispersal Evolution (STRADE)".

## Appendix

### Best-fit Ly$\alpha$ Models for ULLYSES Targets

This appendix provides the observed Ly$\alpha$ profiles and best-fit models with and without resonant scattering for all targets in our sample (see Figure A1). We also include example corner plots showing the marginalized probability distributions for

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

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

**Figure A1.** Observed Ly$\alpha$ emission line (black, solid) and best-fit models with scattering (blue, dashed-dotted) and without scattering (red, dashed) for all targets in our sample. Scattering models are not shown for targets with low-S/N Ly$\alpha$ spectra, for which the MCMC chains did not converge. 2MASS J11432669-7804454 (a) and CS Cha (b) are shown as examples. The complete figure set (40 images) is available in the online Journal.

(The complete figure set (40 images) is available.)

<footnote n="25">²⁵ http://www.astropy.org</footnote>
<figure><img src="image_20.png" /><figcaption>**Figure A2.** Corner plot showing marginalized probability distributions for model parameters sampled for Sz 111, when Lyα scattering was not included. The Gaussian amplitude ($I_0$) is difficult to constrain from the emission line wings alone, but the samplers still converge within our criterion of >30% acceptance fractions.</figcaption></figure>

<table><caption>**Table A1**
Model Parameters and Prior Bounds for Lyα without Scattering</caption><thead><tr><th>Gaussian Amplitude<br>I<sub>0</sub></th><th>Gaussian FWHM<br>FWHM</th><th>Central Wavelength<br>λ<sub>0</sub></th><th>H I Column Density<br>N(H I)<sub>out</sub></th><th>H I Velocity<br>v<sub>out</sub></th></tr></thead><tbody><tr><td>(erg s<sup>-1</sup> cm<sup>-2</sup> Å<sup>-1</sup>)</td><td>(km s<sup>-1</sup>)</td><td>(Å)</td><td>(dex)</td><td>(km s<sup>-1</sup>)</td></tr><tr><td>(10<sup>-16</sup>, 10<sup>-11</sup>)</td><td>(250, 2000)</td><td>(1214.5, 1216.5)</td><td>(18.5, 22.0)</td><td>(-250, 250)</td></tr></tbody></table>

Note. All parameters were sampled over a continuous distribution.
<figure><img src="image_21.png" /><figcaption>**Figure A3.** Corner plot showing marginalized probability distributions for scattering model parameters sampled for Sz 111. Although the column density of H I (log $N_{HI}$) is well-constrained in this case, there are clear degeneracies between the intrinsic Ly$\alpha$ emission line width ($\sigma_i$) and the effective temperature of the scattering medium (log T/K). This is consistent with the results of Li & Gronke (2022) and is likely why the H I column densities reported here do not change when $\sigma_i$ is held fixed.</figcaption></figure>
<table><caption>**Table A2**
Model Parameters and Prior Bounds for Lyα with Scattering</caption><thead><tr><th>Shell Velocity<br>ν<sub>exp</sub></th><th>H I Column<br>Density<br>N(H I)<sub>out</sub></th><th>Effective<br>Temperature<br>log T/K</th><th>Gaussian Width<br>σ<sub>i</sub></th></tr></thead><tbody><tr><td>(km s<sup>-1</sup>)<br>(-495, 495)</td><td>(dex)<br>(15.9, 21.9)</td><td>(2.8, 6.0)</td><td>(km s<sup>-1</sup>)<br>(1, 800)</td></tr><tr><td>Dust Optical<br>Depth<br>τ<sub>d</sub></td><td>Equivalent<br>Width<br>EW<sub>i</sub></td><td>Central Velocity<br>Δ</td><td>Integrated Flux<br>Normalization<br>A</td></tr><tr><td>(0, 5)</td><td>(Å)<br>(1, 6000)</td><td>(km s<sup>-1</sup>)<br>(-100, 100)</td><td>(0.7, 2.5)</td></tr></tbody></table>

**Note.** $ν_{exp}$, $N(HI)$, and $log T/K$ were treated as discrete parameters. All other parameters were sampled over a continuous distribution.

model fits to the Lyα spectra from Sz 111, when scattering was
and was not included (Figure A2, Figure A3). All model
parameters and prior bounds for the MCMC samplers are
provided in Table A1 and Table A2.

ORCID iDs

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

Eleonora Fiorellino https://orcid.org/0000-0002-5261-6216

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

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

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

Kevin France https://orcid.org/0000-0002-1002-3674

Gregory J. Herczeg https://orcid.org/0000-0002-7154-6065

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

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

Péter Ábrahám https://orcid.org/0000-0001-6015-646X

Richard Alexander https://orcid.org/0000-0001-6410-2899

Jerome Bouvier https://orcid.org/0000-0002-7450-6712

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

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

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

Ágnes Kóspál https://orcid.org/0000-0001-7157-6275

Miguel Vioque https://orcid.org/0000-0002-4147-3846

References

Adams, T. F. 1972, ApJ, 174, 439

Ahn, S.-H., Lee, H.-W., & Lee, H. M. 2001, ApJ, 554, 604

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

Alcala, J. M., Natta, A., Manara, C. F., et al. 2014, A&A, 561, A2

Allard, F., Homeier, D., & Freytag, B. 2012, RSPTA, 370, 2765

Andrews, S. M., Wilner, D. J., Espaillat, C., et al. 2011, ApJ, 732, 42

Ansdell, M., Williams, J. P., Marel, N. v. d., et al. 2016, ApJ, 828, 46

Arulanantham, N., France, K., Hoadley, K., et al. 2021, AJ, 162, 185

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Banzatti, A., Pascucci, I., Edwards, S., et al. 2019, ApJ, 870, 76

Bergin, E. A., Du, F., Cleeves, L. I., et al. 2016, ApJ, 831, 101

Bergin, E., Calvet, N., D'Alessio, P., & Herczeg, G. J. 2003, ApJL, 591, L159

Bergner, J. B., Oberg, K. I., Bergin, E. A., et al. 2019, ApJ, 876, 25

Bergner, J. B., Oberg, K. I., Guzman, V. V., et al. 2021, ApJS, 257, 11

Bethell, T. J., & Bergin, E. A. 2011, ApJ, 739, 78

Blakely, D., Francis, L., Johnstone, D., et al. 2022, ApJ, 931, 3

Bohn, A. J., Benisty, M., Perraut, K., et al. 2022, A&A, 658, A183

Bonilha, J. R. M., Ferch, R., Salpeter, E. E., Slater, G., & Noerdlinger, P. D. 1979, ApJ, 233, 649

Bruderer, S. 2013, A&A, 559, A46

Campbell-White, J., Sicilia-Aguilar, A., Manara, C. F., et al. 2021, MNRAS, 507, 3331

Cazzoletti, P., van Dishoeck, E. F., Visser, R., Facchini, S., & Bruderer, S. 2018, A&A, 609, A93

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

Davies, C. L. 2019, MNRAS, 484, 1926

de La Reza, R., Quast, G., Torres, C. A. O., et al. 1986, ESA Special Publication, Vol. 263, New Insights in Astrophysics. Eight Years of UV Astronomy with IUE, ed. E. J. Rolfe & R. Wilson, 107

Dionatos, O., Woitke, P., Gudel, M., et al. 2019, A&A, 625, A66

Ercolano, B., & Pascucci, I. 2017, RSOS, 4, 170114

Espaillat, C., D'Alessio, P., Hernandez, J., et al. 2010, ApJ, 717, 441

Espaillat, C., Ingleby, L., Hernandez, J., et al. 2012, ApJ, 747, 103

Espaillat, C., Muzerolle, J., Najita, J., et al. 2014, in Protostars and Planets VI, ed. H. Beuther et al. (Tucson, AZ: Univ. Arizona Press), 497

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

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

Findeisen, K., Hillenbrand, L., & Soderblom, D. 2011, AJ, 142, 23

Fiorellino, E., Manara, C. F., Nisini, B., et al. 2021, A&A, 650, A43

Folha, D. F. M., & Emerson, J. P. 1999, A&A, 352, 517

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

France, K., Roueff, E., & Abgrall, H. 2017, ApJ, 844, 169

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

France, K., Schindhelm, R. N., Burgh, E. B., et al. 2011, ApJ, 734, 31

France, K., Schindhelm, R. N., Herczeg, G. J., et al. 2012, ApJ, 756, 171

Francis, L., & van der Marel, N. 2020, ApJ, 892, 111

Francis, L., van der Marel, N., Johnstone, D., et al. 2022, AJ, 164, 105

Frasca, A., Boffin, H. M. J., Manara, C. F., et al. 2021, A&A, 656, A138

Furlan, E., Watson, D. M., McClure, M. K., et al. 2009, ApJ, 703, 1964

Gárate, M., Delage, T. N., Stadler, J., et al. 2021, A&A, 655, A18

Gomez de Castro, A. I., Lopez-Santiago, J., Lopez-Martinez, F., et al. 2015, ApJS, 216, 26

Green, J. C., Froning, C. S., Osterman, S., et al. 2012, ApJ, 744, 60

Gronke, M. 2017, A&A, 608, A139

Gronke, M., Bull, P., & Dijkstra, M. 2015, ApJ, 812, 123

Gronke, M., & Dijkstra, M. 2016, ApJ, 826, 14

Gronke, M., Dijkstra, M., McCourt, M., & Peng Oh, S. 2017, A&A, 607, A71

Gronke, M., & Oh, S. P. 2018, MNRAS, 480, L111

Hansen, M., & Oh, S. P. 2006, MNRAS, 367, 979

Hartmann, L., Calvet, N., Gullbring, E., & D'Alessio, P. 1998, ApJ, 495, 385

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

Hashimoto, J., Muto, T., Dong, R., et al. 2021, ApJ, 911, 5

Hendler, N. P., Pinilla, P., Pascucci, I., et al. 2018, MNRAS, 475, L62

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

Hoadley, K., France, K., Alexander, R. D., McJunkin, M., & Schneider, P. C. 2015, ApJ, 812, 41

Huang, J., Andrews, S. M., Dullemond, C. P., et al. 2018, ApJL, 869, L42

Johns-Krull, C. M., & Valenti, J. A. 2001, ApJ, 561, 1060

Kidder, B., Mace, G., Lopez-Valdivia, R., et al. 2021, ApJ, 922, 27

Koenig, X. P., & Leisawitz, D. T. 2014, ApJ, 791, 131

Kudo, T., Hashimoto, J., Muto, T., et al. 2018, ApJL, 868, L5

Kurtovic, N. T., Perez, L. M., Benisty, M., et al. 2018, ApJL, 869, L44

Lawson, W. A., Lyo, A. R., & Muzerolle, J. 2004, MNRAS, 351, L39

Li, Z., & Gronke, M. 2022, MNRAS, 513, 5034

Lima, G. H. R. A., Alencar, S. H. P., Calvet, N., Hartmann, L., & Muzerolle, J. 2010, A&A, 522, A104

Long, F., Herczeg, G. J., Harsono, D., et al. 2019, ApJ, 882, 49

Loomis, R. A., Oberg, K. I., Andrews, S. M., & MacGregor, M. A. 2017, ApJ, 840, 23

Lubow, S. H., & D'Angelo, G. 2006, ApJ, 641, 526

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

Manara, C. F., Mordasini, C., Testi, L., et al. 2019, A&A, 631, L2

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

Manara, C. F., Testi, L., Natta, A., & Alcala, J. M. 2015, A&A, 579, A66

Manara, C. F., Testi, L., Natta, A., et al. 2014, A&A, 568, A18
Mauco, K., Hernandez, J., Calvet, N., et al. 2016, ApJ, 829, 38

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

McJunkin, M., France, K., Schneider, P. C., et al. 2014, ApJ, 780, 150

Miotello, A., Facchini, S., van Dishoeck, E. F., et al. 2019, A&A, 631, A69

Najita, J. R., Strom, S. E., & Muzerolle, J. 2007, MNRAS, 378, 369

Neufeld, D. A. 1990, ApJ, 350, 216

Neufeld, D. A. 1991, ApJL, 370, L85

Orlitova, I., Verhamme, A., Henry, A., et al. 2018, A&A, 616, A60

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

Pascucci, I., & Sterzik, M. 2009, ApJ, 702, 724

Pinilla, P., Benisty, M., Cazzoletti, P., et al. 2019, ApJ, 878, 16

Pinilla, P., Kurtovic, N. T., Benisty, M., et al. 2021, A&A, 649, A122

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

Powell, M. J. D. 1994, in Advances in Optimization and Numerical Analysis, ed. S. Gomez & J.-P. Hennart, Vol. 275 (Dordrecht: Springer), 51

Proffitt, C. R., Roman-Duval, J., Taylor, J. M., et al. 2021, RNAAS, 5, 36

Rebollido, I., Merín, B., Ribas, Á., et al. 2015, A&A, 581, A30

Ribas, Á., Espaillat, C. C., Macias, E., et al. 2017, ApJ, 849, 63

Rodriguez, J. E., Loomis, R., Cabrit, S., et al. 2018, ApJ, 859, 150

Roman-Duval, J., Proffitt, C. R., Taylor, J. M., et al. 2020, RNAAS, 4, 205

Ruaud, M., & Gorti, U. 2019, ApJ, 885, 146

Schindhelm, R., France, K., Burgh, E. B., et al. 2012a, ApJ, 746, 97

Schindhelm, R., France, K., Herczeg, G. J., et al. 2012b, ApJL, 756, L23

Schneider, P. C., Gunther, H. M., & France, K. 2020, Galax, 8, 27

Shakhovskoj, D., Grinin, V., Rostopchina, A., et al. 2006, A&A, 448, 1075

Simon, M., Guilloteau, S., Beck, T. L., et al. 2019, ApJ, 884, 42

Marel, N. v. d., Williams, J. P., Ansdell, M., et al. 2018, ApJ, 854, 177

van Dishoeck, E. F., Jonkheid, B., & van Hemert, M. C. 2006, FaDi, 133, 231

van Zadelhoff, G. J., Aikawa, Y., Hogerheijde, M. R., & van Dishoeck, E. F. 2003, A&A, 397, 789

Verhamme, A., Orlitova, I., Schaerer, D., & Hayes, M. 2015, A&A, 578, A7

Verhamme, A., Orlitova, I., Schaerer, D., et al. 2017, A&A, 597, A13

Verhamme, A., Schaerer, D., Atek, H., & Tapken, C. 2008, A&A, 491, 89

Walsh, C., Nomura, H., & van Dishoeck, E. 2015, A&A, 582, A88

Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868

Xu, Z., Herczeg, G. J., Johns-Krull, C. M., & France, K. 2021, ApJ, 921, 181

Zagaria, F., Clarke, C. J., Rosotti, G. P., & Manara, C. F. 2022, MNRAS, 512, 3538