UV Fluorescence Traces Gas and LyA Evolution in Protoplanetary Disks

Ultraviolet spectra of protoplanetary disks trace distributions of warm gas at radii where rocky planets form. We combine HST-COS observations of H2 and CO emission from 12 classical T Tauri stars to more extensively map inner disk surface layers, where gas temperature distributions allow radially stratified fluorescence from the two species. We calculate empirical emitting radii for each species under the assumption that the line widths are entirely set by Keplerian broadening, demonstrating that the CO fluorescence originates further from the stars (r ~ 20 AU) than the H2 (r ~ 0.8 AU). This is supported by 2-D radiative transfer models, which show that the peak and outer radii of the CO flux distributions generally extend further into the outer disk than the H2. These results also indicate that additional sources of LyA photons remain unaccounted for, requiring more complex models to fully reproduce the molecular gas emission. As a first step, we confirm that the morphologies of the UV-CO bands and LyA radiation fields are significantly correlated and discover that both trace the degree of dust disk evolution. The UV tracers appear to follow the same sequence of disk evolution as forbidden line emission from jets and winds, as the observed LyA profiles transition between dominant red wing and dominant blue wing shapes when the high-velocity optical emission disappears. Our results suggest a scenario where UV radiation fields, disk winds and jets, and molecular gas evolve in harmony with the dust disks throughout their lifetimes.


