Molecular Gas Tracers in Young and Old Protoplanetary Disks

Three sources with detected CO emission and relatively high dust masses were selected from the ALMA survey of Upper Sco by Barenfeld et al. (2016). N₂H⁺ observations of J1608 and J1609 were taken in 2016 (ALMA Cycle 3, PI: Anderson, 2015.1.01199.S) and previously reported by Anderson et al. (2019). A third source, J1614, was chosen to represent a disk with a similar dust mass but higher CO flux (by a factor of 5–20) than the other two. In addition to searching for N₂H⁺ in J1614, we observed HCO⁺, HCN, ¹³CO, and C¹⁸O J = 3–2 line emission from all three sources during ALMA Cycle 6 (PI: Anderson, 2018.1.01623.S).

Band 6 observations of HCO⁺ and HCN were taken on 2018 December 4. Band 7 observations of N₂H⁺ were taken on 2018 December 8. Band 7 observations of ¹³CO and C¹⁸O were taken on 2019 March 8 and 14–15. Observations were taken with 43–47 antennas and maximum baselines of 740.4, 783.5, and 313.7–360.6 m, respectively. Spectral setups included multiple windows to capture line and continuum emission as shown in Table A1. Observed channel widths were 244.141 kHz (∼0.25 km s⁻¹) for spectral line windows and 976.562 or 15,625 kHz for continuum windows. The on-source integration time per source was about 23 minutes for Band 6, 36 minutes for Band 7 observations of N₂H⁺, and 57 minutes for Band 7 observations of ¹³CO and C¹⁸O.

The data were calibrated using the standard ALMA pipeline and analyzed with the Common Astronomy Software Applications (CASA) package (CASA Team et al. 2022). When necessary, calibrated measurement sets were restored using the version of CASA specified in the data set documentation available on the ALMA archive. Otherwise, imaging and analysis were performed using CASA v6.4.0. Pipeline flux and bandpass calibrators were J1427–4206 in Band 6 and J1517–2422 in Band 7. J1625–2527 was used as the pipeline phase calibrator. Self-calibration was tested but ultimately not applied to the data because it did not produce improvements in the signal-to-noise ratio that were significant for this analysis.

Four sources with strong detections of ¹³CO and HCO⁺ emission were selected from the IRAM 30 m survey of Taurus by Guilloteau et al. (2013, 2016). These sources represent a younger age group relative to Upper Sco, including the Class I source IRAS 04302+2247 (Garufi et al. 2021), and included binary/multisystems (see Table 1) to add more diversity to the total sample. Using the SMA (PI: Anderson, 2018B-S046), we searched for N₂H⁺, HCO⁺, and HCN J = 3–2 and CO, ¹³CO, and C¹⁸O J = 2–1 line emission from our sources. Data were taken between 2018 December 16–19, with 28 baselines (15 baselines on December 19) in the "compact" array configuration, providing an angular resolution of 3''–4''. Each source was observed for a total of 5–6 hr. The LO frequencies were 272.489 GHz for the RxA receiver (345 GHz setting) and 225.119 GHz for the RxB receiver (240 GHz setting), covering our selected suite of molecular lines and several additional species including C₂H (not presented here).

The SWARM correlator provides 8 GHz bandwidth per sideband, divided into four equal sized chunks with a uniform spectral resolution of 140 kHz (0.18 km s⁻¹ at 230 GHz and 0.15 km s⁻¹ at 280 GHz). Calibration of visibility phases and amplitudes was performed using periodic observations of quasars 3C111 and 0510+180, at 15 minute intervals. Measurements of Uranus were used to obtain an absolute scale for calibration of the flux densities. 3C279 was used as the bandpass calibrator. All data were phase and amplitude calibrated using the MIR software package. ¹¹ The calibrated SMA data were then exported in uvfits format for subsequent imaging and analysis using the CASA software package v6.4.0. The following data were flagged: 268.0–268.5 GHz in the lower sideband and 276.5–277.0 GHz in the upper sideband for Antenna 3 in the 345 GHz receiver data and 220.8–221.2 GHz in the lower sideband and 229.1–229.5 GHz in the upper sideband for Antenna 2 in the 240 GHz receiver data.

