The Radial Distribution and Excitation of H₂ around Young Stars in the HST-ULLYSES Survey

The fluorescent H₂ lines observed in the CTTS sample can be used to constrain the spatial distribution of H₂ in the circumstellar environment. For our line profile analysis, we focus on the measurement of four progressions ([${v}^{{\prime} }$,${J}^{{\prime} }]$ = [1, 7], [1, 4], [0, 1], and [0, 2]). These progressions are pumped through the (1–2)R(6) λ_lab 1215.73 Å, (1–2)P(5) λ_lab 1216.07 Å, (0–2)R(0) λ_lab 1217.21 Å, and (0–2)R(1) λ_lab 1217.64 Å absorbing transitions, respectively. ¹⁷ The absorbing transitions are within +14 to +487 km s⁻¹ of the Lyα line center. We target the brightest two, unblended emission lines from each progression for analysis here. The H₂ emission line centroids are broadly consistent with the stellar radial velocities within the 15 km s⁻¹ accuracy of the HST/COS wavelength solution, and we will explore multiple emission line components (e.g., Arulanantham et al. 2018) and their kinematic signatures in a future work.

The emission line fluxes are corrected for interstellar reddening using the A_V values given in Table 1. It has been shown that the dust properties of the Taurus star-forming region differ from the Milky Way average extinction curve (Calvet et al. 2004). However, the shape of the extinction curve for star-forming regions and individual targets remains highly uncertain (McJunkin et al. 2016), and the A_V values derived for individual targets carry significant uncertainties (Carvalho & Hillenbrand 2022; Pittman et al. 2022). Both of these issues are exacerbated by the strong wavelength dependence of the dust attenuation curve, where the reddening correction at 1500 Å has a ∼6× larger amplitude relative to 6500 Å for A_V = 1. In light of these uncertainties, we adopt a typical interstellar curve (Cardelli et al. 1989) and an R_V = 3.1. If, on the other hand, we would have adopted a standard interstellar curve with R_V = 5.5 (for A_V = 1), the resultant fluxes would be between 20%–40% of the values reported here.

The total flux from a progression m is given by

$\begin{eqnarray}&&{F}_{m}({{\rm{H}}}_{2})=\frac{1}{2}\sum \left(\frac{{F}_{\mathrm{mn}}}{{B}_{\mathrm{mn}}}\right),\end{eqnarray} \tag{ 2 }$

where F_mn is the reddening-corrected, integrated H₂ emission line flux from rovibrational state m (=$[{v}^{{\prime} }$, ${J}^{{\prime} }$]) in the ${B}^{1}{{\rm{\Sigma }}}_{u}^{+}$ electronic state to n (=[v'', J'']) in the ground electronic state, ${X}^{1}{{\rm{\Sigma }}}_{g}^{+}$. B_mn is the branching ratio between these two states, and 2 is the number of emission lines measured from a given progression. The measurement errors are typically small, so we take the flux error to be the average error of the individual measurements of F_m (H₂). The dominant systematic error on the measured H₂ flux is the correction for interstellar reddening; we do not attempt to account for this uncertainty in the flux and luminosity errors presented below. The total progression luminosity is then L_m (H₂) = (4π d²)F_m (H₂). In Table 2, we present reddening-corrected luminosities for the [1, 4] and [0, 1] progressions. The total H₂ luminosity is the sum of the individual progression luminosities. Because we select a finite number of progressions to analyze, this treatment underestimates the total H₂ fluorescent luminosity by an amount that depends on the temperature and spatial distribution of the absorbing H₂, as well as the shape of the Lyα pumping spectrum. As these four progressions are the brightest in essentially all CTTSs, we consider this to be a 20%–30% underestimate of the total fluorescent H₂ emission (France et al. 2012).

The fluorescent emission lines of H₂ can be used to reconstruct the bright Lyα emission line that is the primary excitation source of this gas (Wood et al. 2002; Herczeg et al. 2004; Schindhelm et al. 2012). In order for this process to produce a high-fidelity reconstruction, a sufficient number of H₂ progressions, covering the majority of the Lyα pumping profile, must be detected at moderate to high S/N (e.g., Schindhelm et al. 2012; Arulanantham et al. 2018). Lyα reconstructions with a small number of low-S/N progressions leave the shape and integrated flux of the Lyα line poorly constrained (Espaillat et al. 2022). Of the 63 accreting sources analyzed here, the majority are only detected in three or four H₂ progressions—54 detections in [1, 4], 41 detections in each [1, 7] and [0, 1], 33 detections in [0, 2], and only three detections in the [2, 12] progression that is often used to constrain the wings of the Lyα line. Therefore, a full Lyα reconstruction cannot be practically applied to the full ULLYSES sample.

