# UV Fluorescence Traces Gas and Ly$\alpha$ Evolution in Protoplanetary Disks

Nicole Arulanantham  Space Telescope Science Institute, 3700 San Martin Drive,
Baltimore, MD 21218, USA  Kevin France  Laboratory for Atmospheric and Space
Physics, University of Colorado, 392 UCB, Boulder, CO 80303, USA  Keri Hoadley
Department of Physics, Mathematics, and Astronomy, California Institute of
Technology, 1200 East California Blvd., Pasadena, CA 91125, USA  Department of
Physics and Astronomy, University of Iowa, Van Allen Hall, Iowa City, IA 52242
P.C. Schneider  Hamburger Sternwarte, Gojenbergsweg 112, 21029 Hamburg,
Germany  Catherine C. Espaillat  Institute for Astrophysical Research,
Department of Astronomy, Boston University, 725 Commonwealth Avenue, Boston,
MA 02215, USA  H.M. Günther  MIT, Kavli Institute for Astrophysics and Space
Research, 77 Massachusetts Ave., Cambridge, MA 02139, USA  Gregory J. Herczeg
Kavli Institute for Astronomy and Astrophysics, Peking University, Yiheyuan 5,
Haidian Qu, 100871, Beijing, China  Alexander Brown  Center for Astrophysics
and Space Astronomy, University of Colorado, 389 UCB, Boulder, Colorado 80309,
USA

###### Abstract

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 $\left(r\sim 20\,\rm{AU}\right)$ than the H2
$\left(r\sim 0.8\,\rm{AU}\right)$. 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 Ly$\alpha$ 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
Ly$\alpha$ 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 Ly$\alpha$ 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.

stars: pre-main sequence, protoplanetary disks, molecules

##  1 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 $\sim$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\sim 300-1500$
K) and ultraviolet emission from hot H2 (UV-H2; $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. 2011b, 2012). The
UV-H2 features originate in a thin surface layer of the gas disk that is
heated to $T>1500$ K by UV radiation, allowing H2 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 $\left(X^{1}\Sigma_{g}^{+}\right)$
into the first and second excited electronic states
$\left(B^{1}\Sigma_{u}^{+},\,C^{1}\Pi_{u}^{+}\right)$ by Ly$\alpha$ photons,
producing fluorescent emission lines as the molecules return to the ground
state (Herczeg et al., 2002, 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$\alpha$-pumping mechanism
as UV-H2 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
$\left(T<600\,\text{ K}\right)$ (Schindhelm et al., 2012a), likely located at
larger radii than hotter UV-H2. 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-H2 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. 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. 2.

We use the molecular gas emission to reconstruct the Ly$\alpha$ radiation
fields that reach the disk surfaces and pump H2 and CO into excited electronic
states.

  3. 3.

We then feed the reconstructed Ly$\alpha$ profiles into 2-D radiative transfer
models that map the full radial distributions of flux from both species.

  4. 4.

However, we discover significant trends in deviations between the models and
data, which demonstrate that the Ly$\alpha$ radiation reaching the molecular
gas is sensitive to the overall state of disk evolution.

  5. 5.

We discuss these results in the context of the disk dust distributions,
presenting statistically significant correlations between Ly$\alpha$
morphologies, UV-H2 and UV-CO fluxes, infrared spectral indices, and the
presence of high-velocity [O I] $\lambda$6300 Å emission.

These results lay the groundwork for more complex models of Ly$\alpha$
irradiation in protoplanetary disks, which will be the subject of future work.

##  2 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 $\Delta v\sim 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^{\prime\prime}\sim 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 $\sim$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 H2 (France et al., 2012; Hoadley et al.,
2015), 12CO and 13CO $A-X$ absorption bands (McJunkin et al., 2013),
Ly$\alpha$ emission (Schindhelm et al., 2012b) and H2 absorption along the
Ly$\alpha$ 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).

Table 1: Sample of CTTS with UV-H2 and UV-CO Emission  Target | $d$ | $A_{V}$ | $M_{\ast}$ | $i_{d}$ |  $n_{13-31}$16 | References aa(1) Bailer-Jones et al. 2018; (2) Gaia Collaboration et al. 2018; (3) France et al. 2017b; (4) Loomis et al. 2017; (5) Bacciotti et al. 2018; (6) Shakhovskoj et al. 2006; (7) Andrews et al. 2011; (8) Lawson et al. 2004; (9) Cabrit et al. 2006; (10) van der Marel et al. 2018; (11) Lawson et al. 1996; (12) Huélamo et al. 2015; (13) Ginski et al. 2018; (14) Guilloteau et al. 2014; (15) Rodriguez et al. 2010; (16) Furlan et al. 2009; (17) Rosenfeld et al. 2012
---|---|---|---|---|---|---
| [pc] | [mag] | $\left[M_{\odot}\right]$ | $\left[{}^{\circ}\right]$ |  |
AA Tau (2011) | 136.7 | 0.34 | 0.8 | 59.1 | -0.51 | 1, 3, 3, 4
AA Tau (2013) | 136.7 | 0.34 | 0.8 | 59.1 | -0.51 | 1, 3, 3, 4
CS Cha | 175.4 | 0.16 | 1.05 | 24.2 | 2.89 | 1, 3, 3, 13
CW Tau | 130.7 |  $\cdots$ | 0.69 | 59 | -0.65 | 1, -, 14, 5
DF Tau | 124.5 | 0.54 | 0.19 | 24 | -1.09 | 1, 3, 3, 6
DM Tau | 144.5 | 0.48 | 0.5 | 35 | 1.3 | 1, 3, 3, 7
LkCa 15 | 158.2 | 0.31 | 0.85 | 49 | 0.62 | 1, 3, 3, 7
RECX-11 | 98.3 | 0.03 | 0.8 | 70 | -0.8 | 1, 3, 3, 8
RECX-15 (2010) | 91.6 | 0.02 | 0.4 | 60 | -0.2 | 1, 3, 3, 8
RECX-15 (2013) | 91.6 | 0.02 | 0.4 | 60 | -0.2 | 1, 3, 3, 8
RY Lupi | 159.1 | 0.10 | 1.71 | 68 | -0.26 | 2, 3, 3, 10
T Cha | 109.6 |  $\cdots$ | 1.5 | 67 | 1.67 | 2, -, 12, 12
UX Tau A | 139.4 | 0.51 | 1.3 | 35 | 1.83 | 1, 3, 3, 11
V4046 Sgr | 72.3 | 0 |  $0.86+0.69$bbV4046 Sgr is a binary system, so the values listed here correspond to the two individual stellar masses (Rosenfeld et al., 2012). | 33 | 0.32 | 1, 3, 17, 15
Table 2: _HST_ -COS Observations  Target | PI:Program IDaaAll the HST-COS data used in this paper can be found in MAST: https://doi.org/10.17909/t9-pxwh-py47 (catalog https://doi.org/10.17909/t9-pxwh-py47). | Observing Modes/Central Wavelengths
---|---|---
AA Tau (2011) | G. Herczeg:11616 | G130M/1291,1327; G160M/1577, 1600, 1623
AA Tau (2013) | K. France:12876 | G130M/1222; G160M/1577, 1589, 1611, 1623
CS Cha | G. Herczeg:11616 | G130M/1291, 1327; G160M/1577, 1600, 1623
CW Tau | K. France:15070 | G130M/1222; G160M/1577, 1589
DF Tau | J. Green:11533 | G130M/1291, 1300, 1309, 1318; G160M/1589, 1600, 1611, 1623
DM Tau | G. Herczeg:11616 | G130M/1291, 1327; G160M/1577, 1600, 1623
LkCa 15 | G. Herczeg:11616 | G130M/1291, 1327; G160M/1577, 1600, 1623
RECX-11 | G. Herczeg:11616 | G130M/1291, 1327; G160M/1577, 1600, 1623
RECX-15 (2010) | G. Herczeg:11616 | G130M/1291, 1327; G160M/1577, 1600, 1623
RECX-15 (2013) | K. France:12876 | G130M/1222; G160M/1577, 1589, 1611, 1623
RY Lupi | C. F. Manara/P. C. Schneider:14469 | G130M/1291, 1318; G160M/1577, 1611
T Cha | A. Brown:15128 | G130M/1291; G160M/1577, 1600
UX Tau A | G. Herczeg:11616 | G130M/1291, 1327; G160M/1577, 1600, 1623
V4046 Sgr | J. Green:11533 | G130M/1291, 1300, 1309, 1318; G160M/1589, 1600, 1611, 1623

The HST-COS spectra of all systems in our sample include strong fluorescent H2
emission lines, with most also showing high signal-to-noise fluorescence from
the 12CO $\left(v^{\prime}-v^{\prime\prime}\right)=\left(14-3\right)$ and
$\left(v^{\prime}-v^{\prime\prime}\right)=\left(14-4\right)$ bands. Two disks
with large dust cavities (LkCa15 and RY Lupi; Francis & van der Marel 2020)
show much weaker UV-CO features than the rest of the group. We include them
here to study combinations of disk geometry and Ly$\alpha$ irradiation that
can produce strong UV-H2 but less prominent UV-CO. Furthermore, five other
disks in our sample have been observed with moderate to high resolution at
sub-mm wavelengths: AA Tau (Loomis et al., 2017), CW Tau (Piétu et al., 2014;
Bacciotti et al., 2018), DM Tau (Kudo et al., 2018), CS Cha, LkCa 15, RY Lupi
(Francis & van der Marel, 2020), T Cha (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 interpreting 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-H2 emission pumped by locally generated Ly$\alpha$
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. 2005, 2006;
Schneider 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).

##  3 Methods

###  3.1 UV-H2

High signal-to-noise emission lines from UV-fluorescent H2 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 H2 $\left(v\geq 2\right)$ in these surface layers
is pumped into excited electronic states by Ly$\alpha$ 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
$\left([v^{\prime},J^{\prime}\right])$ upper level in the excited electronic
state is referred to as a progression.

UV-H2 fluxes from the $\left[v^{\prime},J^{\prime}\right]=\left[1,4\right]$
and $\left[v^{\prime},J^{\prime}\right]=\left[1,7\right]$ progressions (see
Table 3) were measured using an interactive Python tool (SELFiE; Arulanantham
et al. 2018). These features typically have higher signal-to-noise than
emission lines from other progressions, since molecules are pumped into the
upper levels of the transitions by photons with wavelengths near the peak of
the Ly$\alpha$ emission line $\left(1215.73,1216.07\,\rm{\AA}\right)$. This
allows larger fractions of H2 to be pumped into the
$\left[v^{\prime},J^{\prime}\right]=\left[1,4\right]$ and
$\left[v^{\prime},J^{\prime}\right]=\left[1,7\right]$ states (see e.g.,
Herczeg et al. 2004).

Table 3: Measured UV-H2 Emission Lines  Progression | Pumping Wavelength | Line IDaasgn$\left(\Delta J\right)=P$ if $\Delta J=-1$ or $R$ if $\Delta J=+1$ | Rest Wavelength
---|---|---|---
$\left[v^{\prime},J^{\prime}\right]$ | [Å] |  $\left(v^{\prime}-v^{\prime\prime}\right)$sgn($\Delta J$)($J^{\prime\prime}$) | [Å]
$\left[1,7\right]$ | 1215.73 | (1-6)P(8) | 1467.08
|  | (1-7)R(6) | 1500.45
|  | (1-7)P(8) | 1524.65
|  | (1-8)R(6) | 1556.87
$\left[1,4\right]$ | 1216.07 | (1-6)R(3) | 1431.01
|  | (1-6)P(5) | 1446.12
|  | (1-7)R(3) | 1489.57
|  | (1-7)P(5) | 1504.76

