© S. E. van Terwisga et al. 2022

1 Introduction

The evolution of the dust in protoplanetary disks is an essential part of our understanding of wider questions of planet formation, and for the earliest development of planetary systems. The Atacama Large Millimeter/submillimeter Array (ALMA) has revolutionized the study of this field with its high sensitivity to the emission from cold, millimeter-sized dust grains in the mid-planes of these disks (Andrews 2020). One way it has done this is by revealing highly structured emission in individual disks at unprecedented resolutions, which show rings, gaps, and spiral arms (e.g., ALMA Partnership 2015; Andrews et al. 2018a). At the same time, the large number of 12 m-telescopes in the array has made it possible to perform sensitive surveys typically targeting around a hundred disk-bearing stars at once, with sensitivities down to an Earth mass or less of dust. Unbiased and (nearly) complete surveys have now been performed in many of the nearby low-mass star-forming regions (SFRs; e.g., Ansdell et al. 2016; Pascucci et al. 2016; Cieza et al. 2019; Cazzoletti et al. 2019), as well as in several of the Orion clouds (e.g., Mann et al. 2014; van Terwisga et al. 2019, 2020; Ansdell et al. 2020). Each of these is usually treated as a snapshot of one set of conditions, and one age bin. Thus, by comparing the continuum emission from disks in different regions, the evolution of dust can be constrained empirically. The observed properties of protoplanetary disks in a given star-forming region seem to be primarily explained by that region’s age, and by the presence of strong UV radiation fields from nearby young O-type stars.

In environments where no such young and massive stars are present and where the evolution of dust in disks is thought to be dominated by internal processes, one key result has been that if disks are optically thin at millimeter wavelengths, after 1–3 Myr, the average disk does not have more than a few Earth masses of dust left for planet formation. This behavior was seen in surveys of disks at submillimeter wavelengths in the Taurus (Andrews et al. 2013), Lupus (Ansdell et al. 2016, 2018), and Chamaeleon I (Pascucci et al. 2016) SFRs, as well as (more recently) the Corona Australis (Cazzoletti et al. 2019) and ρ Ophiuchus (Cieza et al. 2019; Williams et al. 2019) regions. At later times (6–10 Myr), as exemplified by the survey of Upper Sco disks (Barenfeld et al. 2016), the evolution of dust has proceeded even further and led to even lower masses of dust in the median disk than in the previous regions. A highly similar distribution of disk masses was observed in the ~5 Myr-old λ Ori region (Ansdell et al. 2020).

On the other hand, the presence of O-type stars can have a dramatic effect on protoplanetary disks. In the young (1–3 Myr) Trapezium region, clear evidence of an increasing disk mass with distance from the cluster center has been reported (e.g., Mann et al. 2014; Eisner et al. 2018). Beyond ~2 pc from the Trapezium, however, disk masses are similar to those in Lupus and Taurus (van Terwisga et al. 2019). A strong gradient in disk masses was also seen in σ Ori (Ansdell et al. 2017). This is taken as evidence that protoplanetary disk masses are strongly affected, not only by their age, but by the strong far-ultraviolet radiation of nearby O-stars. In the younger NGC 2024, which contains several massive O-type stars, there is likewise clear evidence that external photoevaporation of protoplanetary disks is important and occurs even as early as 0.5–1 Myr (van Terwisga et al. 2020; Haworth et al. 2021).

Some key questions, however, remain open, despite these surveys’ important results. In particular, it is difficult to say if disks are formed with some universal initial mass distribution, which then evolves similarly in the absence of external photoevaporation, or if different star-forming regions have different initial disk properties or mass loss rates even if no significant external photoevaporation is expected. The strong similarities between the disk mass distributions of approximately equal-age Lupus, Taurus and OMC-2 SFRs, and between the older Upper Sco and σ Ori populations, supports this. However, ρ Oph and CrA are also supposed to be young (comparable in age to Taurus, or, in the case of ρ Oph, even younger) and have significantly lower median disk masses and different disk mass distributions: in CrA and ρ Oph, $M_{\text{dust,median}}={0.3}_{-0.2}^{+0.1}{M}_{\oplus }$ and ${0.8}_{-0.2}^{+0.1}{M}_{\oplus }$, respectively, compared to ${3.3}_{-0.6}^{+0.5}{M}_{\oplus }$ in Taurus. At present, it is not clear to what extent this result may be caused by the underlying stellar age distributions in these regions (Cazzoletti et al. 2019; Galli et al. 2020).

This problem is exacerbated by the relatively small number of regions that were studied: even considering the large uncertainties in the ages, the best-populated bin is the 1–3 Myr old range, which contains a handful of comparable SFRs where no O-type stars are present. An additional complication here is that each region includes a quite wide range of physical conditions, such as stellar densities, but not enough objects to study their impact. To give an example, the Lupus SFR’s Class II disk population (N = 98 Class II disks) is distributed over four main star-forming clouds, with most of the disk-bearing stars located in the Lupus III cloud. However, the local stellar density in Lupus III and, say, Lupus II is quite different: Lupus III is sharply concentrated, while the disk-bearing stars in Lupus II are spread out throughout the cloud. However, the sample sizes are too small to be treated independently in each separate cloud. The physical limits on complete sample sizes in nearby star-forming regions are thus a key problem in our detailed understanding of the evolution of protoplanetary disks.

Within Orion, the question of the initial disk mass distribution is also important. After all, it provides the fundamental constraints behind attempts to solve the proplyd lifetime problem (Scally & Clarke 2001), and the star-formation history of the Orion Nebula Cluster (ONC), an important open problem (e.g., Winter et al. 2019). In the OMC-2 region, disk masses suggest comparable properties to Lupus and Taurus at similar ages.

The large number of Class II Young Stellar Objects (YSOs) in Orion A (at least 2400 were identified with the Spitzer Space Telescope, Megeath et al. 2012; hereafter M12) and the sheer size of the cloud, at about 90 pc in length (Großschedl et al. 2018), means that this region is perfect to study the properties of large samples of disks. In particular, the part of the cloud consisting of L1641 and L1647 is a key region for studying the properties of disks at ~1 Myr ages in giant molecular clouds, away from young O-type stars.

Another reason that Orion A so valuable for observations of protoplanetary disks is that there is a wealth of complementary information to facilitate our analysis. Apart from the previously-mentioned Spitzer-based survey of disks and embedded stars, observations with XMM-Newton at X-ray wavelengths (e.g., Pillitteri et al. 2013) and with Herschel (Furlan et al. 2016) mean that both younger and older YSO populations have been mapped. The embedded disks (of Class 0 and Class I-YSOs) have also been characterized by ALMA at millimeter wavelengths by Tobin et al. (2020). The availability of Gaia data has been essential for understanding the spatial and velocity structure of the cloud’s stellar populations, their ages, and the structure of the Orion A cloud itself in greater detail (Zari et al. 2019; Rezaei Kh. et al. 2020). These observations are complemented by accurate very-long-baseline interferometric and spectroscopic surveys like APOGEE-2 (Kounkel et al. 2017, 2018). Given the known correlation between disk and stellar masses (Pascucci et al. 2016), it is also important that the stellar mass distribution in Orion A has been well-studied, and found to be similar to that of nearby SFRs in this part of the cloud (Luhman 2018).

Recently, Grant et al. (2021) surveyed protoplanetary disks in L1641 with ALMA for the first time. Their sample consists of 104 Class II candidates from M12 that were also detected by Herschel at 70 μm. The authors find a median disk mass of ${11.1}_{-4.6}^{+32.9}{M}_{\oplus }$ for their sample, assuming optically thin emission and the same dust opacity of κ_{v=230 GHz} = 2.3 cm² g⁻¹ used in this paper. As the authors suggest, it is likely that the 70 μm-detection requirement introduces a bias toward more massive disks. Taking this into account, it is plausible that Class II disk masses are lower by some amount, with a median disk mass possibly as low as ${0.02}_{-0.01}^{+0.01}{M}_{\oplus }$, if all unobserved disks in the sample are are treated as nondetections with ALMA and assuming a log-normal disk mass distribution. To resolve this uncertainty, an unbiased survey is needed.

In this article, we present observations from an unbiased survey with ALMA of 873 Class II protoplanetary disks in the L1641 and L1647 clouds, the Survey of Orion Disks with ALMA, or SODA. In Sect. 2, the survey design is described. Section 3 details the observation strategy and data reduction process. In Sect. 4, the first results from this survey are shown, including the detection rates, the unbiased catalog, and the first complete catalog of the most massive disks in L1641 and L1647. We discuss environmental variations in disk masses in Sects. 5.1, 5.1.2 and 5.2.

2 Survey design

The aims of this project require the selection of a large, unbiased sample, with uniform selection criteria. As most other similar studies, this survey relies on a Spitzer survey of YSO colors, which is sensitive to the (optically thick) emission from warm dust in the inner disk, conducted by Megeath et al. (2012). Thus, it should be a good indicator of the presence of a circumstellar disk, without being strongly biased to its overall mass. The catalog also includes transition disks, which are characterized by faint or absent emission from the innermost disk parts of the dust disk, but these are not treated separately here. For consistency, we use the same source identifications as in the original paper, abbreviating as [MGM2012] and the number of the object in their catalog.

Fig. 1 Sample positions on the sky. Disks are marked with white plus signs; disks with masses Mdust > 100 M⊕ are shown as blue circles. The major regions within Orion A (L1647, the three subdivisions of l1641, and the NGC 1980 and ONC clusters) are marked separately. Herschel SPIRE observations at 250, 300, and 500 μm (Soler 2019) form the background. The dashed white circle marks a 4 pc radius around θ1 Ori C.

2.1 Sample selection