INTRODUCTION
The advent of the Atacama Large Millimeter/submillimeter Array (ALMA) has revolutionized studies of circumstellar disks, providing images of planet-forming material at unprecedented spatial resolution and sensitivity. By mapping the dust distributions down to ∼5 AU from the central star, sub-mm observations have revealed both small and large-scale dust substructures in most disks (see e.g., Andrews et al. 2018;Huang et al. 2018a;Long et al. 2018). This includes all observed disks in the Taurus star-forming region with effective emission radii larger than 50 au, implying that disks without detected substructures are simply too compact to resolve the features without even higher resolution observations (Long et al. 2019). The ALMA images are consistent with models of disk evolution, which predict asymmetric dust traps (van der Marel et al. 2013), rings (ALMA Partnership et al. 2015;Pérez et al. 2018;Guzmán et al. 2018;Huang et al. 2018b), and spiral arms (Huang et al. 2018c) as a result of physical mechanisms like planet-disk interactions, dust-trapping, and X-ray photoevaporation (see e.g., Johansen et al. 2004;Bae et al. 2016;Gorti & Hollenbach 2009). However, sensitivity requirements for detecting individual gas species prevent observations of molecular gas at the spatial resolution required to trace key carbon-, oxygen-, and nitrogen-bearing species in regions where rocky planets are expected to form (r < 5 AU). Alternative tracers of warm inner disk gas are required to map the material close to the central star.
Ultraviolet and infrared emission from warm CO (UV-CO, IR-CO; T ∼ 300 − 1500 K) and ultraviolet emission from hot H 2 (UV-H 2 ; T > 1500 K) are readily detected within the inner disk, in regions that are difficult to resolve at sub-mm wavelengths (r < 20 AU; see e.g., Herczeg et al. 2002;Najita et al. 2003;Salyk et al. 2011a;Brown et al. 2013;France et al. 2011bFrance et al. , 2012. The UV-H 2 features originate in a thin surface layer of the gas disk that is heated to T > 1500 K by UV radiation, allowing H 2 to be excited into higher vibrational levels (see e.g., Nomura & Millar 2005;Nomura et al. 2007;Ádámkovics et al. 2016). These vibrationally excited molecules are then "pumped" from the ground electronic state X 1 Σ + g into the first and second excited electronic states B 1 Σ + u , C 1 Π + u by Lyα photons, producing fluorescent emission lines as the molecules return to the ground state (Herczeg et al. 2002(Herczeg et al. , 2004. The features are spectrally resolved in data from HST -STIS and HST -COS, allowing observers to extract both the empirical radial locations of emitting gas (France et al. 2012) and the radial distributions of flux from the spectra (Hoadley et al. 2015;Arulanantham et al. 2018).
The UV-CO emission lines are produced by the same Lyα-pumping mechanism as UV-H 2 fluorescence (France et al. 2011b). However, models of the temperature and column density of emitting gas demonstrate that CO fluorescence originates from a cooler population of molecules (T < 600 K) (Schindhelm et al. 2012a), likely located at larger radii than hotter UV-H 2 . In order to extend our current picture of inner disk gas surface layers to more distant regions, we build on those findings in this work by comparing both UV-H 2 and UV-CO emission lines from disks around T Tauri stars. When both species are included, we tentatively identify a multi-layered sequence of disk evolution that is consistent with similar trends observed at optical and IR wavelengths. Our analysis is structured as follows: 1. We first estimate the purely empirical radial locations of the emitting gas, assuming an ideal model where the features originate in a Keplerian disk.
2. We use the molecular gas emission to reconstruct the Lyα radiation fields that reach the disk surfaces and pump H 2 and CO into excited electronic states.
3. We then feed the reconstructed Lyα profiles into 2-D radiative transfer models that map the full radial distributions of flux from both species.
4. However, we discover significant trends in deviations between the models and data, which demonstrate that the Lyα radiation reaching the molecular gas is sensitive to the overall state of disk evolution.
5. We discuss these results in the context of the disk dust distributions, presenting statistically significant correlations between Lyα morphologies, UV-H 2 and UV-CO fluxes, infrared spectral indices, and the presence of high-velocity [O I] λ6300Å emission.
These results lay the groundwork for more complex models of Lyα irradiation in protoplanetary disks, which will be the subject of future work.

TARGETS AND OBSERVATIONS
Our sample consists of 12 young stars with circumstellar disks that were observed with the Cosmic Origins Spectrograph onboard the Hubble Space Telescope (HST -COS; Green et al. 2012). Two of the targets were observed twice (AA Tau and RECX-15), and we include both sets of data for a total sample size of 14 spectra. We utilize the spectra acquired with the FUV observing modes, which were reduced with the CALCOS pipeline before being aligned and co-added into a final product for analysis (Danforth et al. 2010). The point source resolution of the FUV modes on HST -COS is ∆v ∼ 15 km s −1 (Osterman et al. 2011), which sets the radial velocity accuracy of the spectral features we measure and the minimum spectrally resolved line widths.
Although there are archival HST -COS spectra for many more T Tauri stars, our sample is limited to the objects with distinct CO emission out to J ∼ 20. Basic stellar and disk properties for each target are shown in Table 1, and details of the original observing programs are provided in Table 2. Both the G130M and G160M modes were used for most targets, providing wavelength coverage from ∼1100-1700Å. Various features in the spectra have been presented in previous work, including accretion rate diagnostics (Ardila et al. 2013), fluorescent emission lines from excited electronic states of hot H 2 (France et al. 2012;Hoadley et al. 2015), 12 CO and 13 CO A − X absorption bands (McJunkin et al. 2013), Lyα emission (Schindhelm et al. 2012b) and H 2 absorption along the Lyα profile (Hoadley et al. 2017), the FUV radiation field (Yang et al. 2012;France et al. 2014b), and connections to UV photochemistry (Arulanantham et al. 2020).
The HST -COS spectra of all systems in our sample include strong fluorescent H 2 emission lines, with most also showing high signal-to-noise fluorescence from the  b V4046 Sgr is a binary system, so the values listed here correspond to the two individual stellar masses (Rosenfeld et al. 2012).  (Hendler et al. 2018), and V4046 Sgr (Kastner et al. 2018;Ruíz-Rodríguez et al. 2019). These sub-mm gas and dust images provide additional context for inter-preting spectroscopic data from the disks we survey in this work. Previous long-slit observations of some T Tauri stars with HST -STIS show spatially extended UV-H 2 emission pumped by locally generated Lyα photons, in addition to the on-source emission originating from the inner disk. The features trace cavity walls (T Tau; Saucedo et al. 2003;Walter et al. 2003), shock-excited outflows (RU Lupi, DG Tau; Herczeg et al. 2005Herczeg et al. , 2006Schneider et al. 2013), and photoevaporative winds (GM Aur; Hornbeck et al. 2016). We note that since HST -COS is a slitless instrument, any similar spatially extended emission along its dispersion axis would be mapped to distinguishable velocity shifts in the 1-D spectra (see e.g., France et al. 2012). However, this effect would be hidden if the emission was extended along the cross-dispersion direction and did not exceed the 1.25" aperture radius (see e.g., France et al. 2011a).

UV-H 2
High signal-to-noise emission lines from UV-fluorescent H 2 are detected in all targets in our sample. The features originate in flared surface layers of the gas disk, where temperatures are at least T > 1500 K (Ádámkovics et al. 2014, 2016). Vibrationally excited H 2 (v ≥ 2) in these surface layers is pumped into excited electronic states by Lyα photons, producing a cascade of fluorescent emission lines as the molecules decay back to the ground state. A set of features with the same vibration and rotation ([v , J ]) upper level in the excited electronic state is referred to as a progression.
The interactive line fitting routine approximates each emission line as a Gaussian profile and linear continuum that are convolved with the wavelength-dependent linespread function (LSF) for HST -COS. Integrated line fluxes, FWHMs, and central wavelengths are then extracted from the best-fit Gaussian properties. Assuming the widths of the observed emission lines are dominated by Keplerian rotational broadening, an empirical radius of UV-H 2 emission can be calculated as (1) (Salyk et al. 2011a;France et al. 2012;Brown et al. 2013). A total of 11 emission lines from the [1, 4] and [1,7] progressions were fit with the Gaussian model, and an average FWHM was measured from the best-fit parameters of those features. The average FWHM was then used to estimate the spatial location of the gas with Equation 1.

UV-CO
The same Lyα-pumping mechanism is responsible for UV-CO emission as well (France et al. 2011b). Unlike UV-H 2 , UV-CO is pumped from rotational states in the v = 0 level of the ground electronic state, which allows cooler gas to fluoresce. The UV-CO bands observed with HST -COS are produced by molecules cascading from v = 14 in the A 1 Π electronic state (Schindhelm et al. 2012a).
In contrast to the spectrally resolved UV-H 2 lines, the close spacings of individual rotational states within UV-CO bands ∆v < 15 km s −1 make them impossible to resolve with HST -COS. However, the high J states (J > 10) have energy level spacings that are wide enough to partially separate the rotational structure (see Figure 1). This means that the J > 10 states appear as four to five distinct emission profiles in both UV-CO bands, with each feature consisting of three blended lines. Although the high J features are still not spectrally resolved, the separation between each set is wider than the pile-up of states at the bandhead.
Like the [1,4] and [1,7] progressions of UV-H 2 , CO molecules are pumped into the upper levels of high J states via photons at wavelengths where the Lyα line profile peaks in flux 1215.7 < λ < 1217Å . This produces strong emission lines, even when the temperature is too low to thermally populate the J > 10 energy levels (T ≤ 500 K; France et al. 2011b). The high J states are so prominent in some systems that the signal-to-noise is greater than that in the bandhead (J = 0 − 10). Since the high J states are also more widely separated, we focus this portion of the analysis on those features. However, our approach is still limited by the spectral resolution of the data, making the spatial constraints less robust than the empirical emitting regions derived from the well-resolved UV-H 2 features.
The blended UV-CO features were decomposed by fitting a multi-Gaussian model to the spectra. First, the same interactive line fitting routine used to measure the  UV-H 2 emission lines was used to fit each partially resolved blend of three lines (called "fingers" for simplicity) with a single Gaussian model, convolved with the HST -COS LSF and superimposed on a linear continuum. The central wavelength, Gaussian line width, and amplitude were left as free parameters for each finger. We included as many fingers as possible, excluding only the features that later did not have high enough S/N to be fit with the multi-Gaussian models. The width of each single Gaussian was then used to constrain the widths of the three narrower emission lines composing the corresponding finger, allowing us to obtain multiple "measurements" of the unresolved features. This procedure to fit the fingers is equivalent to fitting multiple UV-H 2 emission lines from a single progression.
Next, each finger was further divided using the multi-Gaussian approach. Each model consisted of three Gaussian emission lines, which were first superimposed and then convolved with the HST -COS LSF. Individual Gaussians were given the same line width, as expected for emission originating from the same Keplerian radius. The amplitude ratios between the three blended features were fixed based on the Lyα fluxes at the appropriate pumping wavelengths and the branching ratios describing the likelihood of each pathway out of the Lyα-pumped upper level. Best-fit models were obtained by varying the amplitude and central wavelength of the transition in the blended triplet with the largest likelihood of transitioning out of its upper level and the single Gaussian line width for all three features. For each target, the radial location of emitting gas was then cal- . Left: Partially resolved "fingers" of UV-CO emission (black curve) were fit with the same interactive fitting routine used on the UV-H2 emission lines. Each finger was fit with a single Gaussian and a linear continuum, with variable amplitude, central wavelength, FWHM, and continuum slope and intercept (teal curve). Right: The observed features were further decomposed using multi-Gaussian models, where the widths of individual transitions within a single finger were required to be identical and the width of the blended lines equal to the FWHM from the single Gaussian model. The radial locations of emitting gas were then calculated from the average decomposed widths of all measured fingers. Observed features which are not included in the Gaussian models shown in the figures were not able to be further divided into individual emission lines. The continuum was assumed to be linear for uniformity across the sample of disks, and uncertainties in the slope and intercept are folded into the errors on the empirical radii. culated from the average decomposed line width from all measured fingers. We note that the linear continuum is fine for most systems in our sample, although it is not a perfect approximation. With this in mind, uncertainties in the continuum slopes and intercepts are also folded into the errors on the empirical radii.
We emphasize that the high J UV-CO emission lines used to derive the empirical radii are not spectrally resolved. While Gaussian decompositions can reproduce the observed features, Weber et al. (2020) point out that such models are still unable to distinguish contributions to the total flux from different emitting regions (e.g., multiple wind scale heights). They note that any emission components originating in winds, rather than the Keplerian disk, will project some velocity signature onto the observed spectra which cannot be untangled with Gaussians alone. Rather than comparing the UV-H 2 and UV-CO empirical radii for each individual system, we instead focus on general trends, which provide a starting point for building more complex models.

RESULTS
Consistent with previous studies (Herczeg et al. 2004;France et al. 2012), we find that the empirical UV-H 2 emitting radius is inside r < 2 AU for all targets (see Table 4). Figure 3 compares the radial locations of UV-CO and UV-H 2 , demonstrating that the UV-CO originates much further from the central star(s) than the UV-H 2 in all seven spectra with strong J > 10 emission.
The UV-CO also appears to trace a different population of gas than the IR-CO emission (Salyk et al. 2011b;Brown et al. 2013), which is produced by material that is roughly radially co-located with the UV-H 2 (France et al. 2012). Although the UV-CO radii are more uncertain than the UV-H 2 locations, even the blended sets of high J lines are narrower than the individual UV-H 2 features, confirming that the emitting gas is further away. For example, the strongest UV-CO blend in the DF Tau spectrum has F W HM = 19.8 km s −1 , compared to the average UV-H 2 line width of F W HM = 63 km s −1 . This further supports our picture of radially stratified fluorescent gas.
We note that spatially resolved HST-STIS observations of UV-H 2 emission lines from DF Tau show blueshifted components at velocity shifts up to ∼ −40 km s −1 (Herczeg et al. 2006). This excess emission originates in a wind, which HST -COS can only distinguish from the on-source disk emission when the source is either oriented along the dispersion direction or the wind extends beyond the 1.25" aperture radius (see e.g., France et al. 2011a). If unresolved wind signatures are present in UV-H 2 emission lines from other targets in this sample, the emitting region boundaries derived here would be biased toward smaller radii. However, Hoadley et al. (2015) found no significant asymmetries between the H 2 line shapes observed with HST -COS from any targets that overlap with this work. Pontop-pidan et al. (2011) also note that sub-Keplerian IR-CO emission from slow, uncollimated winds is likely always present but can be treated as an extension of the disk emission, since the gas maintains its Keplerian rotation as it travels up the wind. In the case of the UV-H 2 and UV-CO features from targets in our sample, the components will only be resolvable with multi-object UV spectroscopy on a future flagship mission (see e.g., France et al. 2017aFrance et al. , 2019 or with UV spectroastrometry (see e.g., Pontoppidan et al. 2008). Figure 3 also includes an estimated CO ice line location for each system, derived from the relationship between stellar luminosity, disk radius, and condensation temperature presented in Long et al. (2018; r ∝ L * /T 4 conden ). We find that the empirical UV-CO emission radius is closer to the star than the derived CO ice line location in all but two systems, where the UV-CO emission comes from just beyond the ice line (DM Tau and RECX-15). However, the regions of disks where temperatures are cool enough for CO freeze-out are likely 2-D surfaces instead of radial cutoffs (see e.g. Bosman et al. 2018;Qi et al. 2019). This thermal structure allows gas-phase CO to survive in surface layers at radii where solids have formed closer to the disk midplane.
We also find that the UV-CO emission in the 2013 spectrum of AA Tau may originate from larger radii (R UV−CO = 23 ± 6 AU) than the same feature in the 2011 spectrum of AA Tau (R UV−CO = 15 ± 6 AU). These empirical R UV−CO values are the medians of all radii where emitting gas is present. The differences in the AA Tau spectra therefore imply that contributions from hot gas in the inner disk appear more prominently in the emission lines from 2011 than in those from 2013, shifting the empirical emitting radius to a smaller value in 2011. This effect was also seen in the UV-H 2 emission lines from the two sets of spectra, with the 2011 spectrum showing an empirical UV-H 2 emitting radius of R UV−H2 = 0.7 ± 0.2 and the 2013 observations showing R UV−H2 = 1.5 ± 0.4.
Changes between the 2011 and 2013 AA Tau UV-H 2 spectra were attributed to extra absorbing material that moved into the line-of-sight after 2011 Hoadley et al. 2015), perhaps caused by an increase in the inner disk scale height (Covey et al. 2021). In the case of the UV-CO, the occulting material masks the contribution from the innermost emitting gas, leading to narrower line profiles that are dominated by more distant molecules. This shifts the empirical radius (R UV−CO ) further out in the disk in 2013, with the dominant contribution coming from closer to the CO ice line. Future HST and JWST observations of AA Tau's mysterious geometry will provide a more detailed picture of interactions between the occulting material, stellar radiation field, and molecular gas disk.

2-D RADIATIVE TRANSFER MODELS OF UV-H 2 AND UV-CO
In order to derive more detailed radial flux distributions from the fluorescent emission lines, we fit 2-D radiative transfer models to the UV-H 2 and UV-CO spectra. We build on the models developed by Hoadley et al. (2015), which create UV-H 2 emission line profiles by distributing Lyα photons to a population of hot gas in a thin surface layer of the disk. This version efficiently accounts for UV-H 2 lines from four different progressions, pumped by photons with four different wavelengths along the red side of the Lyα profile. Since the [1,4] progression has the strongest signal-to-noise across the full sample of disks, we focus on fitting the radiative transfer models to those features. This allows us to extract the most reliable information on the underlying radial distributions of fluorescent gas while maintaining consistency with Hoadley et al. (2015) and minimizing the required computation time to find the best-fit models.
The (v − v ) = (14−3) and (v − v ) = (14−4) bands of UV-CO each consist of lines pumped by >60 different Lyα wavelengths λ = 1214.2 − 1218.8Å . Rather than adding the UV-CO radiative transfer properties to the existing code for UV-H 2 , the models were adapted into an object-oriented framework that can accommodate multiple molecular species. A comparison of UV-H 2 models from the two versions (see Appendix A) demonstrates that both yield similar results that are consistent with those presented in Hoadley et al. (2015).
b Targets with missing values of FWHM CO,(14−4) and R CO did not show UV-CO bands with distinct separations between sets of three high J states (J > 10; see Figure 1), so the multi-Gaussian fitting procedure described in Section 3 could not be applied to their spectra. Since individual UV-CO emission lines are not spectrally resolved, the radii derived here are not as robust as the locations measured from the UV-H2 features.
the stellar mass. The 2-D mass density distribution of gas is then calculated as with radial surface density distribution The characteristic radius (r c ) is the location in the disk where the density distribution begins to fall off exponentially. A power-law with coefficient γ, which describes the radial surface density distribution Σ (r) for a disk with kinematic viscosity ν ∝ r −γ (Lynden-Bell & Pringle 1974), is used to map the distribution interior to r c . Number densities of both H 2 and CO in ground states with vibrational level v and rotational level J ([v, J]) are calculated from the mass density distribution under LTE conditions: where X H2,CO represents the fraction of the total gas density consisting of H 2 or CO, respectively, g [v,J] is the statistical weight of state [v, J], and Z [v,J] is the partition function. Optical depths in each ground state are then calculated as where σ is the cross section for absorbing a Lyα photon at the pumping wavelength. However, overlap between cross sections at adjacent pumping wavelengths can lead to an overestimate of the total optical depth, by allowing the same photon to be absorbed by molecules in multiple rovibrational levels of the ground electronic state.
To correct for this effect, we scale each τ by the fractional optical depth expected for that transition at the temperature (T ) and column density (N ) corresponding to its grid point (McJunkin et al. 2014;Liu & Dalgarno 1996;Wolven et al. 1997), such that  . Average gas emission radii for 12 CO and H2, which are also presented in Table 4. We find that UV-CO (teal circles; calculated in this work) originates at more distant radii than UV-H2 in all systems with strong emission from states with J > 10. In all cases, this emission comes from inside or very close to the The total flux for each transition out of the Lyα-pumped upper level is then calculated as where η is a geometric filling factor describing the fraction of gas exposed to Lyα photons (Herczeg et al. 2004), F * ,Lyα is the relevant Lyα pumping flux at grid point (r, z), s (r, z) is the sightline from the observer to a gas parcel in the disk, and B mn describes the probability of transitioning from [v , J ] → [v , J ]. We do not fit for η directly in this analysis; rather, the value is folded into the total Lyα pumping flux reaching the gas disk. Emission line profiles are calculated by collapsing this 2-D grid in Keplerian velocity space and convolving the models with the HST -COS LSF. We compare the UV-H 2 models to the data in velocity space, but the model UV-CO line profiles are first superimposed in wavelength space to produce a band of emission lines. The observed UV-H 2 emission lines are all shifted to radial velocities of zero before fitting the radiative transfer models to the data, based on the central wavelengths from the Gaussian models. For the UV-CO features, we shift the observed bandheads (J = 0) to radial velocities of zero. The impact of this shift on the modeling results is discussed further in Section 5.4.

Reconstructed Lyα Profiles
Although the Lyα line profile is included in the HST -COS bandpass, the core of the feature is always attenuated by circumstellar and interstellar absorption (see e.g., McJunkin et al. 2014). The line is further contaminated by geocoronal emission, making the highvelocity wings the only portions that are directly observable (McJunkin et al. 2014). However, the line profile can be reconstructed using the emission from fluorescent UV-H 2 lines (Herczeg et al. 2004;Schindhelm et al. 2012b). Progression fluxes, calculated by integrating over all UV-H 2 features from the same upper level i (F H2,i ), are treated as "y" data points. Lyα pumping wavelengths act as corresponding "x" values, such that (x i , y i ) = (λ Lyα,i , F H2,i ). The data points are then fitted with a model profile, assumed to be an outflow-absorbed Gaussian emission line (I Lyα ; Schindhelm et al. 2012b): with τ out calculated as a Voigt profile centered at outflow velocity v out with H I column density N out .
To convert the progression fluxes to intensities that can be fit with the Gaussian model, each y i data point is divided by an equivalent width calculated for Lyαabsorbing H 2 with temperature T H2 and column density N H2 We applied this Lyα reconstruction method to obtain profiles for CW Tau, T Cha, and the 2013 AA Tau spectrum, which had not yet been modeled in the literature. The remaining model profiles were taken from previous work that utilized the same procedure (Schindhelm et al. 2012b;France et al. 2014b;Arulanantham et al. 2018). Uncertainties in the best-fit reconstruction parameters typically make the total reconstructed Lyα fluxes accurate to within ∼20% (Schindhelm et al. 2012b).

UV-H 2 Model Fitting Procedure
UV-H 2 models were fit to emission lines from the [1,4] progression by allowing six parameters to vary (Hoadley et al. 2015): • height of the emitting layer (2 H p < z/r < 7 H p ) • temperature at 1 AU (500 K < T 1 AU < 5000 K) • power law exponent for the temperature distribution (−2.5 < q < 2.5) • power law exponent for the surface density distribution (0 < γ < 2) • radius where the surface density distribution turns over from a power law to an exponential decline (0.1 AU < r char < 20 AU) • total mass of H 2 1 × 10 −5 M < M H2 < 0.1 M Although the model can also generate emission lines from the [1,7], [0,1], and [0,2] progressions, we follow the same approach as Hoadley et al. (2015) and focus on the [1,4] features with the strongest signal-to-noise across the full sample of targets. Since the UV-H 2 line shapes do not change significantly between progressions (see e.g., Hoadley et al. 2015;Arulanantham et al. 2018), we are able to extract sufficient information about the underlying distributions of emitting gas from this single progression.
We obtain a best-fit UV-H 2 model for each disk by minimizing the mean square error (M SE): which describes the goodness-of-fit of the model profiles (ỹ i ) to the observed UV-H 2 emission lines (y i ).
We choose the M SE over test statistics that include the measurement uncertainties, which are estimated by the HST -COS pipeline and are not necessarily derived from a distribution that represents the true errors. This makes it especially difficult to find robust fits for targets with weak but distinct UV-H 2 features, since anomalously small uncertainties in the continuum can bias test statistics (e.g. χ 2 ) toward models that predict nondetections. Uncertainties on the best-fit model parameters, which contribute more to the total error than the measurement errors, were estimated by using a Markov chain Monte Carlo (MCMC) sampler to explore the parameter space defined above (Foreman-Mackey et al. 2013). Models that provide the closest match to the observed UV-H 2 emission lines are representative of the underlying radial distributions of flux from hot, Lyα-pumped H 2 in the inner disk (see Table 5). However, the surface gas masses and scale heights are often degenerate, as are the temperatures at 1 AU and the temperature powerlaw exponents (Hoadley et al. 2015). Furthermore, several targets in our sample have UV-H 2 line profiles with unusually broad wings relative to the narrow widths of the line cores (e.g., DF Tau, RY Lupi; see Herczeg et al. 2006;Banzatti & Pontoppidan 2015;Arulanantham et al. 2018, Hoadley et al. in prep). This makes it difficult for a single model to simultaneously capture both the broad wing emission from material close to the star and the narrow core emission from more distant gas. Both the degenerate model parameters and occasionally non-standard line shapes make it difficult to directly extrapolate uncertainties in the radial flux distributions from the uncertainties in the best-fit model parameters.
Instead, we analyze the radial flux distributions through the 100 models with the smallest M SE values for each target, which represent the top 1% of all models generated during the MCMC resampling. Taken together, this group of models includes combinations of parameters that fit both the broad and narrow emission line components in DF Tau and RY Lupi. For all targets, the sets of 100 best-fit models include the degenerate parameter combinations of z/r, M H2 , T 1 AU , and q. We produce a final radial flux distribution for each system by calculating the median flux value from the set of 100 models at each radial grid point. This method does not necessarily capture a perfect best-fit model, but it allows us to compare the flux distributions while accounting for degenerate model parameters. All 100 models are also overplotted against the observed emission lines in Appendix A, demonstrating that the set collectively accounts for nearly all of the UV-H 2 emission from each target. Although we used a different mass parameterization than Hoadley et al. (2015), the UV-H 2 flux distributions presented here follow the evolutionary trends discovered in that work (see Appendix A).
Finally, we note that the radiative transfer models always produce double-peaked Keplerian line profiles, while the data generally show single-peaked emission lines. Previous work has attributed this to slow, uncollimated winds that fill in the low velocity emission and can be treated as extensions of the disk surfaces (see e.g., Pontoppidan et al. 2011;Hoadley et al. 2015). However, the contributions from multiple scale heights likely impose some velocity signatures on the observed emission lines (see e.g., Weber et al. 2020) that are not accounted for in our modeling approach.

UV-CO Model Fitting Procedure
The main difference in the treatment of UV-H 2 versus UV-CO is the temperature distribution within the fluorescent gas. Herczeg et al. (2004) found that the UV-H 2 features from TW Hya were consistent with emission from a gas layer with T ∼ 2000 − 3000 K and log N (H 2 ) = 18 − 20, so the models presented here allow the temperature of the H 2 distribution to reach as high as the 5000 K dissociation limit (Hoadley et al. 2015). Assuming a CO/H 2 ratio of 10 −4 (France et al. 2014a), the corresponding column densities of CO would be log N (CO) ∼ 14 − 16 in the 2500 K layer.
However, cooling via CO ro-vibrational emission becomes particularly efficient in regions where H I abundances reach the threshold for molecule formation to begin (log N (H) ∼ 21). The gas temperature is then quickly reduced to ∼1000 K in denser surface layers of the disk (Ádámkovics et al. 2014). This implies that hotter layers, where the UV-H 2 emission originates, do not have significant enough populations of CO to either cool the gas or produce detectable UV-CO fluorescence.
These column density and temperature conditions are also consistent with 4.7 µm CO rovibrational emission in surface layers of the inner disk. The rotational temperatures of these features never exceed T ∼ 1800 K and are typically much cooler than the T ∼ 1500 K threshold required for Lyα pumping of H 2 (Salyk et al. 2011b). Since the IR-CO features originate much closer to the star than the UV-CO (R IR−CO ∼ 0.1 − 5 AU; see Figure 3, France et al. 2012), where gas heating via stellar irradiation is highest, the temperature of this warm gas provides a rough upper limit for the temperature distribution of UV-CO. The range of temperatures within the UV-CO emitting region is therefore smaller than what is expected for UV-H 2 , including gas between 20-1800 K to also avoid temperatures where CO is expected to be frozen out of the gas phase (see e.g., Bosman et al. 2018).
Since the UV-CO and UV-H 2 emission lines trace radially separated regions of the gas disk, we expect the Lyα radiation field to look different to each molecular species. To account for this, we add two additional variable parameters to the UV-CO models: the column densities and radial velocities of H I between the UV-H 2 and UV-CO emitting regions. We restrict the H I column densities to range between the values that best fit the UV-H 2 progression fluxes (Schindhelm et al. 2012b) and those estimated for the total ISM absorption (McJunkin et al. 2016). A new Lyα absorption feature is then superimposed on the intrinsic Gaussian emission line generated by the H 2 reconstruction procedure.
In order for Lyα photons to reach warm CO at radii between r = 15 and r = 30 AU, the flared disk geometry forces the UV-CO emitting region to higher values of z/r and larger characteristic radii than the inner disk UV-H 2 . The variable parameters for UV-CO are then: • height of the emitting layer (7 H p < z/r < 11 H p ) • temperature at 1 AU (20 K < T 1 AU < 1800 K) • power law exponent for the temperature distribution (−2.5 < q < 2.5) • power law exponent for the surface density distribution (0 < γ < 2) a 5% of the total UV-H2 flux is contained within this radius.
b Radial location of UV-H2 flux peak.
c 95% of the total UV-H2 flux is contained within this radius.
d Bi-modal posterior distributions are retrieved for T , and q. The reported best-fit values are taken from the mode with the highest peak, and the uncertainties are extended to include the degenerate set of parameters.
• radius where the surface density distribution turns over from a power law to an exponential decline (0.1 < r char < 100 AU) • total mass of CO 10 −20 M < M CO < 10 −5 M • column density of intervening H I (determined individually for each target) • radial velocity of intervening H I (−300 < v < 300 km s −1 ) Best-fit models were again obtained by minimizing the M SE test statistic. Uncertainties on the model parameters were derived by carrying out MCMC resampling, using 1000 walkers to explore the full parameter space described above (Foreman-Mackey et al. 2013; see Table  6). Model fits to the UV-H 2 and UV-CO emission lines from each individual target are provided in Appendix B, along with the corresponding radial flux distributions. We find that the median temperature of fluorescent CO at 10 AU is ∼1100 K, in very good agreement with models of the disk surface that predict similar temperatures in layers where CO column densities begin to increase (Ádámkovics et al. 2014). We also find that 95% of the UV-CO flux is typically contained within r out ∼ 9 AU, demonstrating that the emitting regions extend roughly three times further than the radial distributions of UV-H 2 (r out ∼ 3 AU). However, r out for the UV-CO is only about half the size of the empirical radius presented in Table 4 ( R CO ∼ 20 AU).
The number densities of intervening H I between the UV-H 2 and UV-CO emitting regions were also calculated for each target (see Table 7). Values range between n HI ∼ 0.1 − 3.4 × 10 6 cm −3 , falling well within constraints on the total densities expected to be driven out by inner disk winds (see e.g., Turner et al. 1999). However, we note that the number densities were calculated assuming the difference in r out from the model distributions of UV-H 2 and UV-CO flux. Any systematic uncertainties in that analysis will then propagate through the number density calculations.
Although the 2-D radiative transfer approach provides good fits to the UV-H 2 emission lines, we find that even the best-fit models are unable to fully capture the UV-CO bands. The largest discrepancies between the data and models are seen between λ ∼ 1316.2−1317Å, where the models significantly underpredict the UV-CO flux from AA Tau, CW Tau, RECX-15, and DF Tau. The models also often appear to miss the line centers in states a 5% of the total UV-CO flux is contained within this radius.
b Radial location of UV-CO flux peak.
c 95% of the total UV-CO flux is contained within this radius.
with J > 10, showing emission lines that reproduce either the long or short wavelength sides of the features but not both (see e.g., AA Tau). The deviations are caused by the reconstructed, outflow-absorbed Lyα profiles, which do not appear to match the radiation fields seen by the UV-CO. Figure 4 shows two different best-fit UV-CO models for UX Tau, with the residuals color-coded from blue to red based on the appropriate Lyα pumping wavelengths. The model generated with an outflow absorbed Lyα profile underpredicts the data from transitions with J ∼ 8 − 12, which are pumped at Lyα wavelengths between λ ∼ 1214.5 − 1215.7Å. We also find that the models underpredict the data at longer wavelength transitions pumped by redder Lyα photons, causing the modeled emission lines to miss the observed line centers while still managing to reproduce the spacing between transitions. Since these are the high J features used to derive the empirical radii, the underpredictions are likely responsible for the discrepancies between R CO from the Gaussian fitting procedure (Table 4) and r peak from the radiative transfer models (Table 6). Additional sources of Lyα pumping photons, both from line center and red-shifted wavelengths, must be included to more accurately reproduce the UV-CO emission from large radii (r > 15 au).
This evidence for a missing source of red-shifted Lyα at the disk surface is consistent with previous analyses of UV-H 2 , since the spectra can show significant emission from progressions excited by pumping photons up to ∼1000 km s −1 redward of line center (e.g. DF Tau; Herczeg et al. 2006). Furthermore, even the best-fit UV-H 2 models are not able to fully reproduce both the broad high-velocity wings and narrower line cores in many systems (see e.g., DF Tau, RY Lup, V4046 Sgr). This may also imply that an additional Lyα source is responsible for pumping the more distant gas responsible for the low velocity emission, while the photons generated at the accretion shock reach the gas at smaller radii. Physically, the extra radiation may originate in a disk wind, from which Lyα would appear red-shifted in the frame of the disk. To better capture this missing flux, we added additional low-velocity Gaussian emission components approximating what one might expect to be generated locally in winds or outflows to the model Lyα profiles. Although the excess emission boosted irradiation of the gas by photons near line center (see Figure 4 and Appendix B), it still does not fully account for the deviations between the UV-CO models and data.
A shortcoming of our modeling approach is that we must attempt to simultaneously characterize the strength of Lyα emission at the disk surface, the radial location(s) that receive the strongest irradiation from b H I column densities measured in this work, corresponding to absorption between the UV-H2 and UV-CO emitting regions. f H I number densities were calculated using RCO −RH 2 . Any systematic uncertainties in the estimates for those radii will then propagate through the number density calculations.
Lyα sources beyond the accretion shock (e.g., from disk winds), and the physical structure of the fluorescent gas layer. We find that the degeneracies between sets of model parameters, particularly the gas temperature gradients and Lyα pumping fluxes, make it difficult to converge on accurate radial distributions of UV-CO flux as derived for the UV-H 2 emission (Hoadley et al. 2015).
In the following sections, we instead further analyze the observed Lyα line wings and construct a more detailed picture of interactions between the radiation field and the molecular gas disk. This interpretation lays the ground work for future reconstructions of Lyα profiles from T Tauri stars, which will require more observations of other Lyα sensitive molecules (e.g., HCN, H 2 O, H 2 CO; Bergin et al. 2003;Bethell & Bergin 2011) in disks with available UV spectra.

DISCUSSION
As described in the previous section, the 2-D radiative transfer models of UV-H 2 and UV-CO depend strongly on Lyα irradiation of the gas disk. This is seen most clearly in the best-fit models for the UV-CO bands, which cannot always be reproduced with a reconstructed outflow-absorbed Lyα profile. Since interstellar absorption and geocoronal emission prevent us from directly detecting the Lyα line core, we focus our discussion on what the observed high-velocity wings reveal about the radiation field reaching the molecular gas. We note that the wings were not observed for CW Tau, AA Tau (2013), or RECX-15 (2013), since those spectra were acquired with an orientation of the G130M grating on HST -COS that does not include the Lyα transition (λ central = 1222Å, instead of 1291Å). Figure 5 compares the best-fit model gas temperatures from Table 6 to the ratios of observed UV-CO emission from high J states relative to the bandhead (F J >10 /F J =0−10 ). The relationship between the model gas temperatures at 10 AU and the observed UV-  Figure 4. Best-fit UV-CO models for UX Tau, with residuals color-coded from blue to red based on the appropriate Lyα pumping wavelength. The model acquired with a Lyα radiation field generated exclusively at the central accretion shock (left) underpredicts the data near 1316.5Å, which can be corrected by including an extra Gaussian Lyα profile to capture the missing flux near line center (right). However, it is more difficult to reproduce the UV-CO emission from upper levels pumped by redder Lyα photons, implying that some local source of redshifted photons (e.g., a disk wind) may be responsible for exciting gas at larger radii. We note that the spectral resolution of the data is not high enough to determine whether some portion of the UV-CO emission itself originates in a wind, instead of the disk surface.

Lyα Irradiation Drives the UV-CO Band Shapes
CO flux ratios (hereafter referred to as "band shapes") is an anti-correlation (Spearman rank coefficient = −0.61; p = 0.02), implying that hotter gas does not necessarily lead to more populated high J states.
Since the UV-CO emitting region extends to larger radii than the UV-H 2 , it is possible (or even likely) that Lyα profiles reconstructed from the UV-H 2 emission inside r < 2 AU are not representative of the radiation field at r > 15 AU. To gain a better understanding of the Lyα profiles seen by the UV-CO, we compare the ratios of observed Lyα emission from the red and blue wings of the profiles F to the UV-CO band shapes in Figure 6. Unlike the gas temperatures, the ratios of red wing to blue wing observed Lyα emission show a strong positive correlation with the ratios of UV-CO flux in high J relative to bandhead states (Spearman rank coefficient = 0.699; p = 0.01). The UV-CO band shapes are therefore driven by the shape of the illuminating Lyα profile instead of the temperature distribution within the emitting region.
With this in mind, Figure 7 also compares the ratios of UV-H 2 fluxes from progressions pumped by red ([1, 4] , [1, 7]) and blue Lyα photons ([4, 4]) (France et al. 2012) to the ratios of red to blue Lyα fluxes. Once again, we find a strong correlation between the flux ratios (Spearman rank coefficient = 0.89; p = 5 × 10 −4 ), demonstrating that disks with more blue-wing Lyα emission also show stronger UV-H 2 fluxes in the [4,4] progression (see also McJunkin et al. 2016). Since the shape of the illuminating radiation field clearly has a strong influence on the pumping radiation received at the disk surface, we discuss the mechanisms that dictate the ratios of red to blue Lyα fluxes in the following sections.  Figures 5, 6, 7, and 8 and find that they can be divided into three categories: disks showing more flux in the UV-CO bandhead than in the high J states, systems showing roughly similar bandhead and high J fluxes, and targets with more high J than bandhead emission. As depicted in Table 8, targets that fall within a specific grouping also have similar observed Lyα profiles and infrared spectral indices. For example, disks with more bandhead than high J UV-CO emission also have stronger blue wing Lyα flux and larger SED slopes between 13 and 31 µm that indicate more emission from cold, outer disk dust relative to warm, inner disk material (Furlan et al. 2009; e.g. UX Tau A). We expand on each of the three UV-CO categories in the following subsections. In particular, we investigate the connection between Lyα emission and the structure of the gas and dust disks, which appears analagous to the evolutionary sequence discovered in optical and infrared emission from forbidden lines in disk winds (see e.g., [

Outflow Dominated Profiles
The disks around AA Tau, RECX-15, and DF Tau show the most prominent UV-CO emission from high J states, with relatively suppressed emission from the bandhead. All three disks also have observed Lyα line wings with significantly more red than blue flux. For these systems, the shape of the Lyα radiation field reaching both the UV-CO emitting region and the observer appears generally consistent with the model used to reconstruct the Lyα profile from the UV-H 2 features (see Section 4.2).
Lyα photons from young stellar objects are generated at the accretion shock and in heated funnel flows (see e.g., Alencar et al. 2012), where high opacity material and Stark broadening produce strong, broad emission lines (FWHM∼ 550 − 900 km s −1 ; Schindhelm et al. 2012b). The photons must then travel through optically thick outflowing H I before reaching the molecular gas disk. In the frames of the gas disk and the observer, the outflowing gas removes flux from the blue side of the Lyα line λ ∼ 1214.5 − 1215.7Å . The resulting Lyα radiation field seen by the molecular gas disk then has stronger red than blue flux.
In agreement with this picture, all three targets in this category are driving spatially resolved jets traced by high velocity optical forbidden line emission. AA Tau  This outflowing material must also include H I, which further absorbs the Lyα photons generated near the accretion shock. Figure 9 shows that a handful of targets in our sample have ratios of red wing to blue wing Lyα flux around ∼1: LkCa 15, RECX-11, RY Lupi, T Cha, and V4046 Sgr. However, the UV-CO band shapes from T Cha and V4046 Sgr appear more similar to the disks with stronger blue-wing Lyα emission (CS Cha, DM Tau, and UX Tau A), so we have placed them in that category instead. The remaining three targets span a relatively broad range of bandhead to high J flux ratios.

Intermediate Profiles
RY Lupi and LkCa 15, which host large, dust depleted inner cavities (r ∼ 69 and 76 AU, respectively; see e.g. Francis & van der Marel 2020), have stronger UV-CO emission from the bandhead, while RECX-11 has roughly similar levels of bandhead and high J flux. The SED slopes between 13 and 31 µm for these targets are scattered, with RY Lupi and LkCa 15 having n (13 − 31) = −0.26 and n (13 − 31) = 0.62, respectively (Furlan et al. 2009). Disks in this category seem to represent a middle phase of evolution, where the warm dust producing the 13 µm emission has not necessarily been cleared from the system and is still able to shield UV-fluorescent gas.  None of the targets in this group show high-velocity [O I] 6300Å components from jets (Banzatti et al. 2019), indicating that emission from jets and inner disk outflows is not present. However, accretion is still ongoing in all three disks, with rates ranging from 0.017−1×10 −8 M yr −1 (Ingleby et al. 2011(Ingleby et al. , 2013France et al. 2017b). This is supported by the Hα emission line from RECX-11, which has a wide blue wing (v ∼ 300 km s −1 ) and inverse P Cygni (red) absorption characteristic of freefalling gas accreting onto the star (Ingleby et al. 2011).

Infall Dominated Profiles
Five disks in our sample have stronger integrated UV-CO fluxes in the bandhead (J = 0 − 10) than in the high J states (J > 10): CS Cha, DM Tau, T Cha, UX Tau, and V4046 Sgr. All five systems have steeper disk SEDs than targets in the other two groups, as measured from the slopes between 13 and 31 µm (n (13 − 31) ; Furlan et al. 2009). This indicates that their dust disks produce more emission from cold, outer disk grains than warm, inner disk material, as confirmed by dust-depleted gaps or cavities resolved at sub-mm wavelengths (see e.g., Kudo et al. 2018;Hendler et al. 2018;Ruíz-Rodríguez et al. 2019).
Within this group, 3/5 targets 1 have observed Lyα profiles with stronger blue wing λ < 1214.5Å than red wing λ > 1217Å emission. The strong blue Lyα wings contrast with the model used to reconstruct the Lyα radiation field, which uses outflowing H I to truncate the blue sides of the line profiles (Herczeg et al. 2004 Figure 8. Flux ratios of observed red wing to blue wing Lyα emission versus n (13 − 31) for disks with dust cavities, indicating that disks with more blue Lyα flux show less emission from warm grains relative to emission at longer wavelengths from cold dust. We omit DF Tau from the figure, since it is the only system without dust substructure resolved at infrared or sub-mm wavelengths. This means that the slope of its SED is representative of the degree of grain growth and dust settling towards the midplane, and therefore the disk flaring angle (Furlan et al. 2009). Since the remaining targets do have gaps and cavities in their dust distributions, the slopes of the SEDs instead trace the disk temperatures, and therefore radial locations, where grains have been removed from the disk instead of settled to the midplane (see e.g., Espaillat et al. 2010). DF Tau should instead be compared to other targets without resolved dust substructure, which will be possible once all new data from the ULLYSES program is available. An interactive version of this figure is provided online, providing disk and stellar properties for each target. Schindhelm et al. 2012b). A model with red-shifted ISM absorption can reproduce the observed line wings by suppressing the red Lyα emission (see e.g., McJunkin et al. 2014). However, strong correlations between the Lyα flux ratios, UV-CO band shapes, UV-H 2 flux ratios, and n (13 − 31) (see Figure 9) indicate that the red wing suppression occurs before the Lyα photons enter the ISM. Similar red-absorbed shapes have long been identified in the hydrogen recombination lines (Edwards et al. 1994;Reipurth et al. 1996). The morphologies are attributed to magnetospherically driven mass infall, where high velocity material moving away from both the ob-server and the gas disk absorbs the red sides of emission lines generated near the accretion shock (Calvet & Hartmann 1992;Hartmann et al. 1994). It is possible then that the Lyα profiles with more blue than red emission are also dominated by mass infall. This picture is consistent with the mass accretion rates measured from C IV emission (see e.g., Ardila et al. 2013), which are not necessarily slower in targets with more depleted dust distributions.
Trends in forbidden emission from [O I] and [Ne II] provide some insight into the Lyα profiles for this group of targets. Pascucci et al. (2020) find that CS Cha, DM Tau, T Cha, UX Tau, and V4046 Sgr all have Without these feedback mechanisms, the remaining material in the inner disk may be able to rapidly drain onto the star, as expected from models of inside-out disk dispersal (see e.g., Ercolano et al. 2011). In order to connect the shapes of the observed Lyα profiles to the properties of absorbing H I, we estimate the absorption velocities where the fluxes goes to zero (see Figure 11). In addition to the targets from our original sample, we also include spectra from systems with observed Lyα wings but noisy UV-CO emission lines. These extra targets are 2MASS J16083070-3828268, HT Lup, and MY Lup (Arulanantham et al. 2020), RU Lup (see e.g., Herczeg et al. 2005), RW Aur (see e.g., Woitas et al. 2002), and TW Hya (see e.g., Herczeg et al. 2002).
The absorption velocities were derived by fitting a 2nd degree polynomial to the strongest observed Lyα wing for each target, between the velocities where the wing emission peaks and the geocoronal emission disrupts the observed spectrum. The best-fit polynomial was then extrapolated to find the velocity where the model flux reaches zero. The three disks with more blue than red Lyα emission (DM Tau, UX Tau, CS Cha) have absorption velocities of −339, −260, and −197 km s −1 , respec-tively. The remaining targets, which all have positive absorption velocities and more red than blue Lyα emission, are somewhat scattered, with no clear trend between the ratios of observed Lyα flux ratios and absorption velocities. While mass inflow must be happening in these targets, as evidenced by the high accretion rates, its signature red-shifted absorption may be masked by the ongoing jet and outflow emission. This connection between the absorbing H I and the observed Lyα wings may become more clear when more targets with dominant blue wing emission are observed.

SUMMARY AND CONCLUSIONS
We have presented a joint analysis of Lyα line profiles and fluorescent emission from H 2 and CO molecules in a sample of 12 disks around CTTSs, using spectra from HST -COS. The features originate in surface layers of the gas disk, where optical depths are low enough for substantial amounts of Lyα flux to reach the molecules. We analyze the spectral features both empirically and with 2-D radiative transfer models, finding that 1. The UV-H 2 and UV-CO features originate in radially separated regions of the disk, with UV-H 2 observed closer to the star (r ∼ 0.8 AU) than UV-CO (r ∼ 20 AU) in all targets with partially resolved high J (J > 10) UV-CO states. However, these radii were derived assuming the line widths are set entirely by Keplerian broadening in the disk. Emission from molecules in outflows or winds is not accounted for (see e.g., Weber et al. 2020 Figure 9. Summary of empirical tracers used to categorize the targets in Table 1 (see also Table 8, Figure 10). The "UV-H2" column shows bars representing the total UV-H2 flux from each of three progressions, pumped by blue ([4,4]) and red ([1,4], [1,7]) Lyα photons. This demonstrates that the [4,4] progression fluxes become relatively stronger as the blue Lyα emission increases. CW Tau is not included, since its Lyα wings were not observed with HST -COS. Systems with similar ratios of observed red to blue Lyα emission have similar UV-CO band shapes, UV-H2 fluxes from progressions pumped by blue ([4,4]) and red ([1,7] In the absence of high-velocity outflowing material, infall dominated Lyα profiles are produced as H I absorbs the red sides of the emission lines. UV-H2 and UV-CO emitting regions are represented based on the median empirical radii measured for each group of targets (see Table 8). (c) The infrared spectral indices, n (13 − 31), for disks with dust cavities; 3. However, 2-D radiative transfer models fit to the UV-CO emission lines underpredict the observed flux from progressions pumped by central and red Lyα photons. This causes the radiative transfer models to map the fluorescent gas to smaller radii than the empirical values.
4. It is difficult to map the flux distributions of fluorescent gas beyond the empirically derived radii, since the 2-D radiative transfer modeling results indicate that we are not accounting for some hidden source(s) of Lyα irradiation.
Furthermore, we find that only the three targets with the largest ratios of red to blue Lyα line wing emission (AA Tau, DF Tau, RECX-15) show high velocity [O I] λ6300 emission associated with jets. Taken together, these results imply that jets and outflows turn off at some phase of disk evolution, allowing residual inner disk material to rapidly drain onto the stars. Any existing trends will become more clear over the next few years, as data from the HST Ultraviolet Legacy Library of Young Stars as Essential Standards (ULLYSES) Director's Discretionary program becomes available.

ACKNOWLEDGEMENTS
We are grateful to Zachary Taylor and Klaus Pontoppidan for enjoyable and insightful discussions in various stages of this work. We also thank the anonymous referee for their very detailed and thoughtful comments that greatly improved the clarity of this manuscript. NA was supported in part by NASA Earth and Space Science Fellowship grant 80NSSC17K0531, and KH is supported by the David & Ellen Lee Prize Postdoctoral Fellowship in Experimental Physics at Caltech. This paper also made use of data and financial support from HST program GO-15070, and AB acknowledges support from HST grant HST-GO-15128.001-A to the University of Colorado at Boulder. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/ consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This work utilized the RMACC Summit supercomputer, which is supported by the National Science Foundation (awards ACI-1532235 and ACI-1532236),  UV-H 2 emission lines from almost all targets in the sample presented here were first analyzed by Hoadley et al. 2015, who developed the 2-D radiative transfer models. As described in Section 4, we have adapted the original models into an object-oriented framework that accommodates the UV-CO features as well. In this appendix, we compare our results to the previously published work in Hoadley et al. (2015) and demonstrate that the best-fit model temperatures are consistent within ±1σ error bounds. We also describe the methods of statistical interpretation used here in greater detail.
In this work, uncertainties on the best-fit model parameters were derived using MCMC re-sampling (Foreman-Mackey et al. 2013), which extracts the posterior probability for the model parameters, p (Θ|y i ): where Θ = [z/r, γ, T 1 AU , q, r char , M H2 ]. We have assumed uniform prior distributions for the model parameters (p (Θ)), with boundaries defined in Hoadley et al. 2015 (UV-H 2 ) and Section 4.4 of this work (UV-CO).
We also adopt a Gaussian likelihood function: where y i andỹ i represent the observed and model data points and σ i represents the measurement uncertainties in the data. However, we ignore σ i in the algorithm since it doesn't necessarily follow a Gaussian distribution. We find that the chain converges when 1000 walkers are allowed to probe the parameter space over 150 steps. The MCMC method allows for finer sampling of the posterior distributions than the discrete grid search method used in Hoadley et al. (2015). For example, the temperature distributions in Hoadley et al. (2015) are spaced by 500 K, while the MCMC chain can probe temperature differences down to <1 K. The finer resolution shows distinct double-peaked posterior distributions on the T and q parameters for most targets, which were not possible to capture with the coarser grid search (see Figure 12). Figure 13 compares the best-fit UV-H 2 temperatures from Hoadley et al. (2015) to two temperatures derived in this work. The first temperature, which we call the primary value, corresponds to the strongest peak in the posterior distribution. The secondary value is extracted from the weaker peak. In Figure 12, which shows the marginalized posterior distributions for UX Tau, the primary value is ∼3200 K, and the secondary value is ∼1500 K. CS Cha, RECX-11, and V4046 Sgr do not show significant secondary peaks, so we only plot the primary temperature values in Figure 13.
We find that the temperatures from Hoadley et al. (2015) are consistent with either the primary or secondary temperature value for all targets. The values for AA Tau (2011), DM Tau, RECX-15 (2010), and UX Tau show very close agreement with the secondary peak temperature values, indicating that the primary peak may not have been fully sampled with the grid search method. However, the secondary peak parameters represent equivalent model fits to the data that the MCMC re-sampling simply did not favor as strongly as the primary peak parameters.
To propagate the uncertainties in the best-fit model parameters to error bounds on the radial flux distributions, including those with double-peaked posterior distributions, we generate models from the 100 sets of parameters with the smallest values of the goodness-offit (M SE) test statistic (see Equation 12). This allows us to capture flux distributions corresponding to both the primary and secondary temperature values, without including intermediate temperatures with significantly lower probabilities. The resulting distributions of UV-H 2 flux are shown in Figure 14. We obtain final UV-H 2 flux distributions by calculating the median flux value at each radial grid point, and error bounds are estimated as the minimum and maximum fluxes. The results represent the uncertainties in the fluorescent flux distributions as a result of the uncertainties in the best-fit model parameters. Figure 15 compares the median UV-H 2 flux distributions for the primordial and transitional disks, as classified and presented in the same way as Figure 9 of Hoadley et al. 2015. We find that the flux distributions inside ∼0.1 AU follow the trends discovered in that work, with transitional disks showing no emitting material this close to the star and primordial disks showing strong peaks at these radii. However, we note a discrepancy with the outer radii of the distributions presented here. This is caused by different methods of calculating the partition functions and effective optical depths in the modeling code, which accounts for reducing the column density of fluorescent gas at radii >10 AU.
Instead of seeing UV-H 2 distributions extending to more distant radii in transitional disks than in primordial systems (Hoadley et al. 2015), we observe truncated flux distributions in the transitional disks. This may be attributed to the clearing of warm dust in the inner disk (r < 50 AU), which exposes molecules in surface layers of the disks to intense photodissociating radiation fields that also dissipate the gas in warm winds. By contrast, the extended distributions of fluorescent material in the primordial disks may be indicative of higher total gas masses in systems with larger remaining dust reservoirs. However, as noted by Hoadley et al. (2015), the 2-D radiative transfer models of H 2 in a thin disk surface layer are not particularly sensitive to the mass contained in the underlying bulk gas reservoir. Hoadley 2015 Primary Peak Secondary Peak Figure 13. Comparison of best-fit temperatures at 1 AU from Hoadley et al. (2015) and this work. The primary peak represents the temperature favored by the MCMC chain, while the secondary peak represents a high-likelihood temperature that was favored less than the primary peak.  Figure 15. Comparison of radial flux distributions from primordial versus transitional disks. We find that the adapted radiative transfer models reproduce the trends discovered by Hoadley et al. (2015), where primordial systems have more UV-H2 emission at radii close to the star than transitional disks.  Residuals on the best-fit UV-CO model are color-coded from blue to red, based on the Lyα pumping wavelength required to excite the upper level of the transition. The two epochs of UV-CO spectra from AA Tau (bottom right) are offset by ∼24 km s −1 . Since we do not see significant wavelength shifts between the individual transitions when the data are overplotted (bottom right), this shift is likely caused by the uncertainty in the HST -COS FUV wavelength solution (∼15 km s −1 ). However, the UV-CO and UV-H2 fluxes do change between epochs, and the differences in best-fit models to the two spectra are caused by the uncertainties in the best-fit Lyα profiles required to reproduce the features.  and radial distributions of flux for both species (bottom left) from the disk around DF Tau. Deviations between the UV-H2 models and data point to a two-component line profile, with pieces originating in radially separated populations of hot gas (see e.g., Banzatti & Pontoppidan 2015, Hoadley et al., in prep). Residuals on the best-fit UV-CO model are color-coded from blue to red, based on the Lyα pumping wavelength required to excite the upper level of the transition.