An Evolutionary Study of Volatile Chemistry in Protoplanetary Disks

The sample of protostellar disk candidates consists of five low-mass protostars in the Serpens cluster, described in Table 2. Each was identified as hosting a candidate disk based on nonzero flux at >50 kλ uv distances, corresponding to physical distances <1800 au, in 230 GHz continuum observations (Enoch et al. 2011). Ser-emb 1, 7, and 8 are classified as Class 0 sources, and Ser-emb 15 and 17 as Class I (Enoch et al. 2009). Enoch et al. (2009) estimated that the Class 0 lifetime in Serpens is about 0.2 Myr, assuming a Class I lifetime of 0.54 Myr as found in Evans et al. (2009). A study of the organic molecule emission toward this source sample revealed that Ser-emb 1, Ser-emb 8, and Ser-emb 17 host hot corinos based on the detection of a warm, rich organic chemistry in the protostellar core (Bergner et al. 2019a).

Full observational details can be found in Bergner et al. (2019a). Briefly, targets were observed by ALMA as part of project 2015.1.00964.S from 2016 May to June. Two Band 6 spectral setups (217–233 GHz and 243–262 GHz) were covered. This project makes use of spectral windows covering transitions of CO, C¹⁸O, C₂H, H¹³CN, HC¹⁵N, and DCN. All spectral windows have a channel width of 122 kHz, corresponding to ∼0.14–0.17 km s⁻¹ across the range of frequencies.

The full Class II disk sample is described in Bergner et al. (2019b), and source properties are shown in Table 2. The sample spans a range of physical properties, including age (0.4–14 Myr), stellar mass (0.5–2 M_⊙), and stellar luminosity (0.2–25 L_⊙). In addition to the C₂H, HCN, and C¹⁸O observations presented in Bergner et al. (2019b), here we include previously unpublished observations of DCN, HC¹⁵N, and ¹³CO lines toward a subset of the full disk sample.

The DCN 3–2 transition was observed toward HD 143006, J1604, J1609, J1612, and J1614 as part of ALMA project 2015.1.00964S; and toward CI Tau, DM Tau, and DO Tau as part of project 2016.1.00627.S. The HC¹⁵N 3–2 transition was observed toward CI Tau, DM Tau, DO Tau, LkCa 15, and MWC 480 as part of project 2016.1.00627.S. For sources with DCN or HC¹⁵N observations, we also reanalyze the HCN and H¹³CN 3–2 observations originally presented in Bergner et al. (2019b) to facilitate a consistent analysis of isotopic ratios. Lastly, we include previously unpublished observations of ¹³CO 2–1 toward J1609, J1612, and J1614, for which C¹⁸O emission was not detected. Observations took place from 2016 May to June for project 2015.1.00964S and in 2016 December for project 2016.1.00627.S. In both cases, two Band 6 spectral setups covered similar frequency ranges and spectral resolutions as described for the Serpens observations.

Following pipeline calibration of the ALMA data, we performed self-calibration using the combined continuum emission from spectral windows within a given sideband. Depending on the continuum brightness, up to two rounds of phase self-calibration were performed. Solution intervals ranged from 12.1–42.35 s for the Class 0 sources, 6.05–12.1 s for the Class I sources, 60.5–84.7 s for the Class II Upper Sco sources, and 6.05–48.4 s for the remaining Class II sources. Continuum imaging was performed using Briggs weighting with a robust factor of 0.5. The self-calibration solutions were applied to the spectral line data following continuum subtraction in the uv plane. Calibrated and continuum-subtracted data were CLEANed to a 4σ noise threshold in CASA 5.5.0 using Briggs weighting and the automasking task auto-multithresh in tclean. We applied the standard automasking parameters for short baseline 12 m line data (sidelobethreshold = 2.0, noisethreshold = 4.5, minbeamfrac = 0.3, negativethreshold = 15.0, lownoisethreshold = 1.5). Automasking results were verified by comparing to hand-masked images for select lines including extended and compact emission morphologies, and bright and weak emission strengths. For the Class 0/I sources, images were generated with a channel width of 0.25 km s⁻¹. We used a robust factor of 0.5 for all molecules except C₂H, for which we used a robust factor of 2.0 due to its weak and diffuse emission. For the Class II sources, imaging was performed with a channel width of 0.5 km s⁻¹. A robust factor of 2.0 was used for all disks except for HD 143006 and J1604, whose bright emission permitted robust factors of 1.0 and 0.0, respectively.

3.1.1. Line Detections and Velocity Structures

Figure 1 shows an overview of the continuum and line emission observed toward the five Serpens protostellar disk candidates. As discussed in Bergner et al. (2019a), the continuum emission is more extended in the Class 0 sources (Ser-emb 1, 7, and 8) than in the Class I sources (Ser-emb 15 and 17), reflecting that they are less evolved and more deeply embedded in the envelope.

Zoom In Zoom Out Reset image size

For line observations, the grayscale map represents the intensity integrated across the entire velocity range where emission is detected. The red and blue contours outline the emission originating from red- and blueshifted velocities in the range indicated for each panel. Many velocity ranges were tested for each line in an effort to identify coherent velocity structures in the sources, and those shown represent the clearest structures that could be distinguished. We refer to the total intensity maps and the velocity structure maps as "Full" and "R/B," respectively. Velocity ranges and moment zero rms values for both sets of maps are listed in Table 3. We note that the spatial resolution of these observations (∼0.5''; Table 3) is not high enough to definitively identify the presence of Keplerian disks (see Martín-Doménech et al. 2019). Still, we can identify rotational structures in the compact molecular line emission that are perpendicular to the outflow direction, which we identify as either a disk or the disk-forming region of the rotating inner envelope. The following sections present the line emission velocity structures identified for each source.

Table 3.  Class 0/I Line Observations

Transition	Beam Dim.	Chan. Rmsa	Vel. Rangeb (km s−1)		Mom. 0 Rmsb (mJy beam−1 km s−1)		Fluxc
	('' × '')	(mJy beam−1)	R/B	Full	R/B	Full	(mJy km s−1)
