Early Planet Formation in Embedded Disks (eDisk). I. Overview of the Program and First Results

One approach to distinguish embedded disks and their envelopes is to utilize their kinematic difference; embedded disks are expected to have Keplerian rotation similarly to Class II disks, while envelopes are expected to have rotation conserving their specific angular momenta. As a result, rotational velocities of embedded disks are proportional to r^−0.5, while those of envelopes are proportional to r⁻¹, where r is the cylindrical radius from the rotational axis (e.g., Ohashi et al. 1997a; Simon et al. 2000; Belloche 2013).

This approach was taken to distinguish the embedded disk and its surrounding envelope around the protostar HH 111 (Lee 2011) and also those around six other protostars (Yen et al. 2013) in C¹⁸O 2–1 with the Submillimeter Array (SMA) at angular resolution of ∼1''–4''. In some cases (L1527 IRS and L1448-mm), only rotation proportional to r⁻¹ was found because of insufficient angular resolution. However, L1527 IRS was also observed with the Combined Array for Millimeter-wave Astronomy (CARMA) at an angular resolution of ∼1'' in ¹³CO 2–1, displaying Keplerian rotation instead of that proportional to r⁻¹ (Tobin et al. 2012).

Further observations following up these pioneering works were enabled by ALMA with its high angular resolution and sensitivity. L1527 IRS was observed again as a benchmark at subarcsecond resolutions in C¹⁸O 2–1, clearly demonstrating that its infalling envelope transforms into a Keperian disk at a radius of ∼70 au (Ohashi et al. 2014; Aso et al. 2017). Similar findings have been found in other low-mass Class 0/I protostars (e.g., Yen et al. 2014; Aso et al. 2015; Tobin et al.2016; Yen et al. 2017; Harsono et al. 2018; Maret et al. 2020; Sai et al. 2020). C¹⁸O has often been adopted to investigate kinematics of envelopes and disks, as well as other rarer isotopes, such as C¹⁷O or C³⁴S (e.g., Artur de la Villarmois et al. 2019). Important products from detected Keplerian motions are dynamical masses of protostars and mass accretion rates. Notably, direct measurements of masses for optically invisible protostars are usually difficult because their stellar photospheres are not observable. These physical parameters could shed light on the formation and evolution of protostars (e.g., Yen et al. 2017; Artur de la Villarmois et al. 2019) based on more systematic observations of embedded disks.

The presence of disks around deeply embedded protostars may also affect their chemical structures and, consequently, the molecular line emission observed toward them. For example, such embedded disks have been suggested to be associated with distinct features such as slow accretion shocks that in turn may reveal the presence of the disks (e.g., Sakai et al. 2014; Oya et al. 2016; Sakai et al. 2016; Okoda et al. 2018). In addition, grain growth to millimeter- or centimeter-sized grains in disks may cause the continuum emission to become optically thick, affecting the line emission signatures observed from the embedded disks (e.g., Harsono et al. 2018).

We should, however, note that consistent chemical indicators across larger samples of disks are still lacking. For example, it has been suggested that SO is enhanced at the disk–envelope surface owing to slow accretion shocks in L1527 IRS (Sakai et al. 2014), while in L483 and Elias 29 SO appears to trace the Keplerian rotating disks themselves (Oya et al. 2017, 2019). Oya et al. (2017) also suggested that the emission from CS and complex organic molecules are arising from the disk around L483, while Jacobsen et al. (2019) demonstrated that those species show velocity structures that can be explained better by infalling and rotating material rather than a Keplerian disk. These observations indicate that the molecules tracing the disks may vary among objects, depending on their physical properties (e.g., stellar luminosities) and chemical history in prestellar stage (e.g., Oya et al. 2019).

Although embedded disks have been successfully observed around a number of protostars in the past decade, only a few have shown detectable substructures suggestive of ongoing planet formation to date because higher angular resolutions of 01 or less are required to look for substructures in embedded disks. Sheehan & Eisner (2017) have found a central cavity and ring toward the Class I protostar WL 17, in 3 mm continuum at an angular resolution of 006. Sheehan & Eisner (2018) have also found concentric rings/gaps within the disk of the Class I protostar GY 91, also known as ISO-Oph 54 (see Cieza et al. 2021), in the 0.9 mm continuum emission at an angular resolution of 013. More recently, Sheehan et al. (2020) presented seven protostars, including two Class 0 sources from the Orion Very Large Array/ALMA Nascent Disk and Multiplicity (VANDAM) survey that exhibit ring/gap substructures in their disks, although detailed modeling showed an offset in the inner disk for a few sources, which may suggest that the substructures could be due to close binary formation. In addition, Segura-Cox et al. (2020) have found four annular structures in the Class I protostar Oph IRS 63 in 1.3 mm continuum at an angular resolution of ∼004, although their contrast is relatively low.

