© The Authors 2026

1 Introduction

There is growing evidence that suggests that the planet formation process begins during the embedded protostellar stages (Class 0/I), making the characterization of protostellar disks key to studying both the protostar accretion process and the initial phases of planet formation. The evidence of the early formation of planets comes from high-resolution (~6 au) observations toward more evolved Class II disks (≲1 Myr) that revealed widespread substructures, such as rings and gaps, thought to be linked to planet formation (ALMA Partnership 2015; Andrews et al. 2018a; Long et al. 2018) alongside direct detections of forming planets through continuum emission or gas kinematics (Benisty et al. 2021; Pinte et al. 2018; Izquierdo et al. 2022; Domínguez-Jamett et al. 2025). Another line of evidence of early planet formation in embedded disks is the very low dust masses measured in more evolved Class II disks around stars and brown dwarfs (Manara et al. 2018; Testi et al. 2016; Sanchis et al. 2020). These low values are consistent with models of early pebble growth, migration, and loss through radial drift (Appelgren et al. 2023, 2025), which predict a rapid depletion of the millimeter-sized dust grains in the disk within a few tenths of a megayear, hence also requiring a fast start of the planet formation process. Tentative evidence of a stable or increase in median disk dust masses at around 2–3 Myr is also consistent with early planet formation and subsequent dust regeneration driven by disk-planet interaction (Turrini et al. 2019; Testi et al. 2022; Bernabò et al. 2022; Polychroni et al. 2025). At odds with the results of Appelgren et al. (2025), recent multiwavelength studies have challenged the standard assumption that Class II disks emit optically thin radiation at ~1 mm (Macías et al. 2021; Guerra-Alvarado et al. 2024a; Chung et al. 2024; Painter et al. 2025; Garufi et al. 2025). Accounting for optically thick Class II disks provides a potential alternative scenario that still supports early planetesimal formation since reproducing the observed fluxes and disk sizes would then require the presence of early (≲0.4 Myr) dust traps (Delussu et al. 2024).

Surveys of Class 0/I protostars have progressively improved in resolution, which has helped in the identification and characterization of the dust emission at millimeter wavelengths from these young disks (e.g., Jørgensen et al. 2007; Tobin et al. 2016; Segura-Cox et al. 2018; Maury et al. 2019; Tobin et al. 2020; Ohashi et al. 2023a; Hsieh et al. 2024; Reynolds et al. 2024). Recently, more focus has been placed on determining when dust substructures, similar to those in Class II disks, first appear (Segura-Cox et al. 2020; Nakatani et al. 2020; Ohashi et al. 2023a; Maureira et al. 2024; Guerra-Alvarado et al. 2024b; Hsieh et al. 2025a). Considering only individual studies and surveys with resolutions of ~6–14 au, nearly 60 Class 0/I disks have been observed (Cieza et al. 2021; Ohashi et al. 2023a; Hsieh et al. 2024), and definitive annular substructures have been detected in only five disks: L1489 IRS, ISO-Oph 54, R CrA IRS 2, Ophemb 20, and IRS 63, all of which are Class I disks (Sheehan & Eisner 2018; Segura-Cox et al. 2020; Yamato et al. 2023; Shoshi et al. 2024; Hsieh et al. 2024). Hsieh et al. (2025a) showed that accounting for observational biases, such as inclination and disk size relative to resolution, the detection rate for Class I disks reaches ~60% at bolometric temperatures of ~250 K despite the low apparent numbers. Whether optical depth is what further suppresses observed rates, particularly at earlier times, remains unclear. These results suggest either that planet formation begins during the Class I stage or that many younger disks remain too optically thick at ~1 mm, thus preventing the clear detection of substructures. In fact, Nakatani et al. (2020) reported that a possibly annular substructure in the Class 0/I disk L1527 is detected at 7 mm but not at 1 mm and 3 mm. Similarly, Maureira et al. (2024) found possible evidence of a substructure at 3 mm in a nearly face-on Class 0 disk. Finally, Zamponi et al. (2021) and Guerra-Alvarado et al. (2024b) have found possible evidence of a substructure in two Class 0 disks that are massive and hot and have large scale heights, based on the detection of asymmetries in their profiles. These recent studies support the notion that wavelength-dependent opacity plays a significant role in the observational results. Thus, to build a reliable timeline for planet formation, it is crucial that we quantify the presence and extent of optically thick emission at millimeter wavelengths in young Class 0/I disks.

The implications of optically thick emission extend beyond substructure observability. Dust masses, usually derived assuming optically thin and isothermal emission (Jørgensen et al. 2007; Tobin et al. 2020; Tychoniec et al. 2020), are consequently only lower limits. Optical thickness also biases grain growth estimates, whose sizes may be overestimated if there is significant optically thick emission (Li et al. 2017). Furthermore, young disks may be heated by shocks and accretion (Evans et al. 2015; Zamponi et al. 2021; Xu & Kunz 2021; Maureira et al. 2022; Takakuwa et al. 2024), which can produce very low α values in the presence of optically thick emission, sometimes even below the optically thick limit, without requiring scattering from large grains (Galván-Madrid et al. 2018; Zamponi et al. 2021).

In this study, we present ~7.5 au resolution ALMA continuum observations at 1.3 and 3 mm of Class 0/I disks from the FAUST¹ program (Codella et al. 2021). The high resolution enabled spatially resolved spectral index maps to probe the extent of the optically thick emission. The paper is organized as follows: in Sections 2 and 3, we describe the sample and observations, respectively. In Section 4, we present the continuum maps, spectral index maps, and corresponding radial profiles. In Section 5, we quantify the morphology of the intensity profiles at both wavelengths; reveal correlations between the disk sizes, fluxes, and compare dust and brightness temperatures. In Section 6 we discuss how these results directly imply high fractions of optically thick emission at 1.3 and 3 mm. In this context, we discuss and derive implications for the disk masses, dust temperatures, and grain sizes, including circumbinary disk structures. Section 7 presents a summary and our conclusions.

2 Sample

Protostars in our sample were drawn from the FAUST targets. The only exception is the Class 0 triple system IRAS 16293–2422 (hereafter IRAS 16293) which observations were presented in Zamponi et al. (2021) and Maureira et al. (2020, 2022). We targeted five FAUST fields (RCrA IRS7B, Elias 29, VLA 1623-2417, IRS 63 and L1551 IRS 5) and included in our sample all the Class 0/I sources contained in the respective 3 mm ALMA fields of view. In total, 16 individual disks (7 Class 0,9 Class I) and 3 circumbinary disks (2 Class 0, 1 Class I) were observed at 1.3 mm and 3 mm. In the larger 3 mm field view, we detected 4 additional individual disks (1 Class 0, 3 Class I) and 1 circumbinary disk (Class I). Most sources observed at both wavelengths (14/16) have T_bol < 200 K, making the sample more representative of the early protostellar stages. The sampled molecular clouds and adopted distances are Corona Australis (d = 149 pc, Galli et al. 2020), Taurus (d = 146 pc, Galli et al. 2019) and Ophiuchus. For Ophiuchus, we adopted d = 141 pc (Dzib et al. 2018) for protostars in the eastern part of the cloud corresponding to IRS 63 and IRAS 16293, and d = 137 pc (Ortiz-León et al. 2018) for the protostellar systems toward Ophiuchus A, corresponding to VLA 1623-2417 system (hereafter VLA 1623), Elias 29 and Oph A SM1. Table 1 summarizes the properties of all the sources presented in this work. All listed properties are presented when available in the literature, with the corresponding references. We include a column with the internal luminosity L_int, which is the protostar luminosity without the contribution of UV illumination in the envelope affecting L_bol. L_int is obtained from the relation between the flux at 70 μm and L_int, found in Dunham et al. (2008) through radiative transfer modeling. The flux at 70 μm is taken from the eHOPS catalog (extension of the Herschel Orion protostar survey, or HOPS, Out to 500 ParSecs) updated in Pokhrel et al. (2023)². For the few sources not listed in that catalog (VLA1623 W, OphA SM1, RCrA SMM1C and L1551 IRS 5) the fluxes were obtained from the references provided in Table 1. We also include a column indicating the projected separation to the nearest protostar within the ALMA 3 mm field of view, providing insight into whether the source is relatively isolated or part of a close multiple system. For reference, the last column lists current FAUST papers showing molecular line emission with a resolution of 50 au toward the individual protostars and fields.

3 Observations

The ALMA band 3 and 6 long-baselines observations targeting VLA 1623, RCrA IRS7B, L1551 IRS5, IRS 63, and Elias 29 were taken during 2021 between August and October. The observations were part of the cycle 7 project ID:2019.1.01074.S (PI: Maureira). The most extended configuration C-10 was used for the band 3 observations with baselines ranging from 122 to 16 200 m, while C-8 was used for the band 6 observations with baselines ranging from 470 to 11 500 m. The configuration and frequency setup was chosen to be combined with the observations from FAUST (PI: S. Yamamoto, ID: 2018.1.01205.L). The FAUST observations consisted of two main-array configurations (plus ACA observations for the band 6). Here we only used the most extended configuration for each band, which were taken from October 2018 through March 2020, and corresponded to C-6 for the band 3, while either C-4, C-5 or C-6 for band 6. For the purpose of this project, we refer to the FAUST observations used here as the compact configuration.

The spectral setup in band 3 consisted of four spectral windows centered at frequencies of 93.1, 94.9, 107.0, and 105.1 GHz, with a channel width of 488.281 kHz (~1.5 km s⁻¹) and bandwidth of 1.875 GHz. In band 6, the setup consisted of seven spectral windows centered at frequencies of 217.0, 219.9, 219.5, 218.4, 218.2, 231.6 and 233.7 GHz. The narrower spectral windows between 219.9 GHz and 218.2 GHz were set up for lines following the FAUST setup, with a channel width of 122 kHz (~0.16 km s⁻¹) and bandwidth of 58.6 MHz. The rest of the spectral windows had the same channel width (~0.6 km s⁻¹ at 230 GHz) and bandwidth as the band 3 setup.

We used CASA (CASA Team 2022) to calibrate and image the long baselines observations. Calibration of the raw visibility data was done using the standard pipeline. To create the continuum visibilities, we use the continuum channel ranges provided for the pipeline over all the available spectral windows. We checked that the continuum visibilities were not contaminated by line emission. For this, we produced an cleaned image of the continuum visibilities and afterward appended the model image onto the data visibilities using the ft task in CASA. We then subtracted the model from the data using the CASA task uvsub. Finally, we imaged these residuals checking that only noise was left. For the line rich source L1551 IRS 5 (Bianchi et al. 2020), we also created dirty cubes of the visibilities after line flagging and found no obvious remaining lines. Similarly, for the other line rich sources in this study (IRAS 16293-2422 A/B), the line channels were carefully flagged by hand by inspecting dirty cubes of each spectral window (Maureira et al. 2022). We then averaged in frequency resulting in channel widths of ~10 MHz and ~30 MHz for the band 3 and band 6 data, respectively. These values are conservatively chosen to avoid intensity losses due to bandwidth smearing for sources offset from the phase center.

We performed self-calibration separately for the observations in the extended configuration, which were then combined with the self-calibrated compact configuration observations and self-calibrated again together. Details of the calibration and self-calibration procedures for the FAUST observations can be found in Bianchi et al. (2020) and Maureira et al. (2024). For cleaning and self-calibration we use the CASA tasks tclean and gaincal. The cleaning in the self-calibration iterations was done first with the “multiscale” and later with the “mtmfs” deconvolver parameter, with the parameter nterms=2 for the final phase and amplitude steps. For the parameters “robust” and “uv-taper” we explored values to match the resolution we obtained from the extended configuration observations imaged with a robust of 0.5. This resulted in typical robust values of 0. A manual mask was set and adjusted during the self-calibration process when necessary. Details of the self-calibration and imaging can be found in Appendix A. Table A.1 lists the beam sizes, final rms and imaging parameters of the final continuum images.