In order to provide a uniform estimate for the Lyα luminosity for ULLYSES, we fit F(H₂) versus F(Lyα) relation using the 14 higher-S/N targets reconstructed by Schindhelm et al. (2012) and H₂ fluorescence data from France et al. (2012). These data are shown as the orange points in Figure 4. We extrapolated the ULLYSES measurements of L(H₂) described above onto this relation to estimate the Lyα fluxes for all 63 accreting stars in our sample

$\begin{eqnarray}&&{\mathrm{log}}_{10}(F(\mathrm{Ly}\alpha ))={m}_{\mathrm{Ly}\alpha }({\mathrm{log}}_{10}(F({{\rm{H}}}_{2})))+{b}_{\mathrm{Ly}\alpha },\end{eqnarray} \tag{ 3 }$

where the best-fit coefficients are m_Lyα = 0.791 ± 0.033 and b_Lyα = −1.297 ± 0.400. Propagating the uncertainties on the fit coefficients of the F(H₂) versus F(Lyα) relation indicates a maximum uncertainty on the associated Lyα luminosities of ≲40%. The small dispersion seen in the extrapolated Lyα luminosities is due to the conversion of the flux–flux relationship from France et al. (2012) to luminosity. We caution that while the original Lyα-to-H₂ relationship from France et al. (2012) spans approximately 1.5 orders of magnitude in L(H₂), we are extrapolating that relationship to over three orders of magnitude for the present sample. Figure 5 illustrates a comparison of the estimated Lyα to C iv luminosity between the present work and the France et al. (2012) sample, showing reasonable agreement in the range L(C iv) > 10²⁹ erg s⁻¹ with a larger spread in the lower C iv luminosity region of the ULLYSES sample. The spread in Lyα values at L(C iv) < 10²⁹ erg s⁻¹ is driven by scatter in the relationship between C iv and H₂. Similarly, we find a moderate positive correlation between L(Lyα) and the photometric luminosity (L_*; from Manara et al. 2021, 2022); the Spearman rank correlation coefficient, ρ, and the likelihood that this correlation is consistent with a random distribution, n, for the L(Lyα) versus L_* relation are ρ = 0.498 and n = 2.3 × 10⁻⁴, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Kinematic broadening dominates the observed H₂ line profiles. The thermal broadening of the emission lines is approximately 4.5 km s⁻¹ at the nominal 2500 K H₂ layer (Herczeg et al. 2004; Hoadley et al. 2015); significant additional broadening would require temperatures in excess of the ≈4500 K dissociation temperature of H₂ (Lepp & Shull 1983). If we further assume that any turbulence in the disks is subsonic, then the maximum turbulent velocity will be no larger than a few kilometers per second. Therefore, velocity broadening due to bulk motions and Keplerian rotation dominate the observed line shapes when the FWHM of the emission line is greater than the ≈15–20 km s⁻¹ spectral resolution of COS (see Figure 2; all H₂ lines are well resolved at the spectral resolution of HST/COS; the 1450 Å HST/COS LSF is overplotted as the green dotted line). For the case of H₂ in a circumstellar disk where Keplerian broadening dominates the profile, we employ a simple metric to characterize the average H₂ emitting radius, 〈R_H2〉 (Salyk et al. 2011a; France et al. 2012; Arulanantham et al. 2021), where

$\begin{eqnarray}&&\langle {R}_{{\rm{H}}2}{\rangle }_{m}={{GM}}_{* }{\left(\frac{2\sin (i)}{{\mathrm{FWHM}}_{m}}\right)}^{2}{,}\end{eqnarray} \tag{ 4 }$