We selected all Class II YSOs from Megeath et al. (2012) located in Orion A and below −6° degrees in declination, or above ℓ ~ 209.5° degrees in Galactic longitude. As Fig. 1 shows, this simple selection criterion means that we have a large number (N = 873) of protoplanetary disks in our sample, focused on the L1641 star-forming region. It also ensures we avoid the immediate environment of the ONC, with all of our stars located at least 4 pc away from the bright stars of the Trapezium. This ensures that their bright UV radiation is not actively influencing the evolution of the disks we are interested in. This sample should thus be relatively similar to other nearby star-forming regions, except for its large number of sources. YSOs in the Orion A cloud are young (but with local variations), with typical mean ages of ~1–3 Myr (e.g., Da Rio et al. 2016).

The infrared color criteria used to identify the disk sample can misclassify sources based on their inclination. In particular, edge-on objects might appear to be deeply embedded protostars (or be absent altogether), and face-on embedded sources may appear like more evolved disk-bearing objects (Whitney et al. 2003; Crapsi et al. 2008). To correct for this effect, we refined the classifications from Megeath et al. (2012) to include the observations from the Herschel Orion Protostars Survey (HOPS; Furlan et al. 2016) at longer wavelengths, targeting embedded sources. From the original sample of 882 disk sources, nine were reclassified to be younger (Class I or Flat-spectrum). We mark these sources in Table 2. We thus have an 873-source Class II sample.

We expect fewer re-classifications in the other direction. Two of our presumptive Class II objects, [MGM2012] 512 and [MGM2012] 950, have features indicating they may be younger and more embedded sources than typically associated with a Class II spectral energy distribution (SED). Since the evolution of YSOs is a smooth process, even if the transition from Class I to Class II proceeds rapidly when traced with submillimeter luminosity (Tobin et al. 2020; Tychoniec et al. 2020), some such ambiguous objects are expected to be present in a large sample like the one presented here, but not at elevated rates compared to other Spitzer-based disk surveys, and with a limited statistical impact on our conclusions.

In total, after cross-referencing with the HOPS survey, we have measured fluxes (extracted at the phase center of each pointing, assuming unresolved emission) for 873 Class II disks, plus a total of 18 flat-spectrum, 18 Class I, and 4 Class 0 YSOs within the field of view of our observations. The fluxes of the sources that are not considered to be Class II are listed in Table 2. The larger number of sources than pointings (882) is caused by size of the primary beam of our observations. While no disks were observed outside the primary beam, we could target some additional embedded objects. This also means that the younger YSOs have more uncertain photometry, on the whole, due to the decrease in signal to noise away from the primary beam. In the following, we focus on the results for the Class II objects only.

2.2 Catalog completeness

The completeness of the Spitzer-survey has been thoroughly characterized by Megeath et al. (2016). Due to the presence of strong background emission in some parts of Orion, such as the ONC, Spitzer is less sensitive to fainter stars there (which have lower masses). This can, in turn, introduce a bias in our ALMA survey: lower-mass stars bear fainter and lower-mass disks, although with significant scatter (e.g., Ansdell et al. 2016; Barenfeld et al. 2016; Pascucci et al. 2016). However, Megeath et al. (2016) conclude that the completeness fraction is high in regions of faint nebulosity, and where stellar surface densities are low. This applies in particular to L1641 and L1647, where the highest completeness is reached, of >80%. This is also the region targeted by this survey, and indicates that completeness is not a concern for the area covered by SODA. This is consistent with the results from Großschedl et al. (2019), who added a relatively small number of disk-bearing sources throughout the Orion A cloud, but did not find large differences from the Spitzer-based catalog. Additionally, while in the brightest parts of the ONC the 24 μm band saturates, this is not a concern for the area covered by this survey (Megeath et al. 2012).

Finally, we must consider the possiblity that some of the objects in the sample are not YSOs at all. This is the most difficult category of error to quantify. The primary candidates for such false-positive objects in the sample are background galaxies (Gutermuth et al. 2008), IR nebulae associated with YSOs, but not coinciding with the protostar itself within our search radius, and in-telescope effects like diffraction spikes. Since the color criteria used in this catalog are similar to those used to select disk-bearing stars in nearby SFRs, these contamination rates should also be similar across the literature. We do not expect to detect any of the galaxies possibly contaminating our ALMA observations (e.g., Casey et al. 2014). Likewise, an IR nebula associated with a protostar can lead to an offset from the source position that is large enough to lead to a non-detection. The inclusion of these types of objects will therefore always lead us to underestimate of the median disk luminosity. The galaxy contamination fraction of the catalog is expected to be approximately 6.1 ± 3.1 deg⁻², or 24 ± 12 galaxies in total by Megeath et al. (2012) for the area surveyed here and the color cuts employed. In the analysis of the YSO population in Orion A by Großschedl et al. (2019), the authors infer a similarly low contamination rate for these contaminants in their Class II sample spanning the same area on the sky, of ~2.3%, based on the visual identification of contaminants in near-infrared data from the Visible and Infrared Survey Telescope for Astronomy (VISTA; see also Meingast et al. 2016). In any case, contamination should not be a dominant source of nondetections in this survey.

3 Observations and data reduction

3.1 Observations

We observed the 873 SODA disks using ALMA (Project ID: 2019.1.01813.S, PI: S. van Terwisga) between November 21 and December 2 2019 as part of Cycle 7. All observations were carried out with the most compact 12m-array configuration (C43-1) with a minimum of 41 antennas and under average weather conditions with precipitated water vapor (PWV) = 1.5–3.0 mm. All of our targets were observed in Band 6 including 4 spectral windows centered at frequencies of 224, 226, 240, and 242 GHz, respectively, each with total bandwidth of 1.8 GHz and resolution of 31 MHz. The primary beam for this configuration is ~1.1″, with a Maximum Recoverable Scale (MRS) of 10.9″. This choice of frequencies and array setups maximize the detection of line-free continuum emission from point-like sources in the shortest integration time possible.

We observe each target of our sample with a single ALMA 12-m array pointing for a total integration time of approximately 30 sec on source. However, the large number of targets in our sample makes the simultaneous scheduling of our observations unpractical. As part of our observation strategy, we then divided our SODA sample in 6 scheduling blocks (SBs), each typically containing 150 targets selected by proximity in declination. To guarantee the homogeneity of our SODA sample, all sources within a single SB were observed simultaneously. Moreover, each SB was executed a minimum of two times over the observing run. Pointing, bandpass, and flux calibrations (J0725-0054) were carried out at the beginning of each SB. Regular atmospheric and phase calibrations (J0541-0541) were obtained every 8 min.

3.2 Parallel data reduction at the SURFsara Data Processing Facility

Our project develops a new data reduction framework for large ALMA surveys using high throughput data processing facilities. In standard data reduction schemes, both the calibration and deconvolution of ALMA observations is typically carried out in a sequential manner using personal computers or small server machines. This linear and iterative procedure makes the data processing extremely time consuming and largely inefficient for datasets with a large number of targets, such as ours (see Fig. 2a). Instead, we introduce a new parallel scheme for data reduction at the SURFsara Data Processing facilities in Amsterdam (Fig. 2b)¹,².

The original design and observing strategy of the SODA survey allow a direct parallelization of both calibration and imaging steps at the new SPIDER Data Processing Platform at SURFsara³. The new SODA reduction workflow is described in Fig. 2b. First, each of the 6 SB are independently calibrated using facility-provided scripts in order to recreate their corresponding visibilities. Second, each target field is extracted from the visibility file and individually imaged using identical cleaning parameters (see below). Each imaging process is run using an independent realization of CASA v5.6.1 and distributed over multiple cores and nodes by a Grid-based processing environment. We run independent CASA executions to calibrate all our 6 SBs. More importantly, we employ this new distributed scheme to execute up to 150 deconvolutions per SB simultaneously. Following the results of different performance tests, each of these 150 jobs is executed independently using 3 cores, which provide a total of 24GB of memory per job, in order to optimize the CPU usage, as well as reduction and queuing times. The combination of high performance, flexibility, and number of cores at SURFsara-Spider allows us to reduce the total reduction time by 2 orders of magnitude (~6 × 150 = 900 times) compared to standard linear procedures. Independent experiments show that the use of CASA in parallel mode (i.e., mpicasa) does not introduce additional improvements in the SURFsara architecture. Therefore this method is not considered for our data reduction pipeline.

Originally developed for the reduction of single-field continuum data in SODA, this procedure will be implemented for the analysis of mosaics and spectral cubes in the near future (Hacar et al., in prep.). The reduction of the entire SODA sample required a total of ~3000 CPU hours. After optimization, our new parallel processing allows to run the full calibration and imaging processes of the entire SODA sample in less than a day. These data – containing many independent pointings with relatively simple sources – are particularly suited to this type of approach. A similar reduction strategy has been employed for the data reduction of large radio surveys obtained by interferometers such as the Low-Frequency Array (LOFAR; Mechev et al. 2017; Shimwell et al. 2017) or Very Large Array (VLA; Lacy et al. 2020). To our knowledge, this is the first time this type of parallel data processing is used in the case of ALMA observations.

The optimal use of parallel CASA executions requires the standardization of the imaging process. For all our targets, we combine the four available continuum spectra windows in order to maximize the sensitivity of our maps, but did not further spectrally average them. We deconvolved each of our target fields with the CASA task tclean using standard Briggs weighting and a robust parameter of 0.5. Each of the resulting images is then corrected by their corresponding primary beam response. The individual maps show typical synthesized beam sizes of ~1.4″ × 1.1″, that is, approximately 580 AU × 414 AU at the distance of Orion. The typical noise level of our maps is 0.08 mJy beam⁻¹ which translates into a median sensitivity threshold of 1.5 M_⊕ at the 4σ level for our adopted opacities (see Sect. 4.1.1). This sensitivity is lower than that achieved in nearby low-mass star-forming regions, which achieve similar sensitivities of 0.2 M_⊕ (e.g., Ansdell et al. 2016; Pascucci et al. 2016), but sufficient to detect the majority of sources in such regions assuming similar disk flux distributions. After imaging, we inspected the images for issues, but did not identify any. All our continuum maps are available as png figures and in FITS format on our webpage⁴.