3.1 Imaging for spectral index

To obtain spectral index maps it is important to create images at the different bands that recover as closely as possible the same spatial scales. Thus, we performed additional imaging limiting the baselines so as to cover the same uv-wavelength range in both wavelengths. We use the same tclean setup and scales as during the last self-cal imaging described above, except for the robust parameter. We explored “uv-taper” and “robust” parameters in tclean to closely match the synthesized beams for the two wavelengths. Finally, we use the restoring beam option in tclean to produce images with the same synthesized beam. The resultant resolution and rms of these images are summarized in Table A.2 in the Appendix.

3.2 Archival data

We added ALMA band 3/6 observations toward the Class 0 multiple system IRAS 16293-2422 (project ID: 2017.1.01247.S/PI: Dipierro, and ID: 2016.1.00457.S/PI: Oya) that were observed at similar or slightly lower (~13 au) resolution. The self-calibration and imaging of the data (including spectral index images) is similar to the procedures described above and can be found in Zamponi et al. (2021) and Maureira et al. (2022). The resolution and rms of these maps are listed in Tables A.1 and A.2.

4 Results

4.1 Continuum maps

Figures 1 and 2 show a gallery of the continuum images at 3 mm and 1.3 mm for all the sources in our sample, respectively. The emission is presented with a log stretch to enhance the weaker extended emission. The sources are ordered by increasing bolometric temperature. Multiple protostellar systems with projected separation below 100 au are shown separately in the bottom row. All the sources are presented at the same physical scale.

Emission associated with disks is observed for all targets. There is a clear diversity of dust disk sizes ranging from a few au (Elias 29) up to ~100 au (SMM1C). Overall, the morphology of the disks appears similar at both wavelengths, except for the Class I protostar RCrA IRS7A. A compact disk structure is detected at 1.3 mm for this source. In contrast, the emission at 3 mm from the central region is more extended and shows a non-axisymmetric shape. The 3 mm observations for this source also show narrow features extending up to ~50–150 au from the disk which are absent at 1.3 mm. The morphology, position angle and 1.3–3 mm spectral index of these elongated features (α ≲ −0.2–0.7) supports a non-dust origin (e.g., free-free), possibly associated with material being ejected from the protostar (see Appendix B for further details). More subtle differences between the 1.3 mm and 3 mm emission, such as asymmetries along the disk major axis, are visually apparent in some sources. These asymmetries are more pronounced at 1.3 mm than at 3 mm, suggesting they are due to increasing optical depth toward the center or in more inclined disks (Villenave et al. 2020; Lin et al. 2021; Zamponi et al. 2021; Liu 2021; Ohashi et al. 2023a; Guerra-Alvarado et al. 2024a,b).

The observations also show emission associated with circumbinary disks (CBDs). These systems were previously resolved as close binaries in ALMA continuum observations by Harris et al. (2018), Maureira et al. (2020) and Cruz-Sáenz de Miera et al. (2019) for VLA 1623 A, IRAS 16293 A and L1551 IRS 5, respectively. The Class I RCrA SMM2 is resolved here for the first time into a 14 au separation binary surrounded by a circumbinary disk (CBD). CBDs are only detected in binaries with projected separations below 100 au. The CBDs show clear inner gaps in all sources, except the Class 0 IRAS 16293 A1/A2, along with azimuthal asymmetries in their brightness distribution more clearly observed at 1.3 mm.

In Table C.1, we summarize the center coordinates, position angles and inclinations for all individual protostellar disks derived from the optically thinner 3 mm emission. The values are obtained from a 2D Gaussian fit to the 3 mm images using the CASA task imfit. Inclinations were derived assuming a circular geometry. The derived inclinations range from 16° to 82° with a median of 59°.

4.2 Spectral index maps

We constructed spectral index α maps from the 1.3 mm and 3 mm images built with matching uv-range and synthesized beam (Section 3.1). The spectral index maps were created using $$\alpha =\frac{\ln \left(I_{{\nu }_1}/I_{{\nu }_2}\right)}{\ln \left({\nu }_1/{\nu }_2\right)},$$(1)

where ν₁ and ν₂ correspond to 225³ GHz and 100 GHz, respectively. When examining the maps, some sources showed spectral index maps with systematic gradients along a certain direction, which could be caused by the images at 1.3 mm and 3 mm not being correctly aligned. Small misalignments can arise due to the motion of the source at the different epochs, systematic residual phase errors that are different in different bands or calibrators not being exactly at the same positions at the two frequencies due to opacity effects (Kutkin et al. 2014). To correct for this, we calculated the shift between the images by looking at their cross-correlation⁴, using the scipy function fftconvolve. The detected shifts range from 1 to 2 pixels which is equivalent to 7.5–15 mas in the case of the Corona Australis field, and 10–20 mas in the case of the Ophiuchus field containing IRAS 16293 A/B. Most of the systematic gradients dissipated after the aligning procedure. (See Appendix D for further details.) The statistical error in α for each pixel is found to be typically below 0.1 for all disks and the CBD IRAS 16293 A, reaching a maximum uncertainty of 0.2 toward the edges. For the CBDs around L1551 IRS5 N-S and VLA 1623 Aa–Ab, statistical uncertainties for each pixel⁵ are higher (0.2–0.5). Full maps for the statistical uncertainty of α are in Appendix Figure E.1. In addition to the statistical uncertainty, all α maps have a systematic 1σ error of 0.07 from flux calibration⁶.

Figure 3 shows the resultant spectral index maps for the sample. The upper right corner in each panel shows the minimum and maximum value for α in the map. Most disks display widespread low values (α ≈ 2), with α rising to ≳3 toward the edges. While compact disks without clear radial gradients (e.g., Elias 29) may be affected by free-free contamination leading to the observed low α values, sometimes even below 2, similar low α values are also found in larger disks, extending beyond 30 au (e.g., RCrA SMM1C, RCrA IRS7B, VLA1623 W). This implies that free-free contamination, expected only near the central region, cannot fully account for the low observed values. A previously established exception is the compact RCrA IRS7 A disk, where non-dust emission within 10 au yields a mean α of 0.7.

Figure 3 also shows the 1.3–3 mm spectral index for three CBDs. Overall, the α values toward these structures are higher than those found toward the circumstellar disks, either comparable to the values found toward some of the disk edges (α ≈ 3) or larger. In Figure 4, we compare the distribution of α values in the individual disks versus that toward the CBDs computed using a Gaussian kernel density estimator aggregating all pixels in the spectral index maps for all disk and CBD structures. The distributions clearly peak at different values for α. For the disks, the median of α is 2.1 with 68% of the values in the range 1.8–2.6. For the circumbinary structures, the median value is 3.0, with 68% of the values in the range 2.7–3.2. The individual disk with the largest average α is the Class I IRS 63. Therein, α values other than those toward the ring structure (Segura-Cox et al. 2020) and central region, contribute the most to the second peak of the disk density distribution near α ~ 2.6.

Fig. 1 ALMA 3 mm images of all the sources in our sample at the same physical scale. The emission is presented with a log stretch to enhance the weaker extended emission. The images are organized in two groups. The first one (top three rows) corresponds to all the systems in which the projected separation to the nearest protostellar neighbor is larger than 100 au, while the second one (bottom row) comprises the systems with a protostellar neighbor below 100 au. Unlike the first group, disk-like circumbinary structures are observed for all sources in the second group. For each group the sources are organized by increasing bolometric temperature.

Fig. 2 Same as Figure 1 but for the 1.3 mm observations. Some sources are detected at 3 mm and not observed at 1.3 mm due to the smaller field of view.

4.3 Radial intensity and spectral index profiles

We computed the azimuthally averaged intensity profiles for the larger disks in the sample. This subsample corresponds to those for which the geometric mean of the deconvolved FWHM axes from the 2D Gaussian fit is at least 3 times the size of the synthesized beam in at least one of the wavelengths. In order to compare the intensity at 1.3 mm and 3 mm, the profiles are built from the observations with matching uv range and synthesized beam in both wavelengths (Section 3.1) and after alignment correction (Section 4.2). To build the intensity profile we first deprojected the disks using the center coordinates and inclinations from the 2D Gaussian fit to the 3 mm emission (Table C.1). We then calculated the azimuthally averaged profiles from the deprojected images. For the more inclined disks (RCrA SMM1C, VLA 1623 W, VLA 1623 B) we use only a limited range of azimuthal angles (<10–20 degrees) around the disk major axis. Similarly, for the circumstellar disk L1551 IRS5 N, we also excluded some azimuthal angles to the north, to avoid contamination from emission connected to the circumbinary disk (Figure 2). The uncertainty in each bin is calculated as ${\sigma }_{std}/\sqrt{N}$, where σ_std and N are the standard deviation, and number of beams in each bin, respectively.

Figure 5 shows the resultant disk profiles in logarithmic scale where the intensity is given as brightness temperature T_b calculated using the full Planck function. All the profiles show a behavior that resembles the classic viscous disk solution (Lynden-Bell & Pringle 1974), with an inner power-law and outer exponential cutoff behavior. As discussed in Birnstiel & Andrews (2014), at the outer edges of a disk, where the gas surface density decreases rapidly with radius, the drift speeds are fast for particles of any size. This leads to an early (≲0.1 Myr) sharp outer edge in the dust surface density that should be imprinted in the dust continuum emission. Our observations indeed indicate the presence of a sharp outer edge in the dust emission at 1.3 and 3 mm.

The profiles at 1.3 mm and 3 mm for each disk resemble each other very closely for the majority of the sources. Another way to look at the similarity of the profiles at 1.3 and 3 mm is through the spectral index profile, overlaid on the brightness temperature profiles in Figure 5. The spectral index increases with radius, but the profiles show in more detail that α values in agreement with the optically thick limit of 2 (or below) are observed all the way up to near the location of the disk ‘edge’ marked by the sudden drop of emission. It is only near this ‘edge’ where the 3 mm profile shows brightness temperatures below that of 1.3 mm, resulting in α values increasing above 2. Such behavior near and beyond the ‘edge’ in α is expected when the emission becomes optically thin or at least partially optically thin. The only exception to this behavior is the Class I IRS 63, for which this transition to optically thinner emission occurs at smaller radii than in the rest of the sample. The ring detected in this disk at 1.3 mm (see also Segura-Cox et al. 2020) and for the first time here also at 3 mm, falls within this region where α is already above the optically thick limit. Such regions are clearly almost nonexistent for the other disks in the sample.

Toward the inner part of the disks, the trend is that α can take values below two, as seen in the spectral index maps (Figure 3). As expected, this behavior is related to regions for which the 3 mm brightness temperatures is progressively larger than that of 1.3 mm. Such behavior can be expected when the emission is very optically thick at both wavelengths, and can be due to either positive gradients of temperature for increasing depth into the disk (e.g., Li et al. 2017; Galván-Madrid et al. 2018; Zamponi et al. 2021; Lin et al. 2021) or due to self-scattered emission due to the presence of large grains (Liu 2019; Zhu et al. 2019).

Fig. 3 Maps of 1.3–3 mm spectral index for all sources with observations at 1.3 and 3 mm. The ordering is the same as in Figures 1 and 2. The upper corner of each panel shows the minimum and maximum value within the map, considering only pixels with errors below 0.2 (including both statistical and flux calibration uncertainties).

Fig. 4 Spectral index distribution of all disks versus circumbinary material obtained through a kernel density estimator. The vertical lines at the bottom show the mean value of the spectral index for the individual disk sources and circumbinary structures.

5 Analysis

5.1 Fits to the intensity profiles