In addition to ring structures, spiral structures have also been observed in embedded disks. Lee et al. (2020) have found a pair of symmetric spiral structures in the disk around the Class I protostar HH 111 VLA1 at a resolution of 004 (16 au). The disk is relatively massive (0.33–0.5 M_☉), with Toomre's Q parameter near unity in the outer part of the disk, where the spiral structures are detected, supporting the idea that mass accretion from the envelope is driving the outer disk to be gravitationally unstable. Obviously, it is definitely crucial to observe more embedded sources at angular resolutions higher than 01 to search for substructures, particularly around Class 0 protostars, which show less substructures in their disks to date.

For the eDisk program we focus on sources located in nearby (d < 200 pc) star-forming regions that are relatively bright (L_bol > 0.1 L_⊙). The first criterion is required to spatially resolve observed disks at a given angular resolution (∼004; see below). The second criterion is required to detect disks around sample protostars within a reasonable observation time. We checked how many such protostars there are based on available catalogs of YSOs: the largest catalogs of YSOs in nearby star-forming regions were provided through large survey observations at mid- and far-infrared wavelengths using the Spitzer Space Telescope. Rebull et al. (2010) counted the number of Class 0/I protostars in Taurus based on the near- to mid-IR slopes, α, obtained in the Spitzer observations, finding that there are 21 Class 0/I protostars with L_bol > 0.1 L_⊙, ²¹ excluding potential Taurus members newly identified. Dunham et al. (2015) also counted the number of Class 0/I protostars based on Spitzer observations and found 37 Class 0/I protostars in the regions of Chamaeleon, Corona Australis, Lupus, Musca, and Ophiuchus. Note that 2 of the 37 sources show unreasonably high T_bol (>3000 K) for embedded protostars. Based on these previous works, there would be ∼60 protostars within a distance of 200 pc and with luminosities more than 0.1 L_⊙. Note that protostars in isolated star-forming regions are not included here.

Within a reasonable observation time, it is impossible for us to observe all these protostars from a data volume point of view, as well as the restrictions of the maximum time allocated to a large program in a given LST range. In the Cycle 7 semester 645 hr were available in total for Large Programs, of which we aimed for a required observation time for the eDisk program of ∼100 hr. According to our previous ALMA observations of embedded disks around protostars and also our simulations of observations, we set a target sensitivity for molecular line observations at ∼2.7 mJy beam⁻¹ with a velocity resolution of ∼0.17 km s⁻¹. This sensitivity was able to be achieved with an on-source integration time of ∼2.4 hr, requiring ∼6 hr of observation time, including overhead. Note that this observation time can provide us with a sensitivity of ∼13 μJy beam⁻¹ for the continuum emission, which is sufficiently high to detect continuum emission around protostars. With these estimations of the required observation time, we found that ∼17 sources could be observed within ∼100 hr of observation time.

When 17 sources were selected, we tried to avoid any bias as much as possible; we selected 17 "representative" protostars, which have been relatively well studied before, are located in different star-forming regions (i.e., Chamaeleon, Corona Australis, Lupus, Ophiuchus, Taurus, and also isolated Bok globules), and present a wide distribution of the bolometric luminosities (L_bol) and temperatures (T_bol). Selecting sources spreading out over various star-forming regions can also help to carry out our observations more efficiently within a limited period of the observations. In addition to these 17 sources, two sources, B335 and TMC-1A, for which necessary data have already been obtained with ALMA, were added to the sample. The 19 sample protostellar systems selected for the eDisk program are listed in Table 1.

All of the selected sources have complete spectral energy distributions (SEDs) measured from the near-infrared to submillimeter that are compiled in Appendix A. Based on those SEDs, we rederived their bolometric temperatures (T_bol) and luminosities (L_bol) also utilizing the most recent distances from Gaia measurements listed in Table 1. Note that none of the 19 targets were detected with Gaia directly owing to their embedded natures, but distances to nearby sources, presumably related to the same star-forming complexes, were typically estimated based on Gaia observations and provided the distance estimates quoted in the literature. The final sample composes 12 Class 0 and 7 Class I protostars and represents various evolutionary stages of protostars with a roughly uniform distribution of L_bol and T_bol (Figure 1(a)). For comparison purposes, Figure 1(b) shows the same T_bol–L_bol diagram, but for the 56 protostars with d < 200 pc and L_bol > 0.1 L_⊙ discussed above, except for two sources with unreasonably high T_bol, indicating that the final 19 protostars reasonably represent the nearby and relatively bright protostars.

Zoom In Zoom Out Reset image size

Observations of the eDisk program using the ALMA 12 m array in band 6 were conducted from 2021 April to 2022 July. We aimed to achieve high angular resolutions of ∼004, corresponding to ∼6 au at a distance of 140 pc. Such high angular resolution is necessary for us to investigate substructures of embedded disks and make direct comparisons with previous maps of Class II disks obtained at similar angular resolutions, such as those obtained by DSHARP or ODISEA.