Fig. 2 Data processing workflow for SODA. (a) Standard data reduction. (b) Parallel ALMA data reduction scheme implemented on the SPIDER platform that is part of SURFsara’s Data Processing facilities.

3.3 Flux extraction

To identify disks, we search within a 1.4″ radius around the catalog position for the brightest pixel in the reconstructed image. We use a 4σ detection limit, in order to ensure the number of false positive detections is <1 for our dataset. This false positive rate includes both pure noise and the chance alignment of an IR excess source in the input catalog with a submillimeter galaxy in the background, following Sadavoy et al. (2019). The assumption of unresolved flux is based on observations of nearby star-forming regions at higher resolution, which find median dust emission extents of ~60 AU, while the largest objects are no larger than ~250 AU, and exceedingly rare (≲1% of objects, e.g. van Terwisga et al. 2018; Long et al. 2018; Cieza et al. 2019). To test the validity of this hypothesis, we fit double Gaussians to the emission, using the location of the peak emission as a prior. We identify only 14 objects which have SEDs consistent with Class II objects and best-fit axis ratios inconsistent at the 4σ level with the assumption of unresolved emission. These objects, and the results from the fitting procedure, are shown in Table A.1 and Figs. A.1–A.4 and are marked out in Table 1. Additionally, we discuss two further sources with resolved emission but unresolved central sources in more detail in Sect. 4. Two types of objects are present: sources which seem to be well-fit by a double Gaussian, and which most likely represent objects with some contribution from (resolved) envelope emission, and five sources with residuals suggesting a compact (visual) multiple system ([MGM2012] 129, 754, 1000, 1090, and 1094). We have kept them in the catalog presented here to ensure consistency with regards to the source selection, and because of the small number of such objects. The impact of this flux extraction approach on multiple systems in particular is discussed in the following Section.

4 Results

The SODA survey covers a large number of sources at low resolution. As a result, a wide range of stellar properties and environments are covered. In total, we detect 502 of 873 sources (58 ± 2.6%) in the sample. Comparing this result to nearby low-mass SFRs, we find that the detection rate in SODA is in between Lupus (66 ± 6%) and Cha I (47 ± 5%), cf. Fig. 13. In comparison to the Class I sources in Orion presented by Tobin et al. (2020), which have detection rates >85% in a shallower survey, this is consistent with a continuous mass loss for the circumstellar material.

Figure 3 shows a number of individual detections, both of millimeter-bright disks (in the top two rows) and fainter (but still detectable) disks, and their stellar hosts. The top panel shows, for each source, the ALMA contours on top of near-infrared (J, H, K_s) images from the Vienna Survey in Orion (VISION) using the VISTA telescope (Meingast et al. 2016). This survey has a higher resolution than the input catalog, which is limited by the Spitzer beam. From the figure it is clear that the sample covers a wide range of disk masses. As previous work has shown, there is a correlation – albeit with a large scatter – between the masses of disks and their host stars (e.g., Andrews et al. 2013; Ansdell et al. 2016; Pascucci et al. 2016). Thus, before comparing the SODA sample to other star-forming regions, we should confirm the similarity of the underlying Initial Mass Function (IMF) of L1641 and L1647. L1641 is deficient in B- and late-O-type stars (Hsu et al. 2012, 2013), but similar to other nearby SFRs at lower masses. This means the sample is truly representative of that in other regions, where such high-mass stars are also often absent, and M- and K-type stars dominate the stellar mass distribution.

Multiplicity may also affect the properties of disks. In Fig. 4, a selection of visual multiples with at least one ALMA continuum detection of a disk are shown. Unfortunately, multiplicity is not a well-constrained quantity for the SODA survey sample. Both the input catalog resolution of 1″ and the resolution of our ALMA data prevent us from resolving the most compact binaries (with separations ≲400 AU). For such objects, we effectively consider the system’s average continuum luminosity. The cumulative distribution of total system masses for the large sample of multiple systems in nearby SFRs presented in Zurlo et al. (2020) lacks massive (>50 M_⊕) systems, but the mass distributions are otherwise identical. Thus, we do not expect the way we treat these systems and the marginally resolved ones shown in Appendix A to affect our estimates of median disk luminosities and masses.

In cases where apparent (visual) companions are present, at separations <10” (4000 AU), we have included both sources if both sources have Spitzer photometry indicating a circumstellar disk, as in [MGM2012] 55 and 56 and [MGM2012] 573 and 574 in Fig. 4. Alternatively, if only one source is included, we only include that source. This is for instance the case for [MGM2012] 227 and 975, which are well-resolved, and for the marginally resolved [MGM2012] 129 and 421. If we assume that the trends of disk properties with multiplicity as demonstrated by Harris et al. (2012) and Akeson et al. (2019) in Taurus, and Zurlo et al. (2020, 2021) in Lupus and Ophiuchus also apply in Orion, we do not expect this approach to lead to a (strong) bias. The primary star’s disk in a wide binary (>300 AU in Harris et al. 2012, > 140 AU in Akeson et al. 2019) has a similar luminosity distribution to the disks in single Class II YSOs. Due to the resolution of our data, this separation criterion is always met for (marginally) resolved systems. While a detailed confirmation is beyond the scope of this article, we may therefore expect no large bias will arise from the approach used here. Moreover, not all visual multiple systems are true multiples: they may merely be chance alignments. For objects like [MGM2012] 714’s companion candidate, which does not have an IR excess, adding nondetections for the visual companions leads to a more inconsistent and biased catalog.

Apart from unresolved disks, we also identify two Class II sources associated with extended millimeter-continuum emission. Fig. 5 shows these objects: [MGM2012] 512, which was previously identified by Grant et al. (2021), and [MGM2012] 950, where to our knowledge no extended emission has previously been identified. Both of these sources show striking and, unexpectedly, resolved continuum emission. [MGM2012] 512 is particularly noteworthy, and appears to have an arc of dust close to the star. Such extended emission is not expected to be found in Class II sources, which are expected to have lost their envelopes, although they have counterparts in younger, more embedded objects (Pineda et al. 2020). As Grant et al. (2021) discuss, however, this could be due to either a late accretion event, a perturbed disk (although no suitable candidates can be identified in either 2MASS or VISION), or a younger, still partially embedded protostar viewed from nearly face-on. Face-on embedded protostars are indeed expected to have much bluer spectral indices than those seen at higher inclinations (Whitney et al. 2003; Crapsi et al. 2008). As Grant et al. (2021) discuss, the classification of the source is ambiguous in the literature, and Caratti o Garatti et al. (2012) consider it to be a Class I object. On the other hand, the late accretion scenario has also been proposed to explain the scattered-light and submillimeter molecular line observations of an arc of material accreting onto the Class II SU Aur disk (Ginski et al. 2021).

A misclassification is the most likely explanation for the nature of the [MGM2012] 950 source: to the northwest of the source, there is tentative evidence of an outflow-like structure in the VISION data. However, the Class I [MGM2012] 996 source, which we observed in this survey (but do not include in our statistical analysis due to its classification) has extended continuum emission but a quite different, more centrally concentrated morphology at the same resolution.

Regardless of the origin of the observed extended dust tails seen in these two sources, the advantage of the size of this homogeneous dataset is that it is possible to say that this kind of morphology must be rare, occurring only 0.23 ± 0.16% of the time. Thus, if these structures are the result of younger sources contaminating the sample, this demonstrates that such contamination is not common. This is further supported by the small number of objects with possible resolved emission in Appendix A: the most likely cause for this extended emission is the presence of a (remnant) envelope. In either case, these objects with significant large-scale continuum emission are excellent candidates for follow-up studies of how circumstellar material accretes onto the star and disk, both at millimeter wavelengths and in scattered light.

Fig. 3 Images of selected millimeter-detected sources. For each source, the top panel shows near-infrared (NIR) data (J, H, and Ks-band) form the VISION survey (Meingast et al. 2016) in colors, and 5, 50, and 100σ contour levels. The bottom panel shows the ALMA image at 225 GHz only. White crosshairs mark the positions of detected sources. The top panels show sources with fluxes >4.8 mJy, bottom panels sources with fluxes < 1.6 mJy. The beam is shown in white in the ALMA images, in the lower left.

4.1 Disk luminosities and masses

In this subsection, we discuss how we inferred the millimeter luminosities and masses from the fluxes measured with our survey, as listed in Table 1. Generally, disk masses are inferred by assuming that their millimeter continuum flux is optically thin, with a known opacity and distance. In this case, the mass is found from: $$M_{\text{dust}}=\frac{F_v{d}^2}{k_v{B}_v\left(T_{\text{dust,eff}}\right)}$$(1)

We first discuss the relevant terms and assumptions in this equation. Inferring the disk luminosities only requires assuming one parameter, the distance, which is discussed in Sect. 4.1.2. For disk masses, several additional assumptions must be made. In Sect. 4.1.1, we detail our choices in this regard, and how they may impact the comparison with other surveys in the rest of this article.

4.1.1 Temperature and opacity assumptions for disk masses

For this survey, we assume T_dust,eff = 20 K. This is based on the result in Andrews & Williams (2005) for disks in Taurus, and followed by most other disk surveys. However, this temperature implicitly carries with it assumptions on disk radius: more compact disks will have higher effective temperatures. Likewise, large radiation fields (such as in the ONC) or simply hotter host stars will lead to intrinsically hotter disks. However, this survey’s sample, located far (≥4 pc) away from the massive stars in the Trapezium, should not be truncated by external photo-evaporation, nor be intrinsically hotter than those of other SFRs. Even assuming a higher temperature of for instance 30 K will, however, not have a large effect on the measured dust mass.

