Protostellar and Protoplanetary Disk Masses in the Serpens Region

We present the results from an Atacama Large Millimeter/Submillimeter Array (ALMA) 1.3 mm continuum and $^{12}$CO ($J=2-1$) line survey spread over 10 square degrees in the Serpens star-forming region of 320 young stellar objects, 302 of which are likely members of Serpens (16 Class I, 35 Flat spectrum, 235 Class II, and 16 Class III). From the continuum data, we derive disk dust masses and show that they systematically decline from Class I to Flat spectrum to Class II sources. Grouped by stellar evolutionary state, the disk mass distributions are similar to other young ($<3$ Myr) regions, indicating that the large scale environment of a star-forming region does not strongly affect its overall disk dust mass properties. These comparisons between populations reinforce previous conclusions that disks in the Ophiuchus star-forming region have anomalously low masses at all evolutionary stages. Additionally, we find a single deeply embedded protostar that has not been documented elsewhere in the literature and, from the CO line data, 15 protostellar outflows which we catalog here.


INTRODUCTION
Stars are born from cores within dusty molecular clouds. Though it is widely accepted that planets form in circumstellar disks surrounding these young stars, the effect of disk properties on planet formation is still relatively unknown. Especially uncertain, but important to identify, are the effects of the large cloud scale environment such as stellar density and strong radiation fields, in regulating star, disk, and planet formation.
Most stars, and therefore their accompanying protoplanetary disks, form in clusters (Lada & Lada 2003). If larger than ∼ 200 members, these clusters generally contain at least one OB star (McKee & Williams 1997) and their accompanying strong UV radiation field enhances photoevaporation and may decrease the dust masses of disks nearby (Mann et al. 2014;Ansdell et al. 2017), depleting the reservoir of planet-forming material. However, dynamical interactions in a dense cluster can also reduce disk masses and lifetimes (Bate 2018). To distinguish between photoevaporative (Rosotti et al. 2014) and dynamical effects requires the study of a young, dense star-forming region that lacks massive stars. The Serpens star-forming region (Gutermuth et al. 2008) is one such natural laboratory.
The high visual extinction and low Galactic latitude toward the Serpens and Aquila regions (RA∼18h, Dec∼00 • ) led to prior confusion over membership and distance. With the advent of Gaia (Gaia Collaboration et al. 2016, it has become clear that most young stellar objects (YSOs) in Serpens are clustered into four main groups in a common star-forming region spread over roughly 10 square degrees at a distance ranging from ∼ 400 − 500 pc (Herczeg et al. 2019;Dzib et al. 2010;Ortiz-León et al. 2017;Ortiz-León et al. 2018;Zucker et al. 2019). The tight clustering, associated cloud, and large numbers of protostars all indicate a young, relatively unevolved region.
Multiple studies have observed nearby young starforming regions in an effort to characterize the statistical properties of protoplanetary disk populations. In doing so, these surveys provide constraints on the initial conditions of planet formation. Like Serpens, young star-forming regions (< 3 Myr) have been the primary target of such studies (Lupus, Taurus, Orion, Perseus, Ophiuchus; Andrews et al. 2013;Ansdell et al. 2016; arXiv:2204.08731v1 [astro-ph.SR] 19 Apr 2022 Tychoniec et al. 2020;Tobin et al. 2020;Cieza et al. 2019). These works have predominantly found that (sub-)millimeter luminosities decrease with time due to dust processing whereas the gas and micron-dust may dissipate separately on longer timescale (Michel et al. 2021). However, it is necessary to identify if there are secondary dependencies, in particular the presence of massive stars and high stellar density.
A previous survey of the Serpens region was completed by Law et al. (2017). They used the Submillimeter Array (SMA) to observe 62 disks at a resolution of ∼ 3 . Here we use ALMA to observe the same region to increase the sample size of disks in the region by a factor of more than 5 with greater sensitivity. The primary goal of this project is to analyze the distribution of disk masses within Serpens and compare it to other nearby star forming regions. The layout of the paper is as follows: we describe the survey sample in §2, the ALMA observations in §3, and our results in §4. The mass distributions with comparisons to other regions are depicted in §5. We discuss the implications of the work in §6 and summarizing our findings in §7, while the protostellar outflows seen in the CO line are discussed in Appendix A.

SOURCE SELECTION AND CLASSIFICATION
The Spitzer c2d and Gould Belt surveys identified 1546 candidate YSOs toward the Serpens region (Dunham et al. 2015). We reduced this to a sample size of 320 for our ALMA survey by selecting only those sources with infrared spectral slopes < 1.6, as defined by from 2 to 24 µm, corresponding to Class I, Flat Spectrum, and Class II YSOs (Williams & Cieza 2011); a K-band cutoff of 13 mag to remove very low mass stars, 0.2 M ; and a Gaia DR2 -based distance cutoff of 650 pc.
The visual extinction toward Serpens is high. Over 50% of stars toward this region have A V > 7 mag (Harvey et al. 2007) and Spitzer surveys were restricted to regions of high extinction (A V > 6 mag) to limit confusion from the Galactic background (Dunham et al. 2015). Consequently, the K-band magnitude is significantly extincted which decreases the spectral slope used for YSO classification. In general, the difference between the observed and intrinsic slopes between two wavelengths λ 1 and λ 2 is where A λ is the extinction at a given wavelength and we assume that λ 2 is large enough that A λ2 = 0. From the extinction law in Wang & Chen (2019), when measuring the slope between the 2MASS K band and MIPS 24 µm, ∆α IR ≥ 0.19 for our sample which is sufficient to shift many sources from Class III to Class II, Class II to Flat, and so on. We therefore revised the YSO classification to use the slope between IRAC 3.6 µm and MIPS 24 µm bands rather than the 2MASS K and MIPS 24 µm bands. In this case, ∆α IR = 0.11 for A V = 6 mag, which is small enough to substantially reduce the number of potentially misclassified YSOs. Though Dunham et al. (2015) accounted for extinction in constructing SEDs to calculate α IR in their original identification of YSO candidates, we choose to base our classification system on the 3.6 − 24 µm slope for two main reasons. First, the extinction correction is dependent on an object's intrinsic color, as well as the nature of material along the line of sight and many objects in Serpens have inaccurately determined spectral types, particularly the more embedded ones. Additionally, they assumed a uniform extinction law, though dust properties may vary throughout the region, again especially for the more embedded objects. Our method, using a sightly longer anchor wavelength is less sensitive to YSO temperature and dust properties. Nevertheless, the effect is not large and our classification agrees with Dunham et al. (2015) for 85% of the sample. Figure 1 shows the distribution of our objects in colorspace grouped by class to illustrate how evolutionary classification is derived from colors. We then crosscorrelated our sample with the list of transition disk candidates (TDs) identified by mid-infrared dips in their spectral energy distributions by van der Marel et al. (2016). These disks may have large central gaps or cavities and are demarcated with circles in Figure 1. It is possible that our sample is not complete, particularly for embedded YSOs. However, our sample has similar demographics to other regions presented in Dunham et al. (2015) and can therefore be considered representative.
In the time between proposal, observations, and analysis, the Gaia Early Data Release 3 (EDR3 ) catalog was released. This decreased distance uncertainties and increased the number of sources with parallax measurements. Four sub-regions were identified by visually identifying groups of sources based on their angular coordinates and then confirmed by checking that most sources in these clusters had similar parallaxes. Figure 2 shows the location of each observed source on the sky. Herczeg et al. (2019) refer to these clusters as Serpens Northeast, Serpens Main, Serpens South, and Serpens far-south. The four regions are labeled in Figure 2 of roughly 435 ± 25 pc, 495 ± 90 pc, 455 ± 50 pc, and 410 ± 35 pc, respectively. We calculated the weighted average of the parallaxes of objects that had Gaia EDR3 measurements within a subregion and converted this to an average distance. This averaged distance was assigned to objects without parallax measurements. The stellar density is high in these subregions. The densest is Serpens South where the median nearest neighbor distance is 0.1 pc (Gutermuth et al. 2008, accounting for the revised distance) which is comparable to protostellar core size scales. Serpens also lies near the HII region W40, and the two are estimated to reside at similar distances (Ortiz-León et al. 2018). Previous works have confirmed that the two star-forming regions are physically connected and interacting (Shimoikura et al. 2020). Particularly for objects in Serpens South, star and disk formation may be affected by close proximity to this other stellar nursery.

OBSERVATIONS
We observed the sample of 320 targets in ALMA Cycle 7 program 2019.1.00218.S (PI: van der Marel). The observations were carried out in Band 6 in the C-2 configuration (baselines from 15 to 314 m) with 43 antennas on 2019 December 19 th and 31 st . Each source was ob-served in a single pointing with an integration time of 20 s. The precipitable water vapor ranged from 1.0 to 4.1 mm over the course of the observations and the average system temperature was 100 K. Three wide spectral windows provided a total continuum bandwidth of 5.6 GHz centered on 241.4 GHz and the fourth was a narrow window centered around the 12 CO 2-1 line (230.536 GHz) at 0.32 km s −1 spectral resolution.
We applied standard imaging techniques to the pipeline calibrated visibilities using the Common Astronomy Software Applications package (CASA, Mc-Mullin et al. 2007) tclean task. The Briggs robust weighting parameter was set to 0.5, and the beam size ranged from 1. 2 × 1. 0 to 1. 4 × 1. 0 with a position angle from −60 • to 75 • . The average root-mean-square (rms) uncertainty of the continuum images was 0.45 mJy beam −1 but ranged from 0.13 to 3.2 mJy beam −1 where the highest values were due to clustering (multiple 3σ detections in a single pointing, increasing noise in the image) and low image quality associated with the limited uv coverage in the short integrations. We detected 130 sources in the continuum and show a montage of their images in Figure 3. All are point sources at the observed resolution and the 7 × 7 panels show binary companions in several cases. The same standard imaging techniques were used to produce channel maps of 12 CO once the continuum had been subtracted from the data with uvcontsub from CASA. Each object's channel maps were visually inspected to determine if they contained any emission. These channels were used in the creation of integrated intensity maps for each object. If, for an object, all channels were emission-free, all channels (steps of 0.5 km s −1 ) between -10 and 30 km s −1 (relative to the CO rest frame at 230.536 GHz) were used. Then, the CASA task immoments was used to create moment 0 integrated intensity maps. A mask was applied to these images such that only pixels with values greater than 50 mJy beam −1 channel −1 were included in the maps. The average rms was ∼30 mJy beam −1 channel −1 . In general, our sensitivity is insufficient to detect emission from the disk but we were able to identify outflow emission in several cases that we discuss in Section A.

RESULTS
Continuum flux densities were measured using aperture photometry, as in Ansdell et al. (2016). Circular apertures were centered on the brightest pixel associated with our target. For brighter sources with flux densities > 5 mJy, we used an aperture radius of 2 to ensure most of the emission was captured, and for the fainter sources, a radius of 1. 5. We chose these on the basis of a curve-of-growth analysis that showed this was where the highest signal-to-noise ratio was typically achieved and also that less than ∼10% of an object's flux extended past either of these aperture boundaries. The rms in the measurement was estimated to be the standard deviation of measurements from ten apertures of the same size uniformly distributed in position angle and spaced two times the aperture radius from the central source. In the case of millimeter binaries, apertures that contained any companion sources were removed from flux calculations, though fluxes of the companions are reported later on in Table 2.
Disks were marked as detections if they were within 1 of the source's 2MASS coordinates and if the measured flux density was greater than three times the rms. For a handful of cases which were either highly clustered or where the image was affected by a bright off-center source, we measured fluxes manually using the CASA imview task with an aperture of 1 . For single sources, we verified our measurements with those measured as a point source in the visibility plane through uvmodelfit. Two marginal cases were identified as visual detections, and though they were not counted as such with aperture photometry due to the high degree of clustering in the field, they were confirmed with uvmodelfit. Therefore, these cases use the uvmodelfit results in our tables and subsequent analysis.
The dust masses of the disks were calculated using the simplest assumptions: optically thin emission and constant dust temperature, so as to allow direct comparisons with other regions (essentially, millimeter luminosity). We use a dust opacity coefficient, κ ν = (ν/100GHz) cm 2 g −1 (Beckwith et al. 1990) and, for Class II sources, a dust temperature T dust = 20 K. Then, using the standard formula for a disk with flux density F ν , the dust mass is where B ν is the Planck function, κ ν = 2.41 cm 2 g −1 at our observing frequency, and d is the distance to each   Table 1). Each panel is 7 by 7 , and is scaled from -0.3 mJy to the maximum pixel value of the primary. This causes a few bright secondary sources to be saturated in some fields. White crosses are plotted at the location of maximum intensity for binary companions (defined in §4). A representative beam size is shown in the lower left panel, and the infrared spectral classification for all sources is noted in the bottom right corner of each panel.
source. This approach was the same as that used for Lupus and Ophiuchus Class II disks (e.g., Ansdell et al. 2016;Williams et al. 2019).
Protostellar envelopes keep the disk warmer than in Class II sources (van't Hoff et al. 2020). In a survey of embedded Orion objects, Tobin et al. (2020) derived disk masses using T dust = 40 K for solar luminosity Class I sources. For Flat spectrum sources with a thinner envelope, Tychoniec et al. (2020) used T dust = 30 K. We will compare with these works later and use these characteristic dust temperatures for the flux-to-mass conversion. This gives scaling factors of 1.71 (Class I) and 2.40 (Flat) instead of 4.00 in equation 3. In all cases, non-detections were assigned a mass upper limit corresponding to three times the flux density rms. The median mass limit for the non-detections is 7.6 M ⊕ with a fairly large range (factor of 3) due to the variation in intrinsic noise and, especially, distance across our sample. Unsurprisingly given the sensitivity of our survey, short integration times, and the fact that Class III sources are generally very faint (Lovell et al. 2020;Michel et al. 2021), all 16 Class III sources were undetected. Table 1 lists the basic properties of our target list, including the derived dust masses or upper limits. Flags indicate if the sources have potentially interesting features. O is used for objects that have an associated 12 CO molecular outflow. B is used for objects in binary systems -a 3σ mm detection within 2. 5 (∼1000 AU) of the central source. G is used to indicate sources located at distances outside of the probable range of Serpens members. In the change from Gaia DR2 to EDR3, 13 sources were found to lie at very large (> 675 pc) or small (< 275 pc) distances -well outside the boundaries of the cloud. We exclude these in the analysis of the Serpens mass distributions here but include them in the table of results for completeness. T is used for transition disk candidates as identified by van der Marel et al. (2016). Once these sources had been excluded, our sample comprised 16 Class I, 35 Flat spectrum, 235 Class II, and 16 Class III sources.

Multiplicity
Of the 320 observed sources, four were noted to have binary companions (a 3σ sub-mm detection within 2. 5 of the target source). Because of their small angular separation on the sky and strong sub-mm emission, we assume that they are binary companions at the same distance as our target sources. These objects don't have counterparts in Gaia, but are visible in optical or infrared regimes in published data. Table 2 shows the basic properties of these objects. These objects had not previously been reported as binaries in the literature, and are noted here.
In addition to the four binary sources that were identified, one candidate Class 0 source was also found in the sample. Located in the same pointing as J18275472-0342384, this source is easily detected in the continuum and shows significant 12 CO emission. However, a submm source with these coordinates is not present anywhere in the literature. In optical wavelengths, the Class 0 source seems to be obscured by nebulosity in images of J18275472-0342384. However, Figure 4 reveals a filamentary bridge-like structure of 12 CO connecting the the two objects, suggesting that these objects formed together and that the Class 0 object is just now emerging from its obscuring envelope.

Intra-region Comparison
As disks age and YSOs evolve, disk dust masses are expected to decrease due to a combination of accretion onto the star, photoevaporation, and planet formation. There are 16 Class I, 35 Flat spectrum, and 235 Class II sources that we can use to study evolutionary effects on disk masses in our dataset. We used the Kaplan-Meier estimator within the python lifelines (Davidson-Pilon 2019) package to incorporate the information from the upper limits when creating cumulative mass distributions for each evolutionary group. The distributions for Class I, Flat, and II objects are shown in Figure 5 and demonstrate a clear progression from high to low mass as the inner disk evolves (as characterized by the infrared slope). The uncertainties are large for Class I and Flat object distributions due to the relatively small sample sizes but the detection rates are reasonably high (87% and 63%, respectively). The Class II distribution is better characterized due to the large sample but the detection rate is only 37%. The flattening of the Class II curve below ∼ 10 M ⊕ is due to a small number of detected sources with masses lower than most non-detection's mass upper limits. This is simply a signature of the relatively wide range in mass sensitivity due to the range of noise in the observations and distances to the targets. However, the low detection rate prevents us from a detailed characterization of the sample by, e.g., fitting a gaussian probability distribution function as in Williams et al. (2019).
Above 10 M ⊕ , well beyond the median upper mass limit, the distributions are reasonably complete and can be compared. For this subset, the mean masses are 57 M ⊕ for Class I, 57 M ⊕ for Flat, and 48 M ⊕ for Class II.
The relationship between disk dust mass and stellar mass was also investigated, as it has been previously    noted that the two positively correlate (Andrews et al. 2013). However, only a small fraction of our sources (∼30) had reliable stellar mass measurements. We tried fitting the spectral energy distributions of our sources to various grid models to estimate stellar mass, but the number of degeneracies in our input parameters (spectral type, stellar luminosity, etc.) prevented us from doing so. Nevertheless, we analyzed the M dust − M relation for the few sources in our sample with known M values and plotted our results in Figure 9. Our results showed a slight positive correlation between M dust and  Figure 5).