where M_* is the stellar mass, i is the disk inclination angle, taken primarily from Atacama Large Millimeter/submillimeter Array (ALMA) surveys (Table 1), and FWHM_m is the mean of the Gaussian FWHMs for a given progression m. When disk inclinations were not available (19/63 accreting sources), we assume a disk inclination of 60°, as this is approximately the average inclination found in ALMA disk surveys (Ansdell et al. 2018) and assumed in model studies (Rilinger & Espaillat 2021). We provide correlation coefficients for the full sample and using only the sources with known inclinations in the appropriate subsections below, and note that the "outlier" values of 〈R_H2〉 are not driven by unknown inclinations; 5/6 of the largest H₂ radii presented here are for sources with measured disk inclinations. Table 1 presents the stellar masses and disk inclinations used in this work.

For the H₂ radial distribution analysis, we consider four progressions: the [1, 7] and [1, 4] progressions that are pumped within 100 km s⁻¹ of the Lyα line center and two progressions ([0, 1] and [0, 2]) that are pumped from the wing of the Lyα emission profile (379 and 487 km s⁻¹ from line center, respectively). We analyze multiple progressions because self-absorption can selectively impact the lower-energy level gas (Wood et al. 2002; McJunkin et al. 2016) and that off-star Lyα emission and the local scattering geometry can impact the flux distributions of the different progressions (Walter et al. 2003; Arulanantham et al. 2023). We will return to a discussion of these results in Section 4.

FUV continuum spectra were created for the ULLYSES DR5 sample following the continuum extraction methodology described in France et al. (2014b). We create a grid of 205 unique spectral points between 1138 and 1791 Å, selected by hand to avoid discrete molecular and atomic emission and absorption features, where 0.75 Å (approximately 10 spectral resolution elements) spectral continuum windows can be cleanly measured. We measure the mean and standard deviation of the observed spectra, and these points define the binned flux spectrum and error array (see Figure 3). These binned spectra are then corrected for interstellar reddening and the FUV continuum and H₂ continuum emission (the "1600 Å Bump") are separated as described below.

We measure the FUV continuum emission by fitting the binned spectrum with a second-order polynomial at wavelengths away from the Bump, and residual stellar or airglow emission in the regions Δλ = 1145–1190 Å, 1245–1330 Å, 1395–1401 Å, 1420–1465 Å, 1690–1710 Å, and 1730–1750 Å (see the blue dashed lines in Figure 3). The polynomial fit is extrapolated to the atomic hydrogen ionization limit, 912 Å, to cover the primary photodissociating spectral region for H₂, CO, and N₂. Comparison with archival FUSE observations between 1000 and 1180 Å has shown that this extrapolation is a reasonably good match to the observed short-wavelength flux (France et al. 2014b). The FUV continuum flux, F_FUVCont, is the integral of the polynomial over the 912–1760 Å wavelength region and the FUV continuum luminosity is defined as L(FUV Cont) = 4π d² F_FUVCont.

Figure 6 shows the FUV continuum luminosity as a function of stellar mass. We observe a positive correlation between the continuum and stellar mass (ρ =0.570 and n = 1.0 × 10⁻⁴, as shown in the legend of Figure 6), which is expected in the scenario where the FUV continuum emission is dominated by mass accretion onto the central star (France et al. 2014b). Figure 6 plots the sources for which "good" H₂ radii are measured (see the discussion in Section 4.1), and it is observed that the most luminous sources are also the sources with the largest reddening (A_V ). The attenuation toward FUV-emitting regions of protoplanetary targets has been estimated to be lower than estimates inferred from X-ray, optical, or IR measurements (e.g., McJunkin et al. 2014). Hence, extinction values derived from, e.g., IR observations may overestimate the extinction correction for UV data. This was demonstrated by France et al. (2017), who showed that the largest A_V values derived from IR-based extinctions led to unphysical correlations in FUV luminosity relationships. To avoid biasing the results by the largest, possibly incorrect, A_V values, Figure 6 also shows the 41 stars with A_V < 1 as pink circles. The correlation coefficient and best-fit line are presented for the A_V < 1 sample.

Zoom In Zoom Out Reset image size

The most prominent spectral feature in the binned spectra of some CTTSs is the 1600 Å Bump feature created during the dissociation of H₂ molecules with highly nonthermal rovibrational distributions, which spans ∼30–150 Å and has a peak flux ∼0–4 times the underlying continuum level at 1600 Å. We define the 1600 Å Bump as the excess emission above the FUV continuum in the 1490–1690 Å spectral region, where the Bump spectrum is the binned data minus the FUV continuum fit. An example of the well-defined Bump spectrum of DM Tau is shown in Figure 3. Following previous work, we parameterize the Bump emission as a broad Gaussian plus a second-order polynomial (France et al. 2017; Figure 3, bottom panel). We define the total Bump flux, F_Bump, as the integral of the continuum-subtracted Bump spectrum from 1490 to 1690 Å, and L(1600 Å Bump) = 4π d² F_Bump.