Likewise, the assumptions of optically thin emission and the opacity κ_v assumed in the calculation of disk masses are a source of much unquantified uncertainty. There is now increasing evidence that part of the disk can be optically thick even at the wavelengths used for this survey (e.g., Andrews & Williams 2005; Tripathi et al. 2017; Zhu et al. 2019; Macías et al. 2021, the latter of which find optically thick emission even at 2.1 mm). The net effect of this is that we possibly underestimate our disk masses to some (model-dependent) extent.

Opacities for the millimeter-sized grains can also vary. Here, we used common assumptions on the opacities of the grains, derived from Beckwith et al. (1990): κ_{v=230 GHz} = 2.3 cm² g⁻¹, and a power-law frequency-dependence with ß = 1. In the recent literature, these are commonly-used values (e.g., Ansdell et al. 2016, 2017, 2020; van Terwisga et al. 2019, 2020), but Barenfeld et al. (2016) and Pascucci et al. (2016) both used ß = 0.4 for Band 7 observations. The effect of this, however, is small, especially in the light of other existing uncertainties on disk temperatures, radii, and optical depths; see also van Terwisga et al. (2019).

Fig. 4 Images of selected (candidate) multiple systems. For each source, the top panel shows NIR data (J, H, and Ks-band) form the VISION survey (Meingast et al. 2016) in colors, and 5, 50, and 100σ contour levels. The bottom panel shows the ALMA image at 225 GHz only. White crosshairs mark the source positions of detected sources. Nondetections are marked with grayed-out names. The beam is shown in white in the ALMA images, in the lower left.

Fig. 5 Images of YSOs associated with extended millimeter emission. For each source, the top panel shows NIR data (J, H, and Ks-band) form the VISION survey (Meingast et al. 2016) in colors, and 3, 5, 50, and 100σ contour levels. The bottom panel shows the ALMA image at 225 GHz only. White crosshairs mark the source positions of detected sources. The beam is shown in white in the ALMA images, in the lower left.

Fig. 6 Disk mass distributions inferred with the Kaplan-Meier estimator in Orion A and several nearby star-forming regions. SODA (dark blue) versus Class I and Class 0 disks in Orion (Tobin et al. 2020) (orange and green), and disks in Lupus (Ansdell et al. 2016) (purple), Chamaeleon I (Pascucci et al. 2016) (brown), ρ Oph (Williams et al. 2019) (gray), and Upper Sco (Barenfeld et al. 2016) (yellow).

4.1.2 The distance structure of Orion A

The assumed distance has a potentially important effect on the conversion between fluxes and masses in Eq. (1), and not only between SFRs. Given the size of the Orion A cloud, which spans multiple parsecs projected on the sky, we must be careful to take into account the three-dimensional structure of both the cloud and its young stellar population as best as possible, and literature results suggest that this is indeed quite complex (Kounkel et al. 2017, 2018; Großschedl et al. 2018; Zari et al. 2019). In particular, Großschedl et al. (2018) find quite pronounced distance differences between the southernmost and central parts of Orion A, of some 80 pc.

With the advent of the Gaia satellite, the distances and proper motions of nearby stars are better-known than ever before. However, the generally high optical extinctions mean that not all Class II stars in our sample have high-quality parallaxes in DR2 or EDR3. To solve this problem, we use the same approach as Großschedl et al. (2018): YSOs with well-known distances from DR2 are placed in 16 bins in galactic longitude (oversampled by a factor of two), and the distances of the YSOs in our sample are interpolated based on the distances of these bins. After correcting for the distances structure, we are left with a residual uncertainty in the distance of the sources of ~5% (20 pc), with small variations in galactic longitude. This reflects both the uncertainty in the Gaia data and the intrinsic three-dimensional structure of the cloud.

4.2 Disk mass distribution in SODA

Figure 6 shows the inferred disk mass distribution using the Kaplan-Meier estimator for the full SODA Class II sample of 873 disks. We find a median disk dust mass of ${1.97}_{-0.19}^{+0.15}{M}_{\oplus }$. This figure also compares the observed disk mass distribution to several well-characterized populations of Class II disks in nearby star-forming regions, as well as to the Class I and Class 0 disk mass distributions in Orion A and B combined. Throughout this section and in the following, we use the lifelines Python package (Davidson-Pilon et al. 2020) for fitting the (censored) disk mass distributions. For log-rank correlation tests, we follow Feigelson & Nelson (1985) in subtracting our data from a constant, since lifelines only implements this test for right-censored data.

The large sample size allows us to estimate the properties of this distribution with unprecedented accuracy. Errors in the probability distribution are on average ±0.017, compared to for instance ±0.05 for Taurus and Lupus, an improvement of a factor of 3. Moreover, it is noteworthy how well-sampled even the upper end of the disk mass range is: SODA characterizes even the occurrence of the most massive (M_disk > 100 M_⊕) well, despite their rarity; see also Sect. 4.4.

Typically, in the literature, protoplanetary disk mass distributions are described using the Kaplan-Meier estimator, which is a powerful, nonparametric method for describing censored random variables. However, it is – formally – only useful if the censoring is not dependent on the variable in question. In these types of surveys, that is the case, leading to potential biases. These biases will primarily influence the lower end of the disk mass range probed by a survey, where incompleteness is the largest. Since our survey depth is similar to that of other works, we conclude that our results should be comparable to those of other surveys to approximately the same masses.

Instead of using nonparametric descriptions of the disk mass distribution, however, we may also – if the nature of the distribution is known a priori – use parametric methods to compare regions. As Williams et al. (2019) showed, many protoplanetary disk mass distributions appear to follow log-normal distributions, with highly similar standard deviations, but potentially quite different median disk masses. If we can fit log-normal distributions to our observed disk populations, the mass-dependent censoring problem is less relevant. This is particularly important for samples with fewer sources, as well as for the faint tail of the disk mass distribution. Moreover, fitting a log-normal distribution allows us to directly compare (including statistical errors) the median disk masses of different regions. However, it is important to test if this parametric distribution is indeed a good description of the property we study here.

Interestingly, Fig. 7 shows that SODA sample is not best described by a log-normal distribution, but by a Weibull distribution. Like the log-normal distribution, this distribution has two parameters; its cumulative distribution function is 1 – exp −(x/λ)^k. That the Weibull distribution is indeed a better fit is confirmed from a comparison of the Akaike Information Criterion (AIC): ∆_AIC = 16 in favor of the Weibull distribution. Figure 7 shows the quantile-quantile plots of the empirical distribution of source masses versus a Weibull (top panel) and log-normal (bottom panel) distribution fit to the data. In this figure the axes are linear, not logarithmic, in disk mass. It is clear that the difference is driven by the better fit of the Weibull distribution at masses > 100 M_⊕: the SODA sample has a more massive tail than a log-normal distribution can account for. However, the values for the median dust mass agree to within their respective errors: these are 2.2 ± 0.19 M_⊕ for the Weibull, and ${2.2}_{-0.17}^{+0.17}$ for the log-normal distribution. The differences between the information criteria for these two-parameter distributions are also small. Indeed, neither distribution perfectly describes the intermediate-mass regime: there is more complexity hiding in these data.

Similar, but much more extreme, behavior is seen for the disk masses in the σ Ori cluster (Ansdell et al. 2017). The authors find a steep disk mass gradient in space in that region. As a result the disk mass distribution is likewise not well described by a single log-normal disk mass distribution, but shows an excess of massive sources. This suggests that inhomogeneity in the underlying disk sample may drive the observed behavior away from a single log-normal distribution in the SODA sample as well. In Sect. 5, we investigate possible drivers of this behavior in the SODA sample in greater detail.

4.3 Comparison to literature results

The sample sizes in SODA allow us to split up the sample by the subclouds in Orion A, L1641 and L1647, respectively. Similar data to those presented here (in terms of resolution and frequency) exist for a subset of the SODA sample in L1641, by Grant et al. (2021). That sample was based on Class II objects detected with Herschel in the HOPS fields (Furlan et al. 2016). Surprisingly, the disk mass distribution inferred by these authors has a very high median mass, even including nondetections, which they consider to be the likely consequence of a selection bias. In the most pessimistic case, considering all nondetections in Herschel-bands to be nondetections with ALMA, the disk mass distribution for L1641 could be similar to that of Upper Sco, with a very low detection rate. However, from this sample alone the bias cannot be accurately determined.

In Fig. 8, the SODA data are split up by cloud and compared to the observed sample from Grant et al. (2021). Two results are readily apparent (and confirmed statistically, using a cutoff of p = 0.05 in a log-rank test): L1641 and L1647 have highly similar disk mass distributions, and the Herschel-detected subsample in L1641 (shown in pink) is indeed biased toward higher disk masses. At the same time, it is clear that a significant fraction of sources not detected by Herschel is still detectable by ALMA at the effective sensitivity of these observations, as shown by the difference between the gray distribution and the SODA disk mass distributions.

Fig. 7 Quantile-quantile plots of the empirical distribution (inferred with a Kaplan-Meier estimator) and fitted quantiles using a Weibull distribution (top panel) and log-normal distribution (bottom panel).

Fig. 8 Comparison of SODA disk masses to previous literature. Disk mass (or luminosity) distributions in L1641 (blue) and L1647 (green) from SODA, as well as the observed disk mass distribution in the sample from Grant et al. (2021) (pink) and a scenario in which all of the unobserved sources in Grant et al. (2021) are also nondetections at their sensitivity (gray).

4.4 The most massive disks in Orion A