Class II Regional Comparison
Class II YSOs have much longer lifetimes, on average, than Class I and Flat spectrum protostars (Dunham et al. 2015) so they provide both greater numbers and a longer temporal baseline to study disk evolution. Many millimeter wavelength surveys of Class II "protoplanetary" disks have been carried out both pre-and post-ALMA and generally demonstrate a monotonic decrease in dust mass from ∼ 1 to ∼ 10 Myr (e.g., Ansdell et al. 2017;Villenave et al. 2021). A notable exception is Ophiuchus, which stands out for its systematically low disk masses despite its young age (Williams et al. 2019). A possibly similar region is Corona Australis (Cazzoletti et al. 2019), though new Gaia measurements have expanded the membership considerably and suggest a more evolved state (Galli et al. 2020). However, Serpens has a similar proportion of Class I, Flat and II sources as Ophiuchus and therefore provides a useful test of how unusual Ophiuchus might be.
Rather than compare the Serpens Class II mass distribution with the many others measured in other starforming regions, we restrict our attention here to young regions and Ophiuchus in particular. These regions have an average optical YSO age estimated from pre-main sequence tracks to be 3 Myr and evolutionary states characterized from Spitzer surveys. Figure 6 plots the Kaplan Meier estimator for Class II sources with masses derived from the same Beckwith κ ν prescription and uniform T dust = 20 K in Taurus (Andrews et al. 2013), Lupus (Ansdell et al. 2016), and Ophiuchus (Williams et al. 2019). Although the sensitivity of the Serpens sample is much lower than the other regions, there is a remarkably tight overlap with Lupus and Taurus where the distribution is near complete at masses 10 M ⊕ . This agrees with the results of Law et al. (2017), which saw the same behavior within the same mass range when observing with the SMA, but for a much smaller sample. Given the relatively high proportion of embedded protostars, the Serpens region is likely in a slightly younger evolutionary state than Lupus and Taurus as a whole and more similar to Ophiuchus. Nevertheless, Ophiuchus' Class II mass distribution is the anomaly within these comparisons.