Line fluxes were derived from the new observations described in Sections 2.1–2.2 using the following steps. In addition, N₂H⁺ line fluxes were rederived from the 2015.1.01199.S data set to produce values using the same methodology. Continuum subtraction was performed on the spectral line windows using CASA's uvcontsub. For each source and molecular line, an initial empirical mask was generated based on a smoothed version of the "dirty" image cube (an image cube generated with zero iterations of the CLEAN algorithm). The dirty image cube was smoothed via convolution with a Gaussian beam with dimensions 1.5× the original beam size. For each channel in the smoothed dirty image cube, we selected pixels with flux densities > 3× the standard deviation of a sample of off-source pixels. For the SMA data, the same procedure was used but without smoothing. The largest 2D clump of selected pixels in each channel is taken as the initial mask shape for that channel. We exclude any clumps in spatial regions (or channels) that are too far separated from the central source region (or velocity) to be probable line emission based on visual inspection of the images. For the relatively weak C¹⁸O emission (and N₂H⁺ emission in J1608) where empirically derived masks were overly influenced by noise, the mask generated from the stronger ¹³CO emission was used. The ¹³CO mask was also used in determining the upper limits on fluxes for molecular lines that were not detected. For GG Tau A, RY Tau, and UZ Tau E, empirical masks generated from the strongest emission line, CO J = 2–1 (see Appendix, Figure A2) were used for all other lines.

Images are generated by CASA's tclean using Briggs weighting with a robust parameter of 0.5, a threshold of roughly 3× the standard deviation of the flux density in a sample of line-free channels, and channel widths of 1 km s⁻¹. We provided an input mask to tclean based on the empirical mask generated for each line as described above. For most lines, we used the empirical mask described above as the input mask. For the weak and nondetected lines where the input mask used was that of a brighter line, a channel-summed version of the bright line mask is used. This is done as a precaution to avoid introducing a fake line signature via the tclean mask. To generate the channel-summed mask, the initial empirical bright line mask components in each channel are summed over all channels. This summed mask is applied over the range of channels where line emission is observed. The line emission velocity ranges are −5–15 km s⁻¹ for Upper Sco sources (−15 to 20 km s⁻¹ for HCO⁺ in J1614 given its unusually wide emission line as shown in Figure 1), 4–9 km s⁻¹ for GG Tau A, 2–8 km s⁻¹ for IRAS 04302+2247, 2–10 km s⁻¹ for RY Tau, and 2–13 km s⁻¹ for UZ Tau E. The listed velocity ranges are inclusive for 1 km s⁻¹ channels.

Zoom In Zoom Out Reset image size

Spectra are generated by multiplying the final image by a specific mask using CASA's immath and producing a spectrum of the product using specflux. For the sufficiently bright and largely unresolved line emission in this work, the empirical masks based on the dirty images do not significantly differ from those derived from the clean images when using the same approach and remain appropriate for this step. The black lines in Figure 1 show spectra generated using a channel-summed version of each line's empirical mask over the entire image cube. The blue solid (for detections) and red dotted (for nondetections) lines show the flux collected in each individual channel of the original (non-channel-summed) empirical mask. The blue spectra are integrated over the line velocity range to produce the flux estimates given in Table 2. Noise estimates are determined by applying each line's empirical mask to 30 different line-free regions across the image cube. Because some of the selected line-free regions are off the phase center, the primary beam correction using CASA's impbcor was applied to the image before the noise and line flux estimates were made. The spatial regions sampled are selected to be as close to the source region as possible without sampling line emission to mitigate overestimation of the noise, which is dependent on the distance from the phase center. In our case, the primary beam correction had almost no effect on the line flux estimates. The standard deviation of fluxes from the 30 line-free regions is added in quadrature with 10% of the integrated line flux representing the assumed uncertainty for the flux amplitude calibration (Cortes et al. 2023) and provided as the uncertainties listed in Table 2. Upper limits are estimated as three times this uncertainty value with fluxes above this threshold being considered detections.

Our analysis combines the observations described above with data from the literature. Here we describe the literature data that have been included. Available source distances are from Manara et al. (2022), which are derived by inverting their parallaxes from Gaia Data Release 3 (Gaia Collaboration et al. 2021). The values for AS 209, TW Hya, and V4046 Sgr (which are not included in Manara et al. 2022) are derived in the same manner. For IRAS 04302+2247, the distance is taken to be 161 pc (Garufi et al. 2021). Most sources belong to the Lupus (≲3 Myr), Taurus (∼1–3 Myr), and Upper Sco (∼5–10 Myr) star-forming regions (Manara et al. 2022). AS 209 resides to the northeast of the main Ophiuchus region in Oph N 3a, with an estimated age ∼1 (0.4–2.5) Myr (Andrews et al. 2018). TW Hya and V4046 Sgr are located in the nearby TWA (∼8 Myr) and β Pic MG (∼24 Myr), respectively (Kastner & Principe 2022).