The correlator was configured to observe both 1.3 mm continuum emission and CO (J = 2 − 1) isotopologues and other lines providing additional information about the physical and chemical structures of the gas on the scales of the circumstellar disks surrounding the embedded protostars. Baseband 1 of the correlator was configured in the lower sideband with four spectral windows (spw) to include our primary molecular lines, C¹⁸O (2–1), ¹³CO (2–1), and SO (6₅–5₄), at a velocity resolution of ∼0.17 km s⁻¹. The ¹²CO (2–1) emission was included in baseband 4 in the opposite sideband to study protostellar outflows. The bandwidth of baseband 4 was set to be 937.5 MHz to secure the continuum bandwidth, and as a result, the velocity resolution of the ¹²CO data is 4 times worse than that of the C¹⁸O data. The remaining two basebands are located in the upper and lower sidebands, respectively, with the maximum bandwidth of 1875 MHz. While these two basebands were adopted for the 1.3 mm continuum emission, the central frequency of baseband 2 was adjusted to incorporate as many important molecular lines as possible, such as CH₃OH, SiO, DCN, H₂CO, and c-C₃H₂. Even though the velocity resolution of these observations, ∼1.3 km s⁻¹, is not ideal to study the kinematics from these molecular transitions, their detection can still provide useful constraints on circumstellar structures based on, e.g., the integrated intensities. All the major emission lines observed in our observations are listed in Table 2.

In order to achieve an angular resolution of ∼004, the C43-8 antenna configuration was adopted for our observations, providing long baselines up to ∼12.6 km. In addition to the C43-8 antenna configuration, the C43-5 antenna configuration was also adopted to cover shorter baselines with a minimum of ∼15 m. This additional antenna configuration covering shorter baselines is crucial for us to investigate relatively extended structures up to ∼2''–3'' in the envelopes surrounding our sampled protostars. We should note that a large fraction of the emission from the envelopes around the sampled protostars is still filtered out in our observations. However, this missing flux has less impact on the main objectives of the eDisk program focusing on dust and gas on the scales of the disks in the vicinity of the protostars. The field of view of the 12 m array at the observed frequency is ∼26''. Logs of the observations are provided in Appendix B (see Table 5).

For two sources in our sample, TMC-1A and B335, comparable observations were obtained with ALMA prior to our eDisk program, and we collected those data from the ALMA archive. In addition, due to scheduling constraints, only short-baseline data were obtained for Oph IRS 63 as a part of eDisk. Hence, we relied on ALMA archival observations of Oph IRS 63 for long-baseline data. The information on the ALMA archival data we use in the eDisk program is summarized in Appendix B (see Table 6).

The eDisk data were all calibrated using the ALMA calibration pipeline packaged with the Common Astronomy Software Application (CASA; McMullin et al. 2007) version 6.2.1 and pipeline version 2021.2.0.128. After delivery, the data were restored to a calibrated state running scriptForPI.py in the same CASA version that the calibration pipeline ran with. For four sources (BHR 71 IRS1/IRS2, IRAS 04166+2706, and IRAS 04169+2702), there were execution blocks (EBs) that were given a QA0 status of "semi-pass" (see Table 5 in Appendix B), which in the case of these data sets was because of excessive phase decorrelation during the observations. For those we downloaded the raw data and ran the ALMA calibration pipeline manually to inspect the data sets and to determine whether they could be salvaged using self-calibration. The QA0 semi-pass data that are used in the final images are denoted in the observing logs included in Appendix B.

Outside of standard calibration, we needed to perform additional processing to improve the signal-to-noise ratio (S/N) of the continuum data through self-calibration, to remove the continuum emission from the line data, and to combine data taken in two configurations in order to achieve the balance between flux recovery and angular resolution needed to achieve the science goals of the project. Nearly all of the eDisk targets were sufficiently bright that images created using the standard ALMA calibration were dynamic range limited, meaning that the images were not reaching the thermal noise limit and were instead limited by phase and/or amplitude variations during the science target observations. To both combine the two configurations and self-calibrate, continuum subtract, and image the data, we constructed a set of semi-interactive scripts for each target. Our methodology was based on the DSHARP data reduction procedure (Andrews et al. 2018) with substantial modification. One particular choice we made was to operate on each measurement set individually for calculation and application of phase and amplitude gain solutions and not concatenate our data to a single measurement set. This avoids issues that can be encountered when operating on concatenated measurement sets within various CASA tasks, for example, interpolation of gain solutions between distinct observations that might be far apart in time, spw numbers being offset, and possible issues with combined metadata, to name a few. Many of the issues with concatenated measurement sets can be worked around, but we avoid the known issues and the possible unknown issues by operating on individual measurement sets. The scripts used for each source are provided at http://github.com/jjtobin/edisk.