Class I Regional Comparison
Similarly to Section 5.2, a comparison of the disks around Class I and Flat spectrum protostars in different star forming regions provides insight into the earlier stages of disk evolution. As these evolutionary stages are short (0.44 and 0.35 Myr (Evans et al. 2009), respectively), their numbers are much smaller than those of Class II objects, and Lupus and Taurus do not have sufficient statistics (nor uniform, unbiased surveys) to perform such a comparison. Consequently, we compare Serpens and Ophiuchus with two other regions that are both young and large enough to have many Class I and Flat spectrum sources: Orion and Perseus, surveyed with ALMA by Tobin et al. (2020) and Tychoniec et al. (2020), respectively.
To make the comparison most homogeneous, we recalculated disk masses from the published flux density values using the same κ ν prescription described in S4 and with uniform dust temperatures T dust = 30 K for Flat spectrum and T dust = 40 K for Class I sources. To improve the statistics, especially for Serpens, we then bundled the two categories together in the cumulative mass distribution plots in Figure 7.
The Serpens distribution is very similar to that of Orion and Perseus but, as for the Class II sources, Ophiuchus is atypical. As the detection rate is higher for these brighter sources in our sample, the cumulative distribution is better characterized, extending to cumulative fractions greater than 50%. The median disk dust mass values (including upper limits) are roughly 16, 18, 27, and 2 M ⊕ for Orion, Perseus, Serpens, and Ophiuchus, respectively. We note that these values differ from previously published results due to our recalculation of the masses.
With a detection rate of 70% for the Class I and Flat Spectrum sources, the Serpens distribution is sufficiently complete for the cumulative mass distribution to be fit with an error function and estimate the parameters of a gaussian probability distribution. Indeed, all four distributions in Figure 7 are well fit by this method and the results are shown in Table 3. The means of the distributions are indeed similar for Orion, Perseus and Serpens, whereas Ophiuchus is a factor of 10 lower. Despite this large difference, however, the dispersions are remarkably similar: ∼ 0.7 − 1 in log 10 M dust . In addition to fitting the observed cumulative dust mass distributions to probability distributions, we tested if the differences in distributions was statistically significant. With the lifelines package, we used the log rank test to compare these populations in a non-parametric way that accounted for non-detections. p values, or the probability that the distributions were different simply because of random chance, were larger than 0.05 when comparing Orion, Perseus, and Serpens with each other, so we cannot claim that these three regions' disks come from statistically different distributions. However, when any of these three regions was compared with Ophiuchus, p values were consistently <0.005. Therefore, we can state that the youngest disks in Ophiuchus come from a different parent distribution than those in Orion, Perseus, or Serpens.