The most massive protoplanetary disks in nearby star-forming regions have been a crucial population in the study of disk evolution. These objects are the easiest to observe at high resolution, and may correspond to stars with optically thick material at the largest radii (Andrews et al. 2018a). As such, they allow us to cast light on the formation of, potentially, distant planets, and planetary system architectures that are difficult to observe in main-sequence systems. These massive disks are also relevant since they are the most plausible systems for identifying signatures of gravitational instability. This is especially the case if they are radially extended as well as massive (which is suggested by the disk luminosity-radius relation, see e.g. Andrews et al. 2018b). Such instabilities have now been suggested to occur in a handful of massive nearby disks (e.g., Huang et al. 2018; Booth & Ilee 2020). As such, it is a sample that is also particularly interesting for follow-up observations with ALMA.

Figure 1 shows the locations of all disks with masses (assuming optically thin emission) higher than 100 M_⊕ in dust in Orion A. We find 20 such objects, larger than any other such sample, and the first complete sample of massive protoplanetary disks located within the same star-forming region. The most massive object is [MGM2012]-540, with M_dust = 364 ± 0.5 M_⊕. In L1641, Grant et al. (2021) previously identified only 6 disks with masses more than 100 M_⊕; this survey increases that number to 16 in the same part of the cloud, emphasizing the importance of an unbiased survey for characterizing the upper end of the disk mass distribution.

No clear pattern is visible in their occurrence across the cloud, with the possible exception of a smaller fraction of such sources in L1641-C and L1641-N. Interestingly, these sources also seem to be found only where the cloud is densest, suggesting that they may be younger on average than the full disk sample. However, the relatively small number (2%) of such sources means that it is difficult to identify significant patterns in their occurrence. If this subsample is indeed younger than the median disk in the sample, it may not be representative of the evolution of circumstellar material in the median disk in this star-forming region.

5 Discussion

Until now, we have considered the disk mass (or luminosity) distribution of the SODA sample as a whole. There are several reasons to assume that disk masses may not be homogeneous throughout the cloud, however. Purely observationally, as demonstrated in Fig. 7, there is the massive tail, which may result from the superposition of a lower- and higher-mass population. Testing this hypothesis requires us to study the sample at higher resolution, and to investigate the behavior of disk masses in subsamples. By design, SODA contains a sufficiently large number of sources to make this possible. By forming subsamples of ~100 sources, we can still obtain median disk masses at similar precision to complete surveys of nearby star-forming regions. In particular, in this section we investigate the median disk mass along the length of the cloud (in Sect. 5.1) and in the different clusters of young stellar objects found in L1641 and L1647 (in Sect. 5.1.2). Finally, after investigating the presence of local variations in the median disk mass, we compare the masses of disks in SODA to those of other samples of protoplanetary disks in Orion A and B (Sect. 5.2.2), and to nearby star-forming regions (Sect. 5.2).

Two main effects may lead to spatial variations in the median mass of the protoplanetary disks on the cloud. First, if all disks evolve at similar rates and originate from the same disk mass distribution, we would expect to see mass variations tracing the median age of stars in the cloud, but this is not necessarily constant over the ~100 pc length of Orion A. Indeed, there is evidence for such age variations and a quite significant intrinsic age spread in Orion A (Da Rio et al. 2016). To this end, we compare the disk masses to the ratios of Class II to embedded (Class 0 and I) and Class III YSOs, respectively, which are proxies of the mean stellar age along the line-of-sight. In this section, we use the Pillitteri et al. (2013) catalog for Class III source numbers, and Tobin et al. (2020) for all embedded sources (Class I and 0). Both of these surveys target only L1641. While we show results for L1647 and compare them to the embedded source numbers (derived from Megeath et al. 2012) we are less complete to these kinds of sources there.

A second effect, the existence of variations in the median disk mass between different SFRs, is currently a topic of intense debate as a result of the growing number of surveys of protoplanetary disks in nearby regions. Data from the CrA and ρ Oph disk populations (Cazzoletti et al. 2019; Williams et al. 2019) suggest such variations may indeed exist, but it is not clear at which spatial scale such variations occur, below the largest cloud level.

5.1 The distribution of disk masses across the Orion A cloud

The variation of protoplanetary disk masses as a function of galactic longitude is presented Fig. 9 (left panel), where we binned the survey in seven independent longitudinal bins. Since the “spine” of Orion A is essentially parallel to the galactic longitude (see also Fig. 1), this is a simple measure of how masses change along the cloud. Median disk masses are derived by fitting a log-normal distribution to the mass distribution of objects in each bin, with errors reflecting the errors on the fit. The widths of the distributions are highly similar, and are not shown in the figure but presented in Table 3.

Two results are immediately apparent from Fig. 9. The disk masses in L1647 (corresponding to the first four bins) are statistically indistinguishable, despite spanning a significant distance of tens of parsecs when considering the three-dimensional structure of the cloud (Großschedl et al. 2018). Second, it appears that below about ℓ ~ 212° degrees in galactic longitude, disk masses steadily decrease. Around ℓ = 212° the disk masses seem to have a local maximum. Below ℓ = 212.6°, a linear decrease in median disk mass is marginally favored over a constant median disk mass, with corrected Akaike Information Criteria of AICc = 14.2 versus AICc = 17.6. Moreover, from a log-rank test the disk masses at ℓ = 209.5° and ℓ = 210.0° are inconsistent with those in the ℓ = 214.2° and 213.1° bins (p < 0.05). Overall, however, no large variations in the median disk mass are visible, with contrasts of a factor ~few between the most and least massive bin.

Fig. 9 Trends in disk properties in SODA as a function of galactic longitude. Left: median Class II disk masses (black points) for all longitudinal bins. The location of the Trapezium is marked with a blue star. Right: ratios of embedded (red) and Class III YSOs (purple) to Class II sources in each bin. Round symbols show source fractions derived from Pillitteri et al. (2013) and Tobin et al. (2020) for Class III and embedded sources, respectively; square symbols show the embedded source fraction in L1647, derived from Megeath et al. (2012).

5.1.1 Impact of age gradients

To test how the stellar age along the line-of-sight changes with galactic longitude, the right panel of Fig. 9 shows the fraction of Class II sources relative to embedded (Class 0, I, and flat-spectrum) and to Class III sources in each longitudinal bin. It is clear that, toward the ONC, the number of Class III sources increases sharply, while the embedded source fraction decreases. We also see a decrease in extinction A_Ks toward these latitudes, as listed in the last column of Table 3. This behavior is coincident with the trend in disk masses along the cloud. Together, these observations suggest that the median age along the line of sight increases toward the ONC, as was indeed concluded by Pillitteri et al. (2013). However, if the Class III population is contributed by a population that is old enough, and distinct from the population of Class II and embedded sources, we would not necessarily expect to see any change in disk masses. Figure 9 also shows that the embedded source fraction declines below ℓ = 212°. This cannot be ascribed to only a completeness effect of the Spitzer MIPS 24 μm band, as it only saturates in the very inner part of the ONC.

There is a considerable amount of literature on the physical and age structure of the young stellar populations toward Orion A which can help to shed light on the question if an older population is indeed present toward ℓ = 209.5°, and how it relates to the bulk of the young stellar objects in L1641. Alves & Bouy (2012) and Bouy et al. (2014) suggested a large population of stars associated with NGC 1980 but kinematically similar to L1641 is present in this area, with an age of 5–10 Myr. However, the existence of such a population has been challenged since then: Da Rio et al. (2016) suggest that while older stars are present and have lower extinctions, they are not completely separate from the cloud. On the other hand, an older foreground is not retrieved at all by Fang et al. (2017) and Kounkel et al. (2017), who conclude that the foreground population is indistinguishable both in terms of kinematics and age. A population of stars belonging kinematically to Orion D may be present in this area per Kounkel et al. (2018). A similar group of stars is identified by Zari et al. (2019) and in Chen et al. (2020); these suggested populations have ages in the ~10 Myr range. Combining these results with a three-dimension extinction map, Rezaei Kh. et al. (2020) conclude that the evidence favors an older foreground population along these lines of sight.

Crucially, and regardless of precise three-dimensional structure, the age ranges quoted for the older part of the stellar population in this area of Orion A suggest that we can expect ~10% of these stars to have a near-infrared excess that would lead to these stars being included in this survey (Hernández et al. 2008; Luhman 2022). At the same time, disks at this age are faint at millimeter continuum wavelengths (Barenfeld et al. 2016; Ansdell et al. 2020). Depending on the size of the older population, then, we can explain the behavior observed in Fig. 9: along the line-of-sight toward ℓ = 209°, the median age and the width of the age distribution are larger than at higher latitudes, so the fraction of embedded sources to Class II objects decreases, while the Class III to Class II ratio climbs.

To explore the hypothesis that a mixed population along the line of sight leads to a lower median disk mass and higher Class III to Class II ratio more quantitatively, we artificially construct a sample of disks by randomly sampling (with replacement) disks from both the ℓ > 213.1° area in Orion, and from the Upper Sco catalog published by Barenfeld et al. (2016). This approach thus requires the assumption that disk evolution proceeds similarly in these star-forming regions. For a ratio f = 0.68–0.84 between these populations, the median disk mass of the mixed sample is indistinguishable from ℓ = 209° disks, to within 1σ. This approach is shown in Fig. B.1. Using a Class III to Class II ratio for Upper Sco of 4.2 (Luhman 2022) and ~1 for the southern L1641 sample leads to an N_III/N_II ratio of ~2 for the mixed sample, similar if slightly lower than that observed. However, it is important that the fraction of diskless young stars is determined in different ways in both regions, adding additional uncertainty to this estimate. Ultimately, the limited statistical significance of the median disk mass variations in L1641, the assumptions on the disk mass distribution for 5–10 Myr old disks, and in particular the uncertainties in the (relative) stellar ages of the stars with protoplanetary disks in L1641 prevent us from drawing stronger conclusions.