Dust masses (Table 1) are computed as described in Manara et al. (2022) with the formula used by Ansdell et al. (2016) from Hildebrand (1983):

$$\begin{eqnarray*}&&{M}_{{\rm{dust}}}=\displaystyle \frac{{F}_{\nu }{d}^{2}}{\left.{\kappa }_{\nu }{B}_{\nu }({T}_{{\rm{dust}}}\right)},\end{eqnarray*}$$

where κ_ν = 2.3(ν/230 GHz) cm² g⁻¹ and assuming a dust temperature of 20 K. This assumes optically thin, isothermal dust emission at (sub-)millimeter wavelengths. The submillimeter dust emission may in fact be optically thick (e.g., Zhu et al. 2019; Sierra et al. 2021; Xin et al. 2023), resulting in higher dust masses and complicating our search for trends in this analysis since we would not be tracking the true dust mass. Continuum fluxes from Manara et al. (2022) are used when available. In addition, the continuum flux data from the following works are included: Andrews et al. (2018) for AS 209, van't Hoff et al. (2020) for IRAS 04302+2247, Tsukagoshi et al. (2016) for TW Hya, and Kastner et al. (2018) for V4046 Sgr. The selected continuum fluxes used here are those measured at frequencies of 225–240 GHz, with the exception of Upper Sco for which the continuum fluxes were measured at 340.7 GHz. The uncertainty in the dust mass (see the error bars in Figure 2) is taken as a 10% systematic uncertainty in the flux amplitude calibration (e.g., Andrews et al. 2013, 2018) converted into a mass value using the formula above.

Zoom In Zoom Out Reset image size

The ALMA observations of GG Tau A by Akeson et al. (2019) result in a dust mass ∼ 4 M_⊕ for the circumstellar disk around GG Tau Aa. At their high angular resolution (with a beam of 018 × 016), they are able to separate the circumstellar emission from the wider ring around the GG Tau A multiple system. Given that our line flux observations of GG Tau A were taken at a lower angular resolution and we do not separate out these different components, we use the larger continuum flux value measured by Andrews et al. (2013) for the dust mass (∼277 M_⊕) in our analyses.

Spectral types, stellar luminosities, stellar masses, mass accretion rates, and millimeter dust disk sizes (based on the radius enclosing 68% of the total continuum flux, R_cont,68) are taken from Manara et al. (2022, and references therein) when available. See Table 1 for additional references used for data on AS 209, TW Hya, and V4046 Sgr. Also note that Sz 68, GG Tau A, UZ Tau E, and V4046 Sgr are all binary or multiple-star systems, which may complicate the interpretation of their data and behavior relative to single stars.

Line fluxes are taken from the literature referenced in Table 2. In cases where multiple literature values were available, higher sensitivity observations with stronger signals in the observed spectra are preferred. When multiple analyses were performed on the same or similar data sets, the value derived using the Keplerian mask approach is used. Most fluxes derived from data sets with similar sensitivities are consistent across the different methods used in the literature, but a few varied by up to ∼2× (e.g., HCN in IM Lup and C¹⁸O J = 2–1 in V4046 Sgr). The uncertainty values listed in Table 2 include flux calibration uncertainties of 10%–15% as specified in the literature referenced in the final column of the table (this systematic uncertainty was added in quadrature to previously reported uncertainties if not already included).

The fluxes shown in the analyses in this work (Section 3) are for the J = 3–2 transition for each molecule. When only ¹³CO and C¹⁸O J = 2–1 fluxes are available, they are multiplied by a factor of 2× as a conversion to approximate the J = 3–2 fluxes. The accurate value for this conversion is uncertain. For optically thin emission in local thermodynamic equilibrium (LTE), this conversion would be ∼4–5× for 20–30 K gas. We base our value on the observed ratios of C¹⁸O (3–2)/(2–1) that are available for three sources in our sample. These values are in the range of ∼1–3, not well described by an optically thin LTE analysis. The line fluxes in the following analyses (Section 3) are corrected for differences in source distances by scaling to a uniform distance of 160 pc.