The reduction procedure starts with the construction of a spectrally averaged continuum data set. We identified line-free regions of the individual spw's to include in the continuum data set in two ways: First, we made use of the cont.dat file produced by the ALMA imaging pipeline, and edits to the continuum ranges in that file were made if necessary, based on inspection of the hif_findcont weblog. Second, we manually flagged the spectral lines that are expected to be present in all data sets regardless of whether they were already included in the cont.dat file. We then translated these frequency ranges defined in the kinematic local standard of rest frame (LSRK) to the channel ranges for each EB and flagged those channels. Then, the separate spw's from each EB were spectrally averaged individually. The four high-resolution spw's were averaged over 480 channels, and the three low-resolution spw's were averaged over 60 channels, resulting in ∼30 MHz channels, approximating the low spectral resolution (time domain mode (TDM)) of the ALMA correlator.

We next checked the relative flux density scaling of each individual EB using their azimuthally averaged visibility amplitude data. To do this, we first imaged each EB individually and shifted the brightest source in the image to the phase center using the CASA task fixvis. We then computed the azimuthally averaged visibility amplitudes for each EB, plotted them together, and computed the average scale factors between each EB. The scale factors were estimated from the minimum overlapping uv-distance out to 800 kλ; this outer uv-distance was chosen to avoid a possibility that decorrelation at longer baselines could significantly impact the scale factors. The offsets were typically at the 10%–15% level between EBs, and we scaled each EB to the "best" EB, which was judged by having the least amount of scatter in its data. Thus, our goal was a relative flux normalization of all the EBs, given that we were not improving on the overall flux density scale accuracy. We also note that there is a possibility that some scale factors include intrinsic source variability in addition to calibration differences. While the ALMA Technical Handbook indicates that the accuracy of the flux density scale is expected to be ∼5%, the scaling that we have had to apply to data sets indicates that the flux density scale can have uncertainties up to the rescaling factors that we apply to the data. ²² The scale factor applied to each EB is provided in Appendix C.

There were some cases where this flux normalization could not be performed at this stage owing to decorrelation, which caused the visibility amplitude profile of one or more EBs to have a different behavior as a function of uv-distance than others. In those cases, we did not determine scale factors at this stage, but we completed the self-calibration steps that followed and determined the scale factors after self-calibration. We then reran the whole procedure from the start, applying these scale factors a priori; we referred to this as the "two-pass" method. The "two-pass" method was applied to Oph IRS 43, IRAS 04166+2706, IRAS 04169+2702, and BHR 71 IRS1 and IRS2. The relative scaling of each EB was then applied to the continuum data sets that did not have the brightest source shifted to the phase center because the unshifted data were what we ultimately used in the self-calibration process.

Following flux normalization of each EB, we began the self-calibration process. For the combined processing of the compact and extended configuration data, we first self-calibrated the compact configuration data, and then we self-calibrated the pre-self-calibrated compact configuration data with the extended configuration data. Self-calibration was performed using the data that did not have the phase center shifted. It was not necessary to align our data prior to self-calibration because the time difference between observations was less than 1 yr and all the long-baseline data, where proper motion would have the greatest effect, were taken within a 3-month time period. We create our model for self-calibration using the data from all the measurement sets; however, it is only the calculation and application of the phase and amplitude gain solutions that happen on a per-measurement set basis.

To prepare for self-calibration, we first calculated a ranked list of reference antennas for each EB, equally weighing the amount of flagging on a particular antenna and its geometric position within the array. Following this step, we made an initial image using all the continuum data sets to determine the peak S/N of the image from standard calibration; the creation of all images during the self-calibration process made use of the tclean auto-multithresh algorithm (Kepley et al. 2020) to automatically create clean masks for deconvolution. Following the generation of this initial image, we analyzed the scan lengths within an observation and determined calibration solution intervals that most evenly divide a scan by the number of integrations per scan. This approach will avoid the uneven division of scans that can result in some solutions being calculated from significantly less data. The list of solution intervals to attempt always includes a solution that spans the length of an EB, referred to as "inf_EB," where gaincal is run with the parameter combine="scan,spw." The "inf_EB" solution interval corrects systematic effects, in particular, antenna position offsets that translate to systematic phase offsets between the calibrator and science target. Then, for shorter solution intervals, we are correcting the atmospheric phase variations. We include an "inf" solution, where the "scan" setting is removed from the gaincal parameters and solutions break at scan boundaries, and then for shorter solutions, we divide the median scan length by 2, and all subsequent solution intervals divide the previous interval by 2, until the final solution interval of "int" (one integration) is reached.

Then, for each solution that is to be attempted, we also determine the clean depths in terms N × σ to use for each solution interval in order to progressively clean deeper during each round of self-calibration. We begin with the assumption that we would trust features with intensities > (Peak S/N)/15, which means that for a Peak S/N = 150, we would trust features >10σ. Note that these numbers are empirical, but small changes do not appreciably affect the outcome of self-calibration. So we divide the S/N of the initial image by 15 to obtain our starting clean threshold, and we decreased the thresholds for subsequent solution intervals logarithmically until the clean depth is equivalent to 3σ.