DISCUSSION
With one notable exception, the disk mass distributions of the various young star-forming regions analyzed here share striking similarities and largely overlap one another. We have found that the Class II disks in Serpens have essentially the same distribution as in Taurus and Lupus, at least above 10 M where our survey is near complete. Moreover, the disks around the embedded Class I and Flat spectrum sources have a significantly similar distribution to those in Orion and Perseus. This is despite the differences in environments: the YSO density in Taurus and Lupus is low (∼1-10 pc −3 , Luhman 2018) whereas the Serpens sources are strongly clustered, ∼430 pc −2 Gutermuth et al. 2008).
It is important to note that our observations were made at a resolution of roughly 1 , which could blend possible binaries and lead to a disk dust mass distribution that is biased towards massive objects. However, Serpens is a low-mass star forming region, and as stellar multiplicity correlates with stellar mass (Duchêne & Kraus 2013), we expect this is a minor effect. Enoch et al. (2011) estimated a protostellar multiplicity fraction of ∼ 10% within the cloud which is too small to significantly affect our conclusions. During the review phase of this manuscript, some higher resolution, ∼ 0. 3 data were taken and are currently being analyzed (Tong et al. in prep). These should resolve tighter binaries concealed within our sample.
Orion is a massive star-forming region with several O stars that produce a strong UV radiation field in their immediate vicinity. This is known to enhance photoevaporation and produce lower disk masses at close range (Mann et al. 2014;Ansdell et al. 2017), but most disks are many parcsecs away from these stars and, on the scale of the entire star forming region, there is no discernable effect on the disk population when compared to Serpens. Even though the YSOs in Serpens South are spaced by a typical protostellar core diameter (0.1 pc, Gutermuth et al. 2008), the much smaller disks that are rapidly evolving in their centers are not noticeably affected by the crowding. The notable exception is Ophiuchus, which stands out as anomalous in both the protoplanetary and protostellar disk mass distribution comparisons. It is clearly a young star-forming region as the surrounding molecular cloud has not substantially dispersed, there is a relatively high fraction of Class I and Flat sources (Evans et al. 2009;Dunham et al. 2014), and the Class II YSOs have high luminosities and ages estimated from pre-main sequence tracks of ∼ 0.1 − 1 Myr (Luhman & Rieke 1999). Although there are fewer Class 0 objects in Ophiuchus compared to the other young regions here, such sources are very short-lived and intrinsically rare. The similar proportion of (longer-lived and therefore more common) Class I and Flat objects when compared with Class II objects (Dunham et al. 2015) is a more robust indication that these regions are in similar evolutionary stages and can be equitably compared.
The low disk masses in Ophiuchus were already noted by Williams et al. (2019) for Class II sources and by Tobin et al. (2020) and Tychoniec et al. (2020) for Class I and Flat spectrum sources, respectively. The comparison here with Serpens, which is also a low mass, clustered star-forming region with similar multiple indicators of youth, allows a direct comparison across both protoplanetary and embedded disks and confirms that Ophiuchus is an outlier. The median mass difference is very large for each disk subset -about an order of magnitude. Despite McClure et al. (2010) showing a lower contribution from the envelopes around many Class I and Flat objects, the stark contrast between median masses in Ophiuchus and other regions cannot be ex-plained away by possible source mis-classification. The low masses at all evolutionary stages suggest that the differences are intrinsic and likely began at formation (e.g., see Forbes et al. 2021, for evidence for a global event). Potential avenues to better understand the cause are observational studies of Class 0 sources and cloud-todisk scale numerical simulations of star formation (Bate 2018).