Ser-emb 1
C18O 2–1	0.62 × 0.48	5.4	∣0.75–1.50∣	−1.75–2.00	2.8	6.1	154 ± 16
H13CN 3–2	0.51 × 0.43	5.0	∣1.25–2.00∣	−10.75–7.25	2.6	15.3	233 ± 30
DCN 3–2	0.62 × 0.50	5.2	∣0.75–1.50∣	−2.25–2.50	2.8	6.2	139 ± 15
HC15N 3–2	0.51 × 0.42	5.6	∣0.75–1.50∣	−1.50–2.50	2.9	6.2	55 ± 9
C2H N = 3–2, J = $\tfrac{5}{2}$ – $\tfrac{3}{2}$	0.63 × 0.51	5.0	∣0.00–0.50∣	−1.00–0.75	2.3	3.7	<21
Ser-emb 7
C18O 2–1	0.61 × 0.48	5.5	∣0.00–0.25∣	−3.00–3.00	1.9	7.7	74 ± 9
H13CN 3–2	0.51 × 0.43	5.0	∣0.00–1.25∣	−2.25–1.00	3.3	5.1	27 ± 3
DCN 3–2	0.62 × 0.51	5.2	∣0.00–1.25∣	−2.25–0.75	3.6	5.4	15 ± 2
HC15N 3–2	0.51 × 0.43	5.5	∣0.00–1.25∣	−2.25–0.75	3.6	5.4	<10
C2H N = 3–2, J = $\tfrac{5}{2}$ – $\tfrac{3}{2}$	0.63 × 0.51	4.9	∣0.25–0.50∣	−2.50–0.50	1.8	5.0	13 ± 2
Ser-emb 8
C18O 2–1	0.63 × 0.51	5.3	∣0.75–2.00∣	−3.00–2.00	3.5	7.1	886 ± 240
H13CN 3–2	0.52 × 0.44	5.1	∣3.50–6.00∣	−8.50–5.75	4.6	12.7	426 ± 50
DCN 3–2	0.61 × 0.51	5.6	∣1.50–6.00∣	−3.75–3.75	6.4	7.9	190 ± 20
HC15N 3–2	0.51 × 0.43	5.8	∣1.50–6.00∣	−3.50–4.25	6.8	8.8	212 ± 40
C2H N = 3–2, J = $\tfrac{5}{2}$ – $\tfrac{3}{2}$	0.63 × 0.51	4.8	∣0.00–0.50∣	−3.00–1.00	2.2	5.6	12 ± 4
Ser-emb 15
C18O 2–1	0.63 × 0.48	5.2	∣0.50–1.75∣	−1.50–2.75	3.6	6.5	62 ± 7
H13CN 3–2	0.51 × 0.42	4.9	∣0.25–2.00∣	−1.75–1.50	3.8	5.0	52 ± 6
DCN 3–2	0.63 × 0.50	4.9	∣0.25–2.00∣	−1.25–2.00	3.9	5.3	39 ± 6
HC15N 3–2	0.51 × 0.42	5.7	∣0.25–2.00∣	−0.50–1.50	4.3	4.5	<36
C2H N = 3–2, J = $\tfrac{5}{2}$ – $\tfrac{3}{2}$	0.63 × 0.51	4.9	∣0.00–1.00∣	−1.50–1.00	3.1	4.6	15 ± 4
Ser-emb 17
C18O 2–1	0.63 × 0.50	5.3	∣1.00–3.00∣	−3.00–3.00	4.5	7.5	230 ± 26
H13CN 3–2	0.52 × 0.44	5.0	∣1.00–3.50∣	−4.25–6.25	4.8	10.1	112 ± 19
DCN 3–2	0.63 × 0.53	5.4	∣0.50–3.50∣	−2.25–4.50	5.2	7.5	72 ± 9
HC15N 3–2	0.53 × 0.44	5.5	∣0.50–3.50∣	−2.00–4.25	5.5	8.1	85 ± 10
C2H N = 3–2, J = $\tfrac{5}{2}$ – $\tfrac{3}{2}$	0.63 × 0.51	5.0	∣0.00–0.50∣	−1.25–1.25	2.3	4.7	<35

Notes.

^aFor 0.25 km s⁻¹ channels. ^bVelocity ranges and moment zero rms values correspond to two different sets of moment zero maps: R/B maps include a subset of channels with emission, chosen to reveal velocity structures, and are shown as red/blue contours in Figure 1. Full maps include the full range of channels with emission, shown in grayscale in Figure 1. Fluxes listed in this table correspond to the Full velocity range. ^cExtracted within a mask defined by the DCN emitting region (see Section 3.1.7). 3σ upper limits are reported for non-detections. Uncertainties consist of Gaussian fit uncertainties, added in quadrature with a 10% calibration uncertainty.

Download table as:  ASCIITypeset image

3.1.2. Ser-emb 1

3.1.3. Ser-emb 7

3.1.4. Ser-emb 8

3.1.5. Ser-emb 15

3.1.6. Ser-emb 17

3.1.7. Spectral Extraction

As seen in Figure 1, the structures of these embedded sources can be quite complex. Our aim is to isolate the molecular component of the protostellar disks or disk-forming regions. We therefore extract spectra from a region defined by the DCN emitting region: of all the lines covered here, DCN seems to most consistently trace compact disk-like emission. This is especially apparent in Ser-emb 15, in which DCN traces a flattened, visibly disk-like feature. For Ser-emb 7, 8, 15, and 17 we adopt a 6σ cutoff for the DCN emission mask, and in Ser-emb 1 we choose an 8σ cutoff to isolate the more compact central region. σ values correspond to the rms of the DCN R/B maps listed in Table 3. For Ser-emb 7 we restrict the DCN mask to exclude the northern molecular emission peak (see Section 3.1.3). The angular sizes of the resulting masks are 0.40, 0.18, 0.43, 0.46, and 0.37 square arcseconds for Ser-emb 1, 7, 8, 15, and 17, respectively.

The resulting spectral line profiles are shown in Figure 2. We fit each line as a sum of Gaussians corresponding to different possible velocity components: a "medium" component (Gaussian width constrained to [0.1, 3] km s⁻¹) is included for all lines, and other components are added when needed to reproduce the overall profile: "absorption" (Gaussian width [0.1, 1] km s⁻¹ and amplitude <0), "broad" (Gaussian width [2, 10] km s⁻¹), or "narrow" (Gaussian width [0.1, 0.4]). With this routine we obtain very good fits to the spectra, as can be seen in Figure 2. We use the MCMC package emcee (Foreman-Mackey et al. 2013) to sample posterior distributions of the Gaussian fit parameters. By fitting for the absorption feature, we can recover the original flux density prior to self-absorption from the envelope/cloud, though we note that in the case of C¹⁸O toward Ser-emb 8 and Ser-emb 17 there is some degeneracy between the amplitudes of the emission and absorption features.

Zoom In Zoom Out Reset image size