The average H₂ line emission radius represents the characteristic Keplerian velocity of the gas, assuming the spectral broadening is dominated by the rotational motion of the disk. In a future work, we will present a more detailed line profile decomposition of a subset of bright disks to constrain stratified gas structures and wind contributions (e.g., Herczeg et al. 2006; Gangi et al. 2023), but for the full ULLYSES sample, we assume a single emission line component. The objects included in the analysis are restricted to disk inclinations > 15° (the Doppler motion cannot be confidently measured for face-on disks at the resolution of HST/COS) and do not have any data quality flags (low S/N or spectral blending) in the H₂ emission line catalog.

The final sample of sources that meet these criteria ("good" H₂ radii) includes 52, 38, 39, and 31 disks in the [1, 4], [1, 7], [0, 1], and [0, 2] progressions, respectively. 〈R_H2〉 values range between 0.05 and 2.9 au, with average emitting radii between 0.53 and 0.72 au (median radii between 0.35 and 0.57 au) for the four progressions. The majority of the fluorescent line emission (based on the standard deviation of the radial distributions of the four progressions) arises between 0.1 and 1.4 au. This typical emitting range can be seen in the cluster of radial distributions in Figure 7. The minimum range of H₂ emitting radii (0.05–0.1 au) are comparable to the dust sublimation radius and the inner disk dust rim radii measured by GRAVITY (0.1–0.2 au; GRAVITY Collaboration et al. 2021), suggesting that the molecular gas and dust edges may overlap in some sources.

Zoom In Zoom Out Reset image size

Figure 7 shows the relationship between 〈R_H2〉 and the C iv luminosity (ρ = 0.383 and n = 5.6 × 10⁻³, reducing to ρ = 0.330 and n = 4.9 × 10⁻² if restricted to disks with measured inclinations). Given the variability in mass accretion rates and FUV fluxes over time (Espaillat et al. 2019; Claes et al. 2022; Hinton et al. 2022), we use the L(C iv) measured contemporaneously with the H₂ emission lines to provide a proxy for the "real-time" FUV radiation field of the central star. In Figure 7, we observe a marginal positive correlation between the emitting radius and the strength of the C iv emission. This result may indicate that higher FUV fluxes result in larger equilibrium radii with the heating/dissociating FUV radiation field (e.g., Cazzoletti et al. 2018), suggesting higher levels of H₂ dissociation, or the clearing of H₂-shielding inner disk dust, close to the central star in more FUV-luminous systems (Figure 6). We note that the lowest C iv luminosity systems (log₁₀(L(C iv)) ≲ 28.7) populate the lower left quadrant in Figure 7 (${\mathrm{log}}_{10}\langle {R}_{{\rm{H}}}2\rangle \lesssim 0.3$). Comparing the distribution of L(C iv) with the stellar mass reveals that all stars with log₁₀(L(C iv)) ≲ 28.7 also have M_* < 0.3 M_⊙. This supports the idea that the lowest-mass sources with the lowest L(C iv) show less dissociation in their innermost disks than their more massive counterparts. This would lead the stellar Lyα photons to encounter optically thick gas at smaller radial distances, which manifests as a small 〈R_H2〉 in these systems.

We observe a stronger correlation (ρ = 0.692 and n = 3.2 × 10⁻⁵) between the average H₂ emitting radius and the disk dust mass (Figure 8). The actual distribution is not a straight linear trend; instead, the H₂ emitting radii are roughly constant between ≈0.1–0.4 au for dust masses ≲ 20 M_Earth, with no H₂ radii larger than ≈0.4 au. For dust masses ≳ 20 M_Earth, the H₂ radii range from 0.2 to 2 au. It is established that more massive stars host massive dust disks (Andrews et al. 2013; Pascucci et al. 2016; Ansdell et al. 2017) and that dust disk mass is correlated with mass accretion rate (Manara et al. 2016; Mulders et al. 2017; Manara et al. 2022). The observation of larger gas clearing radii around stars with larger dust masses (and higher FUV continuum luminosities) is consistent with this picture of inner disk structure, as is the scenario presented above where the stars with the lowest masses (and C iv luminosities) were found to host H₂ populations with the smallest emitting radii.