First, we present the disk-integrated line fluxes measured using the method described in Section 2.3 for sources 1–7 in Table 2. Figures 1 and A1 show the spectra and velocity-integrated line intensity (moment 0) maps, respectively, for all seven sources. We detected N₂H⁺ in IRAS 04302+2247 and GG Tau A, in addition to the previously identified N₂H⁺ emission in two of the Upper Sco sources (Anderson et al. 2019). HCO⁺ is detected in all seven sources, while HCN is detected in five. ¹³CO and C¹⁸O are both detected in all three Upper Sco disks. ¹³CO is detected in all four of the Taurus disks, but C¹⁸O is only clearly detected in two. Note that the differences in the detection rates between the Upper Sco and Taurus samples are due to the difference in requested sensitivity of the observations for each survey (see Sections 2.1–2.2). For the Taurus disks, the main isotopologue of CO is also detected in all four disks (see Figure A2). Additional molecules detected but not discussed in this work are CN in all three Upper Sco disks and C₂H and H₂CO in GG Tau A.

With the combined list of newly measured disk-integrated line fluxes and those from the literature, we search for any relationships between the line fluxes (or flux ratios) and other observed stellar and disk properties. The references for the stellar and disk properties used in this analysis are provided in Table 1 and Section 2.4. To provide an appropriate comparison of the flux values, fluxes are scaled to a common distance of 160 pc for these analyses. Fluxes are multiplied by $\tfrac{{d}^{2}}{{(160\,{pc})}^{2}}$ where d is the source distance in parsecs from Table 1. Using the Python linear regression tool linmix ¹² based on the method of Kelly (2007), we test for linear correlations between pairs of parameters. Kelly (2007) produced a Bayesian method that takes into account measurement errors in two variables and nondetections in one variable in linear regression. We apply this method to investigate trends with fluxes and normalized fluxes where we have nondetections in one variable. The method uses Markov Chain Monte Carlo sampling to approximate the posterior distribution. The computed linear correlation coefficient ranges from −1 for a strong negative correlation to 1 for a strong positive correlation, with a value of zero indicating no correlation between the two variables. We report the median correlation coefficient from the linmix runs with one standard deviation as the listed uncertainties. A summary of the computed correlation coefficients is provided in Figure A5.

3.2.1. Line Fluxes versus Dust Mass

3.2.2. Line Fluxes versus Other Stellar and Disk Properties

In addition to the disk dust mass comparison, we also compared the molecular fluxes to other observed source properties from the literature. Through these comparisons we can investigate which source properties are influencing the chemical behavior of the disk and determine disk-integrated molecular fluxes. Data are plotted based on availability, thus not all panels contain the full 20 disks in the sample. Error bars are not included for the following parameters: M_star, L_star, M_acc, or R_cont,68. See Manara et al. (2022) for discussion on the uncertainties involved in acquiring these values.

Anderson et al. (2022) found no trends with stellar mass or luminosity. Our larger sample yields similar results (Figure 3, columns (1)–(2)). As shown in the Appendix, Figure A3, normalizing the molecular line fluxes by M_dust results in a partially decreasing trend with stellar mass and luminosity but the median correlation coefficients are not more than 1–2 standard deviations above zero. This is most noticeable for HCO⁺. The molecular fluxes show generally increasing trends with the mass accretion rate, similar to their behavior with dust mass (Figure 3, column (3)). This is expected given the known relationship between the mass accretion rates and dust masses of protoplanetary disks (Manara et al. 2016; Mulders et al. 2017; Manara et al. 2020). However, given the smaller number of data points and some additional scatter, the linear correlation coefficients are less significantly different from zero. When normalizing the molecular line fluxes by M_dust or the C¹⁸O flux, there are no significant trends with M_acc (see the Appendix, Figure A3).

Zoom In Zoom Out Reset image size

Similar to the behavior with M_dust, molecular fluxes for all four species generally increase with the outer radius of the (sub-)millimeter dust disk (Figure 3, column (4)). This is consistent with the idea that larger disks have more gaseous material and therefore stronger molecular line emission. However, where line emission becomes optically thick this could also be a trend in disk temperature. Perhaps larger disks have more material, pushing emitting regions of optically thick lines to the hotter layers closer to the disk surface. The median correlation coefficients show positive linear correlations between the molecular fluxes and R_cont,68. When normalizing the molecular line fluxes by M_dust or the C¹⁸O flux, there are no significant trends with R_cont,68 (see the Appendix, Figure A3).

3.2.3. Comparison among Line Fluxes

To investigate potential chemical links between molecular species, we compare pairs of molecular line fluxes from different species against each other in Figure 4. Anderson et al. (2022) found a tentative positive correlation among N₂H⁺, HCN, and HCO⁺ fluxes in their set of Lupus sources (8–14 in Table 2). These molecular fluxes were not correlated with their ¹³CO J = 3–2 fluxes. Limited C¹⁸O detections across their sample prevented a comparison with C¹⁸O. In this sample of 20 sources, we find at least a tentative positive correlation in fluxes for all pairs of molecular species. The correlation coefficients computed using linmix (excluding any cases where both line fluxes are upper limits) show strong positive correlations between all combinations of ¹³CO, C¹⁸O, HCO⁺, and N₂H⁺. For HCN, we find strong positive correlations with HCO⁺ and N₂H⁺ and tentative positive correlations with ¹³CO and C¹⁸O. Correlation coefficients are provided in Figure 4. But some of these strong correlations may break down for fainter line fluxes where constraints are lacking due to the observational bias toward brighter disks.