Alternatively, the disk masses in the lowest galactic longitudes might be affected by external photoevaporation. This process is expected to be an important cause for disk mass loss in the immediate vicinity of the Trapezium (Mann et al. 2014; Eisner et al. 2018). However, it seems to be less important even at only ~2 pc distances in the OMC-2 region (van Terwisga et al. 2019). The SODA sample is separated by at least 4pc from the Trapezium. However, it is possible that the proximity of ι Ori (an O9III star) to L1641N does contribute to lowering the masses of the disks there.

5.1.2 Disk masses in YSO clusters

As is clear from the distribution of sources in Fig. 1, Class II YSOs are not spread out uniformly across the Orion A cloud. Even in our sample, which focuses explicitly on the disks located away from the ONC, several regions with higher densities are readily apparent by eye, and such clusters have been identified in this sample systematically by Megeath et al. (2016).

Considering the properties of disks from the perspective of these associations gives us valuable information. Even disregarding the possible foreground contamination in the northern part of our survey area, stars currently in overdensities should be more similar to each other than stars not in overdensities – a result that was previously demonstrated in Megeath et al. (2016). The field population can, after all, have originated in earlier, now dissolved, overdensities, and should thus at least have a larger age dispersion (Kainulainen et al. 2017). Moreover, clusters are less likely to have significant structure along the line-of-sight.

We restrict ourselves to finding clusters in two-dimensional (RA, Dec) space, as a first-order approach that has been successfully used for this dataset before. Even with Gaia EDR3 the fraction of sources with accurate astrometry is fairly low: 284/873, using version 1.0 of the astrometric fidelity classifier presented in Rybizki et al. (2022), which means this is the most practical approach that ensures sufficiently large sample sizes. To search for clusters, we used the OPTICS algorithm (Ankerst et al. 1999). The data distances were scaled, and we required a minimum of five samples in the neighborhood for core points to be considered, as well as a minimum reachability steepness of ξ = 0.05. We additionally set a minimum cluster size of 50 members of Class II or earlier. This choice of parameters allows us to minimize the number of interlopers, while retrieving similar clusters to those reported originally in Megeath et al. (2016), and ensures accurate estimates of the median disk masses of the cluster members. While we use the term cluster, we do not expect any of these to be gravitationally bound to any significant degree. The stellar surface densities of these clusters are too low for that, with between 10–30 YSOs pc⁻² (Megeath et al. 2016).

The resulting 6 clusters and the unclustered field population are shown in Fig. 10. It is clear that sufficiently large overdensities of YSOs are present throughout the cloud. However, toward L1641N, the ratio of cluster members to field stars begins to drop rapidly. We probe a range of surface (over)densities, with Cluster 2 being the least dense. Compared to the clusters found in Megeath et al. (2016), who used a different cluster-finding approach, we see that they do not retrieve cluster 5, while cluster 2 is broken up into two smaller clusters, one of which is (significantly) smaller. Since Megeath et al. (2016) do not require a minimum cluster size, they also find several other small clusters and groups throughout.

In Fig. 11, the median disk masses as inferred from a log-normal fit to the disk mass distribution of each cluster’s members are presented, as well as the average disk mass for all field members. Additionally, Table 4 shows the number of members, median masses, and widths of the log-normal distributions for the clusters and field. The latter is not significantly different from the global average, which is not surprising, since field-disks (508) outnumber cluster-disks by a significant amount, with a median of 59.5 members per cluster.

Figure 11 shows that the differences between disk masses in clusters are not particularly large, like those seen when longitudinally binning the disks along the cloud. However, some variation does exist. Clusters 1 and 2, both located in L1647, have similar median disk masses, and are each statistically indistinguishable from the field. Together, these clusters contain the majority of the Class II sources in L1647 (138 of165).

In L1641, clusters 3, 5, and 6 have median disk masses that are again similar. For all three of these clusters, however, the median disk mass is lower than the field, although this is significant only for cluster 3. Since clusters 5 and 6 are coincident with the part of L1641 where, above, we concluded overall ages may be higher, it is perhaps not surprising that their median disk masses are also lower, especially since we projected the three-dimensional structure of the cloud into two dimensions. It is consistent, in other words, with the cluster approach partially compensating for a (more dispersed, older) foreground population.

Comparing the clusters to the embedded and Class III object populations provides additional insight into what causes the differences between the densest groups of Class II YSOs. Again, we see that the ratio of Class III to Class II sources varies significantly. Clusters 3 and 6 coincide with relatively high numbers of Class III objects, which may provide an explanation for their lower disk masses relative to the field. The differences in N_embedded/N_{Class II} and N_{dass III}/N_{dass II} between clusters 4 and 5 are small. Cluster 4 contains relatively more young, embedded objects. This may be seen as evidence that it is the youngest of the clusters we identified here, and could explain at least in part the relatively high median disk mass of this cluster.

We must stress that the differences between these clusters, which are located far from each other on the sky, are small and not highly significant relative to each other and to the field, even in the most extreme cases. In the surveyed part of Orion A, where we avoid the presence of O-stars, it appears to be the case that disk evolution as traced by millimeter-continuum emission proceeds similarly throughout the cloud when studied with sample sizes comparable to those found for nearby SFRs, and especially when taking into account the fractions of embedded and Class III sources as a proxy for age. If variations in the initial disk mass (or rate of evolution) in the part of Orion A surveyed here do exist, they are not significant with the present analysis, although we cannot exclude their existence completely.

These observations are consistent with the expectations from theoretical work: in environments like the low-mass, low-density clusters we are probing here, the disk evolution models run by Concha-Ramírez et al. (2019) predict no impact from photoevaporation on the disk population. This effect is made even stronger by the lack of O-type stars in this part of Orion A (Hsu et al. 2012, 2013). In the following, we compare the SODA sample to other nearby star-forming regions, in order to see if the observed similarity between different groups of young disk-bearing stars holds at scales larger than GMCs.

Fig. 10 Clusters identified by OPTICS. Cluster members are indicated by colored circles, while field stars are shown as gray plus signs.

Fig. 11 Median disk dust masses in clusters and the field population. Left: median masses for the individual clusters and field population (black) and the median and standard deviation for the clusters together. Right: ratio of the numbers of embedded (red) and Class III YSOs (purple) to Class II sources in each cluster. Round symbols show source fractions derived from Pillitteri et al. (2013) and Tobin et al. (2020) for Class III and embedded sources, respectively; square symbols show the embedded source fraction in L1647, derived from Megeath et al. (2012).

Fig. 12 Comparison of the SODA results to three nearby SFRs. Kaplan–Meier estimators for the cumulative disk mass distribution in L1641S (reduced) (orange), L1641N (blue), the full SODA sample (gray), and three literature regions: Lupus (tan), Cha I (pink), and Upper Sco (purple) (Ansdell et al. 2016; Pascucci et al. 2016; Barenfeld et al. 2016).

5.2 Disk masses: SODA and other star-forming regions

As the previous discussion has shown, the stellar population of the parts of Orion A probed by our survey – from L1647 to NGC1980 – does not have a single age distribution. Surface over-densities reveal some older or younger populations, and there is a trend toward higher ages with decreasing galactic longitude in the cloud as a whole. In each case, this is reflected in the disk mass distribution: older populations have lower median disk masses. In contrast, in those regions where stellar ages appear to be most homogeneous, and in particular in the southern half of L1641S, we do not see significant changes in disk masses. This implies that disk masses in a region like L1647 only depend on the age, and perhaps the global cloud properties on scales of at least several tens of parsecs.

The logical extension of this idea is to compare the masses of disks with similar ages between different SFRs, which we do here, to see if disk properties are the same at similar ages even in quite different SFRs. This requires the careful selection of a subsample of Class II objects from SODA. Because of the complicated age structures of the stellar population in L1641N, we focus primarily on L1641S and L1647. However, the cluster analysis in Sect. 5.1.2 shows that even in L1641S, outliers in age (and disk mass) are present. We therefore define a sample we call L1641S (reduced). This consists of all YSOs with ℓ > 212.8, and all YSOs with 211.8 ≤ ℓ < 212.3. This means that both cluster 3 and cluster 4 are avoided, and presents a sample with the most uniform age- and disk mass distributions, based on information from the number of Class III and embedded sources in the field. For this sample, we adopt an age of 1.5 ± 0.7 Myr, based on Da Rio et al. (2016). At the same time, we define L1641N to consist of all disks in our sample with ℓ < 211.2. This sample has a higher age, and – potentially – a wider age spread, due to the uncertainties in the stellar population structure here.

5.2.1 Nearby star-forming regions

In Fig. 12, the Kaplan-Meier estimators for the disk mass distribution function in L1641S (reduced), L1641N, and the full SODA sample are shown, as well as the disk mass distributions for three nearby star-forming regions with different ages: Lupus (1–3 Myr, Comerón 2008), Cha I (2–3 Myr, Luhman 2008), and Upper Sco (5–10 Myr, Slesnick et al. 2008). Remarkably, the mass distributions for disks in L1641S (reduced) and Lupus are indistinguishable statistically (p < 0.05). On the other hand, the L1641N disks more closely resemble the older Cha I population, albeit marginally (p = 0.04) in a log-rank test, although we caution that the comparison between these populations may not be entirely proper given the uncertainties in the stellar age spread of the L1641N population, which is likely to be larger than that in Chamaeleon I.

The main conclusion that can be drawn from this comparison is that the disks in SODA, when taking into account the age gradient along the spine of Orion A, have similar masses at similar ages as the protoplanetary disks in entirely unrelated, nearby low-mass SFRs. This in turn implies that their initial conditions and subsequent evolution must have been similar. If nothing else, this suggests that the results obtained from high-resolution observations of the disk continuum in these regions can be generalized to those regions of Orion A where no photoionizing O-type stars are found.

5.2.2 Star-forming regions in Orion A and B