SUMMARY
We surveyed 320 protoplanetary disks (16 Class I, 35 Flat spectrum, 235 Class II, and 16 Class III) in the Serpens star-forming region with ALMA at 1.3 mm (Band 6) and analyzed the continuum and 12 CO line data. Our main results are as follows: • We detected 130 sources in the continuum. The derived disk dust masses systematically decline from Class I to Flat spectrum to Class II sources.
• We compared the Class II mass distribution with those of other young star-forming regions with large numbers of YSOs that have been surveyed at high sensitivity. Serpens is very similar to Lupus and Taurus for masses greater than the median detection limit ∼ 10 M ⊕ . This indicates that the YSO stellar density does not greatly affect the disk demographics. However, the Ophiuchus distribution is shifted to much lower masses.
• We compared the mass distribution of Class I and Flat spectrum sources with other young regions where there are near-complete ALMA surveys of their embedded disk population. The Serpens distribution overlaps with those of Perseus and Orion indicating that massive stars do not have a strong effect on the overall disk population. Similar to other studies, we find that Ophiuchus is the outlier with anomalously low masses, and confirm this finding with the log rank test.
• Serpens has many similar properties to Ophiuchus as each are very young, low mass, clustered starforming regions with similar proportions of Class I, Flat spectrum and Class II sources. The direct comparisons shown here reinforce previous conclusions that Ophiuchus disks have anomalously low masses at all evolutionary stages. The mean difference is about an order of magnitude for both Class II and Class I / Flat which suggests that this is a characteristic property of the region that may have been inherited from its initial conditions.
• The survey was too shallow and the resolution too low to unambiguously detect CO line emission from the disks, but we were able to identify several protostellar outflows and catalog them in the Appendix. This provides a relatively unbiased sample that would be interesting to follow up with deeper observations and analysis.