Zoom In Zoom Out Reset image size

Correlations could be driven by physical or chemical relationships between the molecular species, or more likely, a combination of both. The strong correlation between ¹³CO and C¹⁸O is unsurprising because they both track the same underlying CO abundance. In other cases, chemistry can provide explanations for positive correlations between some pairs of molecular species: CO is a precursor to the formation of HCO⁺, while N₂H⁺ and HCN are both carriers of N. However, the chemical relationship between N₂H⁺ and CO, where N₂H⁺ is readily destroyed when CO is abundant (Qi et al. 2013), would suggest they should be anticorrelated. This might be an indicator that physical disk properties, such as gas mass or temperature, are playing a larger role in determining these disk-integrated molecular fluxes. The level of ionization in the disks may also be an important factor, in particular for the fluxes of molecular ions N₂H⁺ and HCO⁺.

The optical thickness of the line emission is another important factor to consider. The molecular line flux correlations generally appear less strong at lower flux values, below ≲1 Jy km s⁻¹ at a distance of 160 pc. This could be a sign of diverging pathways of chemical evolution and gas-phase chemical abundances among the fainter (typically smaller and less massive, see Figures 2–3) disks. Alternatively, the change could represent the transition from optically thin to thick line emission, where optically thick emission is mainly tracking the gas temperature rather than the molecular column density. This would also indicate a correlation in the gas temperatures within the emitting layer of different molecular species. Isotopologue ratios (e.g., ¹³CO/C¹⁸O and HCN/H¹³CN) have revealed that emission from ¹³CO and the dominant isotologues of other molecular carriers are optically thick in a number of protoplanetary disks (e.g., Schwarz et al. 2016; Bergner et al. 2019).

Figure 5 compares the fluxes of all four molecular species (N₂H⁺, HCN, HCO⁺, and C¹⁸O) against each other to give a more complete view of the variations in chemical tracers in each disk. Given the strong correlation between ¹³CO and C¹⁸O, only C¹⁸O is included to simplify the plot. HCO⁺ makes up around 35%–55% of the total molecular flux. It is generally the largest flux contributor of the four molecular species, with the exception of the C¹⁸O-dominant disks Sz 68 and Sz 65 and the HCN-dominant disks J1608 and J1609. The relative HCO⁺ flux is particularly high in J1614, making up more than 80% of the total molecular flux. The spectral shape of the HCO⁺ flux is also an anomaly in this data set (see Figure 1), suggesting special circumstances in this case. Perhaps this disk is displaying an outflow or temporary flare in HCO⁺ emission (Cleeves et al. 2017; Waggoner & Cleeves 2022).

Zoom In Zoom Out Reset image size

The relative HCN flux is highly variable, ranging from less than a few to over 40% of the total molecular flux. It is generally higher in disks where the fraction of C¹⁸O flux is low. The relative C¹⁸O flux also spans a large range of values from less than a few to over 60% of the total molecular flux. The disks with the largest C¹⁸O fractions belong to younger star-forming regions: Taurus, Lupus, and Oph (left side of Figure 5). But overall, the fraction of C¹⁸O is variable across the younger disks. In contrast, the C¹⁸O fractions are all <10% in the older disks (rightmost five disks in Figure 5). Across all sources, N₂H⁺ makes up less than 15% of the total molecular flux. The lowest N₂H⁺ fractions are less than a few percent, but upper limits prevent us from determining the full range and any additional patterns in N₂H⁺ fractions. For at least the cases of Lupus and Upper Sco, disks from the same star-forming region display different relative fractions of molecular fluxes. This may be indicative of a range of gas chemical compositions coexisting within a single star-forming region.