Figure 13 extends this comparison between the SODA sample and Class II disk masses in other nearby star-forming regions, as well as the Class 0 and Class I disk mass distributions in Orion A and B presented by Tobin et al. (2020). The same quantities are also shown in Table 5. Log-normal distributions were fit to all the observed disk populations, and their medians and standard deviations (in log(M_dust/_M⊕) are shown. This time, we plot masses on a logarithmic scale, unlike the linear scales in Figs. 9 and 11, to accommodate the wide range. The left panel focuses on the comparison of the SODA sample to nearby low-mass star-forming regions, where no external photo-evaporation is expected to play an important role. From this panel, we exclude IC348 (Ruíz-Rodríguez et al. 2018), which is dominated by low-mass stars, and thus has lower-than-expected disk masses if not corrected for its bottom-heavy IMF. The right panel of this figure shows the disk masses of a number of previously-studied regions in Orion A and B. Here, only e Ori (Ansdell et al. 2017), which is not well-described by a (single) log-normal distribution due to significant external photoevaporation of disks in the cluster center, is excluded.

Several important conclusions can be drawn from this figure. First, it is clear that the median disk mass of the L1641S (reduced) sample is indistinguishable from that in OMC-2, almost on the opposite side of Orion A from the Trapezium. It also closely resembles not only Lupus (Ansdell et al. 2016), but also Taurus (Andrews et al. 2013), both of which have a similar age (Comerón 2008; Luhman 2004) and similar median disk masses.

Comparison between the median disk masses in SODA and other parts of Orion is complicated by the wide range of environmental conditions in these clouds. While the OMC-2 is similar to the sample studied here, and not affected by O-type stars, disks in the ONC are subject to significant amounts of external photoevaporation due to the proximity of the Trapezium stars. The same is true for the NGC2024 W disks, and recently Haworth et al. (2021) also showed that proplyds are present in NGC2024 E – the youngest sample of disks in this figure. λ Ori, meanwhile, is known to be an older region overall (~5 Myr, Dolan & Mathieu 2001; Mathieu 2008), and indeed less massive than any other population plotted here, regardless of the presence of external photoevaporation.

Fig. 13 Properties of the disk dust mass distributions for disks in nearby star-forming regions and Orion. Median masses (top row) in M⊕ and widths of the log-normal distributions (bottom row), in log(Mdust/M⊕). Left: nearby star-forming regions, sorted by median mass: Lupus (Ansdell et al. 2016), Taurus (Andrews et al. 2013), L1641S (reduced), L1641N, Cha I (Pascucci et al. 2016), ρ Oph (Williams et al. 2019), CrA (Cazzoletti et al. 2019), and Upper Sco (Barenfeld et al. 2016). Right: different populations in Orion A and B: Class 0 and Class I (Tobin et al. 2020) in the full sample and L1641 (in red), L1641S (reduced) and L1641N, OMC-2 (van Terwisga et al. 2019), the ONC (Eisner et al. 2018), NGC2024 E and W (van Terwisga et al. 2020), and λ Ori (Ansdell et al. 2020).

5.2.3 Mass evolution from Class I to Class II in Orion

By combining SODA with the results from the VANDAM survey of Orion by Tobin et al. (2020), we have the most precise picture yet of how protoplanetary disk masses evolve from the earliest embedded phases to Class II (Fig. 13). As these authors also showed, the most embedded Class 0 objects may be lower in mass in L1641 (their sample does not cover L1647). On the other hand, Class I disks do not show significant differences between L1641 and the full sample of such objects in Orion A and B. This suggests that YSOs in L1641 are representative of (nonphotoevaporating) environments in Orion also at younger ages. Thus, we can use these sample to study the mass evolution of disks in a general sense.

From Class I to Class II, the disk mass is reduced by a factor of 13 ± 3, based on these observations. If we assume the evolution from Class I to Class II takes about 1 Myr and continue to use a gas-to-dust mass ratio of 100, that implies an average disk mass loss rate of $\dot{M}~{10}^{-8}{M}_{\odot }{\text{yr}}^{-\text{1}}$ during this period. This number represents the net outcome of a number of different processes: the growth of dust to pebbles and potentially larger bodies which are not detected at these wavelengths, the accretion of material onto the star and its removal through outflows, and the accretion of material from the envelope toward the disk.

Furlan et al. (2016) find values of 10⁻⁶ −5 × 10⁻⁷ M_⊙ yr⁻¹ for the infall rate of new material from the envelope onto the disk and star in the Class I and flat-spectrum sources from SED modeling. However, as these authors suggest, these models may not realistically reproduce the envelope structure, and the stellar masses were not known in these models but are a key parameter for estimates of the infall rate. Thus, they may significantly overestimate the amount of material accreted onto the disk. Indeed, direct observations of inverse P-Cygni profiles associated with infall from envelope species like H₂O and HCO⁺ are remarkably rare in Class I sources (Kristensen et al. 2012; Mottram et al. 2017; van Dishoeck et al. 2021). There is evidence for significant grain growth even in the most embedded phases of protostellar evolution (Harsono et al. 2018), making it difficult to quantify the efficiency of this process from Class I to II. Tychoniec et al. (2020) suggested that planet formation occurring as early as Class 0 may help to explain the observed masses of exoplanetary systems; this would also lead to lower millimeter-dust emission at a rate of approximately 10⁻⁹ M_⊙ yr⁻¹. Interestingly, the rate of mass loss from the disk from Class I to II is of the same order of magnitude as the measured values of the stellar accretion rate for Class II objects (Alcalá et al. 2017) (which in turn are approximately a factor of 10 larger than the mass loss rates via jets; Cabrit 2007; Frank et al. 2014). Combining this observation with the statistically indistinguishable masses of Class I and flat-spectrum disks, we might infer that the observed disk mass at these ages is in (near) steady-state, as also suggested by Tsukamoto et al. (2017). In that case, our closed-box estimate of mass loss describes the evolution from the end of significant envelope accretion, as early as Class I, to the ~Myr-old populations in SODA. However, more detailed studies – in particular observations at higher resolutions of substructures in the youngest disks, and of the interaction between envelope material and the disk – are clearly needed to fully disentangle these competing processes in the (partially) embedded phases and beyond.

5.2.4 Existence of a disk mass – age relation in isolation

The similarity in median disk masses across several distinct populations of Class II YSOs is intriguing. However, it is not without exceptions. Leaving aside IC 342, the two regions that have been suggested in the literature to have low disk masses for their age are CrA and ρ Oph (Williams et al. 2019; Cazzoletti et al. 2019). For CrA, multiple authors have found young ages (in the 1–2 Myr range, Sicilia-Aguilar et al. 2011; Cazzoletti et al. 2019). However, evidence from Gaia now suggests ages of 5–6 Myr for at least a significant fraction of its population (Galli et al. 2020). This suggests that this region's disks should more closely resemble those in Upper Sco in terms of mass, but it is not clear if this older population is large enough to fully explain this region's disk masses.

In ρ Oph, there is similar evidence from Gaia that its YSO population is contaminated to at least some extent by an older, kinematically distinct population, with an age of 5–10 Myr (Grasser et al. 2021). However, again, it is not clear that this population is large enough to fully explain the low observed disk masses. From comparison with L1641N, we can see that it is at least plausible that such a combined population would have median disk masses in the range observed. If this is not the case, however, it suggests that some large-scale property of the ρ Oph SFR must be responsible for the low masses of disks there.

Summarizing, the majority of evidence indicates that disk masses in a given SFR are the same at similar ages (in the absence of external photo-evaporation). Outliers like p Oph and CrA seem to be rare, both in absolute numbers of disks and in terms of different nearby SFRs. If this is the case, it should be possible to use the nearby SFRs as a representative testing ground for theories of disk formation and evolution.

At the same time, this work places strong constraints on the area over which the properties of disks can vary. Compared to previous surveys, the results presented here show for the first time that even within the tens of parsecs spanned by Orion A, no large variations in median disk masses occur, especially after considering the impact of the overall variations in the ages of the sample. This means that if the disks in the outlier regions such as p Oph, CrA, and perhaps L1641N are indeed intrinsically less massive, we must look for the cause at the level of giant molecular clouds. Large-scale cosmic ray ionization trends have been proposed (Kuffmeier et al. 2020), while cloud-scale disk evolution models seem to suggest that turbulence is more important than magnetic fields on the largest scales (e.g., Kuffmeier et al. 2017; Bate 2018).

6 Conclusions

We have conducted an unbiased survey at 225 GHz of 873 Class II disks identified by Spitzer in the southern parts of Orion A, the largest such survey as of yet, by a factor of a few compared to previous work. These observations make it possible to characterize the disk mass distribution of the entire southern part of Orion A as a whole, and enable us to determine the disk mass distribution for clusters of YSOs in the cloud. SODA has unlocked a large, well-characterized sample of protoplanetary disks, which closely resembles smaller nearby star-forming regions. The power of this sample size allows us to distinguish the subtle effects that govern the evolution of protoplanetary disk material at the population level.

Using these data, we have studied the evolution of disk masses, both within the cloud, and in nearby low-mass star-forming regions. We also provide, for the first time, a complete catalog of the most luminous protoplanetary disks at submillimeter wavelengths in L1641 and L1647. We identified an object that may be actively accreting new circumstellar material in sufficient quantities for the accretion stream to be detectable in the 225 GHz continuum, and determine that such objects are rare (~0.1% of all sources); alternatively, it may be a more embedded object seen (nearly) face-on.

Our main conclusions are as follows:

1. Of the 873 Class II disks that make up the SODA sample, we detect 502, or 58% of the full sample. This detection rate lies in between those for Chamaeleon I and Lupus at the same sensitivity;

2. We identify a sample of 20 disks with masses larger than 100 M_⊕, more than doubling the number of such disks in L1641 alone. This is the largest unbiased sample of massive protoplanetary disks in a single cloud to date, and ideally suited to high-resolution follow-up observations;

3. By fitting a log-normal distribution to the observations, we find median disk dust masses of ${2.2}_{-0.2}^{+0.2}{M}_{\odot }$ (full sample) ${2.1}_{-0.2}^{+0.2}{M}_{\odot }$ (L1641) and ${2.6}_{-0.5}^{+0.4}{M}_{\odot }$ (L1647) respectively;

4. Toward the northern part of L1641, the disk masses apparently decrease, even as the fraction of Class III sources increases and embedded sources become rarer. This is consistent with earlier results showing an older foreground population toward this region;

5. We identify six clusters of YSOs with a sufficiently large number of cluster members that their disk properties can be compared statistically. The average disk masses of the YSO clusters we identify in Orion A show differences that are driven primarily by their age;

6. In L1641 S, where we can define a large sample of YSOs with similar ages, we find that disk masses are indistinguishable from those in nearby low-mass star-forming regions in Lupus and Taurus, which have similar ages;

7. Together, these results suggest that disk masses are the same at similar ages (in the absence of external photoevaporation) for many, if not all, nearby SFRs, and that they are similar – and evolve similarly – at scales of ~100 pc in Orion.

The remarkably homogeneous properties of disk samples of the same age are a surprising finding, already hinted at by previous surveys of nearby star-forming regions, which found that disk mass distributions in the Lupus and Taurus star-forming regions are statistically indistinguishable. Now, however, we can show this applies to a larger number of YSOs and YSO clusters, forming in well-separated parts of the same giant cloud. For the first time, the unprecedented size of the SODA disk sample allows us to zoom in on the effects of age gradients and clustering in a single star-forming region. By combining fast ALMA mapping techniques and parallel data processing developed for this study, we open a new window on the evolution of protoplanetary disks.

Acknowledgements

We thank Sierra Grant for insightful discussions on the disk populations in Orion. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2019.1.01813.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 851435). This work made use of the Dutch national e-infrastructure with the support of the SURF Cooperative using using grant no. EINF-259. The authors acknowledge assistance from Allegro, the European ALMA Regional Center node in the Netherlands.