For more extended disks, for which deprojected intensity profiles were derived (Section 4.3), a 2D Gaussian fit can result in obvious residuals that indicate that this type of intensity profile is not a good model. Thus, in order to better quantify the disk properties and compare them across wavelengths, we fit the deprojected radial intensity profiles for these extended disks with a more general form consisting of two independent power-law indices for constraining the inner and outer cutoff behaviors of the profiles. This form is a variation of the classic solution, and has been used to fit the emission of more evolved Class II disks (e.g., Tazzari et al. 2021a). The intensity profile has the form $$I(r)\propto {\left(\frac{r}{r_c}\right)}^{{\gamma }_1}\exp [-{\left(\frac{r}{r_c}\right)}^{{\gamma }_2}],$$(2)

where r is the radius, r_c is the characteristic radius, γ₁ the inner disk power-law index, and γ₂ the power-law index describing the slope of the exponential cutoff. Together with r_c, γ₁, and γ₂, the fourth free parameter is taken as the total flux of the model, F_tot. Each model was convolved with the synthesized beam, and the values for each of the parameters were obtained using the Python module emcee (Foreman-Mackey et al. 2013). We refer to Maureira et al. (2024) for further details of the emcee setup. The resultant fits overlaid in the observations are shown in Appendix Figures F.1 and F.2, and the resulting values for the free parameters are listed in Table F.1.

For the remaining (more compact disks), we use the results of the 2D Gaussian fit to compute model radial profiles. We use the modified self-similar profile in Equation (2) with γ₁ = 0 and γ₂ = 2, resulting in a Gaussian profile. We compute r_c using the value for the deconvolved FWHM (major axis) as FWHM/$(2\sqrt{\ln 2})$. The final parameter F_tot, is computed as the total flux from the 2D Gaussian fit divided by cos(i), which corresponds to the face-on flux, i.e., the same quantity computed from the modified self-similar profile fit to the extended disks. Although for an optically thin disk the total flux does not change with inclination, the cos(i) correction applies for optically thick emission (Tazzari et al. 2021a) which is already suggested by the observed low α values (Figure 3). We discuss the optically thick nature of the emission throughout the article.

5.2 Disk sizes at 1.3 and 3 mm

In the literature, disk sizes are often defined as the radius containing a certain fraction of the total flux. This allows for a better comparison across studies using different prescriptions to fit the disk emission (Miotello et al. 2023). Typically, the reported sizes correspond to R₆₈ or R₉₅ corresponding to the radius containing 68% and 95% of the total flux, respectively. In principle, R₉₅ is a more robust measurement of the disk, but given its dependence on the fainter and rapidly decaying outer emission, its use is more limited in the literature. Tables F.1 and F.2 summarize the resultant R₆₈ and R₉₅ values derived from the modified self-similar profile and the Gaussian fit, respectively.

Figure 6 shows the comparison of the disk sizes at 1.3 and 3 mm using R₆₈, as well as R₉₅. In both cases, the sizes are almost the same for all disks in the sample. On average, the R₆₈ radii are only 7% smaller at 3 mm compared to 1.3 mm. The agreement between R₉₅ at 1.3 and 3 mm is even better, with an average difference of only 3%. We note that similarly small differences (below 9%) were measured for R₆₈ at 0.9, 1.3 and 3 mm, for a sample of Class II disks in Lupus Tazzari et al. (2021a).

5.3 Spectral index versus disk size

Figure 7 shows the mean spectral index value for the individual disks in our sample as a function of R₆₈ at 3 mm. In order to also investigate possible trends with inclination, we exclude disks for which the spectral index map for the disk is not well resolved⁷ (RCrA IRS7A, Elias 29 and IRAS 16293 A1/A2). There is a trend of increasing spectral index with radius, in agreement with what is observed for Class II disks (Tazzari et al. 2021a; Chung et al. 2024). Likewise, there is a trend with inclination such that for a given size, the mean spectral index is lower for more inclined disks. Such a trend is expected when significant optical depth is playing a role in the values for α. This adds scatter to the relation and can help explaining why the largest disk in our sample (i ~ 78°) drops out of the trend of increasing α with size. Similar scatter and drops from the correlation are observed for Class II disks (Tazzari et al. 2021a; Chung et al. 2024).

5.4 Inner slope at 1.3 and 3 mm

Figure 8 shows the resultant slopes for the inner part of the disks from the modified self-similar fit measured at 1.3 mm and 3 mm. The results show a range of power-law index values ranging from ~−1.2 to ~0. The values are similar at both wavelengths, in agreement with the similarity between the profiles. There is a trend for the values measured at 3 mm to be higher than the one measured at 1.3 mm for each source such that on average the index at 1.3 mm is 24% shallower than at 3 mm.

5.5 Comparison of brightness and dust temperature

Figure 9 shows all the disk intensity profiles at 1.3 mm expressed in brightness temperature T_b, using the results from the model radial profile fits (Section 5.1). At each radius, the color of the curves tracks the resulting model spectral index between 1.3 and 3 mm. The median T_b of the fit profiles at 1 au is ~170 K. At 5 au, where the profiles still show α ≈ 2 values, the median is ~100 K. At a distance of 40 au, where most of the disks still show α ≈ 2 values, the median brightness temperature is ~23 K. For comparison, the dashed lines in Figure 9 correspond to minimum and maximum midplane temperature profiles based on the lowest and highest L_bol in the sample. The midplane temperature as a function of radius is calculated as $$T_{mid}^{\mathrm{i}\mathrm{r}\mathrm{r}}={\left(\frac{\phi L_{\ast }}{8\pi r^2{\sigma }_{SB}}\right)}^{0.25},$$(3)

where σ_SB is the Stefan-Boltzmann constant. The above expression is an approximation valid for passively heated disks (Chiang & Goldreich 1997). To provide a lower and upper limit to the sample midplane temperature, we consider the limiting values 0.1 and 29 L_⊙ for the protostar luminosity L_*, and 0.01 and 0.3 for the flaring angle⁸ φ. The values for L_* are taken from the range of L_bol of the sample (Table 1). These lower and upper limits to L_* and φ correspond to H/R between ~0.05 and ~0.15 for a 1 M_⊙ protostar. We observe that the disk T_b profiles as well as median T_b values for the sample cover the same order of magnitude as the predicted range of dust midplane temperatures for irradiated disks, which does not change if we consider instead L_* = L_int. This supports the interpretation of α ≈ 2 as due to optically thick emission. We also note that the observed T_b values for at least two disks (IRAS 16293 B and L1551 IRS 5 N) are above the derived upper limit for the temperature due to irradiation. This is in agreement with a previous study for IRAS 16293 B in which mechanical heating (dissipation due to accretion/viscosity) was needed to explain the high fluxes (Zamponi et al. 2021). Similarly, L1551 IRS 5 is a FU Ori like source (Mundt et al. 1985), for which viscous heating can also be expected. The comparison between the observed T_b and the predicted values due to irradiation for each individual disk is discussed in detail in Section 6.2.3.

Fig. 5 Intensity profiles on logarithmic scales at 1.3 and 3 mm expressed in brightness temperature, calculated using the full Planck function. The shaded area corresponds to the uncertainty in the mean (Section 4.3). The 1.3–3 mm spectral index α calculated from the profiles is shown in magenta with the corresponding value in the right y-axis. The magenta shaded area shows the uncertainty calculated propagating the errors in the intensity profiles. The vertical dotted line represents the geometric mean of the beam size, indicating that emission beyond this point is spatially well resolved.

Fig. 6 Comparison of disk radii, measured as R68 or R95, between 1.3 and 3 mm. The solid line shows the 1:1 ratio, while the dashed lines show a 10% difference.

Fig. 7 Mean spectral index for the Class 0/I disks in the sample as a function of the radius enclosing 68% of the 3 mm flux. The color of the symbols indicates the disk inclination.

Fig. 8 Index of the inner disk power-law γ1 obtained from the fit to the radial intensity profiles (see Section 5.1 for details). The dashed lines mark the index values of −0.5 and −0.75, indicative of passive and viscous heating, respectively.

5.6 Disk size–luminosity relations

Figure 10 shows R₆₈ versus the disk flux at 1.3 and 3 mm, and compares the observed Class 0/I disks with the relations derived for Class II disks in Tazzari et al. (2021a). The relation is shown as a function of the flux of the disk as seen face-on, and all fluxes are rescaled at a common distance of 150 pc. The inclination correction for the face-on flux for optically thick emission (Section 5.1) was also applied by Tazzari et al. (2021a) for the sample of Class II disks, motivated by the non-negligible optically thick fractions found throughout the sample, and an improvement in the tightness of the relation after applying the correction (Tazzari et al. 2021a), which also occurs in our sample. The comparison reveals that for the observed Class 0/I disks there is also a correlation between R₆₈ and the disk flux at 1.3 and 3 mm. Although the slopes appear to be the same, it is clear that the linear relation for the Class 0, and I disks in our sample is shifted from the one derived for the Class II sources, even after considering the observed scatter in the Class II relation Tazzari et al. (2021a). For a given size, the fluxes for the Class 0/I disks in our sample are ~10× higher.

Fig. 9 Intensity profiles at 1.3 mm for all disks resulting from either a Gaussian fit (for the most compact disks) or a modified self-similar profile fit (Section 5.1). The colors of the curves correspond to the resultant 1.3–3 mm spectral index at each radii, except RCrA IRS7A due to the non-dust emission at 3 mm (Appendix B). The circles show the median value for the brightness temperature at certain radii. The dashed lines correspond to the theoretical approximation for the midplane temperature due to irradiation (T ∝ r−0.5) considering the range of bolometric luminosities of the sample (see Section 5.5 for details).

Fig. 10 Size and luminosity relations at 1.3 and 3 mm (circles). Left: R68 vs. flux at 1.3 mm. Right: R68 vs. flux at 3 mm. Fluxes in both panels have a correction for inclination that is valid for optically thick emission. The dashed lines show the linear regression results for Class II sources at the two wavelengths, derived in Tazzari et al. (2021a), which also include the correction for optically thick emission. The dotted lines corresponds to the relation for Class II considering a shift in the fluxes of a factor of 10.

6 Discussion

6.1 Optically thick emission and comparison with Class II disks

The observations reveal that the observed sample of Class 0/I disks, encompassing different molecular clouds, multiplicity, and even dust disk sizes, have significant optically thick emission at 1.3 and 3 mm. This is supported by (a) the low spectral index values for the disks with a median of α = ${2.1}_{-0.2}^{+0.6}$, along with per disk median values all below α ≲2.6, (b) that the median α value for a given size tends to be lower for higher inclinations, and (c) the elevated brightness temperatures with values in agreement (or even above) the predicted temperatures due to irradiation. Results (a) and (b) are very similar to those found toward Class II disks. Low values for the spectral index (~2) have also been measured at different submm/mm wavelengths for large samples of Class II disks across different regions (Tazzari et al. 2021b; Chung et al. 2024; Garufi et al. 2025; Painter et al. 2025). In addition, Tazzari et al. (2021a), Chung et al. (2024) and Garufi et al. (2025) also find that the average per disk spectral index increases with disk size in a similar fashion as our observations of Class 0 and I disks.

However, result (c) shows differences between the Class 0/I sample and observations of Class II disks. The brightness temperatures of Class II disks appear to be lower than those of the Class 0 and I disks in this work. For instance, Andrews et al. (2018b) measures a median of ~5 K at 40 au at 0.98 mm, while the median of our sample at that distance is about 5× higher (~23 K) at 1.3 mm (Figure 9). Likewise, similarly low values $({\mathrm{T}}_b^{1.3mm}\sim 1-10\mathrm{K})$ are observed at 40 au in the DSHARP sample (Huang et al. 2018). These low observed brightness temperatures in the Class II disks have led some studies to conclude that most of the disk emission is optically thin at these wavelengths, even for disks observed at high-resolution (e.g., Huang et al. 2018). However, low brightness temperatures can also be explained if the emission is actually optically thick and dust scattering opacities are significant (Liu 2019; Zhu et al. 2019), or if there are unresolved optically thick structures with a filling factor below 1, for instance optically thick rings with optically thin gaps (Ricci et al. 2012; Andrews et al. 2018b). The latter scenarios are indeed in agreement with multi-wavelength observations of Class II disks mapped with a resolution comparable to the one presented here (Zhu et al. 2019; Carrasco-González et al. 2019; Guidi et al. 2022).