With a set of solution intervals to attempt and their accompanying clean thresholds, the self-calibration process begins. We collect the model creation with tclean, calibration solution generation with gaincal, application of the calibration solutions with applycal, and checking the outcome of a self-calibration iteration with tclean into a single Python function to simplify the process. We always use the combine="spw" parameter in gaincal to use all the available bandwidth for self-calibration, and we use the parameter interp="linearPD," which will scale the phase solutions with frequency to compensate for the fact that the optimal phase solutions will be different between the low- and high-frequency ends of the bandwidth.

At the start of each self-calibration iteration, the self-calibrated data from the previous iteration (corrected data column) are split out into a new measurement set such that each self-calibration solution interval is incremental on the previous iteration. We compare the image S/N before self-calibration and after self-calibration of each iteration is applied to ensure that the S/N is increasing. We do this by running tclean again but using the self-calibration model that was created at the start of a self-calibration iteration as the startmodel for tclean. This approach enables us to effectively determine the improvement in S/N just due to the gain corrections derived from that model. When the image S/N increases within an iteration, we consider that iteration successful, if the S/N decreases, that iteration has failed and we consider the last successful iteration as the final self-calibration iteration. The above steps consider phase-only self-calibration, and we do attempt amplitude self-calibration (solving for amplitude and phase simultaneously) following the phase-only self-calibration using the same criteria for success or failure of amplitude self-calibration. Following the completion of self-calibration for the successful solution intervals, the self-calibrated data for each EB are saved to new measurement sets and are considered the final continuum data.

Once self-calibration is completed on the continuum data, we process the line data. First, we apply the flux normalization to each full measurement set using the scaling factors derived from the continuum. Then, we apply the self-calibration solutions to the full measurement sets for each EB. However, the S/N improvements for the line data are typically not significant because the line data are not dynamic range limited like the continuum. For each EB, we apply all the gain tables for each successful self-calibration iteration, which is necessary because each table is incremental, so applying the full correction requires all the tables. Once self-calibration solutions are applied, we subtract continuum emission from the line data using the uvcontsub task. To subtract the continuum data, we select the line-free regions of the data that will have their continuum data fitted and subtracted. To select the frequency ranges to fit, we select the complement of the frequency ranges that were flagged to produce the continuum data set, that is, the line-free regions of the bandwidth. Following uv continuum subtraction, we have the final spectral line data set, and it is ready for imaging.

We first create the continuum images using tclean and the auto-multithresh algorithm to create the clean masks. Then, we create images using Briggs weighting, for the robust parameters of −2, −1, −0.5, 0, 0.5, 1, and 2, in order to compare the features present with different weighting. We also created images using uv-tapering, to reduce the influence of the long baselines and bring out larger-scale structure, and we created tapered images, with the tapers starting at 1000, 2000, and 3000 kλ, and used robust parameters of 1 and 2. The typical synthesized beam size for the robust =0.5 images was ∼005. The images for most sources use 6000 pixels with a pixel size of 0003; however, some sources require image sizes of ∼14,000 in order to include sources toward the edge of the primary beam.

Then, we create the spectral line images for each line listed in Table 2. We typically created images using Briggs weighting with robust values of 0.5 and 2 and uv-tapering starting at 2000 kλ. The resulting images have synthesized beams of ∼01 and ∼015 for the two robust values and applied tapering, respectively. The tapering was essential for the spectral line images because the surface brightness sensitivity was otherwise too low at the native resolution of the data to detect the spectral lines with high S/N. The spectral line images typically used 4000 pixels and a pixel size of 001. The spectral line images were all created with velocity-defined channels around the spectral lines of interest using their associated rest frequencies (Table 2).

We then save the non-primary-beam-corrected images (flat noise), the primary-beam-corrected images, the final clean mask, and the primary beam response for each individual image. As such, there are ∼52 images created per target when all the combinations of robust parameters, lines, and tapering are considered. The images presented in this overview paper and subsequent papers from the eDisk program are selected based on the balance between resolution and sensitivity, as well as the structures discussed.

Recently it has been pointed out that data imaged with the clean algorithm may suffer systematic errors of the flux scale owing to the so-called Jorsater & van Moorsel (JvM) effect (Jorsater & van Moorsel 1995; see, in particular, Czekala et al. 2021). The reason for this is that the dirty beams in the CLEAN processes deviate from the used Gaussian restoring beams. As discussed in Czekala et al. (2021), the JvM effect is more significant when data taken from different antenna configurations are combined like in the case of the eDisk program. The impact due to the JvM effect to the flux scales of our eDisk maps differs from source to source: maps with more extended diffuse emission are more impacted by the JvM effect. However, it has also been pointed out that the correction method to fix the flux scale issue, implemented by Czekala et al. (2021), may artificially manipulate the noise levels of maps after the correction (Casassus & Cárcamo 2022). As the main discussions of the first set of papers from eDisk are based primarily on the morphology of maps and velocity structure of molecular emission rather than the flux scale of maps, we did not apply any correction to our maps in this overview paper and the related first set of papers.