We wish to isolate line fluxes originating from the disk-like structures, and excluding emission from the parent cloud or outflows. We assume that narrow velocity components (emission or absorption, with Gaussian widths ≲0.5 km s⁻¹) arise from the cloud. We also assume that broad emission components are associated with outflows. This is supported by the line profiles of DCN (which is spatially most closely associated with compact disk-like structures), which do not include a broad emission component. Thus, we assume that the flux coming from the disk-like region corresponds to the medium velocity component. We note that if broad emission does in fact originate from the disk structure, this would increase the fluxes by a factor of ∼2–3 for H¹³CN in Ser-emb 1, 8, and 17 and for HC¹⁵N in Ser-emb 1 and 8. For non-detections of HC¹⁵N and C₂H, we calculate upper limits adopting the source's line width for H¹³CN and C¹⁸O, respectively. The fluxes resulting from this treatment are listed in Table 3. Gaussian fit uncertainties are added in quadrature with a 10% calibration uncertainty, i.e., the estimated absolute calibration uncertainty for ALMA (e.g., Remijan et al. 2020), and propagated through all subsequent analysis.

Figure 3 shows the moment zero maps for the HCN isotopologues observed toward the Class II disks, along with deprojected and azimuthally averaged intensity profiles. ¹³CO moment zero maps for J1609, J1612, and J1614 can be found in Appendix A. Keplerian masking was used to generate the moment zero maps, assuming the star properties listed in Table 2, and disk inclinations and position angles listed in Appendix A (Table A2). Mask outer radii are chosen to encompass the observed emission. A full description of the masking technique can be found in Bergner et al. (2018) and Pegues et al. (2020). Deprojected radial profiles were generated using the disk inclinations and position angles in Table A2.

Zoom In Zoom Out Reset image size

Fluxes are calculated by summing the emission within the moment zero map. We also apply the same Keplerian mask to 200 synthetic image cubes of randomly drawn emission-free channels neighboring the target line. The integrated flux uncertainty is estimated from the standard deviation in the integrated fluxes of these emission-free moment zero maps, added in quadrature with a 10% calibration uncertainty. The moment zero rms varies across the map because each pixel represents the sum of a different number of channels when Keplerian masking is used. We therefore calculate the standard deviation in intensity within each pixel across the 200 emission-free maps, and adopt the median pixel rms as the moment zero map rms. Line observation results for the Class II sources are listed in Appendix A (Table A1).

We consider a line to be detected if the measured flux is >3σ and if emission >3σ is present in the moment zero map. Based on these detection criteria, DCN is detected toward CI Tau, HD 143006, J1604, and J1609; HC¹⁵N toward LkCa 15 and MWC 480; HCN toward HD 143006, J1604, J1609, and J1614; and H¹³CN toward DM Tau, LkCa 15, and MWC 480.

The HCN and H¹³CN lines presented here were imaged with different parameters from Bergner et al. (2019b) in order to be consistent with the DCN and HC¹⁵N images, thus facilitating our analysis of isotopic ratios. The main difference in imaging is in the choice of robust parameter. Here, weak lines are imaged with a robust parameter of 2.0 instead of 1.0, which enabled a detection of H¹³CN toward DM Tau and of HCN toward J1614. We also imaged HCN toward J1604 with a robust of 0.0 since the emission is very bright, allowing us to better resolve the emission morphology in the inner disk. We also note that H¹³CN and HC¹⁵N 3–2 were previously observed toward LkCa 15 and MWC 480 with a similar sensitivity and spatial resolution in Guzmán et al. (2017). We obtain H¹³CN fluxes a factor of ∼1.3 and 1.7 higher compared to Guzmán et al. (2017) for MWC 480 and LkCa 15, respectively, possibly due to our Keplerian masking technique, which reduces the incorporation of emission-free regions into our moment zero images. We find a similar HC¹⁵N flux in MWC 480 but a slightly lower flux in LkCa 15 (by a factor of ∼0.8) compared to Guzmán et al. (2017). This is reflective of the uncertainty inherent in imaging weak molecular line data; our uncertainties are probably underestimated for lines close to the detection threshold.

In disks where both DCN and HCN are detected, their morphologies are qualitatively similar: HD 143006 shows a plateau in intensity toward the disk center, J1604 has a strong emission depletion toward the center, and J1609 is centrally peaked. The inner gap in J1604 is wider in DCN than in HCN, and the intensity is negative within ∼30 au. This negative emission is not seen for channels without line emission, and could be caused by over-subtraction of the continuum where it is coupled to the line emission. That the DCN intensity drops below zero while HCN does not may indicate that it is emitting closer to the midplane. HC¹⁵N and H¹³CN exhibit shallow, extended emission profiles in LkCa 15 and a centrally peaked morphology in MWC 480. HC¹⁵N emission in LkCa 15 is asymmetrical and lacks emission in the northeast region of the disk, whereas H¹³CN emission is more axisymmetric.

With the integrated flux densities $\int {S}_{\nu }d\nu$ from Table 3, we estimate total molecular column densities N_T from:

$$\begin{eqnarray}&&{N}_{T}=\displaystyle \frac{4\pi \int {S}_{\nu }d\nu }{{A}_{{ul}}{\rm{\Omega }}{{hcg}}_{u}}{e}^{{E}_{u}/{T}_{r}}Q({T}_{r}),\end{eqnarray} \tag{ 1 }$$