Zoom In Zoom Out Reset image size

Figure 9 displays the relationship between the WISE W3 – W4 color ¹⁸ and 〈R_H2〉. With the relatively low Spearman correlation coefficient (ρ = 0.326) and high probability of noncorrelation (n = 0.079), we do not find a correlation between the slope of the inner disk dust SED and the average H₂ emitting radius. This is in some tension with the findings of Hoadley et al. (2015), who observed a positive correlation between their modeled radial distribution of H₂ fluorescence and the n_13–31 dust spectral energy distribution derived from Spitzer/IRS spectroscopy. There are two key differences between these studies: first, the Spitzer data employed by Hoadley et al. (2015) are a cleaner metric of warm inner disk dust clearing than the WISE photometric points (Espaillat et al. 2014). The WISE W3 band in particular includes the 10 μm silicate emission feature seen in many CTTS spectra (Olofsson et al. 2010), and the lower angular resolution of WISE increases the possibility of source confusion. Second, Hoadley et al. (2015) employ an emission line forward-modeling analysis of a higher-S/N sample of disk targets; they find correlations between inner disk dust clearing and the maxima of the H₂ radial distributions derived from their models. A follow-on study of the high-S/N ULLYSES sources with Spitzer or JWST spectra would be valuable for a more direct comparison with the Hoadley et al. (2015) results.

Zoom In Zoom Out Reset image size

There are two proposed origins for the 1600 Å Bump emission seen in many protoplanetary disk sources. The first proposes that free electrons (possibly liberated by X-ray emission from the central star) collide with the disk surface population of H₂, creating an electron impact spectrum similar to that observed in the aurorae of giant planets (Bergin et al. 2004; Ingleby et al. 2009). The second route for the formation of the 1600 Å Bump proposes that bright stellar+accretion-generated Lyα emission dissociates water molecules and that the broad continuum Bump is produced during the photodissociation of the nonthermal H₂ created during the breakup (dissociation) of the water molecules (van Harrevelt & van Hemert 2008; France et al. 2017).

Analysis of the ULLYSES DR5 sample finds 24 sources with detectable 1600 Å Bump emission, a detection rate of 38%, roughly consistent with the ≈50% detection rate reported in France et al. (2017). We find a strong correlation between the strength of the Bump and the Lyα field strength. Figure 10 illustrates the relationship between the Lyα and 1600 Å Bump luminosities, where the highly significant correlation (ρ = 0.8 and n = 2.7 × 10⁻⁶) strongly suggests Lyα photons as central to the formation of the H₂ continuum emission feature. For comparison, the correlation between C iv and the 1600 Å Bump luminosity is less pronounced (ρ = 0.5 and n = 1.6 × 10⁻²). We also note that the Lyα luminosities in sources with strong Bump detections are ≳5 × 10³¹ erg s⁻¹, i.e., the Lyα region strongly anchored by the France et al. (2012) Lyα reconstruction sample.

Zoom In Zoom Out Reset image size

Within the sample of 24 1600 Å Bump detections, 14 had WISE W3 – W4 colors without data quality warnings. We plot the Bump luminosity versus the W3 – W4 color in Figure 11. We observe a correlation (ρ = 0.68 and n = 7.6 × 10⁻³) between these quantities, with disks showing the largest depletion of optically thick warm dust (noting the caveats about spectral and spatial contamination in the WISE photometric bands discussed previously) having the largest H₂ continuum emission luminosity. Taken together with the strong correlation with L(Lyα), the emerging picture for the formation of the 1600 Å Bump includes both strong accretion-generated Lyα flux and low optical depth in the inner disk gas and dust environment. Bright Lyα emission alone is not a sufficient condition to produce the H₂ continuum emission, as there are 44 targets in our sample with L(Lyα) > 10³¹ erg s⁻¹ (including, e.g., CVSO-109, Figure 3; Espaillat et al. 2022) yet only 21 of those are detected in the 1600 Å Bump feature. The second required condition is that those Lyα photons are able to reach the molecular reservoirs in the inner disk. As the gas and dust opacities decline with time, Lyα photons can more readily reach the molecular populations of the inner disk, where they dissociate water molecules and the resultant nonthermal H₂ gas.

Zoom In Zoom Out Reset image size