We present representative maps of the continuum images obtained in the eDisk program on the same spatial scale in Figure 2. The same continuum images but enlarged are also presented in Figure 3. Robust parameters used to create these continuum images and their angular resolutions are shown in Table 3. The rms noise level of these continuum maps is ∼25 μJy beam⁻¹ on average. Continuum emission is clearly detected at all source positions, and 4 out of the 19 sources also show continuum emission arising from companions; 3 of them, i.e., Ced 110 IRS4, R CrA IRS 7B, and R CrA IRAS 32, are newly found to be binaries in our observations. These newly found binaries are named based on the IAU convention: Ced 110 IRS4A and Ced 110 IRS4B, R CrA IRS 7B-a and R CrA IRS 7B-b, and R CrA IRAS 32A and R CrA IRAS 32B.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

All the continuum emission images are spatially resolved, exhibiting elongated structures. These elongated structures strongly suggest that the dust continuum traces disks around the targeted sources. Note that for some sources envelope contamination may still play a role for both dust and gas emission; to assess that, more detailed radiative transfer modeling will be needed on a source-by-source basis. The extent of the continuum emission, measured through 2D Gaussian fittings and listed in Table 3, varies significantly from source to source; L1489 IRS has the largest emission with an angular extension of ∼39 × 13, while the companion of Oph IRS 43 (VLA 2) is the most compact, marginally resolved (2σ) with an angular extent ∼0019 × 0016. We note that some of the continuum emissions were not well fitted with a 2D Gaussian, suggesting that the extents of the emissions provided in Table 3 for those sources may not reflect their actual extents. In such a sense, the extents of the continuum emissions provided here should be considered as a quick reference of their extents measured by a simple procedure. More accurate measurements of the extents of the dust emissions through more sophisticated methods will be presented in future papers. A comparison between the extents of the continuum emission from the different sources and their bolometric temperatures is shown in Figure 4. In agreement with the systematic studies of protostars in Orion with the Very Large Array (VLA; Tobin et al. 2020), no clear correlation between size and evolutionary stage is seen—although the two largest disks in our sample are both associated with Class I protostars (L1489 IRS and IRAS 04302+2247). Most of the disks in the eDisk sample are smaller than about 100 au, but we caution that these are estimates of the sizes of the dust disks and that the gas disks may be larger as is often inferred for more evolved disks (see, e.g., Miotello et al. 2022).

Zoom In Zoom Out Reset image size

From a comparison between the major and minor axes of the continuum emission, the inclination of the system can also be estimated on the assumption that the continuum emission is a vertically thin disk. We find that 7 out of 19 sources, excluding companions, have inclination angles of more than 70° (see Table 3), roughly consistent with the expectation from a random distribution. We caution that the estimate of the inclination is a lower limit to the true inclination because the disks have most likely some vertical thickness. Still, these are different from most Class II disks that have been observed at high angular resolution to date, which tend to be more face-on (Andrews et al. 2018; Cieza et al. 2021). Figure 5 shows a comparison between the inclination angle and the bolometric temperature, demonstrating no correlation between them. This suggests that the disk inclination has less impact on our classification of the 19 protostellar systems based on the bolometric temperature even though it was also suggested that higher inclination angles of disks could misclassify protostars as Class 0 sources (Tamura et al. 1996; Nakazato et al. 2003; Tobin et al. 2008; Fischer et al. 2017). We also note that some of the continuum structures show relatively high brightness temperatures of more than 100 K. Such high brightness temperature may suggest that the detected continuum emission is at least partially optically thick. Basic physical parameters of the continuum emission are summarized in Table 3. We note that statistical errors are not given in this table because the actual uncertainties of these measurements are more likely determined by systematic errors, such as the flux calibration error and the gain calibration error.

Zoom In Zoom Out Reset image size

When we carefully inspect the continuum maps, it is found that some of them show ring-like structures that can be visually identified by the naked eye. The most obvious source is L1489 IRS. Oph IRS 63 also weakly shows ring-like structures. More details of the continuum maps obtained in the eDisk program will be presented in subsequent papers. Ring-like structures in continuum emission around these two sources have been suggested by previous studies as well (Segura-Cox et al. 2020; Ohashi et al. 2022), although our map of L1489 IRS seems to show ring-like structures more clearly than in the previous work (Ohashi et al. 2022). The continuum emission around IRAS 04169+2702 seems to show a bean-like structure with brighter emission on its southwest side, which could also be a ring-like structure. We should stress that a more detailed and systematic analysis of continuum emission is definitely required to identify their substructures, including those that cannot be visually identified by the naked eye, as will be discussed in a forthcoming paper. Nevertheless, the first analysis of the observations suggests that continuum emission around our Class 0/I sample does not show very sharp ring or spiral structures that are often seen in Class II continuum disks (e.g., DSHARP, ODISEA). This striking difference will be discussed in more detail in Section 6.1. On the other hand, several disk continuum images show brightness asymmetries, particularly along their minor axes. IRAS 04302+2247, GSS 30 IRS3, IRAS 16253–2429, IRAS 16544–1604, R CrA IRS 7B-a, R CrA IRAS 32A, and R CrA IRAS 32B show clear asymmetry that can be identified by eye, while asymmetry can also be identified with more careful inspection of continuum intensity distributions for some additional sources, such as L1527 IRS. For these sources showing brightness asymmetries, the geometric centers of the continuum emission derived from single 2D Gaussian fittings listed in Table 3 are shifted from their actual continuum peak positions. The possible origin of the brightness asymmetry is discussed in Section 6.1.