The total disk gas mass is an essential parameter regulating planet formation mechanisms and the number of planets that can form in a given system. Disk masses are dominated by H₂ gas, which is unobservable at cold temperatures throughout the bulk of the disk. As such, indirect tracers, most often submillimeter dust and CO isotopologue emission, are traditionally used to provide disk gas mass estimates. While both are readily observable and useful for a first look, the accuracy of the resulting gas mass estimates ultimately depends on the accuracy of the gas/dust or CO/H₂ conversion factor. The dust may evolve separately from the gas over time, resulting in uncertainty in the gas/dust ratio. CO/H₂ can also vary spatially and temporally within the disk and among different systems (e.g., Schwarz et al. 2016; Zhang et al. 2019, 2021). Measuring CO/H₂ values in protoplanetary disks requires at least one additional measured quantity beyond CO isotopologue emission, preferably one that is tied to the H₂ content of the gas. Hydrogen deuteride (HD) measurements enabled by the Herschel Space Observatory (Bergin et al. 2013; McClure et al. 2016) revealed CO/H₂ up to 100× below the typically assumed interstellar value ∼10⁻⁴ in a few protoplanetary disks. Since the end of Herschel's mission, there have not been any observatories equipped for the necessary HD measurements. Meanwhile, developments have been made in using N₂H⁺ in combination with CO to provide constraints on both the CO/H₂ and H₂ mass (Anderson et al. 2019, 2022; Trapmanet al. 2022). N₂H⁺ is readily observable by ground-based radio observatories including ALMA and the SMA.

We use our data to compare disk gas mass estimates based on measurements of (sub-)millimeter dust, CO gas, and from combinations of molecular gas tracers. In Figure 6, we plot an initial look at how disk gas masses vary with age across our selected sample. The top panel shows gas masses as equal to 100 × M_dust. Dust masses are computed using (sub-)millimeter continuum emission as described in Section 2.4. More than half of this sample lie above 0.01 M_⊙, with the remaining sources extending to values 1–2 orders of magnitude lower. The Lupus sources represent values across the full range, whereas the Taurus sources and Upper Sco sources are each grouped around a narrower range of disk dust mass values. Previous comparisons across star-forming regions find decreasing median dust masses with age (e.g., Ansdell et al. 2017; Cieza et al. 2019; Villenave et al. 2021). The sources in this sample were not selected to be representative of the population in their respective star-forming regions. Selection is often biased toward brighter sources, therefore this comparison is not representative of typical sources across the star-forming regions included in this work.

Zoom In Zoom Out Reset image size

The second panel in Figure 6 shows gas masses based on CO isotopologue observations. These masses were computed using the fits functions provided by Miotello et al. (2016, 2017) from their large grid of models. Masses were based on C¹⁸O fluxes from either the J = 3–2 or J = 2–1 transitions using the appropriate functions for each from Miotello et al. (2016, 2017). Where fluxes for both transitions are available, the one with a higher signal-to-noise ratio was used. For sources where C¹⁸O has not been detected, ¹³CO fluxes and appropriate functions are used instead. With the exception of the young source IRAS 04302+2247, the CO-based gas masses are lower than the dust-based values by roughly 1–2 orders of magnitude. While IRAS 04302+2247 is part of the Taurus star-forming region and that age is used here, it may represent an earlier evolutionary stage as it is classified as a Class I source (Lucas & Roche 1997; Garufi et al. 2021). The distinction between gas masses derived from (sub-)millimeter dust and those derived from CO isotopologues is well documented in large-population ALMA studies of multiple star-forming regions (e.g., Ansdell et al. 2016; Long et al. 2017; Miotello et al. 2017). As discussed in these previous works, the differences between the gas mass estimates derived from these two methods could be the result of gas/dust ratios and/or CO/H₂ abundances decreasing over time as the disks evolve.

Anderson et al. (2019, 2022) and Trapmanet al. (2022) have shown that by using a combination of N₂H⁺ and CO observations, we can place constraints on both the total disk gas mass and CO/H₂. In this way, we can distinguish between the effects of depleted gas/dust ratios versus depleted CO/H₂. The extensive modeling analysis needed to constrain CO/H₂ and H₂ masses for this sample of disks is beyond the scope of this work. But to give a rough idea of how this might affect our inferred gas masses, we perform a simple scaling of the CO-based total H₂ gas masses from the middle panel of Figure 6. These values assume a CO/H₂ abundance ∼ 10⁻⁴. We increase these values by a factor corresponding to the CO/H₂ abundances derived from the observed N₂H⁺/C¹⁸O flux ratios for each source. For CO/H₂, we used the modeled relation between the N₂H⁺ (3–2)/C¹⁸O (2–1) flux ratio and the CO/H₂ abundance for the TW Hya model with a cosmic-ray ionization rate of 10⁻¹⁸ s⁻¹ from Trapmanet al. (2022). In reality, this relation will depend on the individual physical/chemical structure of each disk (Trapmanet al. 2022) and further exploration of how much this relation varies over a wider parameter space is needed.