In the case of larger samples of Class II disks, where the resolution of the observations does not allow the presence of substructures to be probed, it has been shown that the size–luminosity relation for Class II disks can be explained if less than half of the flux is optically thick (Andrews et al. 2018b). Andrews et al. (2018b) discuss that higher fractions made the disks too luminous for a given size. Since we have shown that the Class 0 and I disks in our sample are brighter for a given size that Class II disks, we discuss whether this difference can be explained with a higher fraction of the Class 0/I disk’s area being optically thick, in line with the resolved dust emission and spectral index maps presented in this work. To test this idea, we followed Andrews et al. (2018b) and compared the observed location of the Class 0/I disks in the size-luminosity relation with the predicted location while assuming the disks are fully optically thick. To do a simple first comparison, we followed Andrews et al. (2018b) and assumed a temperature radial profile given by irradiation with a power law in radius r and protostar luminosity L_* = L_int. The results in Andrews et al. (2018b) can be explained using T = T₀(r/r₀)^−0.5(L_*/L_⊙)^0.25 with T₀ = 30 K at r₀ = 10 au for the temperature parametrization. This comes from the approximation of midplane temperature in passive disks where T ∝ (φL_*/r²)^0.25, assuming a constant flaring angle, φ (Chiang & Goldreich 1997). We then assumed a radial intensity profile given by $$I_v=\{\begin{array}{ll}{\mathcal{F}}_{\text{thick}}\cdot B_{\nu }(T)& r\le R_0\\ 0& r>R_0\end{array},$$(4)

where ℱ_thick is a constant that takes values between zero and one and thus reflects the fraction of the disk flux that is optically thick. Following Andrews et al. (2018b), R₀ is taken as the characteristic radius r_c, as beyond this radius typically less than 10–20% of the total flux is emitted. In cases where a larger percentage of the flux is outside r_c, we set R₀ = R₉₀. For single protostars we use directly L_* = L_int with the values in Table 1. For close multiples for which we only have L_int values for the entire system, we assume that each member has the same L_*, and equal to the one measured for the entire system. Thus, for those disks, the predicted optically thick disk luminosities are a very conservative upper limit. As a check, we showed that if the sample of Class 0/I disks were to have ℱ_thick = 0.3 the predicted location in the size luminosity relation would match that of the Class II disks (see Figure G.1 in the Appendix), thus supporting the methodology and the need for a higher optically thick fraction to be able to reach the higher luminosities observed for the younger disks.

Figure 11 top panels show the predicted location of the Class 0/I disks if we set ℱ_thick = 1, i.e., fully optically thick disks. The predicted location is marked with empty squares, while the observed location is marked with filled circles. The predicted location more or less matches the observed location for some of the disks, but there are several disks that remain too luminous. To zoom into the difference between predicted and observed luminosity, the bottom panels in Figure 11 show the ratio between the observed disk flux and the one predicted using ℱ_thick = 1 for each disk as a function of size. Observed fluxes that are at least 2× higher than what is predicted by fully optically thick emission are observed in 6 out of 13 disks at 1.3 mm, and 7 out of 15 disks at 3 mm, thus ~46% at both wavelengths. This percentage is really a lower limit considering that for several disks, those in closer multiple systems, the ratio is only a lower limit given that we have assigned the luminosity of the system to each member. The agreement between the predicted and observed luminosity seems to be size-dependent such that the luminosities are matched for the larger disks (R₆₈ ≳40–50 au), while discrepancies are more apparent on the more compact disks. We checked if a flaring might be the reason for the discrepancies and find that while a non-constant flaring does help, it cannot not remove all the discrepancies (see Appendix G for details). Assuming L_* = L_bol along with a radius-dependent flaring still does not completely resolve the discrepancies for IRAS 16293B, Elias 29, and L1551 IRS5N disks (see Figure G.3), which is a conservative number of sources given the lower limits due to the assumed luminosities in multiple systems.

A factor of two or more is difficult to explain by possible errors in L_* or flaring as it would require at least a factor of 2^1/0.25 = 16 higher than L_* or flaring, assuming Rayleigh Jeans and the temperature prescription in Equation (3). An extreme case is the single Class I disk Elias 29, the smallest disk in our sample. At 1.3 mm, this disk is from 2 to 5× brighter than its fully optically thick prediction, considering all above temperature parametrizations. The results suggest that the temperatures in these disks are higher than what is expected from irradiation alone, and additional heating such as viscous heating might need to be considered. As this type of heating is expected to dominate at smaller radii, it would also naturally explain why smaller disks show the higher discrepancies between the predicted and observed fluxes.

Fig. 11 Comparison between the observed disk luminosity and the predicted values assuming the disks are fully optically thick (ℱthick = 1). Top: same as Figure 10. Filled circles are the observations and squares are the predicted luminosity for each disk if the emission is fully optically thick. Bottom: size versus the ratio of the observed disk luminosity to the predicted one assuming the entire disk is optically thick. Lower limits for the ratio are due to the assumed internal luminosity being an upper limit.

6.2 Optical depths at 1.3 and 3 mm

Here we attempt to quantify and discuss the optical depth without the need for L_*. Using the modeled disks intensity profiles⁹ and resultant spectral index (Figure 9), we can compute a profile for the optical depth at 3 mm τ_{3 mm}. For this, we use a modified black-body for the intensity profile I_ν ~ B_ν(T)(1 − e^−τ_ν). For the temperature, we assume it to be equal to the observed brightness temperature until r(α = 2). For larger radii we extrapolated this value assuming a r^−0.5 profile. We find that all disks, except IRS 63, have at least 70% (and up to near 100%) of their 3 mm flux emitted from regions with τ_{3 mm} ≳ 1. The median extent of the τ_{3 mm} ≳ 1 region is the disk characteristic radius r_c. These values are just a lower limit at 1.3 mm. Assuming β ~ 1 (dust emissivity index) the sample has a median of ~87% of the disk flux with τ_1.3mm above 1, in line with the previous discussion for which we assumed a temperature profile based on L_int (Section 6.1). These results are summarized in Appendix Figure H.1.

6.2.1 Disk masses and dynamical state

To provide estimates of the disk masses, we use the above derived optical depth profiles at 3 mm, extrapolating inward from the maximum derived near r(α = 2) considering a r^p profile with p = 0 and p = −1.5, as a lower and higher mass limiting cases, respectively. Then, the mass of the disk is calculated as $M_{\text{disk\hspace{0.17em}}}=2\pi {\int }_0^R\mathrm{\Sigma }(r)rdr$, where Σ(r) = τ_{3 mm}(r)/κ_{3 mm} and κ_{3 mm} is the dust opacity, and R = r_c. We used an opacity at 3 mm of ~ 1 cm² g⁻¹ (Beckwith et al. 1990; Birnstiel et al. 2018). Figure 12 shows the estimated range of disk masses in solids and for the gas assuming a gas-to-dust ratio of 100. The mass in solids spans four orders of magnitude, from ~2 M_⊕ (Elias 29, R₉₅ ≈ 2 au) up to a few times 10³ M_⊕ (RCrA SMM1C, R₉₅ ≈ 100 au). Most of the disk dust masses are in the range 30–900 M_⊕, corresponding to gas masses in the range 0.01–0.3 M_⊙. The values are in agreement with the estimates in Tychoniec et al. (2020) for Class 0/I disks in Perseus using observations at 9 mm assuming optically thin conditions, with an opacity value consistent with the one we used at 3 mm if extrapolated using β ~ 1. The derived range of Class 0/I disk masses in Perseus was found to be comparable or above the masses in observed exoplanets (Tychoniec et al. 2020), and thus a similar conclusion can be applied to the masses derived here. The results show that these relatively high masses, required to explain exoplanetary systems, are not unique to a particular cloud, and that it is important to consider a high fraction of optically thick emission if masses are estimated for young disks from continuum ALMA observations, even at 3 mm. Compared with numerical simulations, the estimated masses also agree well with the median disk gas mass of ~0.04–0.06 M_⊙ found in the synthetic populations in Bate (2018) and Lebreuilly et al. (2021).

We caution that an important source of uncertainty in this mass estimation, as well as others in the literature, is the value for the dust opacity at a particular wavelength, which depends on the dust composition, morphology and distribution of grain sizes (Testi et al. 2014; Miotello et al. 2023; Birnstiel 2024; Liu et al. 2024). As a reference for some of the most significant variations, the estimated masses would increase by a factor of 2–3 if we adopted DSHARP opacities with a maximum grain size of a_max ~ 1 mm Birnstiel et al. (2018), Ricci et al. (2010) opacities with a_max ≲ 100 μm, or DIANA opacities with a_max ≲ 300 μm (Woitke et al. 2016). Conversely, the masses would decrease by approximately a factor of 2 if we used Ricci et al. (2010) opacities with a_max ~ 1 mm. Another poorly constrained source of uncertainty in observations of young Class 0/I disks are possible changes of the gas-to-dust ratio (Lebreuilly et al. 2020; Miotello et al. 2023; Ohashi et al. 2023b).

The derived mass range in Figure 12, considering the above uncertainties, also agrees better with Class 0/I disk mass estimates in Orion from single-wavelength (λ ≲ 1.3 mm) observations that are based on marginally gravitationally unstable disk models (Xu 2022). In contrast, fully passive models give lower masses and no unstable disks using the same observations (Sheehan et al. 2022). To explore the dynamical state of our disks, Figure 12 compares the derived disk masses with protostellar masses from the literature (Table 1). When no measurement was available, we assumed M_* = 1 M_⊙ (empty symbols). For close binaries for which only the combined protostar mass is reported in the literature, we assume equal mass for each component. Figure 12 shows that seven disks, including some in close binaries, approach or exceed the gravitational instability limit M_disk/M_* ~ 0.1 (Kratter & Lodato 2016). In the case of IRAS 16293 B, this matches the results in Zamponi et al. (2021) where a model of a hot and massive unstable disk was able to reproduce the 1.3–3 mm fluxes, similar to the case of the Class 0/I edge-on disks L1527 IRS using observation from 0.9 mm to 7 mm (Ohashi et al. 2022b). Disks with M_disk/M_* ≳ 0.1 also appear in simulations (Bate 2018; Lebreuilly et al. 2021), though less frequently when magnetic fields and non-ideal MHD are included (Lebreuilly et al. 2021). One might then wonder why these or many other young disks do not show for instance spiral arms features. One possibility, discussed in Xu (2022), is that the high optical depths limit the detectability of such features, as the observations cannot see all the way through the disk, plus the perturbations in density might be only of order unity. In addition, the disks might be only marginally gravitationally unstable with Toomre Q parameter ~1–2 (Toomre 1964; Kratter & Lodato 2016). In this case, spiral perturbations do not grow exponentially. Using our surface density and temperature profiles, we compute Q profiles assuming flat or r^−1.5 optical depth profiles, thus following the same assumptions employed for the mass constraints¹⁰ (Figure 13). Even under conservative assumptions, several disks reach Q ~ 1–2 or less at least in the outer parts, supporting the scenario that gravitational instabilities are playing a role during the protostellar stages.

Fig. 12 Estimated disk masses versus the protostar mass. The right y-axis corresponds to the dust mass, while the left y-axis corresponds to gas mass calculated assuming a gas-to-dust ratio of 100. The solid line represents a criterion for gravitational instability such that disks above the line are unstable. Empty symbols mark sources with no estimates in the literature of the stellar mass for which we assume M* = 1 M⊙. For close binaries for which only the combined protostar mass is reported in the literature we assigned half that mass to each component.

6.2.2 Dust substructure and molecular line detection

The only clear substructure in our sample is the ring in the Class I disk IRS 63 (Figure 2), previously identified by Segura-Cox et al. (2020), which uniquely shows α values above the optically thick limit well before the disk edge in our sample. Another potential substructure appears in the Class 0 disk OphA SM1 intensity profile at 3 mm, showing a deviation in the radial profile near 30 au (Figure 5), which could be due to an annular substructure (Maureira et al. 2024). Notably, the IRS 63 ring lies beyond the optically thick region, while the SM1 deviation also occurs where τ_3mm < 1 (Maureira et al. 2024). These results support the idea that annular substructures can emerge as early as the Class 0 stage but are often hidden by optically thick emission. In line with this, Hsieh et al. (2025a) find that substructure detection sharply increases to ~60% at T_bol ~ 200–400 K (Class I) for features directly visible in the 1.3 mm images, while the detection rate drops to zero for Class 0 disks. The latter show slightly elevated average brightness-to-dust temperature ratios, suggesting higher optical depths. In Class 0 disks, large scale heights may further obscure gaps; for example, bumps in the Class 0 NGC 1333 IRAS4A profile can be explained with a gap in a hot and highly flared disk (Guerra-Alvarado et al. 2024b).

The higher optical depths at 1.3 and 3 mm also complicates kinematic and chemical studies. Optically thick dust can cause molecular lines above the disk to appear in absorption (e.g., Pineda et al. 2012; Oya et al. 2018; Sahu et al. 2019; Garufi et al. 2021), or remain undetected as part or all molecular emission can be masked within the optically thick dust, leading to under-estimated abundances. VLA observations already demonstrate that longer wavelengths can reveal species missed in ALMA ~1 mm observations (López-Sepulcre et al. 2017; De Simone et al. 2020). Thus, abundance measurements at disk scales must account for high dust optical depths (Hsieh et al. 2025b). Future high-sensitivity facilities such as SKA¹¹ and ngVLA¹² will be crucial to advance studies of young disk kinematics and chemistry.

Fig. 13 Profiles of the Toomre Q parameter as a function of radius. The left panel shows the case where we extrapolate the optical depth in order to remain constant inward of the radius for which α ≈ 2, while the right panel assumes the optical depth increases as r−1.5 therein, thus providing an upper and lower limit for Q. The curves are colored by the resultant dust optical depth at 3 mm. The horizontal dotted line corresponds to Q = 1.7, a reference value below which the disk can develop gravitational instabilities. (See Section 6.2.1 for more details.

6.2.3 Slope of the temperature radial profile

Given that the disk emission is optically thick, the radial profile of the brightness temperature can, in principle, trace the temperature radial profile and thus inform us about the dominant heating mechanism. According to our analysis in Section 5.5, the power-law index of the brightness temperature with radius spans from approximately -1.2 to 0 (Figure 8). These values extend both above and well below theoretical expectations for passive (−0.5) and viscous (−0.75) heating. Moreover, the slopes at 3 mm tend to be steeper than at 1.3 mm, consistent with lower optical depths at longer wavelengths, regardless of whether the disk is marginally or fully optically thick.

The variation in slope across the sample, including several disks with nearly flat profiles, may indicate significant differences in absolute optical depth. Flatter profiles likely correspond to very high optical depths. To understand this, we can consider two idealized cases under the Rayleigh-Jeans and modified blackbody approximations: optically thin emission for which T_b ∝ r^q+p with q and p the power-law index of the temperature and surface density as a function of radius, and optically thick emission for which T_b ∝ r^q. It follows that when emission goes from optically thick to thin, the brightness profile becomes steeper. This is consistent with the Class I disk IRS 63, which shows the steepest slope (~ −1) and has low spectral index values (α ≲ 2) confined only to the inner disk, suggesting reduced optical depth in outer regions (Figures 5 and 8). However, in extremely optically thick disks, the temperature observed at each wavelength corresponds to the τ ~ 1 surface, which can vary in shape and altitude across wavelengths. In such cases, we may not probe the midplane at all, preventing us from probing the radial temperature profile. For instance, the edge-on Class 0 disk HH212 shows flat brightness profiles that steepen with increasing wavelength (Lin et al. 2021), due to the τ ~ 1 surface becoming less flat at longer wavelengths. Their modeling inferred a physical temperature profile with q ~ −0.7, despite the observed flatness. Similarly, the nearly face-on Class 0 disk IRAS 16293B shows flat profiles. Zamponi et al. (2021) demonstrated that even in such orientation, the τ ~ 1 surfaces at 1.3 mm and 3 mm can lie above the midplane and remain relatively flat in the inner regions. Based on these studies, detailed modeling and multi-wavelength observations are necessary to constrain the true temperature power-law index. Still, some disks in our sample already show 3 mm slopes near −0.75 (e.g., RCrA IRS7B a, RCrA SMM1C) or −0.5 (VLA 1623B), suggesting that q may vary among disks even at similar evolutionary stages.

6.3 Dust temperatures and iceline locations

Given that the high optical depths allowed us to measure dust temperatures, we estimated where different icelines would be located in our Class 0/I disk sample. When performing this estimation and based on the discussion in Sect. 6.2.3, it was important to keep in mind that the high optical depths in some disks likely prevent us from tracing the midplane temperatures, and thus the radial location discussed here is only a first approximation.

To obtain the location of an iceline at temperature T_ice, we used the $T_b^{1.3mm}$ profiles obtained from the model intensity profiles fit (Section 5.1). We find the radius r_ice such that $T_b^{1.3mm}=T_{\text{ice}}$ provided that α ≤ 2. When α > 2, the value for r_ice is taken as a lower limit. When T_ice > max(T_b) throughout the profile, we set r_ice = 1 au as an upper limit. In Figure 14, the top panel shows the computed values for r_ice for difference icelines for all the disks in the sample with 1.3 mm observations¹³. The adopted values for T_ice should be considered approximate, as they are sensitive to variations in ice composition and local density conditions (Martín-Doménech et al. 2014; Fayolle et al. 2016; Potapov et al. 2019; van ’t Hoff & Bergner 2024).

We find that the median location¹⁴ of the water iceline (T_ice=150 K) is r_ice ~ 3 au. There is significant spread across the sample, which can be partially associated to the range of disk sizes in the sample from ~1 to ~100 au (symbol sizes). Only upper limits are found for several disks (IRS7B b, SMM1C, CXO 34, VLA 1623 Aa/Ab), as well as values as high as 8 au (L1551 IRS 5 S), and lower limits of 13 au and 23 au, for L1551 IRS5 N and IRAS 16293 B, respectively. Figure 14 bottom panel shows the location of the icelines as a function of the ratio r_ice/r_c, in order to show what fraction of the disk size we can expect T_d > T_ice. There is a group of disks for which the water iceline is located only in the inner region of the disks (≲0.1 r_c, IRS 63, IRS7B a and b, SMM1C, VLA1623 B and CXO 34), and another group of disks have values ≳0.3r_c (IRS 7A, Elias 29, IRAS 16293 B, L1551 IRS 5 N and S, VLA1623 Aa and Ab). Despite the larger disk size fraction in the latter disks with dust temperatures above 100–150 K, emission from complex organic molecules (COMs) have only been detected for the largest hotter disks: IRAS 16293 B (Chandler et al. 2005; Pineda et al. 2012; Jørgensen et al. 2016), and L1551 IRS 5 N/S (Bianchi et al. 2020). We note that with the exception of IRAS 16293-2422, all sources in the sample are part of FAUST and therefore have dedicated observations aimed at revealing the emission of COMs. For the compact sources RCrA IRS7A (R₉₅ ~ 7 au) and Elias 29 (R₉₅ ~ 2 au), further inspection of the ALMA spectral setups from FAUST reveals no bright, compact sources associated with COM emission at a resolution of 50 au. This could be due to the high dust optical depths, but also due to the compact disk size in these sources, as both make the detection of complex molecules more difficult. Overall, COMs have been detected in our sample only in sources for which the size of the region with dust temperature above 100 K is computed to be at least 10 au. In the case of Elias 29, Oya et al. (2019) also proposed that the lack of COMs toward this source is due to a more elevated envelope temperature (>20 K) during the prestellar phase which could prevent efficient depletion of the parent species of COMs (e.g., CO). Thus, the lack of COM detection in Elias 29 and IRS 7A could also be related to an intrinsic chemical difference originated in the environment (Oya et al. 2025).

The other, colder icelines in Figure 14, follow similar trends. The median r_ice values for the N₂ (T_ice=20 K) and CO (T_ice=25 K) icelines are 13 au and 9 au respectively. The location of these two icelines in most of the disks is located beyond the characteristic disk radius, which is consistent with the extensive use of C¹⁸O for tracing keplerian rotation in young disks.

Fig. 14 Location of different icelines in the observed Class 0/I disks derived from the fit intensity profiles at 1.3 mm. Top: location plotted as distance from the center of the disk. Bottom panel: location plotted as the radius of the iceline over the characteristic disk radius rc. The ‘◃’ and ‘▹’ symbols correspond to upper and lower limits, respectively. The sizes of the symbols are proportional to log10 rc. The median value for each iceline (ignoring upper/lower limits) is marked with a diamond symbol. The numbers of disks considered for the median is indicated next to diamond marker at the top panel. (See Section 6.3 for more details.

6.4 Disk sizes at 1.3 and 3 mm and grain growth

The high optical depth in the observed disks prevents us from directly using the resolved spectral index to infer the dust emissivity index β, a proxy for the maximum grain sizes, throughout the disk. However, we note that the outer edges of several disks in our sample show a steep α_{1.3–3 mm} gradient (Figure 3). This rapid increase, with α values changing from ~2 to ~3, is typically associated with the location of the falloff in the intensity profiles (Figure 5) and hence related to the transition from optically thick to thin emission. From this, we can conclude that the optically thin dust in the outer disk region is in agreement with β ≈ α − 2 ~ 1 (or higher considering non R–J effects at lower temperatures). Such values do not rule out the presence of grains that are a few tens of micrometers up to millimeters in size, depending on the adopted opacity (Testi et al. 2014; Birnstiel 2024; Zamponi et al. 2024; Carrasco-González et al. 2019). Observations at longer wavelengths with the VLA, and in the future ngVLA/SKA, are required to overcome the high optical depth, thus allowing us to reveal the grain size distribution in these young disks. Indeed, Radley et al. (2025) recently derived the presence of mm-sized grains in the Class 0 and I disks in the VLA 1623 A/B and W system (Figures 1 and 2), by combining ALMA high-resolution observations at 1.3 and 3 mm, with VLA observations at 7 mm, 1.4 cm and 3 cm.

An independent constraint to the grain sizes can be obtained by looking at how the disk size changes with wavelength. In this scenario, the observed sharp outer disk edge in the dust emission profile is not tracing a drop of the dust surface density, but instead a sudden change in the dust opacity, related to the presence of maximum grain sizes of about a few hundred micrometers interior to the intensity profile falloff (Rosotti et al. 2019). This is because the dust opacity at a given wavelength increases sharply with maximum grain size at around a_max ~ λ/2π. Thus, the observed disk edge or ‘knee’ at 1.3 mm and 3 mm would be due to grain sizes larger than a few hundred micrometers interior to the knee and smaller beyond that location. In that scenario, the disk radius is also smaller for longer wavelengths, as the sudden increase in the opacity moves with wavelength (Tazzari et al. 2021a). In Figure 10, we showed that both R₆₈ and R₉₅ are comparable to within 4–7% for the Class 0/I disks studied here. Tazzari et al. (2021a) finds also small differences (below 9%) for R₆₈ between 0.9, 1.3 and 3 mm for Class II disks in Lupus. They concluded that the small difference could be reproduced by considering constant β values in the range 0.5–1 over a substantial part of the disks with possible larger grains only in the very inner regions, i.e., not strong variations in the grain sizes with radius for most of the disk emitting area. According to the models presented in Tazzari et al. (2021a) (see their Figure 8), the similarity in sizes between 1.3 and 3 mm, the high fraction of optically thick emission and the associated typical spectral index of ~2 obtained for our Class 0/I disk sample can be all reproduced by also considering β values in the range 0.5–1, similar to the conclusions for the Class II disks.

6.5 Class 0/I circumbinary disks

Figure 1 shows that CBDs are detected only in systems with projected separations below 100 au (14–54 au in our sample), while wider multiples (≥103 au) lack CBDs. This matches simulations showing that the fraction of CBDs with ages up to 0.1 Myr drops rapidly from about 50% for binaries with semimajor axes a in the range 10–100 au, to almost zero for a ≳ 100 au (Bate 2019; Elsender et al. 2023). The same trend holds when considering other ALMA high-resolution Class 0/I surveys (eDisk; Ohashi et al. 2023a, CAMPOS; Hsieh et al. 2024) as well a sample of multiples in Perseus (Reynolds et al. 2024). Although we do not have a measurement of a for the majority of these close Class 0/I multiples (see Maureira et al. 2020 and Hernández Garnica et al. 2024 for exceptions), the lack of CBDs above ~100 au projected separation strongly supports a real drop around the predicted threshold of a ~ 100 au (Elsender et al. 2023).

6.5.1 Class 0/I CBDs optical depths and masses

For the three CBDs for which we have observations at 1.3 and 3 mm (VLA 1623 A, IRAS 16293 A, and L1551 IRS 5) we find that their mean spectral index are in the range 2.6–3.1 (Figure 4). Such a difference between the spectral index in CBDs versus circumstellar disks implies a change in the optical depth in these populations, with lower values for CBDs. Maureira et al. (2022) studied in detail the 1.3–3 mm spectral index in IRAS 16293 A CBD. By assuming a constant value of β throughout the CBD, the optical depth at 1.3 mm is estimated in the range ~0.4–1, corresponding to 0.2–0.3 at 3 mm. This is indeed much lower than the typical optical depths at 3 mm in the circumstellar disks discussed above (see Figure 13) but still comparable with the values observed at the edges of the circumstellar disk. We estimated the optical depth for the CBDs in the sample by comparing the predicted temperatures due to irradiation with the observed intensities at 1.3 and 3 mm. We performed this comparison over a cut oriented along the CBD’s major axis, and passing through the mid-point between the individual CSDs. Details about the intensity cuts and the temperature profile can be found in Appendix I. The resultant optical depths are below 0.5 for L1551 IRS 5 CBD and VLA 1623 CBD at both wavelengths, while for IRAS 16293 A CBD the optical depth at 1.3 mm is higher, reaching values close to 1, in agreement with previous constraints (Maureira et al. 2022). The intensity, temperature and optical depth cut profiles are summarized in Figure I.1.

Having estimated the optical depth in these CBDs, we can measure the amount of material in young Class 0/I CBDs. Such measurements are important given the detections of circumbinary planets around binary stars (Doyle et al. 2011; Martin 2018; Standing et al. 2023). For VLA1623A and L1551 IRS5 CBDs, we measure the 1.3 mm integrated flux in the CBDs within a contour above 3σ, and use the optically thin approximation, adjusting for the temperature and β based on the above discussion. Details can be found in Appendix I. The rate of mass underestimation due to the assumption of optically thin emission is less than 25% for these optically thinner structures, which is less than the typical uncertainties due to the choice of dust opacity. The resultant dust masses are 94 M_⊕ for VLA1623A CBD, and 78 M_⊕ for L1551 IRS 5 CBD, this corresponds to ~0.03 M_⊙, and ~0.02 M_⊙ in gas (or ~26 M_Jup), assuming a gas-to-dust ratio of 100. These masses are comparable to, or higher than the ones derived for the CSDs in these systems (Figure 12). As shown in Maureira et al. (2022), the circumbinary disk around IRAS 16293 A shows additional spots of bright emission surrounding the CSDs, which are consistent with localized shock heating. Taking into consideration this complex temperature structure, Maureira et al. (2022) derived a range of masses for this CBD of 0.02–0.04 M_⊙ for the gas (~20–40 M_Jup) or 67–133 M_⊕ in solids¹⁵. Overall, the gas mass range of ~20–40 M_Jup for these Class 0/I CBDs is at least comparable to that of the currently detected circumbinary planets with masses between 0.1–13 M_Jup¹⁶.

6.5.2 Dust emissivity and grain sizes in the Class 0/I CBDs

We used a modified black body to derive β profiles along the major axis of the CBDs, using the optical depth and temperature profiles derived in the previous section and the full Planck function. The results are shown in Figure 15, where we binned the β profiles in bins of 20 au. The mean β_CBD for VLA 1623 A and IRAS 16293 A (Class 0) are 1.37 ± 0.15 and 1.51 ± 0.01, respectively. For these binaries, the derived values can be taken as a conservative lower limit given that the internal luminosities of the binaries therein also include a third protostellar source nearby (see Table 1). For the Class I L1551 IRS5, the mean β_CBD is 1.1 ± 0.07. The results for the CBDs suggests a possible trend of higher β values for the Class 0 sources. The measurements of β for the Class 0 CBDs suggests either ISM-like grains or grains that have grown up to a few hundred μm, as similar β are expected in both scenarios (Agurto-Gangas et al. 2019; Ricci et al. 2010). On the other hand, the Class I circumbinary disk L1551 IRS 5 could host larger mm-sized grains, based on the estimated lower β values. This is in agreement with recent findings of low β ≲ 1 values in the dusty and wide outflow cavity of this source, extending up to 2000 au (Sabatini et al. 2025). These grains might then have been lifted from the circumbinary disk. Finally, we note that the β_CBD for VLA 1623 A and L1551 IRS5 are also in agreement with the values found for the envelope around these sources corresponding to β_env ~ 1.43 ± 0.05 and β_env ~ 0.94 ± 0.04, respectively (Cacciapuoti et al. 2025).

Fig. 15 Derived β1.3–3 mm for the CBDs in the sample. The values are calculated for a cut centered at the midpoint between the CSDs, oriented along the CBD major axis. The cuts are displayed in the top panels. The dotted lines correspond to the mean value for each source. Positive offsets correspond to the eastern side of the cut. (See Section 6.5.2 for details.

7 Summary and conclusions

We have presented and analyzed 1.3 and 3 mm ALMA continuum observations with a median resolution of 7.5 au toward a sample of 16 Class 0 and I disks in single and multiple systems and across several molecular clouds. The sample is more representative of the early protostellar stages, with 14 of the sources having bolometric temperatures lower than 200 K. In addition, about a third of the disks analyzed in this sample (Elias 29, L1551 IRS 5 N/S and IRAS 16293 A1/A2) have bolometric luminosities greater than 20 L⊙. Such luminosities are higher than those typically derived for Class II pre-main-sequence stars, but they may be indicative of protostars undergoing higher accretion rates, as predicted for early stages (Dunham et al. 2014) and confirmed by direct measurements in some Class I protostars (e.g., Fiorellino et al. 2023). The observations revealed individual disks with a variety of sizes ranging from R₉₅ ≲ 3 au to 100 au. The observations also revealed four CBD structures reaching ~100 au scales. The observations allowed us to build resolved spectral index maps to investigate the optical depth in the disk structures. Our results and conclusions are summarized as follows:

Class 0/I disks. We find that the disks in our sample are mostly optically thick at 1.3 and 3 mm. The elevated optical depth results in most of the estimated disk masses being in the range 30–900 M_⊕, corresponding to gas¹⁷ masses in the range 0.01–0.3 M_⊙. Such values allow for a fraction of the disks to be marginally gravitationally unstable, with derived Toomre Q values ≲ 2 for the outer disk regions. Some disks appear to have temperatures that are in excess of what is expected from irradiation alone. The median location for the water iceline (T_ice = 150) K is ~3 au, with lower limits for the hotter disks from 10 to 20 au. This suggests the need to account for other types of heating, such as accretion or viscous heating for the hotter disks.

Comparison with Class II disks. Class 0/I disks exhibit low spectral indices (α_CSD ~ 2) similar to Class II disks and follow the same trend of increasing α with disk size. As for Class II disks, Class 0/I disk sizes at 1.3 and 3 mm differ by less than 10%, indicating a minimal wavelength-dependent structure in this wavelength range. A clear correlation between disk size and millimeter luminosity is observed, with a slope comparable to that of Class II disks. However, Class 0/I disks are ten times brighter at a given size, which is consistent with their higher brightness temperatures. This luminosity offset can be explained with a higher fraction (~1) of optically thick emission for the younger disks, though some cases may also require higher temperatures. Consequently, assuming similar dust properties, Class 0/I disks likely host higher dust masses than their more evolved Class II counterparts.

Class 0/I circumbinary disks. Circumbinary disks are detected only around binaries with projected separations of ≤100 au, consistent with other observations from the literature and numerical simulations indicating a sharp decline in the occurrence of CBD for larger separations. These CBDs typically exhibit optical depths ≲1 at 1.3 and 3 mm, which accounts for their higher spectral indices (α_CBD ~ 3) compared to individual disks. The estimated CBD masses range from 80 to 130 M_Earth and are comparable to or exceed those of known circumbinary planets. The CBD dust emissivity indices, β, derived from 1.3 and 3 mm data fall between 1–1.5, with the Class I CBD L1551 IRS 5 showing β ~1, suggesting larger millimeter-sized grains, while the two Class 0 CBDs (VLA1623 and IRAS 16293 A) exhibit slightly higher β values (~1.4–1.5), which can be explained with ISM-like grain populations.

The derived high optical depths and high optical depth fractions at both 1.3 and 3 mm toward the sample of individual Class 0/I disks strongly support the possibility that annular substructures related to the planet formation process may emerge as early as the Class 0 stage, but they remain difficult to detect with current observations (e.g., Maureira et al. 2024; Hsieh et al. 2025b). Optically thick dust also hinders the estimation of grain sizes within these disks and poses challenges for studying the disk kinematics and chemistry at the location of the optically thick dust. In this study, we could only obtain an estimation of ${\beta }_{\text{disk}}^{\text{edge}}\gtrsim 1$ at the disk edge based on ${\alpha }_{\text{disk\hspace{0.17em}}}^{\text{edge\hspace{0.17em}}}\gtrsim 3$ therein. Observations at longer wavelengths are necessary to overcome these issues (e.g., De Simone et al. 2020; Ohashi et al. 2022b; Radley et al. 2025), and thus future observations with SKAO and ngVLA along with more sensitive observations with ALMA to reach wider and fainter populations will be key for advancing our understanding of early disk and planet formation and evolution.

Acknowledgements

M.J.M. is grateful to Munan Gong, Kedron Silsbee, Matthew Bate, Til Birnstiel, and Joaquin Zamponi for valuable conversations that contributed to the improvement of this work. This work was supported by the Max-Planck Society. D.J. is supported by NRC Canada and by an NSERC Discovery Grant. N.C. acknowledges support from the European Research Council (ERC) under the European Union Horizon Europe programme (grant agreement No. 101042275, project Stellar-MADE). L.P., E.B., G.S., and Cl.Co. acknowledge the PRIN-MUR 2020 BEYOND-2p (Astrochemistry beyond the second period elements, Prot. 2020AFB3FX), the project ASI-Astrobiologia 2023 MIGLIORA (Modeling Chemical Complexity, F83C23000800005), the INAF-GO 2024 fundings ICES, the INAF-GO 2023 fundings PROTO-SKA (Exploiting ALMA data to study planet forming disks: preparing the advent of SKA, C13C23000770005), the INAF Mini-Grant 2022 “Chemical Origins” (PI: L. Podio), the INAF Mini-Grant 2023 TRIESTE (“TRacing the chemIcal hEritage of our originS: from proTostars to planEts”; PI: G. Sabatini) and the National Recovery and Resilience Plan (NRRP), Mission 4, Component 2, Investment 1.1, Call for tender No. 104 published on 2.2.2022 by the Italian Ministry of University and Research (MUR), funded by the European Union – NextGenerationEU– Project Title 2022JC2Y93 Chemical Origins: linking the fossil composition of the Solar System with the chemistry of protoplanetary disks – CUP J53D23001600006 – Grant Assignment Decree No. 962 adopted on 30.06.2023 by the Italian Ministry of Ministry of University and Research (MUR). EB aknowledges contribution of the Next Generation EU funds within the National Recovery and Resilience Plan (PNRR), Mission 4 – Education and Research, Component 2 – From Research to Business (M4C2), Investment Line 3.1 – Strengthening and creation of Research Infrastructures, Project IR0000034 – “STILES – Strengthening the Italian Leadership in ELT and SKA”.

References

Appendix A Self-calibration and imaging

We self-calibrated the long-baselines observations iteratively using CASA and the taks tclean and gaincal. We use the ’multiscale’ deconvolver with 3 or 4 scales depending on the source for the first cleaning (pre-selfcal) as well as the ones performed after each phase-only self-calibration solution is applied. For both phase and amplitude self-calibration, we image with a robust parameter of 0.5. We started with phase-only self-calibration, performing seven iterations with decreasing solint intervals from ’inf’ (with combine=scan) down to ’int’. Then, for each source and band, we kept the solution with the shortest solint interval that would still result in improved S/N and visible reduction of artifacts. In most cases this interval was 30 s, but in somes cases smaller solint still resulted in improvements (Elias 29 and VLA 1623 both in band 6). After applying these solutions, we performed another round of cleaning using the ’mtmfs’ deconvolver with nterms= 2 which generally resulted in further improved image fidelity. After this process of phase-only self-calibration the S/N improved from 10% up to 40% in the band 3, and from 15% up to a factor of 3.5 in the band 6. We then continued with amplitude self-calibration also using the’mtmfs’ deconvolver with nterms= 2. We performed two steps with solint values of inf and 600s. In each step the phase solution was updated using the shortest successful solint value from the phase-only process. We kept the amp+phase selfcal solution only if the S/N remained the same or improved, there were no artifacts or changes to the flux scale (with no other apparent change in the image). Only a few observations met this criterion. For all of them, the solution with solint=inf was used and the improvement in S/N over the phase-only solution was only marginal (1-5%).

We then self-calibrated the long baselines and FAUST compact configuration continuum observations together. The self-calibration for the FAUST observations consisted of rounds of solint intervals with’int’ for phase-only solutions, followed by amplitude+phase solutions computed per-scan. For the combined self-calibration, we averaged both datasets in time and frequency. We averaged in frequency up to the maximum that would not result in a reduction of the intensity larger than ~1% for sources away from the phase center¹⁸. The final channel widths was 30 MHz or 60 MHz, depending on where additional sources appeared in the field of view. Due to a PI error in the input of the coordinates for the long-baselines observations, for all targets except IRS 63 there was an offset between the phase center of the compact configuration and the extended configuration. The shift varied from 0.2" up to 6". To account for this, we used the ’mosaic’ gridder option in tclean with mosweight=False in all the imaging (including self-calibration) for these targets. Similar to the long-baselines only procedure, we used 3 or 4 scales for cleaning. We explored robust and uv-taper values to match the resolution we obtained from the extended configuration observations. This resulted in typical robust values of 0. We performed six successive steps of phase-only self-calibration. For this first round, we used the final image after self-calibration of the extended configuration as the model with “inf" as the solint interval. This aligned the two datasets to the peak position of the extended configuration and thus corrects for differences in the peak positions across datasets arising from proper motion or astrometric errors of the observations. The resultant model was then used for the next steps reducing the solint intervals from 60s to 9s, always keeping the last solution that resulted in improved S/N and image fidelity. After these solutions are applied we performed another round of cleaning using the “mtmfs" deconvolver with nterms=2. After this process of phase-only self-calibration the S/N improved only marginally (few percent) in some cases and up to 20% in others. The improvements were similar in both bands. We then continued with amplitude self-calibration also using the “mtmfs" deconvolver with nterms= 2. We first explored three consecutive steps for amplitude self-calibration. In all of them we use solint= inf. In the first and second one we combine all spectral windows and use solnorm=True and False, respectively. For the third and final we keep solnorm=False and attempt a solution per spectral window. By making solnorm=True in the first round we are normalizing the solution to an average to correct for the temporal fluctuations without changing the flux scale. In half of the cases, this process failed at the first step, with very obvious artifacts in the resultant image. For the other half, we kept the solutions of the third step. For the ones that failed, we tried amplitude self-calibration using only the two last steps (i.e., keeping solnorm=False), which produced acceptable results for two more cases for which we kept the solutions per spectral window. For the rest, amplitude was not applied as the resultant images always showed significant new artifacts, reduced S/N or changes in the overall flux scale. For all the cases in which amplitude self-calibration was applied the S/N remained the same or improved from a few percent up to 10%. In some cases (e.g., RCrA IRS7B field band 6) the final image has very weak background emission which contributes to the noise but that, after comparing with the compact configuration only data, it correlates with real extended emission. This was also observed in some fields of the ALMA-IMF Large Program (Ginsburg et al. 2022). We extended the mask in the cases in which this extended emission was fairly simple and present at the source position and performed a final clean. Table A.1 summarizes the final rms and beam sizes of the final images.

Appendix B RCrA IRS 7A non-dust emission at 3 mm

The continuum emission toward RCrA IRS 7A at 3 mm in this study shows two elongated jet-like features to the south-east (SE) and west (W) of the protostar (Figure 1). The peak intensity of the SE feature is 86.7 μmJy beam⁻¹ (~5.7 K), while for the W feature is 182 μmJy beam⁻¹ (~9.8 K) measured in the image using CARTA (Comrie et al. 2021). Such features are not detected at 1.3 mm (Figure 2). Using the 3σ value for the undetected 1.3 mm intensity (σ ≈ 52 μmJy beam⁻¹), the upper limit for the 1.3 - 3mm spectral index for the SE and W features are α ≲ 0.7 and α ≲ −0.2, respectively. Such a low and even negative spectral index is in agreement with free-free emission from a radio jet (Anglada et al. 2018).

We also inspected the fluxes measured with the most extended configuration of the FAUST observations (~0.3" or 45 au. The compact source in those observation shows some tails of extended emission in the SE and W directions at both wavelenghts (Sabatini et al. 2024). The derived spectral index from the integrated flux in the SE and W features therein are α ~ 0.3 and α ~ 1.6, respectively. While the values for the SE feature are low and in agreement with free-free emission, higher value in the lower resolution observations toward the W feature are possibly due to more extended emission from dust emission being detected in this more compact configuration, which becomes weaker in the higher-resolution observations.

Toward the disk, we derive a mean spectral index of ~0.6, which also supports significant contribution from non-dust emission in the disk (Figure 3), in agreement with previous 3.0-3.7 cm observations with the VLA (Liu et al. 2014) as well as the detection in the VLA Sky Survey (VLASS) at 2-4 GHz with a resolution of 2-3" (Lacy et al. 2020).

Appendix C Disk center, position angle, and inclinations

Table C.1 summarizes the coordinates for the center of the disk from a 2D Gaussian fit to the 3 mm emission. The respective position angles and inclination derived assuming circular geometry.

Appendix D Alignment between 1.3 mm and 3 mm

We use the cross-correlation between the 1.3 mm and 3 mm images in order to correct for small offsets which led to visible gradients in the spectral index maps for sources in the field of Corona Australis and IRAS 16293 A/B in Ophiuchus. The cross-correlation was done individually for the sources in these fields to account for both systematic shifts in the field as well as shifts due to individual astrometric motions. Table D.1 summarizes the resultant shifts of the 1.3 mm image with respect to the 3 mm image. Figure D.1 shows the spectral index before and after applying the corresponding shifts. The most visually obvious gradients (IRS7B, IRS7A and CXO 24) are removed after the alignment. In the case of IRAS 16293 B some asymmetry remains, which is expected given that the source has an azimuthal asymmetry at 1.3 mm and 3 mm (Figures 2 and 1). RCrA SMM1C shows no major change, while the circumbinary emission toward IRAS 16293 A shows a smoother distribution overall. We note that none of the results in this work are significantly affected by this small correction.

Fig. D.1 Spectral index maps before and after applying a shift to the 1.3 mm image to correct for systematic gradients observed in the fields of Corona Australis and IRAS 16293 A/B (Ophiuchus). The left and right images correspond to before and after the correction, respectively. The upper-right corner contains the minimum and maximum values of the spectral index in the map, considering only the pixels with errors up to 0.2.

Appendix E Spectral index error maps

Figure E.1 shows maps of the statistical uncertainty of the 1.3-3 mm spectral index maps in Figure 3.

Fig. E.1 Spectral index (1.3-3 mm) statistical uncertainty maps corresponding to the spectral index maps in Figure 3.

Appendix F Fit to intensity profiles of extended disks

Figures F.1 and F.2 show the fit to the intensity profiles and residuals for the observations at 1.3 mm and 3 mm, respectively.

Fig. F.1 Fits to the 1.3 mm intensity profiles and residuals. The orange curves at the top of each panel show 100 profiles drawn randomly from the posterior distributions. The observed profile is shown with a black solid line. The resultant median residual is shown at the bottom of each panel (solid black line) as well as 1σ dispersion around the median residual (orange shaded area). The gray shaded area corresponds to the 1σ uncertainty of the intensity profile in all panels.

Fig. F.2 Same as Figure F.1 but for the 3 mm intensity profiles.

Appendix G Other predictions for the size-luminosity relation

To investigate further the issue of the disks being more luminous than fully optically thick passive heating predicts, we performed another test pushing the temperature to higher values by taking into account that the flaring angle φ also changes as a function of the star luminosity and radius of the disk (Chiang & Goldreich 1997). Following Chiang & Goldreich (1997) to compute φ as a function of r and L_*, we arrive to a slightly different parametrization for the temperature such that T = T₀(r/r₀)^−0.43(L_*/L_⊙)^0.29. For a 1 M_⊙¹⁹ and 1 L_⊙ star we get that T₀ ≈ 50 K at r₀ = 10 au. With this new parametrization we get higher temperatures thus getting closer to the observed location of the Class 0 and I disks in the size-luminosity relation (Figure G.2). In Figure G.3 we repeat this exercise but considering L_* = L_bol, instead of L_* = L_int. In both cases the discrepancy decreases but still persist for at least several sources.

Fig. G.1 Comparison between the observed disk luminosity and the predicted values assuming ℱthick = 0.3 as in Andrews et al. (2018b) for Class II disks. Same as Figure 10. Filled circles are the observations and squares are the predicted luminosity for each disk if ℱthick = 0.3 of the flux if optically thick.

Fig. G.2 Comparison between the observed disks luminosity and the predicted values assuming the disks are fully optically thick, considering a nonconstant flaring angle φ and L* = Lint. Top: Same as Figure 10. Filled circles are the observations and squares are the predicted luminosity for each disk if the emission is fully optically thick. Bottom: Size versus the ratio of the observed disk luminosity to the predicted one assuming the entire disk is optically thick. Lower limits for the ratio are due to the assumed internal luminosity being an upper limit.

Fig. G.3 Comparison between the observed disks luminosity and the predicted values assuming the disks are fully optically thick, considering a nonconstant flaring angle φ and L* = Lbol. Top: Same as Figure 10. Filled circles are the observations and squares are the predicted luminosity for each disk if the emission is fully optically thick. Bottom: Size versus the ratio of the observed disk luminosity to the predicted one assuming the entire disk is optically thick. Lower limits for the ratio are due to the assumed luminosity being an upper limit.

Appendix H Optically thick fractions at 1.3 and 3 mm

Figure H.1 shows the fractions of the total flux and fraction of the disk radius containing emission with τ > 1 (left panel) at 1.3 and 3 mm.

Fig. H.1 Fractions of the total flux and fractions of the disk characteristic radius rc containing emission with τ > 1 at 1.3 and 3 mm. Each symbol corresponds to an individual disk.

Appendix I Optical depth and masses in the circumbinary disks

To estimate the optical depth in the circumbinary disk structures, we performed intensity cuts along the direction of the major axis of the circumbinary disk structures, centered at the midpoint between the circumstellar disks. The intensity at each location is averaged over a beam. For the P.A. we used 36°, 48°, and 165°, for VLA1623A, IRAS 16293A, and L1551 IRS 5, respectively. We masked out emission from the individual CBDs at the center. In the case of IRAS 16293A we also masked out the emission coming from hot dust spots features which are due to shocks, and thus their temperature cannot be modeled by only considering irradiation. The resultant intensity profiles are shown in Figure I.1. To infer the optical depth we assumed that that the intensity is given by a modified black body and assumed a temperature profile from irradiation following the parametrization discussed in discussed in Appendix G. The luminosities were taken from Table 1. The resultant optical depths are shown in the bottom panels of Figure I.1. We note that using other parametrization more typical of envelopes (e.g., Galametz et al. 2019) results in higher temperatures and thus even lower optical depths. However, assuming a colder temperature parametrization such as that in Andrews et al. (2018b), also discussed in Section 6.1 in this work, results in temperatures in some regions that go below the observed brightness temperatures. In order to improve the derivation of the temperature and optical depth, in a more model independent way, more wavelengths are required.

To estimate the mass in the VLA1623A and L1551 IRS5 CBDs, we measured the 1.3 mm integrated flux in the CBDs and used the optically thin approximation $$M=\frac{d^2{S}_{\nu }}{B_{\nu }\left(T_d\right){\kappa }_{\nu }},$$(I.1)

where S_ν is the integrated flux density, B_ν is the Planck function, T_d the dust temperature, and d is the distance. Given the low optical depths for these CBDs (τ_1.3mm ≲ 0.5 Appendix I), the fraction of mass underestimation due to optically thin assumptions is less than 25%, which is lower than the typical uncertainties due to the choice of dust opacity. For the temperature, we calculated the mean value along the major axis weighted by the 1.3 mm intensity (see Figure I.1). The obtained values are T_d = 17 K and T_d = 44 K for VLA1623 and L1551 IRS5, respectively. To be consistent with the CSD mass measurements in this work, we use the same value for κ_{3 mm} = 1 cm² g⁻¹ and extrapolate this value to 1.3 mm (223 GHz) according to the mean β values measured for these CBDs (Section 6.5.2). These are β = 1.1 for L1551 IRS5 and β = 1.37 for VLA 1623. For IRAS 16293 we used the mass values derived in Maureira et al. (2022) and corrected for the opacity as above using β = 1.5. This yielded values that are 3.7× smaller than those in Maureira et al. (2022).

Fig. I.1 Intensity profiles in a cut along the major axis direction for the three CBDs in our sample observed at 1.3 and 3mm. The zero-offset corresponds to the mid-point between the circumstellar disks. Positive offsets track positive R.A. offsets. The red dotted line is the predicted temperature profile due to irradiation and the bottom panels show the resultant optical depth at each wavelength based on this temperature profile. The emission from the circumstellar disks or temperature substructures in the case of IRAS 162923A are masked out from the central region.

All Tables

All Figures

Fig. 1 ALMA 3 mm images of all the sources in our sample at the same physical scale. The emission is presented with a log stretch to enhance the weaker extended emission. The images are organized in two groups. The first one (top three rows) corresponds to all the systems in which the projected separation to the nearest protostellar neighbor is larger than 100 au, while the second one (bottom row) comprises the systems with a protostellar neighbor below 100 au. Unlike the first group, disk-like circumbinary structures are observed for all sources in the second group. For each group the sources are organized by increasing bolometric temperature.

In the text

	Fig. 2 Same as Figure 1 but for the 1.3 mm observations. Some sources are detected at 3 mm and not observed at 1.3 mm due to the smaller field of view.
In the text

	Fig. 3 Maps of 1.3–3 mm spectral index for all sources with observations at 1.3 and 3 mm. The ordering is the same as in Figures 1 and 2. The upper corner of each panel shows the minimum and maximum value within the map, considering only pixels with errors below 0.2 (including both statistical and flux calibration uncertainties).
In the text

	Fig. 4 Spectral index distribution of all disks versus circumbinary material obtained through a kernel density estimator. The vertical lines at the bottom show the mean value of the spectral index for the individual disk sources and circumbinary structures.
In the text

Fig. 5 Intensity profiles on logarithmic scales at 1.3 and 3 mm expressed in brightness temperature, calculated using the full Planck function. The shaded area corresponds to the uncertainty in the mean (Section 4.3). The 1.3–3 mm spectral index α calculated from the profiles is shown in magenta with the corresponding value in the right y-axis. The magenta shaded area shows the uncertainty calculated propagating the errors in the intensity profiles. The vertical dotted line represents the geometric mean of the beam size, indicating that emission beyond this point is spatially well resolved.

In the text

	Fig. 6 Comparison of disk radii, measured as R68 or R95, between 1.3 and 3 mm. The solid line shows the 1:1 ratio, while the dashed lines show a 10% difference.
In the text

	Fig. 7 Mean spectral index for the Class 0/I disks in the sample as a function of the radius enclosing 68% of the 3 mm flux. The color of the symbols indicates the disk inclination.
In the text

	Fig. 8 Index of the inner disk power-law γ1 obtained from the fit to the radial intensity profiles (see Section 5.1 for details). The dashed lines mark the index values of −0.5 and −0.75, indicative of passive and viscous heating, respectively.
In the text

Fig. 9 Intensity profiles at 1.3 mm for all disks resulting from either a Gaussian fit (for the most compact disks) or a modified self-similar profile fit (Section 5.1). The colors of the curves correspond to the resultant 1.3–3 mm spectral index at each radii, except RCrA IRS7A due to the non-dust emission at 3 mm (Appendix B). The circles show the median value for the brightness temperature at certain radii. The dashed lines correspond to the theoretical approximation for the midplane temperature due to irradiation (T ∝ r−0.5) considering the range of bolometric luminosities of the sample (see Section 5.5 for details).

In the text

Fig. 10 Size and luminosity relations at 1.3 and 3 mm (circles). Left: R68 vs. flux at 1.3 mm. Right: R68 vs. flux at 3 mm. Fluxes in both panels have a correction for inclination that is valid for optically thick emission. The dashed lines show the linear regression results for Class II sources at the two wavelengths, derived in Tazzari et al. (2021a), which also include the correction for optically thick emission. The dotted lines corresponds to the relation for Class II considering a shift in the fluxes of a factor of 10.

In the text

Fig. 11 Comparison between the observed disk luminosity and the predicted values assuming the disks are fully optically thick (ℱthick = 1). Top: same as Figure 10. Filled circles are the observations and squares are the predicted luminosity for each disk if the emission is fully optically thick. Bottom: size versus the ratio of the observed disk luminosity to the predicted one assuming the entire disk is optically thick. Lower limits for the ratio are due to the assumed internal luminosity being an upper limit.

In the text

Fig. 12 Estimated disk masses versus the protostar mass. The right y-axis corresponds to the dust mass, while the left y-axis corresponds to gas mass calculated assuming a gas-to-dust ratio of 100. The solid line represents a criterion for gravitational instability such that disks above the line are unstable. Empty symbols mark sources with no estimates in the literature of the stellar mass for which we assume M* = 1 M⊙. For close binaries for which only the combined protostar mass is reported in the literature we assigned half that mass to each component.

In the text

Fig. 13 Profiles of the Toomre Q parameter as a function of radius. The left panel shows the case where we extrapolate the optical depth in order to remain constant inward of the radius for which α ≈ 2, while the right panel assumes the optical depth increases as r−1.5 therein, thus providing an upper and lower limit for Q. The curves are colored by the resultant dust optical depth at 3 mm. The horizontal dotted line corresponds to Q = 1.7, a reference value below which the disk can develop gravitational instabilities. (See Section 6.2.1 for more details.

In the text

Fig. 14 Location of different icelines in the observed Class 0/I disks derived from the fit intensity profiles at 1.3 mm. Top: location plotted as distance from the center of the disk. Bottom panel: location plotted as the radius of the iceline over the characteristic disk radius rc. The ‘◃’ and ‘▹’ symbols correspond to upper and lower limits, respectively. The sizes of the symbols are proportional to log10 rc. The median value for each iceline (ignoring upper/lower limits) is marked with a diamond symbol. The numbers of disks considered for the median is indicated next to diamond marker at the top panel. (See Section 6.3 for more details.

In the text

	Fig. 15 Derived β1.3–3 mm for the CBDs in the sample. The values are calculated for a cut centered at the midpoint between the CSDs, oriented along the CBD major axis. The cuts are displayed in the top panels. The dotted lines correspond to the mean value for each source. Positive offsets correspond to the eastern side of the cut. (See Section 6.5.2 for details.
In the text

Fig. D.1 Spectral index maps before and after applying a shift to the 1.3 mm image to correct for systematic gradients observed in the fields of Corona Australis and IRAS 16293 A/B (Ophiuchus). The left and right images correspond to before and after the correction, respectively. The upper-right corner contains the minimum and maximum values of the spectral index in the map, considering only the pixels with errors up to 0.2.

In the text

	Fig. E.1 Spectral index (1.3-3 mm) statistical uncertainty maps corresponding to the spectral index maps in Figure 3.
In the text

Fig. F.1 Fits to the 1.3 mm intensity profiles and residuals. The orange curves at the top of each panel show 100 profiles drawn randomly from the posterior distributions. The observed profile is shown with a black solid line. The resultant median residual is shown at the bottom of each panel (solid black line) as well as 1σ dispersion around the median residual (orange shaded area). The gray shaded area corresponds to the 1σ uncertainty of the intensity profile in all panels.

In the text

	Fig. F.2 Same as Figure F.1 but for the 3 mm intensity profiles.
In the text

	Fig. G.1 Comparison between the observed disk luminosity and the predicted values assuming ℱthick = 0.3 as in Andrews et al. (2018b) for Class II disks. Same as Figure 10. Filled circles are the observations and squares are the predicted luminosity for each disk if ℱthick = 0.3 of the flux if optically thick.
In the text

Fig. G.2 Comparison between the observed disks luminosity and the predicted values assuming the disks are fully optically thick, considering a nonconstant flaring angle φ and L* = Lint. Top: Same as Figure 10. Filled circles are the observations and squares are the predicted luminosity for each disk if the emission is fully optically thick. Bottom: Size versus the ratio of the observed disk luminosity to the predicted one assuming the entire disk is optically thick. Lower limits for the ratio are due to the assumed internal luminosity being an upper limit.

In the text

Fig. G.3 Comparison between the observed disks luminosity and the predicted values assuming the disks are fully optically thick, considering a nonconstant flaring angle φ and L* = Lbol. Top: Same as Figure 10. Filled circles are the observations and squares are the predicted luminosity for each disk if the emission is fully optically thick. Bottom: Size versus the ratio of the observed disk luminosity to the predicted one assuming the entire disk is optically thick. Lower limits for the ratio are due to the assumed luminosity being an upper limit.

In the text

	Fig. H.1 Fractions of the total flux and fractions of the disk characteristic radius rc containing emission with τ > 1 at 1.3 and 3 mm. Each symbol corresponds to an individual disk.
In the text

Fig. I.1 Intensity profiles in a cut along the major axis direction for the three CBDs in our sample observed at 1.3 and 3mm. The zero-offset corresponds to the mid-point between the circumstellar disks. Positive offsets track positive R.A. offsets. The red dotted line is the predicted temperature profile due to irradiation and the bottom panels show the resultant optical depth at each wavelength based on this temperature profile. The emission from the circumstellar disks or temperature substructures in the case of IRAS 162923A are masked out from the central region.

In the text