The interactive line fitting routine approximates each emission line as a
Gaussian profile and linear continuum that are convolved with the wavelength-
dependent line-spread 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-H2 emission can
be calculated as

$\left<R_{H_{2}}\right>=GM_{\ast}\left(\frac{2\sin i_{d}}{FWHM}\right)^{2}$
(1)

(Salyk et al., 2011a; France et al., 2012; Brown et al., 2013). A total of 11
emission lines from the $\left[1,4\right]$ and $\left[1,7\right]$ 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.

###  3.2 UV-CO

The same Ly$\alpha$-pumping mechanism is responsible for UV-CO emission as
well (France et al., 2011b). Unlike UV-H2, 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^{\prime}=14$ in the $A^{1}\Pi$ electronic state
(Schindhelm et al., 2012a).

In contrast to the spectrally resolved UV-H2 lines, the close spacings of
individual rotational states within the
$\left(v^{\prime}-v^{\prime\prime}\right)=(14-4)$ and
$\left(v^{\prime}-v^{\prime\prime}\right)=(14-3)$ UV-CO bands $\left(\Delta
v<15\,\rm{km\,s}^{-1}\right)$ make them impossible to resolve with _HST_ -COS.
However, the high $J^{\prime\prime}$ states $\left(J^{\prime\prime}>10\right)$
have energy level spacings that are wide enough to partially separate the
rotational structure (see Figure 1). This means that the $J^{\prime\prime}>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^{\prime\prime}$ 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-H2, CO molecules are pumped into
the upper levels of high $J^{\prime\prime}$ states via photons at wavelengths
where the Ly$\alpha$ line profile peaks in flux
$\left(1215.7<\lambda<1217\,\rm{\AA}\right)$. This produces strong emission
lines, even when the temperature is too low to thermally populate the
$J^{\prime\prime}>10$ energy levels (T $\leq 500$ K; France et al. 2011b). The
high $J^{\prime\prime}$ states are so prominent in some systems that the
signal-to-noise is greater than that in the bandhead
$\left(J^{\prime\prime}=0-10\right)$. Since the high $J^{\prime\prime}$ 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-H2 features.

Figure 1: Left: HST-COS spectrum of DF Tau, including a single UV-H2 emission
line from the $\left[v^{\prime},J^{\prime}\right]=[1,4]$ progression and the
full $\left(v^{\prime}-v^{\prime\prime}\right)=(14-3)$ band of UV-CO. The
wavelength of the UV-H2 transition is marked with a red dashed line, and lab
wavelengths of individual UV-CO transitions are denoted with blue dashed-
dotted lines. Right: Zoomed in view of the high $J^{\prime\prime}$ region from
the same spectrum. Each blend of three transitions can be decomposed into
individual Gaussian emission line components, to empirically estimate the
radii from where the UV-CO originates (see Equation 1).

Figure 2: 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.

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-H2 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-H2 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$\alpha$ fluxes at the appropriate pumping
wavelengths and the branching ratios describing the likelihood of each pathway
out of the Ly$\alpha$-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 calculated 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^{\prime\prime}$ 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-H2 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.

##  4 Results

Consistent with previous studies (Herczeg et al., 2004; France et al., 2012),
we find that the empirical UV-H2 emitting radius is inside $r<2$ AU for all
targets (see Table 4). Figure 3 compares the radial locations of UV-CO and
UV-H2, demonstrating that the UV-CO originates much further from the central
star(s) than the UV-H2 in all seven spectra with strong $J^{\prime\prime}>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-H2
(France et al., 2012). Although the UV-CO radii are more uncertain than the
UV-H2 locations, even the blended sets of high $J^{\prime\prime}$ lines are
narrower than the individual UV-H2 features, confirming that the emitting gas
is further away. For example, the strongest UV-CO blend in the DF Tau spectrum
has $FWHM=19.8$ km s-1, compared to the average UV-H2 line width of $FWHM=63$
km s-1. This further supports our picture of radially stratified fluorescent
gas.

We note that spatially resolved HST-STIS observations of UV-H2 emission lines
from DF Tau show blueshifted components at velocity shifts up to $\sim-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-H2 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 H2 line shapes observed with HST-COS from any targets that overlap
with this work. Pontoppidan 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-H2 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. 2017a, 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\propto\sqrt{L_{\ast}/T_{\rm{conden}}^{4}}$). 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.