Based on this approximation, the gas masses for a subset of our sample were increased by factors ∼5–135 relative to the CO-based values. Some disks have N₂H⁺/C¹⁸O flux ratios that are consistent with a CO/H₂ of 10⁻⁴, including disks where only an upper limit on this flux ratio is available, so their gas mass estimates were not changed from the CO-based values. The gas mass value for IRAS 04302+2247 was also not increased given its already high value. The final result is a larger spread in the range of gas masses for each age group (Figure 6, bottom panel). The additional constraints on CO/H₂ increase the gas mass estimates for some disks to values that are closer to the dust-based masses (Figure 6, top panel). The mass estimates based on HD emission (Bergin et al. 2013; McClure et al. 2016; Trapmanet al. 2022) also suggest that the CO-based values are underestimating the total amount of H₂ gas present. Our conclusions are limited by the size of and selection biases present in the current sample. Ultimately, more data are needed to explore the evolution of disk gas masses, especially data from older sources (ages > 5 Myr).

This data set includes 20 disks spread across various birth environments and ages. Within this sample, we find examples of disks with very different disk-integrated molecular flux ratios, even for disks residing in the same star-forming region. In Figures 2–3 color is used to distinguish sources from different star-forming regions or moving groups. In this case, any apparent disparities among disks from different star-forming regions are likely caused by differences in the sample selection criteria used for the different sets of observations. As a result, the selected Taurus and individual sources (AS 209, TW Hya, and V4046 Sgr) tend to have larger and brighter disks relative to those from Lupus and Upper Sco. When normalized by M_dust or C¹⁸O, the molecular fluxes fall within generally the same range within 1–2 orders of magnitude regardless of star-forming region—with a few exceptions. There are two sources from Lupus (Sz 65 and Sz 68, both appearing as upper limits in Figures A3–A4) for which constraints on the normalized N₂H⁺ and HCN values fall well below the rest of the sample. This could be an indication that these two sources have a distinct disk gas chemical composition from the others.

For the comparisons made in Figures 2–4, we see no obvious distinctions between the younger sources (the AS 209, Taurus, and Lupus disks; 1–3 Myr old) and the older sources (the Upper Sco, TW Hya, and V4046 Sgr disks; >5 Myr old). Although, for at least 4/5 of the older sources, their N₂H⁺, HCN, and HCO⁺ fluxes normalized by M_dust or C¹⁸O are among the highest values in the sample. There are hints of more uniformity among these disks (Figure 5). The fifth source, J1614, stands out by having a relatively large amount of HCO⁺ emission with an unusual line shape. However, even when excluding HCO⁺, J1614 still differs from the other >5 Myr old sources due to its low HCN/C¹⁸O flux ratio.

It should be noted that different age ranges are not uniformly sampled in this work. The older populations may suffer the greatest selection bias, often toward brighter objects that still have molecular gas. Overall, aged populations are more poorly sampled in multimolecular observations. Furthermore, this comparison may be muddled by uncertainty in the individual source ages and potential age variations within the star-forming regions (e.g., Krolikowski et al. 2021). A more uniform unbiased survey of disks across different age brackets is needed to investigate fully the trends in the molecular fluxes with disk age.

Sample bias and incompleteness largely affect this current analysis. The selected sources tend to be among the largest and brightest in their respective star-forming regions. The star-forming regions included in this comparison are also not equally sampled, with fewer sources particularly for the older populations. Furthermore, because the sample is compiled from different observing programs, the sensitivity is not the same for all sources or molecular lines. Additional uncertainty is introduced when trying to make comparisons using different transitions of CO isotopologues (see Section 2.4). Current ALMA Large Programs in progress, namely AGE-PRO: The ALMA Survey of Gas Evolution in PROtoplanetary Disks and The ALMA Disk-Exoplanet C/Onnection (ALMA DECO), will provide a more uniform sampling and analysis of multiple molecular lines across different star-forming regions.

Within our sample are four systems that contain multiple stars. Given the current data set, there are no obvious flux trends that set these systems apart from the single star systems. Although it should be noted that the data sets of stellar and disk properties are not complete for our sample, especially for the stellar multiple systems. In addition, using spatially resolved observations that can distinguish between different disk components in these complex systems (e.g., disks around individual stars versus disks around binaries) may be necessary to make more meaningful comparisons.