The C¹⁸O 2–1 emission is also detected toward all 19 sources in our sample. Here we show the C¹⁸O emission detected in one of our sources, R CrA IRS 7B, as an example of what can be seen toward the eDisk targets.

Lindberg et al. (2014) have carried out ALMA observations of IRS 7B in band 7 at an angular resolution of ∼04 and detected continuum emission and several molecular lines, including C¹⁷O 3–2. Both continuum and C¹⁷O emission have extended structures elongated in the southeast-to-northwest direction. The C¹⁷O emission also shows a clear velocity gradient along the elongation, which was interpreted as Keplerian rotation (Lindberg et al. 2014). Although a systemic velocity of 6.2 km s⁻¹ was adopted for the analysis of the velocity structures of C¹⁷O in Lindberg et al. (2014), we adopt ∼6.0 km s⁻¹ as the systemic velocity of IRS 7B in this paper (see more details below).

C¹⁸O 2–1 was detected at more than 3σ in the LSR velocity range between −3.15 and 14.22 km s⁻¹ at an angular resolution of ∼013 × 011 as shown in Figure 6. Note that the velocity resolution of the C¹⁸O data shown here is ∼0.33 km s⁻¹, which is double the original velocity resolution to achieve a higher sensitivity. Since the systemic velocity of this system is ∼6 km s⁻¹ (see Section 6.2), LSR velocities smaller and larger than 6 km s⁻¹ correspond to blue- and redshifted, respectively. At higher blueshifted velocities (V_LSR = −3.15 to 2.86 km s⁻¹) and redshifted velocities (V_LSR = 8.20−14.22 km s⁻¹), C¹⁸O is mainly detected on the southeast side and the northwest side of IRS 7B-a within ∼05, respectively. At lower blueshifted (V_LSR = 3.19−4.20 km s⁻¹) and redshifted velocities (V_LSR = 7.20−7.87 km s⁻¹), more extended C¹⁸O emission is detected within the field of 2'' × 2'', demonstrating that our observations can also trace extended emission as well. At velocities close to the systemic velocity, absorption of C¹⁸O is detected. In particular, very deep absorption is seen at the position of IRS 7B-a at V_LSR between 5.20 and 5.87 km s⁻¹. The extent of the deep absorption is almost the same as that of the continuum emission associated with IRS 7B-a, suggesting that the deep absorption is due to cold, foreground C¹⁸O emission against the bright continuum emission. Compact and weak C¹⁸O emission is detected at the position of IRS 7B-b at several channels of redshifted velocities, such as V_LSR = 9.21 and 11.21 km s⁻¹, and also absorption at V_LSR = 5.20 and 5.53 km s⁻¹. The absorption at the position of IRS 7B-b may also be due to cold foreground C¹⁸O against the continuum emission associated with IRS 7B-b.

Zoom In Zoom Out Reset image size

To investigate the overall distribution and velocity structures of C¹⁸O around IRS 7B-a, moment 0 and 1 maps of IRS 7B are shown in Figure 7(b). For comparison, the continuum map of IRS 7B is also shown in Figure 7(a). The moment 0 map shows two emission components with semi-elliptic shapes, one located to the southeastern side of IRS 7B-a and the other located on the northwestern side of IRS 7B-a, forming a disk-like structure divided into two along its minor axis as a whole. The C¹⁸O disk-like structure is basically similar to that seen in the continuum associated with IRS 7B-a, although it extends slightly larger than the continuum emission, suggesting that C¹⁸O probably traces gas associated with the protoplanetary disk around IRS 7B-a, but further investigation of the kinematics of C¹⁸O is required for us to have a firm conclusion (see below). Although both continuum and C¹⁸O show similar disk-like structures, their intensity distributions are quite different, e.g., C¹⁸O is stronger in outer parts of the disk-like structure and it gets weaker toward the position of IRS 7B-a, where the continuum emission shows a peak. This discrepancy would be most likely due to absorption of the warm continuum emission by colder C¹⁸O gas.

Zoom In Zoom Out Reset image size

The moment 1 map presented in Figure 7(b) shows a velocity structure of the C¹⁸O emission having blueshifted velocities with respect to the systemic velocity of ∼6 km s⁻¹ on the southeastern side and redshifted velocities on the northwestern side. This clear velocity gradient along the major axis of the disk-like structure is naturally interpreted as rotation. The nature of the velocity structures can be investigated further with a position–velocity (PV) diagram cutting along the major axis, as discussed in Section 6.2. The diagram shows that the velocity of C¹⁸O increases toward the position of IRS 7B-a, a clear sign of differential rotation. Further discussions on the nature of the rotation are presented in Section 6.2.