Table 4: Measured UV-H2 and UV-CO FWHMs and Average Emitting Radii  Target | $\rm{FWHM}_{\rm{H_{2}},[1,4]}$ | $\left<R_{[1,4]}\right>$ | $\rm{FWHM}_{\rm{H_{2}},[1,7]}$ | $\left<R_{[1,7]}\right>$ |  $\rm{FWHM}_{\rm{CO},(14-4)}$aaIndividual rotational states within the UV-CO bands are too close to each other in wavelength space to spectrally resolve with HST-COS. Instead of measuring FWHMs directly, we fit multi-Gaussian profiles to states with $J^{\prime\prime}>10$, where the separation between rotational states becomes wider. The FWHMs quoted here correspond to the best-fit parameters from the multi-Gaussian models. |  $\left<R_{CO}\right>$bbTargets with missing values of $\rm{FWHM}_{\rm{CO},(14-4)}$ and $\left<R_{CO}\right>$ did not show UV-CO bands with distinct separations between sets of three high $J^{\prime\prime}$ states ($J^{\prime\prime}>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.
---|---|---|---|---|---|---
| [km s-1] | [AU] | [km s-1] | [AU] | [km s-1] | [AU]
AA Tau (2011) | $61\pm 9$ | $0.7\pm 0.2$ | $60\pm 3$ | $0.75\pm 0.07$ | $13\pm 3$ | $15\pm 6$
AA Tau (2013) | $43\pm 6$ | $1.5\pm 0.4$ | $33\pm 7$ | $2\pm 1$ | $11\pm 1$ | $23\pm 6$
CS Cha | $30\pm 11$ | $3\pm 2$ | $20\pm 6$ | $7\pm 4$ |  $\cdots$ |  $\cdots$
CW Tau | $59\pm 3$ | $0.51\pm 0.06$ | $53\pm 9$ | $0.6\pm 0.2$ |  $\cdots$ |  $\cdots$
DF Tau | $63\pm 3$ | $0.17\pm 0.02$ | $64\pm 6$ | $0.17\pm 0.03$ | $6.4\pm 0.6$ | $16\pm 3$
DM Tau | $32\pm 5$ | $0.6\pm 0.2$ | $27\pm 4$ | $0.8\pm 0.3$ | $3.3\pm 0.8$ | $53\pm 26$
LkCa 15 | $59\pm 9$ | $0.5\pm 0.1$ | $50\pm 1$ | $0.68\pm 0.04$ |  $\cdots$ |  $\cdots$
RECX-11 | $56\pm 4$ | $0.8\pm 0.1$ | $52\pm 2$ | $0.92\pm 0.08$ | $14\pm 1$ | $13\pm 2$
RECX-15 (2010) | $35\pm 4$ | $0.9\pm 0.2$ | $41\pm 1$ | $0.63\pm 0.02$ | $6\pm 3$ | $26\pm 21$
RECX-15 (2013) | $47\pm 2$ | $0.48\pm 0.04$ | $38\pm 4$ | $0.7\pm 0.1$ | $14\pm 4$ | $6\pm 3$
RY Lupi | $51\pm 10$ | $1.9\pm 0.8$ | $48\pm 5$ | $2.1\pm 0.4$ |  $\cdots$ |  $\cdots$
T Cha | $58\pm 2$ | $1.4\pm 0.1$ | $54\pm 10$ | $1.7\pm 0.6$ |  $\cdots$ |  $\cdots$
UX Tau A | $31\pm 5$ | $1.5\pm 0.5$ | $30\pm 3$ | $1.7\pm 0.4$ | $7.7\pm 0.5$ | $25\pm 3$
V4046 Sgr | $51\pm 3$ | $0.63\pm 0.07$ | $47\pm 1$ | $0.74\pm 0.03$ | $13\pm 3$ | $6\pm 3$
Median | 51 km s-1 | 0.75 AU | 48 km s-1 | 0.78 AU | 9 km s-1 | 20 AU
Figure 3: Average gas emission radii for 12CO 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^{\prime\prime}>10$. In all cases, this emission
comes from inside or very close to the CO ice line (open circles; estimated
from Long et al. 2018). Dashed black lines connect the UV-CO and CO ice line
radii for the same target. By contrast, IR-CO emission (yellow circles; Salyk
et al. 2011b; Brown et al. 2013) appears to trace a population of gas at
radial locations closer to the UV-H2 (France et al., 2012).

We also find that the UV-CO emission in the 2013 spectrum of AA Tau may
originate from larger radii $\left(R_{\rm{UV-CO}}=23\pm 6\text{ AU}\right)$
than the same feature in the 2011 spectrum of AA Tau $\left(R_{\rm{UV-
CO}}=15\pm 6\text{ AU}\right)$. These empirical $R_{\rm{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-H2 emission lines from the two sets of
spectra, with the 2011 spectrum showing an empirical UV-H2 emitting radius of
$R_{\rm{UV-H_{2}}}=0.7\pm 0.2$ and the 2013 observations showing $R_{\rm{UV-
H_{2}}}=1.5\pm 0.4$.

Changes between the 2011 and 2013 AA Tau UV-H2 spectra were attributed to
extra absorbing material that moved into the line-of-sight after 2011
(Schneider et al., 2015; 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 $\left(R_{\rm{UV-CO}}\right)$
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.

##  5 2-D Radiative Transfer Models of UV-H2 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-H2
and UV-CO spectra. We build on the models developed by Hoadley et al. (2015),
which create UV-H2 emission line profiles by distributing Ly$\alpha$ photons
to a population of hot gas in a thin surface layer of the disk. This version
efficiently accounts for UV-H2 lines from four different progressions, pumped
by photons with four different wavelengths along the red side of the
Ly$\alpha$ 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 $\left(v^{\prime}-v^{\prime\prime}\right)=(14-3)$ and
$\left(v^{\prime}-v^{\prime\prime}\right)=(14-4)$ bands of UV-CO each consist
of lines pumped by $>$60 different Ly$\alpha$ wavelengths
$\left(\lambda=1214.2-1218.8\text{ \AA}\right)$. Rather than adding the UV-CO
radiative transfer properties to the existing code for UV-H2, the models were
adapted into an object-oriented framework that can accommodate multiple
molecular species. A comparison of UV-H2 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).

###  5.1 Description of the Models

Each radiative transfer model is built on the same 2-D disk structure
described by Hoadley et al. (2015), which starts with a radial temperature
gradient and pressure scale height. The temperature gradient is normalized to
$T_{1AU}$ and decays with power law exponent $q$:

$T\left(r\right)=T_{1AU}\left(\frac{r}{1\,\text{AU}}\right)^{-q}$ (2)

The pressure scale height follows as:

$H_{p}=\sqrt{\frac{kT\left(r\right)}{\mu
m_{H}}\times\frac{r^{3}}{GM_{\ast}}},$ (3)

where $k$ and $G$ are the Boltzmann and gravitational constants, $\mu=2.33$ is
the mean molecular weight of the gas, $m_{H}$ is the mass of a single hydrogen
atom, and $M_{\ast}$ is the stellar mass. The 2-D mass density distribution of
gas is then calculated as

$\rho\left(r,z\right)=\frac{\Sigma\left(r\right)}{\sqrt{2\pi}H_{p}}\times\exp\left[-0.5\left(\frac{z}{H_{p}}\right)^{2}\right]$
(4)

with radial surface density distribution

$\Sigma\left(r\right)=\Sigma_{c}\left(\frac{r}{r_{c}}\right)^{-\gamma}\exp\left[-\left(\frac{r}{r_{c}}\right)^{2-\gamma}\right].$
(5)

The characteristic radius $\left(r_{c}\right)$ is the location in the disk
where the density distribution begins to fall off exponentially. A power-law
with coefficient $\gamma$, which describes the radial surface density
distribution $\Sigma\left(r\right)$ for a disk with kinematic viscosity
$\nu\propto r^{-\gamma}$ (Lynden-Bell & Pringle, 1974), is used to map the
distribution interior to $r_{c}$.

Number densities of both H2 and CO in ground states with vibrational level $v$
and rotational level $J$ $\left(\left[v,J\right]\right)$ are calculated from
the mass density distribution under LTE conditions:

$n_{[v,J]}\left(r,z\right)=X_{H_{2},CO}\times\frac{\rho\left(r,z\right)}{\mu
m_{H}}\times
g_{[v,J]}\times\frac{\exp{\left[-E_{[}v,J]/kT\left(r\right)\right]}}{Z_{[v,J]}\left(T\right)},$
(6)

where $X_{H_{2},CO}$ represents the fraction of the total gas density
consisting of H2 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

$\tau_{\lambda}\left(r,z\right)=\sum_{z}^{z-H_{p}}z\sigma
n_{[v,J]}\left(r,z\right),$ (7)

where $\sigma$ is the cross section for absorbing a Ly$\alpha$ 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 $\tau$ by the fractional optical depth expected for that transition at
the temperature $\left(T\right)$ and column density $\left(N\right)$
corresponding to its grid point (McJunkin et al., 2014; Liu & Dalgarno, 1996;
Wolven et al., 1997), such that

$\tau_{eff}\left(r,z\right)=\tau_{\lambda}\left(r,z\right)\times\frac{\tau_{\lambda}\left(r,z\right)}{\tau_{all}\left(T,N\right)}.$
(8)

The total flux for each transition out of the Ly$\alpha$-pumped upper level is
then calculated as

$F_{\lambda}\left(r,z\right)=\eta
F_{\ast,Ly\alpha}\left(\frac{R_{\ast}}{r}\right)^{2}\left(\frac{r\times\cos
i_{disk}}{s\left(r,z\right)}\right)^{2}B_{mn}\sum_{\tau_{eff}}\left(1-\text{e}^{-\tau_{eff}}\right),$
(9)

where $\eta$ is a geometric filling factor describing the fraction of gas
exposed to Ly$\alpha$ photons (Herczeg et al., 2004), $F_{\ast,\rm{Ly}\alpha}$
is the relevant Ly$\alpha$ pumping flux at grid point $\left(r,z\right)$,
$s\left(r,z\right)$ is the sightline from the observer to a gas parcel in the
disk, and $B_{mn}$ describes the probability of transitioning from
$\left[v^{\prime},J^{\prime}\right]\rightarrow\left[v^{\prime\prime},J^{\prime\prime}\right]$.
We do not fit for $\eta$ directly in this analysis; rather, the value is
folded into the total Ly$\alpha$ 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-H2
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-H2 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 $\left(J^{\prime\prime}=0\right)$ to radial velocities of
zero. The impact of this shift on the modeling results is discussed further in
Section 5.4.

###  5.2 Reconstructed Ly$\alpha$ Profiles

Although the Ly$\alpha$ 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 high-velocity 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-H2 lines (Herczeg et
al., 2004; Schindhelm et al., 2012b). Progression fluxes, calculated by
integrating over all UV-H2 features from the same upper level $i$
$\left(F_{H_{2},i}\right)$, are treated as “y” data points. Ly$\alpha$ pumping
wavelengths act as corresponding “x” values, such that
$\left(x_{i},y_{i}\right)=\left(\lambda_{Ly\alpha,i},F_{H_{2},i}\right)$. The
data points are then fitted with a model profile, assumed to be an outflow-
absorbed Gaussian emission line ($I_{\rm{Ly}\alpha}$; Schindhelm et al.
2012b):

$I_{\rm{abs}}\left(\lambda\right)=I_{\rm{Ly}\alpha}\left(\lambda\right)\times\exp^{-\tau_{\rm{out}}}\\\
$ (10)

with $\tau_{\rm{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$\alpha$-absorbing H2 with temperature $T_{H_{2}}$ and column
density $N_{H_{2}}$

$W_{\lambda_{\rm{pump}}}=\int_{\lambda_{\rm{pump}}}\left(1-\exp^{-\tau_{\lambda_{\rm{pump}}}\left[\rm{T_{H_{2}}},\rm{N_{H_{2}}}\right]}\right).$
(11)

We applied this Ly$\alpha$ 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$\alpha$
fluxes accurate to within $\sim$20% (Schindhelm et al., 2012b).

###  5.3 UV-H2 Model Fitting Procedure

UV-H2 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 $\left(2\,H_{p}<z/r<7\,H_{p}\right)$

  * •

temperature at 1 AU $\left(500\,\rm{K}<T_{1\,\rm{AU}}<5000\,\rm{K}\right)$

  * •

power law exponent for the temperature distribution $\left(-2.5<q<2.5\right)$

  * •

power law exponent for the surface density distribution
$\left(0<\gamma<2\right)$

  * •

radius where the surface density distribution turns over from a power law to
an exponential decline $\left(0.1\,\rm{AU}<r_{\rm{char}}<20\,\rm{AU}\right)$

  * •

total mass of H2 $\left(1\times
10^{-5}\,M_{\odot}<M_{\rm{H_{2}}}<0.1\,M_{\odot}\right)$

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-H2 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-H2 model for each disk by minimizing the mean square
error ($MSE$):

$MSE=\frac{1}{N}\sum\limits_{i=0}^{N}\left(y_{i}-\tilde{y_{i}}\right)^{2},$
(12)

which describes the goodness-of-fit of the model profiles ($\tilde{y_{i}}$) to
the observed UV-H2 emission lines $\left(y_{i}\right)$. We choose the $MSE$
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-H2
features, since anomalously small uncertainties in the continuum can bias test
statistics (e.g. $\chi^{2}$) toward models that predict non-detections.
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-H2 emission lines are
representative of the underlying radial distributions of flux from hot,
Ly$\alpha$-pumped H2 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 power-law exponents (Hoadley et al., 2015). Furthermore,
several targets in our sample have UV-H2 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 $MSE$ 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_{H_{2}}$, $T_{1\,\text{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-H2
emission from each target. Although we used a different mass parameterization
than Hoadley et al. (2015), the UV-H2 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.

Table 5: Best-Fit Parameters for UV-H2 Radiative Transfer Models  Target | $T_{1\text{ AU}}$ | $q$ | $r_{char}$ |  $r_{in}$aa5% of the total UV-H2 flux is contained within this radius. |  $r_{peak}$bbRadial location of UV-H2 flux peak. |  $r_{out}$cc95% of the total UV-H2 flux is contained within this radius.
---|---|---|---|---|---|---
| [K] |  | [AU] | [AU] | [AU] | [AU]
AA Tau (2011) | $2700^{+900}_{-1200}$ | $0.3^{+0.6}_{-0.26}$ | $14\pm 6$ | $0.2^{+0.2}_{-0.1}$ | $0.34^{+0.01}_{-0.1}$ | $2^{+5}_{-1}$
AA Tau (2013) | $2400^{+1200}_{-900}$ | $0.45^{+0.55}_{-0.4}$ | $17^{+3}_{-6}$ | $0.5\pm 0.4$ | $0.7^{+0.3}_{-0.06}$ | $4^{+8}_{-2}$
CS Cha | $3600\pm 900$ | $0.15^{+0.7}_{-0.04}$ | $19.3^{+0.6}_{-8}$ | $0.2^{+0.2}_{-0.15}$ | $0.4\pm 0.1$ | $3^{+6}_{-1}$
CW Tau | $3200^{+800}_{-1200}$ | $0.6^{+0.4}_{-0.5}$ | $7.3^{+5}_{-7}$ | $0.1^{+0.1}_{-0.06}$ | $0.26^{+0.03}_{-0.1}$ | $2^{+3}_{-1}$
DF Tau | $1300^{+1200}_{-700}$ | $-0.3^{+0.5}_{-1.7}$ | $17^{+3}_{-5}$ | $0.04^{+0.02}_{-0}$ | $0.04^{+0.03}_{-0}$ | $2^{+5}_{-1}$
DM Tau | $3300\pm 800$ | $0.7^{+0.5}_{-1}$ | $10\pm 6$ | $0.19\pm 0.08$ | $0.26\pm 0.05$ | $3^{+2}_{-1}$
LkCa 15 | $2400^{+800}_{-1300}$ | $0.14^{+1.1}_{-0.1}$ | $15\pm 5$ | $0.2^{+0.2}_{-0.15}$ | $0.3^{+0.3}_{-0.07}$ | $2^{+8}_{-1}$
RECX-11 | $2400^{+600}_{-800}$ | $0.3^{+0.5}_{-1.1}$ | $6\pm 5$ | $0.1^{+0.2}_{-0.05}$ | $0.25^{+0.02}_{-0.09}$ | $5^{+6}_{-3}$
RECX-15 (2010) | $4100^{+800}_{-1100}$ | $0.9^{+0.4}_{-0.7}$ | $15\pm 5$ | $0.06^{+0.1}_{-0.02}$ | $0.2^{+0.02}_{-0.1}$ | $5^{+6}_{-3}$
RECX-15 (2013) | $2100\pm 1200$ | $0.2^{+0.7}_{-1}$ | $5^{+5}_{-4.8}$ | $0.06^{+0.1}_{-0.02}$ | $0.15^{+0.05}_{-0.07}$ | $4^{+5}_{-2}$
RY Lupi | $1900^{+1200}_{-800}$ | $0.3^{+0.5}_{-0.7}$ | $11\pm 6$ | $0.3\pm 0.2$ | $0.5^{+0.05}_{-0.15}$ | $11^{+6}_{-5}$
T Cha | $4000^{+900}_{-1300}$ | $-0.4^{+1.8}_{-1}$ | $12^{+7}_{-5}$ | $0.7^{+0.2}_{-0.6}$ | $1^{+0}_{-0.05}$ | $3^{+4}_{-1}$
UX Tau AddBi-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. | $3200^{+700}_{-1400}$ | $-0.7^{+1.8}_{-0.6}$ | $15\pm 5$ | $0.5^{+0.1}_{-0.2}$ | $0.8^{+0}_{-0.2}$ | $3\pm 1$
V4046 Sgr | $2200^{+1000}_{-700}$ | $0.5^{+0.8}_{-0.3}$ | $17^{+3}_{-7}$ | $0.3^{+0.1}_{-0.2}$ | $0.4^{+0}_{-0.05}$ | $2^{+3}_{-1}$
Median | 2550 K | 0.3 | 14.5 AU | 0.195 AU | 0.3 AU | 3.3 AU

###  5.4 UV-CO Model Fitting Procedure

The main difference in the treatment of UV-H2 versus UV-CO is the temperature
distribution within the fluorescent gas. Herczeg et al. (2004) found that the
UV-H2 features from TW Hya were consistent with emission from a gas layer with
$T\sim 2000-3000$ K and $\log N\left(\rm{H_{2}}\right)=18-20$, so the models
presented here allow the temperature of the H2 distribution to reach as high
as the 5000 K dissociation limit (Hoadley et al., 2015). Assuming a CO/H2
ratio of $10^{-4}$ (France et al., 2014a), the corresponding column densities
of CO would be $\log N\left(\rm{CO}\right)\sim 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 $\left(\log N\left(\rm{H}\right)\sim 21\right)$. The gas temperature is
then quickly reduced to $\sim$1000 K in denser surface layers of the disk
(Ádámkovics et al., 2014). This implies that hotter layers, where the UV-H2
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
$\mu$m CO rovibrational emission in surface layers of the inner disk. The
rotational temperatures of these features never exceed $T\sim 1800$ K and are
typically much cooler than the $T\sim 1500$ K threshold required for
Ly$\alpha$ pumping of H2 (Salyk et al., 2011b). Since the IR-CO features
originate much closer to the star than the UV-CO $\left(R_{IR-CO}\sim
0.1-5\text{ AU; see Figure \ref{Gaussian_Radii_Comp},
{\cite[citet]{\@@bibref{Authors
Phrase1YearPhrase2}{France2012}{\@@citephrase{}}{\@@citephrase{}}}}}\right)$,
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-H2, 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-H2 emission lines trace radially separated regions of
the gas disk, we expect the Ly$\alpha$ 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-H2 and UV-CO emitting regions. We restrict the H I column
densities to range between the values that best fit the UV-H2 progression
fluxes (Schindhelm et al., 2012b) and those estimated for the total ISM
absorption (McJunkin et al., 2016). A new Ly$\alpha$ absorption feature is
then superimposed on the intrinsic Gaussian emission line generated by the H2
reconstruction procedure.

In order for Ly$\alpha$ 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-H2. The
variable parameters for UV-CO are then:

  * •

height of the emitting layer $\left(7\,H_{p}<z/r<11\,H_{p}\right)$

  * •

temperature at 1 AU $\left(20\,\rm{K}<T_{1\,\rm{AU}}<1800\,\rm{K}\right)$

  * •

power law exponent for the temperature distribution $\left(-2.5<q<2.5\right)$

  * •

power law exponent for the surface density distribution
$\left(0<\gamma<2\right)$

  * •

radius where the surface density distribution turns over from a power law to
an exponential decline $\left(0.1<r_{\rm{char}}<100\,\rm{AU}\right)$

  * •

total mass of CO
$\left(10^{-20}\,M_{\odot}<M_{\rm{CO}}<10^{-5}\,M_{\odot}\right)$

  * •

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 $MSE$ 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-H2 and
UV-CO emission lines from each individual target are provided in Appendix B,
along with the corresponding radial flux distributions.

Table 6: Best-Fit Parameters for UV-CO Radiative Transfer Models  Target | $T_{10\text{ AU}}$ | $q$ | $r_{char}$ |  $r_{in}$aa5% of the total UV-CO flux is contained within this radius. |  $r_{peak}$bbRadial location of UV-CO flux peak. |  $r_{out}$cc95% of the total UV-CO flux is contained within this radius.
---|---|---|---|---|---|---
| [K] |  | [AU] | [AU] | [AU] | [AU]
AA Tau (2011) | $950^{+90}_{-400}$ | $-0.8^{+0.2}_{-0.3}$ | $25^{+20}_{-14}$ | $0.16_{-0.12}^{+0.84}$ | $1.8^{+0.5}_{-0}$ | $11^{+4}_{-3}$
AA Tau (2013) | $880\pm 110$ | $-0.7^{+0.35}_{-0.2}$ | $47^{+50}_{-22}$ | $0.1^{+0.2}_{-0.06}$ | $1.6^{+0.7}_{-1.0}$ | $10^{+2}_{-4}$
CS Cha | $1300\pm 300$ | $-0.6^{+0.4}_{-0.3}$ | $50^{+10}_{-30}$ | $0.04^{+0.2}_{-0.01}$ | $0.4\pm 0.1$ | $9^{+30}_{-4}$
CW Tau | $1000\pm 100$ | $0.15^{+0.03}_{-0.05}$ | $99^{+1}_{-5}$ | $0.04^{+0.06}_{-0}$ | $0.04^{+0.06}_{-0}$ | $7\pm 2$
DF Tau | $1100^{+150}_{-800}$ | $-0.038\pm 0.1$ | $38^{+50}_{-10}$ | $0.04^{+0.02}_{-0}$ | $0.04^{+0.02}_{-0}$ | $1.4\pm 0.4$
DM Tau | $1700^{+100}_{-600}$ | $-1.0^{+0.5}_{-0.8}$ | $40^{+20}_{-30}$ | $0.07^{+1}_{-0.03}$ | $1^{+0.6}_{-0.9}$ | $10^{+14}_{-2}$
LkCa 15 | $400^{+400}_{-350}$ | $-0.6^{+3}_{-0.4}$ | $69^{+30}_{-20}$ | $0.2^{+6}_{-0.15}$ | $2^{+5}_{-1}$ | $14^{+20}_{-5}$
RECX-11 | $1700\pm 100$ | $-1.1^{+0.7}_{-0.4}$ | $9^{+10}_{-4}$ | $0.2^{+0.2}_{-0.15}$ | $1.4\pm 1$ | $7\pm 1$
RECX-15 (2010) | $560^{+500}_{-200}$ | $-0.46\pm 0.1$ | $62^{+30}_{-10}$ | $0.05\pm 0.02$ | $0.65^{+0.3}_{-0.05}$ | $12^{+6}_{-5}$
RECX-15 (2013) | $540^{+500}_{-200}$ | $-0.4^{+0.2}_{-0.1}$ | $67^{+30}_{-10}$ | $0.043^{+0.006}_{0.004}$ | $0.39^{+0.4}_{-0.08}$ | $9^{+3}_{-5}$
RY Lupi | $1600^{+200}_{-500}$ | $-0.7^{+0.65}_{-0.3}$ | $7^{+35}_{-6}$ | $0.04^{+0.2}_{-0}$ | $0.06^{+0.2}_{-0.02}$ | $6^{+6}_{-1}$
T Cha | $1100^{+500}_{-700}$ | $-0.96^{+0.4}_{-1.3}$ | $63\pm 30$ | $0.04^{+0.3}_{-0}$ | $0.1^{+0.3}_{-0.06}$ | $6^{+8}_{-2}$
UX Tau A | $1600^{+200}_{-500}$ | $-1.2^{+1}_{-0.6}$ | $24\pm 20$ | $0.8^{+8}_{-0.7}$ | $3^{+6}_{-1}$ | $12^{+50}_{-2}$
V4046 Sgr | $1550\pm 80$ | $-0.9\pm 0.1$ | $34^{+5}_{-9}$ | $0.08^{+0.6}_{-0.04}$ | $0.75^{+0.25}_{-0.06}$ | $7^{+2}_{-1}$
Median | 1100 K | -0.7 | 44 AU | 0.06 AU | 0.7 AU | 9 AU

We find that the median temperature of fluorescent CO at 10 AU is $\sim$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}\sim 9$ AU, demonstrating that the emitting regions extend
roughly three times further than the radial distributions of UV-H2
$\left(r_{out}\sim 3\,\rm{AU}\right)$. However, $r_{out}$ for the UV-CO is
only about half the size of the empirical radius presented in Table 4
$\left(\left<R_{CO}\right>\sim 20\,\rm{AU}\right)$.

The number densities of intervening H I between the UV-H2 and UV-CO emitting
regions were also calculated for each target (see Table 7). Values range
between $n_{\rm{HI}}\sim 0.1-3.4\times 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-H2 and UV-CO flux. Any systematic uncertainties in that
analysis will then propagate through the number density calculations.

Table 7: H I Column Densities at Various Disk Checkpoints  Target |  $N\left(\rm{H\,I}\right)$, UV-H2aaH I column densities from Schindhelm et al. 2012b; Arulanantham et al. 2020, corresponding to absorption between the accretion shock and UV-H2 emitting region. |  $N\left(\rm{H\,I}\right)$, UV-CObbH I column densities measured in this work, corresponding to absorption between the UV-H2 and UV-CO emitting regions. | $\Delta N\left(\rm{H\,I}\right)\left(\rm{CO}-\rm{H}_{2}\right)$ |  $N\left(\rm{H\,I}\right)$, totalccH I column densities from McJunkin et al. 2014; These values are dominated by ISM absorption. |  $R_{\rm{CO}}-R_{\rm{H}_{2}}$eeDifference in outer radii of UV-CO and UV-H2 model flux distributions, which contain 95% of the fluorescent emission. DF Tau and RY Lupi both have negative values, implying that the UV-H2 emitting region is further from the star than the UV-CO. However, this is likely caused by underfitting the high-velocity UV-H2 line wings, which would drop emission at radii close to the star from the model flux distributions. |  $n_{\rm{HI}}$ffH I number densities were calculated using $R_{\rm{CO}}-R_{\rm{H}_{2}}$. Any systematic uncertainties in the estimates for those radii will then propagate through the number density calculations.
---|---|---|---|---|---|---
| [dex] | [dex] | [dex] | [dex] | [AU] | [106 cm-3]
AA Tau (2011) | $19.05\pm 0.15$ | $20.37^{+0.36}_{-0.46}$ | 1.3 | $20.73^{+0.03}_{-0.13}$ | 9 | 1.7
AA Tau (2013) | $18.29\pm 0.3$ | $19.93^{+0.8}_{-0.4}$ | 1.6 |  $\cdots$ | 6 | 0.9
CS Cha | $19.05\pm 0.15$ | $19.4\pm 0.4$ | 0.35 | $20.4^{+0.03}_{-0.15}$ | 6 | 0.15
CW Tau | $19.2\pm 0.8$ | $20.29^{+0.08}_{-0.07}$ | 1.09 |  $\cdots$ | 5 | 2.4
DF Tau | $19.5\pm 0.2$ | $19.57^{+0.34}_{-1}$ | 0.07 | $20.93\pm 0.03$ | -0.6 | 0.62
DM Tau | $18.59\pm 0.03$ | $19.4^{+0.5}_{-0.6}$ | 0.8 | $20.88^{+0.05}_{-0.09}$ | 7 | 0.20
LkCa 15 | $18.9\pm 0.3$ | $19.69^{+0.4}_{-0.7}$ | 0.8 | $20.68^{+0.15}_{-0.09}$ | 12 | 0.23
RECX-11 | $18.89\pm 0.4$ | $19.18^{+0.38}_{-0.15}$ | 0.3 | $19.7^{+0.03}_{-0.2}$ | 2 | 0.25
RECX-15 (2010) | $18.98\pm 0.22$ | $19.25^{+0.12}_{-0.25}$ | 0.27 | $19.58^{+0.05}_{-0.18}$ | 7 | 0.08
RECX-15 (2013) | $18.98\pm 0.22$ | $19.23^{+0.13}_{-0.2}$ | 0.25 |  $\cdots$ | 5 | 0.1
RY Lupi | $19.1^{+0.4}_{-0.2}$ | $19.1^{+0.3}_{-0.1}$ | 0 | 20.34ddH I column density estimated from France et al. 2017b | -5 | 0
T Cha | $18.9\pm 0.6$ | $20.20^{+0.7}_{-0.5}$ | 1.3 |  $\cdots$ | 3 | 3.4
UX Tau A | $18.26\pm 0.12$ | $19.4^{+0.5}_{-0.7}$ | 1.1 | $20.9^{+0.05}_{-0.17}$ | 9 | 0.17
V4046 Sgr | $18.81\pm 0.22$ | $19.83^{+0.01}_{-0.05}$ | 1.02 | $19.85^{+0.03}_{-0.28}$ | 5 | 0.82
Median | 18.94 | 19.49 | 0.8 | 20.54 | 6 AU |  $\mathbf{2.4\times 10^{5}}$ cm-3

Although the 2-D radiative transfer approach provides good fits to the UV-H2
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 $\lambda\sim 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 with
$J^{\prime\prime}>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$\alpha$
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$\alpha$
pumping wavelengths. The model generated with an outflow absorbed Ly$\alpha$
profile underpredicts the data from transitions with $J^{\prime\prime}\sim
8-12$, which are pumped at Ly$\alpha$ wavelengths between $\lambda\sim
1214.5-1215.7$ Å. We also find that the models underpredict the data at longer
wavelength transitions pumped by redder Ly$\alpha$ 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^{\prime\prime}$ features used to derive the empirical radii, the
underpredictions are likely responsible for the discrepancies between
$\left<R_{CO}\right>$ from the Gaussian fitting procedure (Table 4) and
$r_{peak}$ from the radiative transfer models (Table 6). Additional sources of
Ly$\alpha$ pumping photons, both from line center and red-shifted wavelengths,
must be included to more accurately reproduce the UV-CO emission from large
radii $\left(r>15\,\text{au}\right)$.

This evidence for a missing source of red-shifted Ly$\alpha$ at the disk
surface is consistent with previous analyses of UV-H2, since the spectra can
show significant emission from progressions excited by pumping photons up to
$\sim$1000 km s-1 redward of line center (e.g. DF Tau; Herczeg et al. 2006).
Furthermore, even the best-fit UV-H2 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$\alpha$ 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$\alpha$ 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$\alpha$ 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$\alpha$ emission at the disk
surface, the radial location(s) that receive the strongest irradiation from
Ly$\alpha$ 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$\alpha$ pumping fluxes, make it difficult to converge on accurate
radial distributions of UV-CO flux as derived for the UV-H2 emission (Hoadley
et al., 2015). In the following sections, we instead further analyze the
observed Ly$\alpha$ 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$\alpha$
profiles from T Tauri stars, which will require more observations of other
Ly$\alpha$ sensitive molecules (e.g., HCN, H2O, H2CO; Bergin et al. 2003;
Bethell & Bergin 2011) in disks with available UV spectra.

Figure 4: Best-fit UV-CO models for UX Tau, with residuals color-coded from
blue to red based on the appropriate Ly$\alpha$ pumping wavelength. The model
acquired with a Ly$\alpha$ 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$\alpha$ 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$\alpha$
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.

##  6 Discussion

As described in the previous section, the 2-D radiative transfer models of
UV-H2 and UV-CO depend strongly on Ly$\alpha$ 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$\alpha$
profile. Since interstellar absorption and geocoronal emission prevent us from
directly detecting the Ly$\alpha$ 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$\alpha$ transition ($\lambda_{central}=1222$ Å, instead of 1291 Å).

###  6.1 Ly$\alpha$ Irradiation Drives the UV-CO Band Shapes

Figure 5 compares the best-fit model gas temperatures from Table 6 to the
ratios of observed UV-CO emission from high $J^{\prime\prime}$ states relative
to the bandhead
$\left(F_{J^{\prime\prime}>10}/F_{J^{\prime\prime}=0-10}\right)$. The
relationship between the model gas temperatures at 10 AU and the observed UV-
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^{\prime\prime}$ states.

Since the UV-CO emitting region extends to larger radii than the UV-H2, it is
possible (or even likely) that Ly$\alpha$ profiles reconstructed from the
UV-H2 emission inside $r<2$ AU are not representative of the radiation field
at $r>15$ AU. To gain a better understanding of the Ly$\alpha$ profiles seen
by the UV-CO, we compare the ratios of observed Ly$\alpha$ emission from the
red and blue wings of the profiles
$\left(F_{1217-1220\,\rm{\AA}}/F_{1210-1214.5\,\rm{\AA}}\right)$ to the UV-CO
band shapes in Figure 6. Unlike the gas temperatures, the ratios of red wing
to blue wing observed Ly$\alpha$ emission show a strong positive correlation
with the ratios of UV-CO flux in high $J^{\prime\prime}$ 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$\alpha$ profile
instead of the temperature distribution within the emitting region.

With this in mind, Figure 7 also compares the ratios of UV-H2 fluxes from
progressions pumped by red $\left(\left[1,4\right],\left[1,7\right]\right)$
and blue Ly$\alpha$ photons $\left(\left[4,4\right]\right)$ (France et al.,
2012) to the ratios of red to blue Ly$\alpha$ fluxes. Once again, we find a
strong correlation between the flux ratios (Spearman rank coefficient $=0.89$;
$p=5\times 10^{-4}$), demonstrating that disks with more blue-wing Ly$\alpha$
emission also show stronger UV-H2 fluxes in the $\left[4,4\right]$ 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$\alpha$ fluxes in the following sections.

Figure 5: UV-CO model temperatures at 10 AU versus observed flux ratios of
$J^{\prime\prime}>10$ to $J^{\prime\prime}=0-10$ UV-CO emission. The
significant anti-correlation demonstrates that higher gas temperatures do not
necessarily lead to increased high $J^{\prime\prime}$ emission. CS Cha is
denoted with a cross (see also Figures 6, 7, and 8), because it does not have
high-resolution [O I] $\lambda$6300 spectra that reveal the presence or
absence of a high-velocity component. Uncertainties on the measured flux
ratios are smaller than the marker sizes, so vertical error bars are not
shown. An interactive version of this figure is provided online, providing
disk and stellar properties for each target. Figure 6: Flux ratios of observed
red wing (1217-1220 Å) to blue wing (1210-1214.5 Å) Ly$\alpha$ emission versus
ratios of UV-CO $J^{\prime\prime}>10$ to $J^{\prime\prime}=0-10$ emission. The
log-linear relationship implies that the shape of the Ly$\alpha$ profile is
the dominant factor in the $J^{\prime\prime}$ state distribution of UV-CO flux
(see also Table 7 of Schindhelm et al. 2012a). An interactive version of this
figure is provided online, providing disk and stellar properties for each
target. Figure 7: Flux ratios of observed red wing to blue wing Ly$\alpha$
emission versus flux ratios of red
$\left(\left[1,4\right],\left[1,7\right]\right)$ to blue $\left[4,4\right]$
UV-H2 progressions demonstrate the same trend seen between Ly$\alpha$ and UV-
CO in Figure 6. An interactive version of this figure is provided online,
providing disk and stellar properties for each target. Figure 8: Flux ratios
of observed red wing to blue wing Ly$\alpha$ emission versus
$n\left(13-31\right)$ for disks with dust cavities, indicating that disks with
more blue Ly$\alpha$ 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.

###  6.2 Three Categories of UV-CO Band Shapes and Ly$\alpha$ Morphologies

Figure 9 provides the observed Ly$\alpha$ wings, UV-CO
$\left(v^{\prime}-v^{\prime\prime}\right)=\left(14-3\right)$ bands, UV-H2
progression fluxes, infrared spectral indices, and [O I] $\lambda$6300 line
shapes for all targets in our sample. We order the targets based on where they
fall in 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^{\prime\prime}$ states, systems showing roughly similar bandhead and high
$J^{\prime\prime}$ fluxes, and targets with more high $J^{\prime\prime}$ than
bandhead emission. As depicted in Table 8, targets that fall within a specific
grouping also have similar observed Ly$\alpha$ profiles and infrared spectral
indices. For example, disks with more bandhead than high $J^{\prime\prime}$
UV-CO emission also have stronger blue wing Ly$\alpha$ flux and larger SED
slopes between 13 and 31 $\mu$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$\alpha$
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., [O I], [Ne II]; Banzatti et al. 2019;
Pascucci et al. 2020). Figure 10 illustrates this sequence, demonstrating how
the Ly$\alpha$ profile might shift as the disk changes.

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$\alpha$ photons. This demonstrates that the [4,4] progression fluxes become relatively stronger as the blue Ly$\alpha$ emission increases. CW Tau is not included, since its Ly$\alpha$ wings were not observed with HST-COS. Systems with similar ratios of observed red to blue Ly$\alpha$ emission have similar UV-CO band shapes, UV-H2 fluxes from progressions pumped by blue ([4,4]) and red ([1,7], [1,4]) Ly$\alpha$ photons, infrared spectral indices, and [O I] $\lambda$6300 emission line components. Outflow-dominated targets are enclosed in the purple, dotted box, intermediate targets in the orange, solid box, and infall-dominated targets in the green, dashed box. The UV-H2 [4,4] flux from RY Lup is unusually high, as driven by the (4-9)P(5) 1526.55 Å feature (Arulanantham et al., 2018). Table 8: Group Properties Based on Empirical Tracers Discussed in the Text  |  $<R_{H_{2}}>$aaGiven the larger uncertainties in the 2-D radiative transfer models, we report the median empirical radii for both the UV-H2 and UV-CO features from each group. Radii are calculated from the Gaussian emission line widths, under the assumption that all emitting material is in Keplerian rotation. | $<R_{CO}>$ | UV-CO | Ly$\alpha$ | Dust Evolution | Intervening H IddNumber densities of intervening H I were estimated from the column densities and distances for each target presented in Table 7. A median value was then calculated for each group. The number of targets in any individual group is small and the scatter is $\sim$40%, so we urge the reader to interpret these results with caution. | Targets
---|---|---|---|---|---|---|---
| [AU] | [AU] | $\frac{F\left(J^{\prime\prime}>10\right)}{F\left(J^{\prime\prime}<10\right)}$ |  $\frac{F\left(\text{red}\right)}{F\left(\text{blue}\right)}$bb$F\left(\text{red}\right)/F\left(\text{blue}\right)=F\left(\lambda\sim 1217-1220\,\rm{\AA}\right)/F\left(\lambda\sim 1210-1214.5\,\rm{\AA}\right)$ |  $n\left(13-31\right)$cc$n\left(13-31\right)=\frac{\log\left(\lambda_{31}F_{\lambda_{31}}\right)-\log\left(\lambda_{13}F_{\lambda_{13}}\right)}{\log\left(\lambda_{31}\right)-\log\left(\lambda_{13}\right)}$; Furlan et al. (2009) |  $10^{6}$ cm-3 |
Outflow Dominated |  $\sim$0.7 |  $\sim$16 | $>\sim 1$ | $>\sim 1$ | $<\sim-0.5$ |  $\sim$0.51 | AA Tau, DF Tau, RECX-15
Intermediate |  $\sim$0.65 |  $\sim$13 | $\sim 1$ | $\sim 1$ | $\sim-0.5-0.5$ |  $\sim$0.23 | RECX-11, LkCa 15, RY Lupi
Infall Dominated |  $\sim$1.4 |  $\sim$25 | $<\sim 1$ | $<\sim 1$ | $>\sim 0.5$ |  $\sim$0.20 | CS Cha, DM Tau, UX Tau,
|  |  |  |  |  |  | T Cha, V4046 Sgr
Figure 10: Illustration of Ly$\alpha$ evolution over the disk lifetime. High-
velocity jets are present in targets with outflow dominated profiles, where
outflowing H I absorbs the blue side of the emission line and fluorescent
emission is seen from as close as $r<0.1$ AU. When outflowing mass loss rates
decline, intermediate, symmetric Ly$\alpha$ profiles are observed and the
UV-H2 and UV-CO emitting regions shift to larger radii. At this stage, any
asymmetries in the observed Ly$\alpha$ wings are present at low velocities
that are mostly hidden by geocoronal emission. However, this does not
necessarily represent an intermediate phase of disk evolution. In the absence
of high-velocity outflowing material, infall dominated Ly$\alpha$ 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).

####  6.2.1 Outflow Dominated Profiles

The disks around AA Tau, RECX-15, and DF Tau show the most prominent UV-CO
emission from high $J^{\prime\prime}$ states, with relatively suppressed
emission from the bandhead. All three disks also have observed Ly$\alpha$ line
wings with significantly more red than blue flux. For these systems, the shape
of the Ly$\alpha$ radiation field reaching both the UV-CO emitting region and
the observer appears generally consistent with the model used to reconstruct
the Ly$\alpha$ profile from the UV-H2 features (see Section 4.2).

Ly$\alpha$ 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$\sim 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$\alpha$ line $\left(\lambda\sim
1214.5-1215.7\,\text{\AA}\right)$. The resulting Ly$\alpha$ 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 shows a sequence of knots with [S II] $\lambda$6717 emission
at $v=-200$ km s-1, which are consistent with model predictions of jet
velocities (Cox et al., 2013). Prominent jets are also resolved in optical [O
I] 6300 Å observations of RECX-15 (Woitke et al., 2011) and DF Tau (Uvarova et
al., 2020). At high spectral resolution, signatures of high-velocity jets
($|v|\sim$ few hundred km s-1) originating near the accretion shock are
detected in [O I] 6300 Å spectra from all three targets (see e.g., Banzatti et
al. 2019). The [O I] features also contain low-velocity ($|v|<50$ km s-1)
blue-shifted components, which originate near the base of MHD winds (Banzatti
et al., 2019). This outflowing material must also include H I, which further
absorbs the Ly$\alpha$ photons generated near the accretion shock.

####  6.2.2 Intermediate Profiles

Figure 9 shows that a handful of targets in our sample have ratios of red wing
to blue wing Ly$\alpha$ flux around $\sim$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$\alpha$ 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^{\prime\prime}$ flux ratios.

RY Lupi and LkCa 15, which host large, dust depleted inner cavities ($r\sim
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^{\prime\prime}$ flux. The SED slopes between 13
and 31 $\mu$m for these targets are scattered, with RY Lupi and LkCa 15 having
$n\left(13-31\right)=-0.26$ and $n\left(13-31\right)=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 $\mu$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\times 10^{-8}$ $M_{\odot}$ yr-1
(Ingleby et al., 2011, 2013; France et al., 2017b). This is supported by the
H$\alpha$ emission line from RECX-11, which has a wide blue wing ($v\sim 300$
km s-1) and inverse P Cygni (red) absorption characteristic of free-falling
gas accreting onto the star (Ingleby et al., 2011).

####  6.2.3 Infall Dominated Profiles

Five disks in our sample have stronger integrated UV-CO fluxes in the bandhead
$\left(J^{\prime\prime}=0-10\right)$ than in the high $J^{\prime\prime}$
states $\left(J^{\prime\prime}>10\right)$: 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 $\mu$m
$\left(n\left(13-31\right);\,\text{{\cite[citet]{\@@bibref{Authors
Phrase1YearPhrase2}{Furlan2009}{\@@citephrase{}}{\@@citephrase{}}}}}\right)$.
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 targets111T Cha and V4046 Sgr both display fairly
symmetric blue and red wing Ly$\alpha$ profiles. However, the shapes of the
UV-CO bands observed from these two systems appear more similar to the disks
with stronger blue wing Ly$\alpha$ emission, so we place them in this category
instead. have observed Ly$\alpha$ profiles with stronger blue wing
$\left(\lambda<1214.5\,\rm{\AA}\right)$ than red wing
$\left(\lambda>1217\,\rm{\AA}\right)$ emission. The strong blue Ly$\alpha$
wings contrast with the model used to reconstruct the Ly$\alpha$ radiation
field, which uses outflowing H I to truncate the blue sides of the line
profiles (Herczeg et al., 2004; Schindhelm et al., 2012b). A model with red-
shifted ISM absorption can reproduce the observed line wings by suppressing
the red Ly$\alpha$ emission (see e.g., McJunkin et al. 2014). However, strong
correlations between the Ly$\alpha$ flux ratios, UV-CO band shapes, UV-H2 flux
ratios, and $n\left(13-31\right)$ (see Figure 9) indicate that the red wing
suppression occurs before the Ly$\alpha$ 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 observer 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$\alpha$ 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$\alpha$ 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 larger ratios of
low velocity [Ne II] 12.8 $\mu$m to [O I] 6300 Å luminosities than targets in
the category with blue-suppressed Ly$\alpha$. That work demonstrates that low-
velocity [Ne II] emission strengthens as the $n\left(13-31\right)$ spectral
index increases and the low-velocity [O I] emitting region migrates to larger
radii (see e.g., Simon et al. 2016; Banzatti et al. 2019; Pascucci et al.
2020). Furthermore, none of the targets in this group show high velocity [O I]
emission, implying that the jets and other high velocity outflows have turned
off. 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$\alpha$ 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$\alpha$
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$\alpha$ 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$\alpha$ emission (DM Tau, UX Tau, CS Cha) have absorption
velocities of $-339$, $-260$, and $-197$ km s-1, respectively. The remaining
targets, which all have positive absorption velocities and more red than blue
Ly$\alpha$ emission, are somewhat scattered, with no clear trend between the
ratios of observed Ly$\alpha$ 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$\alpha$ wings may become more clear when more targets with
dominant blue wing emission are observed.

Figure 11: Upper left: Extrapolated H I absorption velocity for a target with
an outflow dominated Ly$\alpha$ profile. Lower left: Same as upper left, for a
target with an infall dominated Ly$\alpha$ profile. These velocities trace the
opacity of intervening H I between the accretion shock and observer, rather
than the motion of the material. Right: Ratio of observed red-wing to blue-
wing Ly$\alpha$ fluxes versus extrapolated H I absorption velocities. The
three targets with negative absorption velocities do not have high velocity [O
I] $\lambda$6300 Å components, implying that emission from inner disk jets and
outflows is not present (see e.g., Banzatti et al. 2019). Targets with
intermediate Ly$\alpha$ flux ratios
$\left(F_{Ly\alpha,red}/F_{Ly\alpha,blue}\sim 1\right)$ also do not show jet
emission.

##  7 Summary and Conclusions

We have presented a joint analysis of Ly$\alpha$ line profiles and fluorescent
emission from H2 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$\alpha$ flux to reach the molecules. We analyze the spectral features both
empirically and with 2-D radiative transfer models, finding that

  1. 1.

The UV-H2 and UV-CO features originate in radially separated regions of the
disk, with UV-H2 observed closer to the star $\left(r\sim 0.8\,\rm{AU}\right)$
than UV-CO $\left(r\sim 20\,\rm{AU}\right)$ in all targets with partially
resolved high $J^{\prime\prime}$ $\left(J^{\prime\prime}>10\right)$ 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), although
we note that slow, uncollimated winds can maintain Keplerian rotation and
behave as an extension of emission from the disks (see e.g., Pontoppidan et
al. 2011).

  2. 2.

The ratios of red to blue Ly$\alpha$ line wing emission fall into three
categories, which correlate significantly with:

    1. (a)

The ratios of high $J^{\prime\prime}$ to bandhead
$\left(J^{\prime\prime}>10/J^{\prime\prime}\leq 10\right)$ UV-CO flux;

    2. (b)

The ratios of UV-H2 flux from progressions pumped by red to blue Ly$\alpha$
wavelengths ([4,4]/([1,7]+[1,4]);

    3. (c)

The infrared spectral indices, $n\left(13-31\right)$, for disks with dust
cavities;

  3. 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$\alpha$ photons. This causes the radiative transfer models to map the
fluorescent gas to smaller radii than the empirical values.

  4. 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$\alpha$
irradiation.

Furthermore, we find that only the three targets with the largest ratios of
red to blue Ly$\alpha$ line wing emission (AA Tau, DF Tau, RECX-15) show high
velocity [O I] $\lambda$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.

##  8 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), the University of Colorado Boulder, and
Colorado State University. The Summit supercomputer is a joint effort of the
University of Colorado Boulder and Colorado State University. This research
made use of Astropy,222http://www.astropy.org a community-developed core
Python package for Astronomy (Astropy Collaboration et al., 2013; The Astropy
Collaboration et al., 2018).

##  Appendix A
Verifying UV-H2 Modeling Results with Hoadley et al. (2015)

UV-H2 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 $\pm 1\sigma$ 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\left(\Theta|y_{i}\right)$:

$p\left(\Theta|y_{i}\right)\propto
p\left(\Theta\right)p\left(y_{i}|\Theta\right),$ (13)

where $\Theta=[z/r,\gamma,T_{1\,\text{AU}},q,r_{\text{char}},M_{H_{2}}]$. We
have assumed uniform prior distributions for the model parameters
$\left(p\left(\Theta\right)\right)$, with boundaries defined in Hoadley et al.
2015 (UV-H2) and Section 4.4 of this work (UV-CO). We also adopt a Gaussian
likelihood function:

$p\left(y_{i}|\Theta\right)=\exp\left[-\sum\limits_{i=0}^{N}\left(\frac{y_{i}-\tilde{y_{i}}}{\sigma_{i}}\right)^{2}\right],$
(14)

where $y_{i}$ and $\tilde{y_{i}}$ represent the observed and model data points
and $\sigma_{i}$ represents the measurement uncertainties in the data.
However, we ignore $\sigma_{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 12: Corner plot created from the MCMC chain for UX Tau. Most targets in
our sample show double-peaked posterior distributions for $T_{1\,\text{AU}}$
and $q$ and degeneracies between $z/r$ and $M_{H_{2}}$ (Hoadley et al., 2015).

Figure 13 compares the best-fit UV-H2 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 $\sim$3200 K, and the secondary value is $\sim$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.

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.

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-of-fit $\left(MSE\right)$ 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-H2 flux are shown in Figure 14. We obtain
final UV-H2 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 14: Radial distributions of UV-H2 flux from the 100 best-fit models to
the UX Tau spectrum (pink). The black dots represent the median flux values at
each radial grid point.

Figure 15 compares the median UV-H2 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 $\sim$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-H2 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
$\left(r<50\,\rm{AU}\right)$, 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 H2 in a thin disk surface
layer are not particularly sensitive to the mass contained in the underlying
bulk gas reservoir.

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.

##  Appendix B
2-D Radiative Transfer Modeling Results for UV-H2 and UV-CO

Figure 16: Best-fit UV-CO 2-D radiative transfer (top left) and Gaussian
models (top right), 100 UV-H2 radiative transfer models with the smallest MSEs
(middle left), comparison of UV-H2-/UV-CO-based and observed Ly$\alpha$
profiles (middle right) and radial distributions of flux for both species
(bottom left) from the disk around AA Tau in 2011. Residuals on the best-fit
UV-CO model are color-coded from blue to red, based on the Ly$\alpha$ pumping
wavelength required to excite the upper level of the transition. The largest
deviations between the models and data correspond to states with upper levels
pumped by central and redder Ly$\alpha$ photons, causing the models to
underpredict the observed flux. This effect is not caused by radial velocity
effects, but rather missing Ly$\alpha$ irradiation, as supported by the model
reproducing the observed spacing between features while still missing the
amplitudes of the emission lines. The full figure set, showing UV-CO and UV-H2
radiative transfer model results, empirical UV-CO models (when available), and
reconstructed and observed Ly$\alpha$ profiles for all targets, is included in
the online journal.

Figure 17: Best-fit UV-CO 2-D radiative transfer (top left) and Gaussian
models (top right), 100 UV-H2 radiative transfer models with the smallest MSEs
(middle left), comparison of UV-H2-/UV-CO-based Ly$\alpha$ profiles (middle
right) and radial distributions of flux for both species (bottom left) from
the disk around AA Tau in 2013. Residuals on the best-fit UV-CO model are
color-coded from blue to red, based on the Ly$\alpha$ 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 $\sim$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 ($\sim$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$\alpha$ profiles required to reproduce the features.

Figure 18: Best-fit UV-CO model (top left), 100 UV-H2 radiative transfer
models with the smallest MSEs (top right), comparison of UV-H2-/UV-CO-based
and observed Ly$\alpha$ profiles (bottom left) and radial distributions of
flux for both species (bottom right) from the disk around CS Cha. Residuals on
the best-fit UV-CO model are color-coded from blue to red, based on the
Ly$\alpha$ pumping wavelength required to excite the upper level of the
transition.

Figure 19: Best-fit UV-CO model (top left), 100 UV-H2 radiative transfer
models with the smallest MSEs (top right), comparison of UV-H2-/UV-CO-based
Ly$\alpha$ profiles (bottom left) and radial distributions of flux for both
species (bottom right) from the disk around CW Tau. Residuals on the best-fit
UV-CO model are color-coded from blue to red, based on the Ly$\alpha$ pumping
wavelength required to excite the upper level of the transition.

Figure 20: Best-fit UV-CO 2-D radiative transfer (top left) and Gaussian
models (top right), 100 UV-H2 radiative transfer models with the smallest MSEs
(middle left), comparison of UV-H2-/UV-CO-based and observed Ly$\alpha$
profiles (middle right) 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$\alpha$ pumping wavelength
required to excite the upper level of the transition.

Figure 21: Best-fit UV-CO 2-D radiative transfer (top left) and Gaussian
models (top right), 100 UV-H2 radiative transfer models with the smallest MSEs
(middle left), comparison of UV-H2-/UV-CO-based and observed Ly$\alpha$
profiles (middle right) and radial distributions of flux for both species
(bottom left) from the disk around DM Tau. Residuals on the best-fit UV-CO
model are color-coded from blue to red, based on the Ly$\alpha$ pumping
wavelength required to excite the upper level of the transition.

Figure 22: Best-fit UV-CO model (top left), 100 UV-H2 radiative transfer
models with the smallest MSEs (top right), comparison of UV-H2-/UV-CO based
and observed Ly$\alpha$ profiles (bottom left) and radial distributions of
flux for both species (bottom right) from the disk around LkCa 15. Residuals
on the best-fit UV-CO model are color-coded from blue to red, based on the
Ly$\alpha$ pumping wavelength required to excite the upper level of the
transition.

Figure 23: Best-fit UV-CO 2-D radiative transfer (top left) and Gaussian
models (top right), 100 UV-H2 radiative transfer models with the smallest MSEs
(middle left), comparison of UV-H2-/UV-CO-based and observed Ly$\alpha$
profiles (middle right) and radial distributions of flux for both species
(bottom left) from the disk around RECX-11. Residuals on the best-fit UV-CO
model are color-coded from blue to red, based on the Ly$\alpha$ pumping
wavelength required to excite the upper level of the transition.

Figure 24: Best-fit UV-CO 2-D radiative transfer (top left) and Gaussian
models (top right), 100 UV-H2 radiative transfer models with the smallest MSEs
(middle left), comparison of UV-H2-/UV-CO-based and observed Ly$\alpha$
profiles (middle right) and radial distributions of flux for both species
(bottom left) from the disk around RECX-15 (spectra from 2010). Residuals on
the best-fit UV-CO model are color-coded from blue to red, based on the
Ly$\alpha$ pumping wavelength required to excite the upper level of the
transition.

Figure 25: Best-fit UV-CO 2-D radiative transfer (top left) and Gaussian
models (top right), 100 UV-H2 radiative transfer models with the smallest MSEs
(top right), comparison of UV-H2-/UV-CO-based and observed Ly$\alpha$ profiles
(bottom left) and radial distributions of flux for both species (bottom right)
from the disk around RECX-15 (spectra from 2013). Residuals on the best-fit
UV-CO model are color-coded from blue to red, based on the Ly$\alpha$ pumping
wavelength required to excite the upper level of the transition.

Figure 26: Best-fit UV-CO model (top left), 100 UV-H2 radiative transfer
models with the smallest MSEs (top right), comparison of UV-H2-/UV-CO-based
and observed Ly$\alpha$ profiles (bottom left) and radial distributions of
flux for both species (bottom right) from the disk around RY Lupi. Residuals
on the best-fit UV-CO model are color-coded from blue to red, based on the
Ly$\alpha$ pumping wavelength required to excite the upper level of the
transition.

Figure 27: Best-fit UV-CO model (top left), 100 UV-H2 radiative transfer
models with the smallest MSEs (top right), comparison of UV-H2-/UV-CO-based
and observed Ly$\alpha$ profiles (bottom left) and radial distributions of
flux for both species (bottom right) from the disk around T Cha. Residuals on
the best-fit UV-CO model are color-coded from blue to red, based on the
Ly$\alpha$ pumping wavelength required to excite the upper level of the
transition.

Figure 28: Best-fit UV-CO 2-D radiative transfer (top left) and Gaussian
models (top right), 100 UV-H2 radiative transfer models with the smallest MSEs
(middle left), comparison of UV-H2-/UV-CO-based and observed Ly$\alpha$
profiles (middle right) and radial distributions of flux for both species
(bottom left) from the disk around UX Tau. The best-fit UV-H2 model is unable
to simultaneously capture emission from the high-velocity wings of the line
profiles ($v>75$ km s-1) and the narrow line cores, so the true flux
distribution likely extends closer to the central star than the inner radius
shown here $\left(r<0.4\,\rm{AU}\right)$. Residuals on the best-fit UV-CO
model are color-coded from blue to red, based on the Ly$\alpha$ pumping
wavelength required to excite the upper level of the transition.

Figure 29: Best-fit UV-CO 2-D radiative transfer (top left) and Gaussian
models (top right), 100 UV-H2 radiative transfer models with the smallest MSEs
(middle left), comparison of UV-H2-/UV-CO-based and observed Ly$\alpha$
profiles (middle right) and radial distributions of flux for both species
(bottom left) from the disk around V4046 Sgr. Residuals on the best-fit UV-CO
model are color-coded from blue to red, based on the Ly$\alpha$ pumping
wavelength required to excite the upper level of the transition.

## References

  * Ádámkovics et al. (2014) Ádámkovics, M., Glassgold, A. E., & Najita, J. R. 2014, ApJ, 786, 135
  * Ádámkovics et al. (2016) Ádámkovics, M., Najita, J. R., & Glassgold, A. E. 2016, ApJ, 817, 82
  * Alencar et al. (2012) Alencar, S. H. P., et al. 2012, A&A, 541, A116
  * ALMA Partnership et al. (2015) ALMA Partnership et al. 2015, ApJ, 808, L3
  * Andrews et al. (2018) Andrews, S. M., et al. 2018, ApJ, 869, L41
  * Andrews et al. (2011) Andrews, S. M., Wilner, D. J., Espaillat, C., Hughes, A. M., Dullemond, C. P., McClure, M. K., Qi, C., & Brown, J. M. 2011, ApJ, 732, 42
  * Ardila et al. (2013) Ardila, D. R., et al. 2013, ApJS, 207, 1
  * Arulanantham et al. (2020) Arulanantham, N., et al. 2020, arXiv e-prints, arXiv:2002.09058
  * Arulanantham et al. (2018) —. 2018, ApJ, 855, 98
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al. 2013, A&A, 558, A33
  * Bacciotti et al. (2018) Bacciotti, F., et al. 2018, ApJ, 865, L12
  * Bae et al. (2016) Bae, J., Zhu, Z., & Hartmann, L. 2016, ApJ, 819, 134
  * Bailer-Jones et al. (2018) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Mantelet, G., & Andrae, R. 2018, AJ, 156, 58
  * Banzatti et al. (2019) Banzatti, A., Pascucci, I., Edwards, S., Fang, M., Gorti, U., & Flock, M. 2019, ApJ, 870, 76
  * Banzatti & Pontoppidan (2015) Banzatti, A., & Pontoppidan, K. M. 2015, ApJ, 809, 167
  * Bergin et al. (2003) Bergin, E., Calvet, N., D’Alessio, P., & Herczeg, G. J. 2003, ApJ, 591, L159
  * Bethell & Bergin (2011) Bethell, T. J., & Bergin, E. A. 2011, ApJ, 739, 78
  * Bosman et al. (2018) Bosman, A. D., Walsh, C., & van Dishoeck, E. F. 2018, A&A, 618, A182
  * Brown et al. (2013) Brown, J. M., Pontoppidan, K. M., van Dishoeck, E. F., Herczeg, G. J., Blake, G. A., & Smette, A. 2013, ApJ, 770, 94
  * Cabrit et al. (2006) Cabrit, S., Pety, J., Pesenti, N., & Dougados, C. 2006, A&A, 452, 897
  * Calvet & Hartmann (1992) Calvet, N., & Hartmann, L. 1992, ApJ, 386, 239
  * Covey et al. (2021) Covey, K. R., Larson, K. A., Herczeg, G. J., & Manara, C. F. 2021, AJ, 161, 61
  * Cox et al. (2013) Cox, A. W., Grady, C. A., Hammel, H. B., Hornbeck, J., Russell, R. W., Sitko, M. L., & Woodgate, B. E. 2013, ApJ, 762, 40
  * Danforth et al. (2010) Danforth, C. W., Keeney, B. A., Stocke, J. T., Shull, J. M., & Yao, Y. 2010, ApJ, 720, 976
  * Edwards et al. (1994) Edwards, S., Hartigan, P., Ghandour, L., & Andrulis, C. 1994, AJ, 108, 1056
  * Ercolano et al. (2011) Ercolano, B., Clarke, C. J., & Hall, A. C. 2011, MNRAS, 410, 671
  * Espaillat et al. (2010) Espaillat, C., et al. 2010, ApJ, 717, 441
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * France et al. (2017a) France, K., et al. 2017a, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10397, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, 1039713
  * France et al. (2014a) France, K., Herczeg, G. J., McJunkin, M., & Penton, S. V. 2014a, ApJ, 794, 160
  * France et al. (2011a) France, K., et al. 2011a, ApJ, 743, 186
  * France et al. (2019) —. 2019, BAAS, 51, 167
  * France et al. (2017b) France, K., Roueff, E., & Abgrall, H. 2017b, ApJ, 844, 169
  * France et al. (2014b) France, K., Schindhelm, E., Bergin, E. A., Roueff, E., & Abgrall, H. 2014b, ApJ, 784, 127
  * France et al. (2011b) France, K., et al. 2011b, ApJ, 734, 31
  * France et al. (2012) —. 2012, ApJ, 756, 171
  * Francis & van der Marel (2020) Francis, L., & van der Marel, N. 2020, arXiv e-prints, arXiv:2003.00079
  * Furlan et al. (2009) Furlan, E., et al. 2009, ApJ, 703, 1964
  * Gaia Collaboration et al. (2018) Gaia Collaboration et al. 2018, A&A, 616, A1
  * Ginski et al. (2018) Ginski, C., et al. 2018, A&A, 616, A79
  * Gorti & Hollenbach (2009) Gorti, U., & Hollenbach, D. 2009, ApJ, 690, 1539
  * Green et al. (2012) Green, J. C., et al. 2012, ApJ, 744, 60
  * Guilloteau et al. (2014) Guilloteau, S., Simon, M., Piétu, V., Di Folco, E., Dutrey, A., Prato, L., & Chapillon, E. 2014, A&A, 567, A117
  * Guzmán et al. (2018) Guzmán, V. V., et al. 2018, ApJ, 869, L48
  * Hartmann et al. (1994) Hartmann, L., Hewett, R., & Calvet, N. 1994, ApJ, 426, 669
  * Hendler et al. (2018) Hendler, N. P., et al. 2018, MNRAS, 475, L62
  * Herczeg et al. (2002) Herczeg, G. J., Linsky, J. L., Valenti, J. A., Johns-Krull, C. M., & Wood, B. E. 2002, ApJ, 572, 310
  * Herczeg et al. (2006) Herczeg, G. J., Linsky, J. L., Walter, F. M., Gahm, G. F., & Johns-Krull, C. M. 2006, ApJS, 165, 256
  * Herczeg et al. (2005) Herczeg, G. J., et al. 2005, AJ, 129, 2777
  * Herczeg et al. (2004) Herczeg, G. J., Wood, B. E., Linsky, J. L., Valenti, J. A., & Johns-Krull, C. M. 2004, ApJ, 607, 369
  * Hoadley et al. (2015) Hoadley, K., France, K., Alexander, R. D., McJunkin, M., & Schneider, P. C. 2015, ApJ, 812, 41
  * Hoadley et al. (2017) Hoadley, K., France, K., Arulanantham, N., Loyd, R. O. P., & Kruczek, N. 2017, The Astrophysical Journal, 846, 6
  * Hornbeck et al. (2016) Hornbeck, J. B., et al. 2016, ApJ, 829, 65
  * Huang et al. (2018a) Huang, J., et al. 2018a, ApJ, 852, 122
  * Huang et al. (2018b) —. 2018b, ApJ, 869, L42
  * Huang et al. (2018c) —. 2018c, ApJ, 869, L43
  * Huélamo et al. (2015) Huélamo, N., de Gregorio-Monsalvo, I., Macias, E., Pinte, C., Ireland, M., Tuthill, P., & Lacour, S. 2015, A&A, 575, L5
  * Ingleby et al. (2011) Ingleby, L., et al. 2011, ApJ, 743, 105
  * Ingleby et al. (2013) —. 2013, ApJ, 767, 112
  * Johansen et al. (2004) Johansen, A., Andersen, A. C., & Brandenburg, A. 2004, A&A, 417, 361
  * Kastner et al. (2018) Kastner, J. H., et al. 2018, ApJ, 863, 106
  * Kudo et al. (2018) Kudo, T., Hashimoto, J., Muto, T., Liu, H. B., Dong, R., Hasegawa, Y., Tsukagoshi, T., & Konishi, M. 2018, ApJ, 868, L5
  * Lawson et al. (1996) Lawson, W. A., Feigelson, E., & Huenemoerder, D. P. 1996, MNRAS, 280, 1071
  * Lawson et al. (2004) Lawson, W. A., Lyo, A.-R., & Muzerolle, J. 2004, MNRAS, 351, L39
  * Liu & Dalgarno (1996) Liu, W., & Dalgarno, A. 1996, ApJ, 467, 446
  * Long et al. (2019) Long, F., et al. 2019, ApJ, 882, 49
  * Long et al. (2018) —. 2018, ApJ, 869, 17
  * Loomis et al. (2017) Loomis, R. A., Öberg, K. I., Andrews, S. M., & MacGregor, M. A. 2017, ApJ, 840, 23
  * Lynden-Bell & Pringle (1974) Lynden-Bell, D., & Pringle, J. E. 1974, MNRAS, 168, 603
  * McJunkin et al. (2013) McJunkin, M., France, K., Burgh, E. B., Herczeg, G. J., Schindhelm, E., Brown, J. M., & Brown, A. 2013, ApJ, 766, 12
  * McJunkin et al. (2016) McJunkin, M., France, K., Schindhelm, E., Herczeg, G., Schneider, P. C., & Brown, A. 2016, ApJ, 828, 69
  * McJunkin et al. (2014) McJunkin, M., France, K., Schneider, P. C., Herczeg, G. J., Brown, A., Hillenbrand, L., Schindhelm, E., & Edwards, S. 2014, ApJ, 780, 150
  * Najita et al. (2003) Najita, J., Carr, J. S., & Mathieu, R. D. 2003, ApJ, 589, 931
  * Nomura et al. (2007) Nomura, H., Aikawa, Y., Tsujimoto, M., Nakagawa, Y., & Millar, T. J. 2007, ApJ, 661, 334
  * Nomura & Millar (2005) Nomura, H., & Millar, T. J. 2005, A&A, 438, 923
  * Osterman et al. (2011) Osterman, S., et al. 2011, Ap&SS, 335, 257
  * Pascucci et al. (2020) Pascucci, I., et al. 2020, arXiv e-prints, arXiv:2009.09114
  * Pérez et al. (2018) Pérez, L. M., et al. 2018, ApJ, 869, L50
  * Piétu et al. (2014) Piétu, V., Guilloteau, S., Di Folco, E., Dutrey, A., & Boehler, Y. 2014, A&A, 564, A95
  * Pontoppidan et al. (2011) Pontoppidan, K. M., Blake, G. A., & Smette, A. 2011, ApJ, 733, 84
  * Pontoppidan et al. (2008) Pontoppidan, K. M., Blake, G. A., van Dishoeck, E. F., Smette, A., Ireland, M. J., & Brown, J. 2008, ApJ, 684, 1323
  * Qi et al. (2019) Qi, C., et al. 2019, ApJ, 882, 160
  * Reipurth et al. (1996) Reipurth, B., Pedrosa, A., & Lago, M. T. V. T. 1996, A&AS, 120, 229
  * Rodriguez et al. (2010) Rodriguez, D. R., Kastner, J. H., Wilner, D., & Qi, C. 2010, ApJ, 720, 1684
  * Rosenfeld et al. (2012) Rosenfeld, K. A., Andrews, S. M., Wilner, D. J., & Stempels, H. C. 2012, ApJ, 759, 119
  * Ruíz-Rodríguez et al. (2019) Ruíz-Rodríguez, D., Kastner, J. H., Dong, R., Principe, D. A., Andrews, S. M., & Wilner, D. J. 2019, AJ, 157, 237
  * Salyk et al. (2011a) Salyk, C., Blake, G. A., Boogert, A. C. A., & Brown, J. M. 2011a, ApJ, 743, 112
  * Salyk et al. (2011b) —. 2011b, ApJ, 743, 112
  * Saucedo et al. (2003) Saucedo, J., Calvet, N., Hartmann, L., & Raymond, J. 2003, ApJ, 591, 275
  * Schindhelm et al. (2012a) Schindhelm, E., France, K., Burgh, E. B., Herczeg, G. J., Green, J. C., Brown, A., Brown, J. M., & Valenti, J. A. 2012a, ApJ, 746, 97
  * Schindhelm et al. (2012b) Schindhelm, E., et al. 2012b, ApJ, 756, L23
  * Schneider et al. (2013) Schneider, P. C., Eislöffel, J., Güdel, M., Günther, H. M., Herczeg, G., Robrade, J., & Schmitt, J. H. M. M. 2013, A&A, 557, A110
  * Schneider et al. (2015) Schneider, P. C., France, K., Günther, H. M., Herczeg, G., Robrade, J., Bouvier, J., McJunkin, M., & Schmitt, J. H. M. M. 2015, A&A, 584, A51
  * Shakhovskoj et al. (2006) Shakhovskoj, D., et al. 2006, A&A, 448, 1075
  * Simon et al. (2016) Simon, M. N., Pascucci, I., Edwards, S., Feng, W., Gorti, U., Hollenbach, D., Rigliaco, E., & Keane, J. T. 2016, ApJ, 831, 169
  * The Astropy Collaboration et al. (2018) The Astropy Collaboration et al. 2018, ArXiv e-prints
  * Turner et al. (1999) Turner, N. J., Bodenheimer, P., & Różyczka, M. 1999, ApJ, 524, 129
  * Uvarova et al. (2020) Uvarova, A. V., Günther, H. M., Principe, D. A., & Schneider, P. C. 2020, AJ, 160, 39
  * van der Marel et al. (2013) van der Marel, N., et al. 2013, Science, 340, 1199
  * van der Marel et al. (2018) —. 2018, ApJ, 854, 177
  * Walter et al. (2003) Walter, F. M., et al. 2003, AJ, 126, 3076
  * Weber et al. (2020) Weber, M. L., Ercolano, B., Picogna, G., Hartmann, L., & Rodenkirch, P. J. 2020, MNRAS, 496, 223
  * Woitas et al. (2002) Woitas, J., Ray, T. P., Bacciotti, F., Davis, C. J., & Eislöffel, J. 2002, ApJ, 580, 336
  * Woitke et al. (2011) Woitke, P., et al. 2011, A&A, 534, A44
  * Wolven et al. (1997) Wolven, B. C., Feldman, P. D., Strobel, D. F., & McGrath, M. A. 1997, ApJ, 475, 835
  * Yang et al. (2012) Yang, H., et al. 2012, ApJ, 744, 121