There are many relevant parameters that could affect the observed molecular line fluxes that have not been explored in this work. This is largely due to limitations in the data readily available for this combined sample. Additional stellar properties such as UV and X-ray spectra, stellar activity, and magnetic field strength could be explored for connections in the future. Fluxes will also depend on the individual density and temperature distribution of each disk. Determining whether these physical conditions are relatively consistent across populations (or subgroups within populations) or are unique to each individual source is an important step in understanding disk evolution and planet-forming environments.

Previous works have explored surveys of disk-integrated molecular line fluxes for different samples of disks and molecular species. Pegues et al. (2023) find that molecular line fluxes generally increase with continuum fluxes across their sample of T Tauri and Herbig disks (stellar masses of 0.15–2.5 M_⊙), consistent with our trends in line flux with M_dust. Our results agree with Pegues et al. (2021) in finding no correlation between the flux ratio of HCN/C¹⁸O and M_star (Figure A4). Pegues et al. (2023) also show that HCO⁺/C¹⁸O flux ratios (when both species are detected) are relatively flat with continuum flux and stellar luminosity for their T Tauri sample, similar to our findings.

While we find positive correlations between the disk-integrated line fluxes of all our molecular species pairs, this is not always the case. van Terwisga et al. (2019) found that CN fluxes do not strongly correlate with ¹³CO in disks from the Lupus star-forming region. Furthermore, the CN fluxes do not strongly correlate with the submillimeter continuum fluxes. However, for a subset of the sample (Miotello et al. 2019) show a strong correlation between CN and C₂H fluxes and between both species and the stellar luminosity. The abundances of CN and C₂H are largely affected by photochemistry, which might explain their correlation.

With a diverse set of 14 disks across a range of ages, stellar masses, and stellar luminosities (including some overlap with our sample), Bergner et al. (2019) find no correlation between disk-integrated HCN and C¹⁸O fluxes when both are normalized to the continuum flux. Here we find only a tentative trend between the HCN and C¹⁸O fluxes relative to strong positive correlations between other molecular species. In contrast, Pegues et al. (2021) report a significant positive correlation between HCN and C¹⁸O fluxes, although the correlation is not as strong as those between other molecular species pairs (C₂H and HCN, H₂CO and HCN). In addition, they find a positive correlation between C¹⁸O and M_star, whereas we find no correlation. The sample of Pegues et al. (2021) includes disks around five low-mass stars (0.14–0.23 M_⊙) alongside solar-type T Tauri and Herbig Ae disks. In comparison, our sample is more concentrated toward stellar masses ∼ 0.5–1 M_⊙ and does not include such low-mass stars (<0.25 M_⊙). Current trends are determined with limited numbers of disks and line detections, thus more sensitive observations of larger and more diverse samples of disks are ultimately needed to explore these discrepancies as well as to verify existing trends.

Molecular abundances are sensitive to various physical and chemical properties of protoplanetary disks (Öberg et al. 2023). Molecular emission is further sensitive to excitation conditions, predominately how the abundance distribution aligns with the disk temperature and density structures. Physical/chemical disk models have explored how the molecular fluxes of some of our observed species depend on various stellar and disk parameters. Miotello et al. (2016) show that disk-integrated ¹³CO and C¹⁸O fluxes are sensitive to the disk mass and outer radial extent. The modeling of Boyden & Eisner (2023) shows HCO⁺ J = 4–3 emission increases with the outer disk radius, but decreases with increasing stellar mass and disk dust mass. Fedele & Favre (2020) find HCN J = 4–3 emission is insensitive to changes in stellar and physical disk parameters, including the total stellar luminosity, but sensitive to elemental ratios of C/O (particularly as C/H is reduced) and N/H. But note that the J = 3–2 HCO⁺ and HCN transitions observed here have lower upper state energies and likely trace a colder and therefore deeper vertical layer of the disk. Anderson et al. (2022) show that ¹³CO, HCO⁺, HCN, and N₂H⁺ J = 3–2 emission generally increase together for increasing gas/dust ratios, but diverge with volatile depletion (C/H and O/H). As C and O are removed, ¹³CO and HCO⁺ emission decreases, while HCN and N₂H⁺ emission increases.

While our observed trends agree with some of the model results (e.g., increasing line emission with disk radius), direct comparison with existing model grids is difficult. All model results depend on the assumptions made about the interplay of different disk parameters, for example, the relationship between gas and dust components of the disk and between disk mass and outer radius. The assumed relationships may not be reflective of our observed sample of disks. While beyond the scope of this work, a model grid exploring N₂H⁺, HCN, HCO⁺, C¹⁸O, and ¹³CO J = 3–2 emission across the parameter space relevant to the star and disk properties of this observed sample would aid in the interpretation of the trends we find.