When we compare our continuum maps from Figure 2 with those toward Class II sources from, e.g., the DSHARP and ODISEA programs, the absence of clear ring or spiral structures identifiable by the naked eye with few exceptions is remarkable. The noticeable difference between protostars and Class II sources could suggest that substructures are not well developed in the embedded phase yet, and they could be quickly formed when Class 0/I protostars evolve into Class II sources. Interestingly, the two sources, L1489 IRS and Oph IRS 63, showing visually identified ring-like structures and another source, IRAS 04169+2702, possibly showing a ring-like structure in the eDisk sample are all Class I protostars. In addition, two protostars previously found to have ring-like structures in their continuum maps are also Class I protostars (WL 17 and GY 91). Note that substructures found in dust continuum around GY 91 are less sharp than those seen around Class II sources (see Sheehan & Eisner 2018; Cieza et al. 2021), while WL 17 does show a sharp ring-like structure (Sheehan & Eisner 2017). Although seven sources from the Orion VANDAM survey exhibit relatively clear ring structures, these ring structures could be due to close binary formation (Sheehan et al. 2020).

We should note, however, that a reason behind the lack of obvious sharp substructures in continuum emission around the Class 0/I protostars could also be that the continuum emission at 225 GHz may be optically thick. In fact, some of the continuum emission observed in the eDisk program may be optically thick as indicated by their higher brightness temperatures (see Table 3). To address this possibility, it will be necessary to observe continuum emission at lower frequencies, where dust continuum emission is expected to be optically thinner. One of the sources in the eDisk sample, L1527 IRS, has been observed in 7 mm continuum with VLA, finding no clear gaps (Sheehan et al. 2022).

We should also note that some of the disk structures seem to be more inclined (see Table 3), particularly as compared with Class II disks observed at high angular resolution to date, which could make substructures more difficult to observe in disks around our sample. In addition, some of the continuum emissions are also more compact, and our angular resolution may not be high enough to detect substructures. More detailed analyses and discussions of substructures in the continuum emission, including the visibility analysis, will be presented in forthcoming papers.

In contrast to substructures, several of the target sources show brightness asymmetries, particularly along their minor axes (see Figure 3). A possible explanation for these brightness asymmetries along the minor axes might come from a geometrical effect, namely, the far side of a highly inclined disk with a finite vertical thickness is brighter than the near side. Such a near-side−far-side brightness asymmetry is an indication that the dust has yet to completely settle onto the disk midplane. In more edge-on systems such as HH 212, the nonsettled dust produces a dark lane near the midplane sandwiched between two brighter features (Lee et al. 2017; Lin et al. 2021). In contrast to embedded disks, Class II disks rarely show any asymmetry along their minor axes (e.g., Andrews et al. 2018; Long et al. 2019). The difference could suggest that dust is not quite settled in the embedded phase but quickly settles as embedded disks evolve into Class II disks.

In addition to investigating substructures in continuum emission, another important purpose of the eDisk program is to identify Keplerian rotation around our target protostars. Here we discuss the case of R CrA IRS 7B-a as an example showing how we identify Keplerian rotation.

One of the methods to identify Keplerian rotation, which is also adopted by the eDisk program, is fitting to PV diagrams. Our fitting process includes two independent methods. One method uses the ridge in 1D intensity profiles along the positional axis (the velocity axis) at a given velocity (position). The other method uses the outer edge in the 1D intensity profiles. We call the former the ridge method, while we call the latter the edge method. In each method, PV diagrams are fitted with a power-law function. The function is a single or double power-law function, which is defined with a power-law index, p_in. The systemic velocity, v_sys, is also included as a free parameter in the fitting process. More details of the fitting method and process are explained in Appendix D.

As was shown in Section 5.2, the C¹⁸O emission in R CrA IRS 7B-a shows a velocity gradient along its disk-like elongated structure centered at the position of R CrA IRS 7B-a, suggestive of rotation. The C¹⁸O PV diagram cutting along the major axis of the dust emission around R CrA IRS 7B-a presented in Figure 8 shows that the rotation can be described as differential rotation, and it is fitted in both the edge and ridge methods described above. The fitting results are presented in Table 4. For demonstration purposes, only fittings with a single power-law function are presented here. Both edge and ridge fittings with a single power law show that p_in is close to 0.5, suggesting that the rotation can be explained as Keplerian rotation. The dynamical mass of the central protostar was estimated to be ∼3.2 M_☉ with the edge method, while it was estimated to be ∼2.1 M_☉ with the ridge method. As explained in Appendix D, the actual dynamical mass of the central star is considered to be between ∼2.1 and ∼3.2 M_☉. Discussions on further details of the Keplerian disk around R CrA IRS 7B-a are beyond the scope of this overview paper, and they will be presented in the subsequent paper.

Zoom In Zoom Out Reset image size