where A_ul is the Einstein coefficient for the transition, Ω is the angular size of the emitting region, g_u is the upper state degeneracy, E_u is the upper state energy, T_r is the rotational temperature, and Q(T_r) is the temperature-dependent molecular partition function. All spectral line parameters can be found in Table 4. For Ω, we adopt the angular size of the masks used for spectral extraction (listed in Section 3.1.7). The rotational temperatures are not known and cannot be determined from these data alone, so we calculate column densities adopting rotational temperatures of 20, 30, and 50 K. This is informed by the lukewarm temperature profile derived for the embedded disk L1527 (∼20–50 K midplane temperatures; van 't Hoff et al. 2018), though the effect of a much warmer (100 K) rotational temperature is discussed in Section 6.

Table 4.  Spectral Line Data

Molecule	Transition	Frequency	Eup	log(Aul)	gu	Q(30 K)	Referencesa
		(GHz)	(K)	(s−1)
CO	2–1	230.538	16.6	−6.16	5	11	(1)
13CO	2–1	220.399	15.9	−6.22	10	23	(2)
C18O	2–1	219.560	15.8	−6.22	5	12	(2), (3)
C2Hb	N = 3–2, J = $\tfrac{5}{2}$ – $\tfrac{5}{2}$, F = 3–3	262.209	25.2	−5.40	7	59	(4), (5), (6)
	N = 3–2, J = $\tfrac{7}{2}$ – $\tfrac{5}{2}$, F = 4–3 and F = 3–2c	262.004 and 262.006	25.2	−4.28	16	59	(4), (5), (6)
HCN	3–2	265.886	25.5	−3.08	21	43	(7)
H13CN	3–2	259.012	24.9	−3.11	21	44	(8), (9), (10)
DCN	3–2	217.239	20.9	−3.34	21	53	(11)
HC15N	3–2	258.157	24.8	−3.12	7	15	(8), (12)

Notes.

^aAll line parameters are taken from the CDMS catalog (Müller et al. 2001, 2005) based on data from the following: (1) Winnewisser et al. (1997), (2) Klapper et al. (2001), (3) Winnewisser et al. (1985), (4) Müller et al. (2000), (5) Padovani et al. (2009), (6) Sastry et al. (1981), (7) Ahrens et al. (2002), (8) Fuchs et al. (2004), (9) Cazzoli & Puzzarini (2005), (10) Maiwald et al. (2000), (11) Brünken et al. (2004), (12) Cazzoli et al. (2005). ^bThe J = $\tfrac{5}{2}$ – $\tfrac{5}{2}$ transition was observed toward the Class 0/I sources, and the J = $\tfrac{7}{2}$ – $\tfrac{5}{2}$ complex toward the Class II sources. ^cThe hyperfine components are blended, so the line parameters refer to the total J = $\tfrac{7}{2}$ – $\tfrac{5}{2}$ line. The total degeneracy is found by adding the degeneracies of each hyperfine component, and the total Einstein coefficient is found from the degeneracy-weighted sum of Einstein coefficients of each hyperfine component.

Download table as:  ASCIITypeset image

The resulting column densities are shown in Table 5. We find that the choice of rotational temperature has only modest effects on the column density. For subsequent analysis, we adopt the 30 K case as our fiducial value, and incorporate the variation in column density resulting from a 20 or 50 K rotational temperature into the uncertainties.

Column density estimates for C₂H, HCN, and C¹⁸O in the Class II sources are derived using the observations presented in Bergner et al. (2019b), with a few exceptions. First, for CI Tau, DM Tau, DO Tau, LkCa 15, and MWC 480, we do not have ALMA observations of HCN 3–2, so we previously used HCN fluxes measured with the Submillimeter Array for these disks. However, for the comparison presented here, we opt to solve for the HCN column density using H¹³CN as a proxy since HCN suffers from optical depth effects (Bergner et al. 2019b). Second, because C¹⁸O 2–1 is not detected toward J1609, J1612, or J1614, we previously used CO to estimate C¹⁸O fluxes. However, here we use ¹³CO instead, as CO is likely optically thick in these sources. These two changes allow for an improved comparison of the chemistry across our sample.

To estimate disk-averaged column densities, integrated fluxes are extracted from moment zero maps using an elliptical mask. The mask is defined by the maximum radius containing emission at a 4σ level, projected along the same position angle and inclination used to generate the Keplerian mask (Table A2). This treatment allows us to determine the disk-integrated flux and the emission region Ω self-consistently, while avoiding severe beam dilution; note though that these fluxes differ slightly from those integrated across the entire moment zero map, as in Section 3.2. For non-detected lines, we derive upper limits using the continuum major axis (fit at the 3–5σ level) as the mask radius.

Equation (1) was used along with the spectral line parameters listed in Table 4 to estimate column densities from disk-integrated fluxes. As for the Class 0/I sources, we solve for column densities assuming rotational temperatures of 20, 30, and 50 K. The integrated fluxes, emitting regions Ω, and resulting column densities are listed in Table B1 in Appendix B. As in the case of the Serpens sources, the column densities are not very sensitive to the choice in rotational temperature.

As demonstrated in Bergner et al. (2019b), the C₂H and HCN lines covered here are generally optically thick in Class II disks. We therefore adjust the Class II HCN and C₂H column densities following Goldsmith & Langer (1999):

$$\begin{eqnarray}&&{N}_{T}={N}_{T,\mathrm{thick}}\times \displaystyle \frac{\tau }{1-{e}^{-\tau }}.\end{eqnarray} \tag{ 2 }$$

This correction excludes HCN column densities derived from H¹³CN, as these are already likely to be optically thin.

For the disks fit in Bergner et al. (2019b), the average disk-averaged optical depth for the C₂H N = 3–2, J = $\tfrac{7}{2}$ – $\tfrac{5}{2}$ and HCN J = 3–2 lines are 2.1 ± 0.4 and 5.8 ± 1.8, respectively. Uncertainties represent the standard deviations across the disk sample, and are propagated throughout all analyses involving τ-corrected column densities. Note that these values represent the total line optical depth of all hyperfine components, τ_total. In principle, the flux from each hyperfine component should be corrected individually for optical depth, however, in the disk-integrated spectra used in this work, the hyperfine components are blended. Because of this, the optical depth corrections are approximate. For C₂H, the F = 4–3 and F = 3–2 hyperfine components have similar relative intensities (0.572 and 0.428, respectively), and the estimated τ_total is not too high. We therefore assume that the opacity of the blended hyperfine lines can be described by the weighted opacity of the individual components, leading to τ_blend = 0.51 × τ_total. For the C₂H N = 3–2, J = $\tfrac{7}{2}$ – $\tfrac{5}{2}$ τ_total of 2.1 ± 0.4, this results in a τ_blend of 1.1 ± 0.2 for use in Equation (2). In the case of HCN, most opacity (92.6%) in the J = 3–2 transition is contained in the overlapping F = 2–1, 3–2, 4–3, and 2–3 components at 265.886 GHz. These components are sufficiently overlapping that they may be treated as a single line. We assume that the τ_total of 5.8 ± 1.8 can be used directly in Equation (2), though this may over-correct the contribution from the satellite lines (F = 3–3 and F = 2–2), which are less optically thick than the main complex.

We now combine the Class 0/I and Class II column density results to explore correlations between HCN, C₂H, and C¹⁸O across all stages. The size scales that the column densities are extracted from correspond to ∼150 au for the Class 0/I sources, and ∼150–600 au for the Class II sources. For sources with H¹³CN observations (i.e., all Class 0/I and some Class II sources), we convert from H¹³CN to HCN column densities using the local ISM ¹²C/¹³C ratio of 68 ± 15 (Milam et al. 2005). We note that ¹³C fractionation is not expected to be important for HCN in protostellar or disk environments (e.g., Roberts et al. 2002; Hily-Blant et al. 2019). For J1609, J1612, and J1614, C¹⁸O column densities are found from ¹³CO assuming the local ISM ¹⁶O/¹⁸O ratio of 560 ± 25 (Wilson & Rood 1994) and the same ¹²C/¹³C ratio as before.

We also estimate disk-averaged H₂ column densities for our source sample to enable a comparison of molecular abundances. Following Hildebrand (1983) and Andrews & Williams (2005), we convert the 260 GHz continuum flux to a disk mass by:

$$\begin{eqnarray}&&{M}_{d}=\displaystyle \frac{{d}^{2}{F}_{\nu }}{{\kappa }_{v}\zeta {B}_{\nu }(T)},\end{eqnarray} \tag{ 3 }$$

where M_d is the disk mass, d is the distance, κ_v is the dust opacity, ζ is the dust-to-gas ratio, and B_ν(T) is the Planck function given a dust temperature T. We adopt a κ_v of 2.62 cm² g⁻¹ following the power-law scaling of Beckwith et al. (1990) with a power-law index β = 1. We assume a dust-to-gas ratio of 1:100 and a uniform dust temperature of 20 K (Andrews & Williams 2005). We acknowledge that the dust opacity, dust-to-gas ratio, dust size distribution, and dust temperature may vary across our source sample, and further insight into the physical structures of these sources is needed for a more robust abundance comparison. The effects of these uncertainties are discussed further in Section 4.5.

We convert the disk mass M_d to the total number of H₂ molecules ${{ \mathcal N }}_{{\rm{H}}2}$ assuming a typical gas molecular weight of 2.37× the mass of atomic hydrogen. Finally, we convert to an H₂ column density by dividing ${{ \mathcal N }}_{{\rm{H}}2}$ by the area of the mask used to extract the continuum flux. For the Class 0/I sources this is the same DCN emission mask used for the spectral line extraction. Note that for Ser-emb 7, this does not coincide with the continuum peak position. For the Class II sources, we integrate the continuum flux within an ellipse fit to the 260 GHz continuum at a 3–5σ threshold, depending on the continuum brightness. While the Class II continuum emission is generally more compact than the molecular line emission (typically by a factor of ∼2–5 in angular area; see Tables B1 and B2), we assume that the H₂ column densities derived in this way are representative of the disk average. Continuum fluxes, source areas Ω, and H₂ columns are listed in Table B2 in Appendix B.

Figure 4 shows comparisons of the C₂H, HCN, and C¹⁸O column densities (top) and column density ratios with respect to H (bottom), which we will refer to as abundances. For each panel, we also include the Spearman ρ correlation coefficient calculated for the Class 0/I sources and the Class II sources. This measures the degree to which the column densities/abundances of each molecule pair are monotonically related. Possible ρ values range from −1 (strongly anticorrelated) to 1 (strongly correlated), with 0 indicating no correlation. We use the scipy.stats.spearmanr module to calculate Spearman ρ coefficients and their p values. We consider a p value <5% to be statistically significant. Note that this treatment excludes upper limits.

Zoom In Zoom Out Reset image size

The observed correlations are similar when considering column densities or abundances, and our discussion will mainly focus on the abundances. Statistically significant positive correlations are seen between the HCN and C₂H abundances in the Class II sources, and between the HCN and C¹⁸O abundances in the Class 0/I sources. In all other cases, abundance correlations are weakly positive and statistically insignificant.

There are no clear abundance jumps/discontinuities between the Class 0/I and Class II sources. Within the protostellar stage, abundances are typically higher in the Class 0 sources than the Class I sources. Within the Class II sources, there are no clear trends between <4 and >4 Myr sources. Chemical implications of these trends are discussed in Section 6.2.

It is important to note that our sample contains just two Class 0 and three Class I sources, and a larger sample size is needed for a more robust assessment of correlation and scatter within these evolutionary stages. Also, better constraints for sources with non-detections are needed to improve this analysis. In particular, the C₂H transition observed toward the Class 0/I sources is over an order of magnitude weaker than the transition observed toward the Class II sources, and the upper limits for non-detections are not very constraining.

With these observations, we can also evaluate how the C¹⁸O abundance changes with evolutionary stage. Here, for the Class II disks, we re-derived column densities and abundances using C¹⁸O and continuum fluxes extracted from a circle defined by the C¹⁸O beam major axis. By focusing on the inner-most disk regions accessible at this resolution, we aim to reduce the impact of CO freeze-out on our observations. We also exclude J1609, J1612, and J1614 from this analysis, as their C¹⁸O abundances are estimated from ¹³CO, which may be optically thick. Class II source ages and uncertainties are taken from Table 2. For the Class 0 and Class I sources, we adopt ages of 0.1 ± 0.1 and 0.35 ± 0.15 Myr, respectively. This is based on the estimated Class 0 and Class I lifetimes of 0.2 and 0.5 Myr, respectively (Enoch et al. 2009; Evans et al. 2009).

Figure 5 shows the resulting C¹⁸O abundances as a function of the estimated source ages. We emphasize that the C¹⁸O optical depth is a major uncertainty in assessing both the absolute C¹⁸O/H abundances and the relative trends between evolutionary stages. Still, we can use constraints from the literature to estimate how much our observations are affected. C¹⁸O emission has been found to be optically thick in several embedded disks (van 't Hoff et al. 2018; Zhang et al. 2020). The C¹⁸O 2–1 optical depths in the embedded disks in Zhang et al. (2020) are likely around ∼4–8, based on the reported C¹⁸O and ¹³C¹⁸O fluxes and assuming a ¹²C/¹³C ratio around 70. We therefore expect that moderate optical depth effects may affect the protostellar sources in our sample. C¹⁸O has also been found to be optically thick within the CO snow lines of Class II disks, with opacities ranging from ∼8 in TW Hya (Zhang et al. 2017) to several tens in HD 163296 (Booth et al. 2019). It is therefore probable that C¹⁸O optical depth effects range from moderate to severe for the Class II disks in our sample.

Zoom In Zoom Out Reset image size

In light of this, it is likely not appropriate to compare our measured C¹⁸O/H abundances to the ISM abundance of ∼2.5 × 10⁻⁷ (assuming CO/H = 1.4 × 10⁻⁴ and ¹⁶O/¹⁸O = 560). We therefore include in Figure 5 the apparent ISM abundance given an optical depth of 8 and 40 as more realistic benchmarks. Based on the previously mentioned literature measurements, we estimate that the τ = 8 case is the most appropriate comparison for the Class 0/I sources. The Class II disks likely span a range of optical depths between the τ = 8 and τ = 40 cases. As a sanity check, we compare the depletion factors inferred from our C¹⁸O/H abundances to the depletion factors found from detailed modeling in Zhang et al. (2019). We recover similar CO depletion factors for HD 163296 and DM Tau in the τ = 40 scenario (i.e., ∼1–10× depleted), and for IM Lup in the τ = 8 scenario (∼10–100× depleted). This further supports that this is a reasonable range of optical depths to assume for the Class II disks in our sample. Still, we stress the need for observations of a rarer CO isotopologue to better constrain the magnitude of optical depth effects in our observations.

As seen in Figure 5, the Class 0 C¹⁸O/H abundances in our sample are consistent with ISM levels if C¹⁸O is moderately optically thick (τ = 8). On the other hand, even assuming τ = 8, the Class I C¹⁸O/H abundances are depleted in CO relative to ISM levels by about an order of magnitude. The Class II sources are likely depleted in CO by factors of a few up to 100, depending on the optical depths in individual sources. We cannot distinguish if there are any trends in CO depletion with age within the Class II stage given the optical depth uncertainties. However, even in the τ = 40 case we see that most disks exhibit around an order of magnitude in CO depletion relative to ISM levels. The exceptions, V4046 Sgr and J1604, show abundances that are comparable to ISM levels if the C¹⁸O optical depths are large.

It is important to consider possible effects of CO freeze-out on our derived abundances. Although we extract fluxes from the inner beam of the Class II disks to minimize freeze-out effects, in some sources the midplane CO snow line radius may be smaller than the beam. Still, while freeze-out can lower the gas-phase CO abundance in the dense midplane region, detailed models show that freeze-out around the midplane cannot explain the low CO abundances in the emissive surface layers, which our observations are sensitive to (Favre et al. 2013; Zhang et al. 2019). Moreover, for the disks around the Herbig Ae stars HD 163296 and MWC 480, which have a warmer temperature structure and midplane snow lines likely beyond the radius of our extraction mask, we measure C¹⁸O/H abundances comparable to those of other Class II sources. Thus, CO freeze-out cannot fully explain the low CO abundances observed in our sample. For the embedded disks, we expect that the temperature structures are warmer (van 't Hoff et al. 2018), and CO freeze-out should not impact the inner protostellar core. Still, it is important to recognize that episodic accretion events may play an important role in setting the CO freeze-out behavior in protostars (Jørgensen et al. 2015). Improved constraints on the protostellar physical structures and evolutionary histories are needed to better understand whether accretion variations are likely to impact our conclusions.

There are several additional uncertainties that may affect the abundances we derive for all molecules. Optically thick continuum emission could be problematic for recovering accurate molecular line fluxes (Cleeves et al. 2016; Harsono et al. 2018). Additionally, our assumption of a constant dust opacity for all sources may introduce bias into this comparison, as κ depends on the grain size distribution and could therefore vary systematically with source age. However, it is not clear that there is a more appropriate assumption, as dust opacities are not well constrained for any source age. Extrapolating κ_{1.3 mm} = 1.0 cm²g⁻¹ derived for dense protostellar cores in Ossenkopf & Henning (1994) (assuming β = 1) yields a κ_{1.1 mm} about 2.5× smaller than the value that we used. If the true opacity is smaller than our adopted value, this would result in an underestimated H₂ column and overestimated molecular abundances for the Class 0/I sources. Similarly, if the dust-to-gas ratio increases with evolutionary stage, this would result in an overestimated H₂ column and in turn underestimated molecular abundances in the older sources.

Another uncertainty in deriving abundances is the temperature of the gas and dust emitting in the protostellar sources. For instance, gas in the inner protostellar regions may be warmer than 30 K. If we instead adopt an excitation temperature of 100 K in the Class 0/I sources, this increases the column densities of all molecules by a factor of ∼2. Additionally, because the disk cools and becomes more optically thick from the Class 0 to Class I stage, the adopted 20 K characteristic dust temperature may result in mass overestimates for the younger protostellar sources or underestimates for the older protostellar sources (Jørgensen et al. 2009; Dunham et al. 2014b). Adopted dust temperatures of 15 or 30 K instead of 20 K results in ∼40%–50% variations in the derived dust masses (and therefore H₂ columns).

As discussed in Section 1, isotopic fractionation occurs under specific physical conditions, and fractionation patterns are often used to infer the conditions in which molecules formed. We aim to explore deuterium and ¹⁵N fractionation using our observations of the HCN isotopologues. Because the major isotopologue HCN is not covered for the Class 0/I sources, we convert from H¹³CN to HCN fluxes assuming the local interstellar ¹²C/¹³C ratio of 68 ± 15 (Milam et al. 2005). As mentioned in Section 4.3, ¹³C fractionation is not expected to be important in protostellar or disk environments (e.g., Roberts et al. 2002; Hily-Blant et al. 2019). HCN/DCN and HCN/HC¹⁵N ratios are then found from the ratio of column densities determined in Section 3.1.7.

The HCN isotopologue ratios are listed in Table 6 for an assumed rotational temperature of 30 K. HCN/HC¹⁵N ratios range from 88 in Ser-emb 17–287 in Ser-emb 1. For comparison, the ¹⁴N/¹⁵N ratio in the local ISM is 450 ± 22, similar to the value of 441 ± 5 measured in the Solar wind (Wilson & Rood 1994; Marty et al. 2011). Thus, we see some degree of ¹⁵N enhancement in HCN for all of the Serpens sources. HCN/DCN ratios ratios are comparatively more consistent across the Serpens sources, ranging only from 52 in Ser-emb 15–87 in Ser-emb 8. Relative to the local ISM H/D ratio of 4.3 ± 0.4 × 10⁴ and the protosolar value of 5.0 ± 0.9 × 10⁴ (Geiss & Gloeckler 2003; Linsky et al. 2006), these sources show significant deuterium fractionation.

As noted previously, the choice of rotational temperature has a small impact on our results. This is particularly true for column density ratios: changing the rotational temperature from 30–50 K or 30–20 K results in HCN/DCN variations around 6% and HCN/HC¹⁵N variations around 1%. This is small compared to other sources of uncertainty (i.e., the ¹²C/¹³C ratio, Gaussian line fitting, and telescope calibration). In addition to these known sources of uncertainty, we emphasize that this treatment assumes that the HCN isotopologues emit co-spatially. Higher spatial resolution and multiline observations are needed to confirm if all isotopologues are indeed characterized by the same emitting area and rotational temperature.

We now estimate the HCN/DCN and HCN/HC¹⁵N ratios in the Class II disks with detections. In some disks H¹³CN was observed instead of HCN, so like for the Class 0/I disks we convert from H¹³CN to HCN fluxes assuming a ¹²C/¹³C ratio of 68 ± 15 (Milam et al. 2005). When performing this analysis, it is important that the fluxes are extracted from the same spatial region. The CASA task imsmooth is used to smooth the HCN and H¹³CN images to the same resolution as the DCN images. We then use the same mask to extract fluxes of each major and minor isotopologue pair (i.e., HCN or H¹³CN, and HC¹⁵N or DCN). We define the extraction mask as pixels with >3σ emission of both the major and minor isotopologue in order to avoid diluting the intensity of the minor isotopologue. Note that the fluxes used in this analysis differ from the disk-integrated values presented previously. We solve for the major:minor isotopologue column density ratio according to:

$$\begin{eqnarray}&&\displaystyle \frac{{N}_{T}(a)}{{N}_{T}(b)}=\displaystyle \frac{\int {S}_{\nu }d\nu (a)}{\int {S}_{\nu }d\nu (b)}\displaystyle \frac{Q({T}_{r})(a)}{Q({T}_{r})(b)}\displaystyle \frac{{e}^{{E}_{u}(a)/{T}_{r}}}{{e}^{{E}_{u}(b)/{T}_{r}}}\displaystyle \frac{{A}_{{ul}}(b){g}_{u}(b)}{{A}_{{ul}}(a){g}_{u}(a)},\end{eqnarray} \tag{ 4 }$$

where all parameters are the same as described for Equation (1). As before, we adopt a rotational temperature of 30 K. Variations in the isotopologue ratios derived using 20 or 50 K rotational temperatures are incorporated into the uncertainties. As described in Section 4.2, HCN is likely optically thick in the Class II disks. We therefore apply an opacity correction assuming τ = 5.8 ± 1.8 (see Section 4.2) to the HCN fluxes used in this analysis.

Table 6 lists the derived HCN isotopologue column density ratios. The HCN/HC¹⁵N ratio could be measured for two sources in the sample. In LkCa 15 and MWC 480, the HCN/HC¹⁵N ratios of 100 and 106, respectively, are low compared to the local ISM and protosolar values around 450 (Wilson & Rood 1994; Marty et al. 2011), and within the range measured for the Serpens sources (Section 5.1). The ratios we derive are similar to those found in Guzmán et al. (2017). In DM Tau, H¹³CN is detected but HC¹⁵N is not detected, so we place a lower limit on HCN/HC¹⁵N of 115. Neither HC¹⁵N nor H¹³CN was detected toward CI Tau or DO Tau, and so we cannot place any useful limits on the HCN/HC¹⁵N ratio in those sources.

The HCN/DCN ratio could be measured for three sources: HD 143006, J1604, and J1609. For these sources, the HCN/DCN ratio is ∼10 assuming HCN is optically thin, or 54–75 for an HCN opacity of 5.8. We can also place limits on the HCN/DCN ratio for CI Tau of <33, for DM Tau of >70, and for J1614 of >17 (assuming τ_HCN = 5.8). Note that CI Tau is detected in DCN but not H¹³CN (Section 3.2), resulting in an HCN/DCN upper limit. In DO Tau and J1612, neither DCN nor the major isotopologue (HCN or H¹³CN) is detected, so no limits can be placed.

For comparison, the Class II HCN/DCN ratios measured from (optically thin) H¹³CN/DCN flux ratios in Huang et al. (2017) range from ∼14–194, assuming a 30 K excitation temperature. Both the opacity-uncorrected and corrected values that we find are consistent with this range. Regardless of the optical depth, the ratios we measure are very low compared to the local ISM and presolar values around 5 × 10⁴ (Geiss & Gloeckler 2003; Linsky et al. 2006). Compared to the Class 0/I HCN/DCN ratios around 70, the ratios in the Class II sources are either low in the optically thin case or comparable in the optically thick case.

Our observations cover multiple HCN isotopologues in both the Class 0/I and Class II sources, enabling a comparison of HCN/HC¹⁵N and HCN/DCN ratios throughout the disk lifetime. Figure 6 shows the resulting column density ratios derived in Sections 5.1 and 5.2. For comparison, we also show literature measurements of other evolutionary stages: HCN/HC¹⁵N in the pre- and protostellar sources TMC-1 (Ikeda et al. 2002), the dark clouds L183 and L1544 (Hily-Blant et al. 2013), and the protostellar envelopes of IRAS 16293A, R CrA IRS7B, and OMC-3 MMS6 (Wampfler et al. 2014); the Class II disks AS 209, LkCa 15, MWC 480, HD 163296, V4046 Sgr, (Guzmán et al. 2017), and TW Hya (Hily-Blant et al. 2019); and the comets Hale–Bopp and Holmes (Bockelée-Morvan et al. 2008). For HCN/DCN we include measurements in the pre- and protostellar sources TMC-1 (Turner 2001), a sample of infrared dark clouds (Feng et al. 2019), and a sample of protostellar envelopes (Jørgensen et al. 2004); the Class II disks AS 209, LkCa 15, MWC 480, HD 163296, V4046 Sgr, (Huang et al. 2017), and TW Hya (Qi et al. 2008); and the comets Hale–Bopp (Meier et al. 1998) and Hyakutake (Bockelée-Morvan et al. 2008). For the Class II measurements taken from Guzmán et al. (2017) and Huang et al. (2017), we re-derive column density ratios assuming a 30 K rotational temperature to enable a consistent comparison with our results.

Zoom In Zoom Out Reset image size

The HCN/HC¹⁵N ratios in the Class 0/I sources range from ∼90–290, and decrease with evolutionary stage as traced by bolometric temperature. Most Class II disks have an HCN/HC¹⁵N ratio between ∼80 and 150. The exception, TW Hya, has a disk-averaged ratio above 200, though the ratio in the inner disk is lower (121 ± 11; Hily-Blant et al. 2019). The HCN/HC¹⁵N ratios in the Class II disks are comparable to the Class 0/I sources Ser-emb 8 and 17 (and the lower limit in Ser-emb 15), and low compared to the other pre- and protostellar sources. Cometary HCN/HC¹⁵N ratios are consistent with (if on the high end of) the values measured in Class II disks.

Compared to the HCN/HC¹⁵N ratio, the HCN/DCN ratio is quite flat across our sample of Class 0/I disk candidates, ranging from just 52–87. These ratios are comparable to measurements in pre- and protostellar sources. In the Class II sources, the new measurements we derive have large uncertainties due to potential optical depth effects (Section 5.2). Even so, it seems that the HCN/DCN ratio is slightly lower in young Class II sources compared to the protostellar disk sources. The HCN/DCN ratio again increases slightly in >3 Myr Class II sources, beginning with DM Tau. Cometary HCN/DCN ratios are high compared to most Class II disks. The exception is V4046 Sgr, whose anomalously high HCN/DCN ratio is consistent with the upper limit measured toward Hale–Bopp.

In Section 4.4 we presented evolutionary patterns in C¹⁸O/H abundances. While there are a number of important uncertainties in our derived abundances, namely, the C¹⁸O optical depths, we can still identify several trends from the existing data. We emphasize that more rigorous abundance retrievals and/or observations of rarer CO isotopologues are needed to confirm our findings.

The Class 0 C¹⁸O/H abundances are not distinguishable from ISM levels with our data, assuming that C¹⁸O optical depths up to a factor of ∼10 are possible in these sources. It is also important to note that the disk cannot always be isolated in the Class 0 sources: the morphology of C¹⁸O emission in the Class 0 sources (Figure 1) is more suggestive of envelope emission than protostellar core/disk emission. Observations of a rarer isotopologue would help to isolate emission from the inner disk-forming region (e.g., Zhang et al. 2020).

The C¹⁸O/H abundances in the Class I sources are difficult to reconcile with an ISM abundance, even assuming moderate C¹⁸O optical depths. It is likely that these sources are depleted in CO by about an order of magnitude compared to ISM levels. We note that Zhang et al. (2020) recently found evidence for ISM-like CO abundances in three protostellar disks, which may signify that a range of depletion timescales and outcomes are possible. Indeed, chemical models by Drozdovskaya et al. (2014) demonstrate that differences in infall physics can strongly influence the chemical processing that occurs in embedded disks. The sources in Zhang et al. (2020) are among the most massive protostellar disks known, and their higher CO abundances may reflect a different physical/chemical environment compared to the Serpens sources. Still, the results from this work and Zhang et al. (2020) both indicate that CO depletion is a rapid process taking place on a ∼0.5–1 Myr timescale.

The decrease in C¹⁸O/H abundances from the Class 0 to Class I sources is either an artifact of envelope material clouding our view of the Class 0 core/disk regions, or reflects an active depletion mechanism. The drop in CO abundances also coincides with an increase in ¹⁵N fractionation of HCN, which is likely due to enhanced photochemistry during envelope clearing (see Section 6.3). It is possible that envelope clearing also facilitates CO depletion, e.g., through increased exposure to X-rays which drive CO destruction by reactions with He⁺ (e.g., Favre et al. 2013). We also note that single-dish studies show that CO abundances on envelope scales increase from Class 0 to Class I objects (Jørgensen et al. 2002, 2005), converse to what we see here. This implies that the chemistry in the inner protostellar regions is distinct from the envelope, and highlights the importance of understanding chemical evolution during infall.

The Class II sources have a wide range of probable C¹⁸O optical depths, from factors of a few to several tens (Zhang et al. 2017; Booth et al. 2019), making it difficult to interpret depletion trends with age within the Class II stage. Still, even with a conservative estimate of C¹⁸O optical depth, τ = 40, the majority of disks exhibit C¹⁸O/H abundances about an order of magnitude below ISM levels. The similar depletion levels in the Class II sources compared to the Class I sources suggests that the low CO abundances seen widely in Class II disks (e.g., Dutrey et al. 2003; Favre et al. 2013; Cleeves et al. 2016; Zhang et al. 2017) are set very early in the disk lifetime. Indeed, Cleeves et al. (2016) suggest depletion in an earlier evolutionary stage as a possible explanation for CO depletion in the young disk IM Lup (∼0.5 Myr; Andrews et al. 2018).

While most disks have apparent C¹⁸O/H abundances around 10⁻¹⁰–10⁻⁹, two sources, V4046 Sgr and J1604, have notably higher abundances closer to 10⁻⁸. Interestingly, both are old (∼13 Myr) transition disks. The high C¹⁸O abundances in these disks could be explained by a distinct chemistry taking place at the heavily irradiated inner disk edge, or could be an artifact of different dust properties skewing the derived gas masses in (old) transition disks compared to full disks.

In summary, our results are consistent with a scenario in which CO is depleted relative to ISM levels on a ∼0.5–1 Myr timescale. We do not see evidence for significantly higher depletion levels in the Class II stage compared to the Class I stage. This suggests that depletion levels are largely set in the embedded stages, and further depletion during the Class II stage is not especially efficient. An important consideration if CO depletion takes place early is that models of physical and chemical mechanisms for CO depletion in disks generally begin with ISM CO abundances and assume physical properties of Class II disks (e.g., Reboussin et al. 2015; Eistrup et al. 2018; Krijt et al. 2018; Schwarz et al. 2018). Our results underscore the need to test whether these (or other) mechanisms could deplete CO in the context of protostellar rather than protoplanetary disk physics and chemistry.

We now assemble evidence from the C₂H, HCN, and C¹⁸O emission morphologies (Figure 1) and abundance correlations (Figure 4) to gain insight into the chemical and physical factors driving their formation. In particular, we wish to explore how the chemistries of the three molecules relate to one another. As described in Section 1, models predict that CO will be removed from the disk atmosphere over its lifetime, as some combination of physical and chemical processes sequester volatiles in the midplane (e.g., Meijerink et al. 2009; Krijt et al. 2016; Schwarz et al. 2018). In an increasingly O-poor (i.e., high C/O) gas in the disk atmosphere, HCN and C₂H formation are predicted to be enhanced (Du et al. 2015). The expectation would therefore be that CO abundances decrease and HCN and C₂H abundances increase as the disk evolves. While our results do not fit this evolutionary scenario, the underlying chemistry (i.e., efficient C₂H and HCN formation in O-depleted gas) can still explain most of our results.

6.2.1. Class 0/I Sources

6.2.2. Class II Sources

As shown in Figure 6, the protostellar sources show a decreasing HCN/HC¹⁵N ratio with evolutionary stage as traced by bolometric temperature. This indicates active HC¹⁵N fractionation in the protostellar stage. This trend is readily explained by photodissociation-driven fractionation turning on as the envelope clears and the disk is increasingly exposed to UV from the protostar. If temperature-driven fractionation were responsible, we would expect to see a similar trend in the evolution of DCN fractionation. We note that a photochemical fractionation scenario requires gas-phase HCN formation via N₂ dissociation, implying that HCN observed in the protostellar cores is at least in part formed in situ in the gas. Photochemical fractionation is also thought to be important in Class II disks based on radially resolved HCN/HC¹⁵N measurements (Guzmán et al. 2017; Hily-Blant et al. 2019). The agreement between HCN/HC¹⁵N ratios in the most evolved protostellar sources and the Class II sources supports that the same mechanism is at play. Still, our inferences are based on disk-averaged properties, and higher-resolution observations that reveal radial variations in the HCN/HC¹⁵N ratio would help to confirm whether and how photochemical fractionation takes place in these sources.

Unlike HC¹⁵N, deuterium fractionation in the protostellar stage is either inefficient or proceeds with comparable efficiency across all Class 0 and Class I sources, resulting in flat HCN/DCN ratios. This likely reflects that deuterium fractionation is typically temperature driven rather than photodissociation driven, and thus should not exhibit strong trends with envelope clearing. With the data available to date, we see a tentative trend in which HCN/DCN ratios decrease slightly from Class I to young Class II disks, followed by a modest increase in HCN/DCN ratios in >3 Myr Class II disks (with the exception of the anomalously high increase in V4046 Sgr). These trends may trace evolution of the physical structures of Class II disks, however a spatially resolved analysis is needed to better understand the dependence of D fractionation patterns on disk physical properties.