Appendix A Potentially resolved sources

Fig. A.1 Sources with evidence of resolved emission from double Gaussian fitting. Data, model, and residuals are shown per source from left to right. Contours indicate the +4σ (solid black) and -4σ (dashed white) contours. The beam is shown in white in the middle image.

Fig. A.2 Sources with evidence of resolved emission from double Gaussian fitting, as above.

Fig. A.3 Sources with evidence of resolved emission from double Gaussian fitting, as above.

Fig. A.4 Sources with evidence of resolved emission from double Gaussian fitting, as above.

Appendix B The impact of a mixed population of YSOs on the median disk mass

Fig. B.1 Properties of a mixture of two populations of Class II objects with different ages. Left: Median disk masses of distributions sampling N = 105 disks with replacement from ℓ > 213.1° and Upper Sco (Barenfeld et al. 2016) with ratio f (black), and the median disk masses and uncertainties for Upper Sco (red), SODA between ℓ < 210° (blue), and between ℓ > 213.1° (green). Right: Kaplan-Meier estimators for the disk mass distributions resulting from a draw of 500 samples with f = 0.77 (dark gray), in SODA between ℓ > 213.1° (green), ℓ < 210° (blue), and in Upper Sco (red).

References

All Tables

All Figures

Fig. 1 Sample positions on the sky. Disks are marked with white plus signs; disks with masses Mdust > 100 M⊕ are shown as blue circles. The major regions within Orion A (L1647, the three subdivisions of l1641, and the NGC 1980 and ONC clusters) are marked separately. Herschel SPIRE observations at 250, 300, and 500 μm (Soler 2019) form the background. The dashed white circle marks a 4 pc radius around θ1 Ori C.

In the text

	Fig. 2 Data processing workflow for SODA. (a) Standard data reduction. (b) Parallel ALMA data reduction scheme implemented on the SPIDER platform that is part of SURFsara’s Data Processing facilities.
In the text

Fig. 3 Images of selected millimeter-detected sources. For each source, the top panel shows near-infrared (NIR) data (J, H, and Ks-band) form the VISION survey (Meingast et al. 2016) in colors, and 5, 50, and 100σ contour levels. The bottom panel shows the ALMA image at 225 GHz only. White crosshairs mark the positions of detected sources. The top panels show sources with fluxes >4.8 mJy, bottom panels sources with fluxes < 1.6 mJy. The beam is shown in white in the ALMA images, in the lower left.

In the text

Fig. 4 Images of selected (candidate) multiple systems. For each source, the top panel shows NIR data (J, H, and Ks-band) form the VISION survey (Meingast et al. 2016) in colors, and 5, 50, and 100σ contour levels. The bottom panel shows the ALMA image at 225 GHz only. White crosshairs mark the source positions of detected sources. Nondetections are marked with grayed-out names. The beam is shown in white in the ALMA images, in the lower left.

In the text

Fig. 5 Images of YSOs associated with extended millimeter emission. For each source, the top panel shows NIR data (J, H, and Ks-band) form the VISION survey (Meingast et al. 2016) in colors, and 3, 5, 50, and 100σ contour levels. The bottom panel shows the ALMA image at 225 GHz only. White crosshairs mark the source positions of detected sources. The beam is shown in white in the ALMA images, in the lower left.

In the text

Fig. 6 Disk mass distributions inferred with the Kaplan-Meier estimator in Orion A and several nearby star-forming regions. SODA (dark blue) versus Class I and Class 0 disks in Orion (Tobin et al. 2020) (orange and green), and disks in Lupus (Ansdell et al. 2016) (purple), Chamaeleon I (Pascucci et al. 2016) (brown), ρ Oph (Williams et al. 2019) (gray), and Upper Sco (Barenfeld et al. 2016) (yellow).

In the text

	Fig. 7 Quantile-quantile plots of the empirical distribution (inferred with a Kaplan-Meier estimator) and fitted quantiles using a Weibull distribution (top panel) and log-normal distribution (bottom panel).
In the text

	Fig. 8 Comparison of SODA disk masses to previous literature. Disk mass (or luminosity) distributions in L1641 (blue) and L1647 (green) from SODA, as well as the observed disk mass distribution in the sample from Grant et al. (2021) (pink) and a scenario in which all of the unobserved sources in Grant et al. (2021) are also nondetections at their sensitivity (gray).
In the text

Fig. 9 Trends in disk properties in SODA as a function of galactic longitude. Left: median Class II disk masses (black points) for all longitudinal bins. The location of the Trapezium is marked with a blue star. Right: ratios of embedded (red) and Class III YSOs (purple) to Class II sources in each bin. Round symbols show source fractions derived from Pillitteri et al. (2013) and Tobin et al. (2020) for Class III and embedded sources, respectively; square symbols show the embedded source fraction in L1647, derived from Megeath et al. (2012).

In the text

	Fig. 10 Clusters identified by OPTICS. Cluster members are indicated by colored circles, while field stars are shown as gray plus signs.
In the text

Fig. 11 Median disk dust masses in clusters and the field population. Left: median masses for the individual clusters and field population (black) and the median and standard deviation for the clusters together. Right: ratio of the numbers of embedded (red) and Class III YSOs (purple) to Class II sources in each cluster. Round symbols show source fractions derived from Pillitteri et al. (2013) and Tobin et al. (2020) for Class III and embedded sources, respectively; square symbols show the embedded source fraction in L1647, derived from Megeath et al. (2012).

In the text

	Fig. 12 Comparison of the SODA results to three nearby SFRs. Kaplan–Meier estimators for the cumulative disk mass distribution in L1641S (reduced) (orange), L1641N (blue), the full SODA sample (gray), and three literature regions: Lupus (tan), Cha I (pink), and Upper Sco (purple) (Ansdell et al. 2016; Pascucci et al. 2016; Barenfeld et al. 2016).
In the text

Fig. 13 Properties of the disk dust mass distributions for disks in nearby star-forming regions and Orion. Median masses (top row) in M⊕ and widths of the log-normal distributions (bottom row), in log(Mdust/M⊕). Left: nearby star-forming regions, sorted by median mass: Lupus (Ansdell et al. 2016), Taurus (Andrews et al. 2013), L1641S (reduced), L1641N, Cha I (Pascucci et al. 2016), ρ Oph (Williams et al. 2019), CrA (Cazzoletti et al. 2019), and Upper Sco (Barenfeld et al. 2016). Right: different populations in Orion A and B: Class 0 and Class I (Tobin et al. 2020) in the full sample and L1641 (in red), L1641S (reduced) and L1641N, OMC-2 (van Terwisga et al. 2019), the ONC (Eisner et al. 2018), NGC2024 E and W (van Terwisga et al. 2020), and λ Ori (Ansdell et al. 2020).

In the text

	Fig. A.1 Sources with evidence of resolved emission from double Gaussian fitting. Data, model, and residuals are shown per source from left to right. Contours indicate the +4σ (solid black) and -4σ (dashed white) contours. The beam is shown in white in the middle image.
In the text

	Fig. A.2 Sources with evidence of resolved emission from double Gaussian fitting, as above.
In the text

	Fig. A.3 Sources with evidence of resolved emission from double Gaussian fitting, as above.
In the text

	Fig. A.4 Sources with evidence of resolved emission from double Gaussian fitting, as above.
In the text

Fig. B.1 Properties of a mixture of two populations of Class II objects with different ages. Left: Median disk masses of distributions sampling N = 105 disks with replacement from ℓ > 213.1° and Upper Sco (Barenfeld et al. 2016) with ratio f (black), and the median disk masses and uncertainties for Upper Sco (red), SODA between ℓ < 210° (blue), and between ℓ > 213.1° (green). Right: Kaplan-Meier estimators for the disk mass distributions resulting from a draw of 500 samples with f = 0.77 (dark gray), in SODA between ℓ > 213.1° (green), ℓ < 210° (blue), and in Upper Sco (red).

In the text