APPENDIX
A. PROTOSTELLAR OUTFLOWS Young stars are the drivers of protostellar outflows, or bipolar jets ejecting mass from a YSO (Arce et al. 2007). These outflows sweep up material nearby, and may therefore influence the formation and evolution of planet-forming disks. Additionally, as outflows are more energetic than the ambient medium and therefore move at distinct velocities, they can be identified separately from surrounding cloud material by manually inspecting channel maps. Fifteen molecular outflows originating from the sources in our sample were identified in 12 CO, and shown in Figure 8. Four of the outflows were associated with Class I sources, seven with Flats, and four with Class IIs. Therefore, roughly 20% of Class I or Flat sources seem to be producing protostellar outflows, while only ∼ 2% of Class II sources do, showing that outflows in Serpens seem to dissipate over time.
Generally speaking, within our sample, outflows tended to emanate from young sources with strong millimeter detections and inferred dust masses of hundreds of M ⊕ . However, sources J18291118-0204307, J18292072-0137173, J18293766-0152049, J18302712-0210564, and J18320048-0207440 are not detected in the continuum and have disk dust mass upper limits on the order of only a few M ⊕ . They clearly are visible in 12 CO and have the characteristic collimation of an outflow, as well as emission in the specific velocity range of other outflows from detected sources in our survey. Deeper observations will be required to learn more about their properties.  . Stellar mass as a function of disk dust mass. A correlation is seen between the two quantities, but with a high degree of scatter, similarly to what was found in Law et al. (2017, Figure 5).